1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Upsilon(2S) -> Upsilon(1S) pi+pi-
  9  class CLEOII_1998_I467642 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_1998_I467642);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(),"UFS");
 24      book(_h_charged_ex   ,1,1,1);
 25      book(_h_charged_inc  ,2,1,1);
 26      book(_h_neutral_ex   ,3,1,1);
 27      book( _h_lept_ctheta ,4,1,1);
 28      book(_h_lept_phi     ,4,1,2);
 29      book(_h_piStar_ctheta,5,1,1);
 30      book(_h_piStar_phi   ,5,1,2);
 31      book(_h_pi_ctheta    ,6,1,1);
 32      book(_h_pi_phi       ,6,1,2);
 33
 34    }
 35
 36    void findDecayProducts(const Particle & mother,
 37			   unsigned int & nstable,
 38			   Particles& pip, Particles& pim,
 39			   Particles& pi0, Particles & onium) {
 40      for(const Particle & p : mother.children()) {
 41        int id = p.pid();
 42      	if ( id == PID::PIMINUS) {
 43	  pim.push_back(p);
 44	  ++nstable;
 45	}
 46       	else if (id == PID::PIPLUS) {
 47       	  pip.push_back(p);
 48       	  ++nstable;
 49       	}
 50       	else if (id == PID::PI0) {
 51       	  pi0.push_back(p);
 52       	  ++nstable;
 53       	}
 54	else if (abs(id)%1000==443 || abs(id)%1000==553) {
 55	  onium.push_back(p);
 56	  ++nstable;
 57	}
 58	else if ( !p.children().empty() ) {
 59	  findDecayProducts(p,nstable,pip,pim,pi0,onium);
 60	}
 61	else
 62	  ++nstable;
 63      }
 64    }
 65
 66    void findLeptons(const Particle & mother,
 67		     unsigned int & nstable,
 68		     Particles& lp, Particles& lm) {
 69      for(const Particle & p : mother.children()) {
 70        int id = p.pid();
 71      	if ( id == 11 || id == 13 ) {
 72	  lm.push_back(p);
 73	  ++nstable;
 74	}
 75       	else if (id == -11 || id==-13) {
 76       	  lp.push_back(p);
 77       	  ++nstable;
 78       	}
 79	else if ( !p.children().empty() ) {
 80	  findLeptons(p,nstable,lp,lm);
 81	}
 82	else
 83	  ++nstable;
 84      }
 85    }
 86
 87
 88    /// Perform the per-event analysis
 89    void analyze(const Event& event) {
 90      // loop over unstable particles
 91      for(const Particle& ups : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==100553)) {
 92	unsigned int nstable(0);
 93	Particles pip, pim, pi0, onium;
 94	findDecayProducts(ups,nstable,pip,pim,pi0,onium);
 95	// check for onium
 96	if(onium.size() !=1 || onium[0].pid()!=553 || nstable !=3) continue;
 97	//pi+pi-
 98	if( pip.size()==1 && pim.size() ==1) {
 99	  FourMomentum q = pip[0].momentum()+pim[0].momentum();
100	  _h_charged_ex->fill(q.mass()/GeV);
101	  _h_charged_inc->fill(q.mass()/GeV);
102	  // leptons from Upsilon(1S) decay
103	  nstable = 0;
104	  Particles lp, lm;
105	  findLeptons(onium[0],nstable,lp,lm);
106	  if(nstable==2&&lp.size()==1&&lm.size()==1) {
107	    _h_lept_ctheta->fill(cos(lp[0].momentum().polarAngle()));
108	    _h_lept_phi   ->fill(lp[0].momentum().azimuthalAngle());
109	  }
110	  // pions in rest frame
111	  LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(q.betaVec());
112	  FourMomentum ppi = boost.transform(pip[0].momentum());
113	  Vector3 axis1 = q.p3().unit();
114	  Vector3 axis2 = Vector3(0.,0.,1.).cross(axis1).unit();
115	  Vector3 axis3 = axis1.cross(axis2).unit();
116	  _h_piStar_ctheta->fill(axis1.dot(ppi.p3().unit()));
117	  double phi = atan2(axis3.dot(ppi.p3().unit()),axis2.dot(ppi.p3().unit()));
118	  if(phi<0.) phi+=2.*M_PI;
119	  _h_piStar_phi->fill(phi);
120	  // pions in lab frame
121	  _h_pi_ctheta->fill(cos(q.polarAngle()));
122	  _h_pi_phi   ->fill(q.azimuthalAngle());
123	}
124	// 2pi0
125	else if (pi0.size()==2) {
126	  FourMomentum q = pi0[0].momentum()+pi0[1].momentum();
127	  _h_neutral_ex->fill(q.mass()/GeV);
128	}
129      }
130    }
131
132
133    /// Normalise histograms etc., after the run
134    void finalize() {
135      normalize(_h_charged_ex);
136      normalize(_h_charged_inc);
137      normalize(_h_neutral_ex);
138      normalize( _h_lept_ctheta);
139      normalize(_h_lept_phi);
140      normalize(_h_piStar_ctheta);
141      normalize(_h_piStar_phi);
142      normalize(_h_pi_ctheta);
143      normalize(_h_pi_phi);
144    }
145
146    /// @}
147
148
149    /// @name Histograms
150    /// @{
151    Histo1DPtr _h_charged_ex,_h_charged_inc,_h_neutral_ex;
152    Histo1DPtr _h_lept_ctheta,_h_lept_phi,_h_piStar_ctheta,_h_piStar_phi,_h_pi_ctheta,_h_pi_phi;
153    /// @}
154
155
156  };
157
158
159  RIVET_DECLARE_PLUGIN(CLEOII_1998_I467642);
160
161}