1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief chi_b(2S) decays
  9  class CLEOII_1992_I32611 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_1992_I32611);
 14
 15
 16    /// @name Analysis methods
 17    ///@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(),"UFS");
 23      // Book histograms
 24      // averages
 25      _h_N_aver  = {Profile1DPtr(),Profile1DPtr(),Profile1DPtr()};
 26      _h_R2_aver = {Profile1DPtr(),Profile1DPtr(),Profile1DPtr()};
 27      book(_h_N_aver [0],1,3,1);
 28      book(_h_N_aver [1],1,2,1);
 29      book(_h_N_aver [2],1,1,1);
 30      book(_h_R2_aver[0],1,3,2);
 31      book(_h_R2_aver[1],1,2,2);
 32      book(_h_R2_aver[2],1,1,2);
 33      // dists
 34      _h_N  = {Histo1DPtr(),Histo1DPtr(),Histo1DPtr()};
 35      _h_R2 = {Histo1DPtr(),Histo1DPtr(),Histo1DPtr()};
 36      book(_h_N [0],2,1,1);
 37      book(_h_N [1],2,1,2);
 38      book(_h_N [2],2,1,3);
 39      book(_h_R2[0],3,1,1);
 40      book(_h_R2[1],3,1,2);
 41      book(_h_R2[2],3,1,3);
 42    }
 43
 44    void findDecayProducts(Particle parent, Particles & children, unsigned int & nCharged) {
 45      for(const Particle & p: parent.children()) {
 46	if(p.children().empty()) {
 47	  if(isCharged(p)) ++nCharged;
 48	  children.push_back(p);
 49	}
 50	else
 51	  findDecayProducts(p,children,nCharged);
 52      }
 53    }
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      Particles chib = apply<UnstableParticles>(event,"UFS").particles(Cuts::pid==110551 or
 58								       Cuts::pid==120553 or
 59								       Cuts::pid==100555);
 60      for(const Particle & p : chib) {
 61	unsigned int iHist = (p.pid()%10)/2;
 62	unsigned int nCharged(0);
 63	Particles children;
 64	findDecayProducts(p,children,nCharged);
 65	// ncharged
 66	_h_N     [iHist]->fill(nCharged);
 67	_h_N_aver[iHist]->fill(0.5,nCharged);
 68	// R_2
 69	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 70	vector<FourMomentum> mom;
 71	mom.reserve(children.size());
 72	for(const Particle & p2: children) {
 73	  mom.push_back(boost.transform(p2.momentum()));
 74	}
 75	// compute R2
 76	double H0(0.),H2(0.);
 77	for(const FourMomentum & p1:mom) {
 78	  double mod1 = p1.p3().mod();
 79	  Vector3 axis = p1.p3().unit();
 80	  for(const FourMomentum & p2:mom) {
 81	    double mod2 = p2.p3().mod();
 82	    double cTheta = axis.dot(p2.p3().unit());
 83	    H0 += mod1*mod2;
 84	    H2 += mod1*mod2*0.5*(3.*sqr(cTheta)-1.);
 85	  }
 86	}
 87	double R2=H2/H0;
 88	_h_R2     [iHist]->fill(R2);
 89	_h_R2_aver[iHist]->fill(0.5,R2);
 90      }
 91    }
 92
 93
 94    /// Normalise histograms etc., after the run
 95    void finalize() {
 96      for(unsigned int ix=0;ix<3;++ix) {
 97	normalize( _h_N [ix]);
 98	normalize( _h_R2[ix]);
 99      }
100    }
101
102    ///@}
103
104
105    /// @name Histograms
106    ///@{
107    vector<Histo1DPtr> _h_N,_h_R2;
108    vector<Profile1DPtr> _h_N_aver,_h_R2_aver;
109    ///@}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(CLEOII_1992_I32611);
116
117}