1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Upsilon(1S) -> gamma pi+pi- K+K-
 10  class BABAR_2018_I1667191 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2018_I1667191);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==200553 or
 24						Cuts::abspid==100553);
 25      declare(ufs, "UFS");
 26      DecayedParticles UPS(ufs);
 27      UPS.addStable( 553);
 28      declare(UPS, "UPS");
 29      // histograms
 30      for(unsigned int ix=0;ix<2;++ix)
 31	book(_h_mass[ix],1+ix,1,1);
 32      book(_h_pi,3,1,1);
 33      for(unsigned int ix=0;ix<3;++ix)
 34	for(unsigned int iy=0;iy<2;++iy)
 35	  book(_h_angle[ix][iy],4+ix,1,1+iy);
 36    }
 37
 38
 39   /// Recursively walk the decay tree to find decay products of @a p
 40    void findDecayProducts(Particle mother, Particles& gamma,
 41			   Particles & pip, Particles & pim,
 42			   Particles & Kp , Particles & Km,unsigned int & nstable) {
 43      for(const Particle & p: mother.children()) {
 44	if     (p.pid()== 211) pip.push_back(p);
 45	else if(p.pid()==-211) pim.push_back(p);
 46	else if(p.pid()== 321) Kp .push_back(p);
 47	else if(p.pid()==-321) Km .push_back(p);
 48	else if(p.pid()==22) gamma.push_back(p);
 49	else if(p.children().empty())
 50	  nstable+=1;
 51	else {
 52	  findDecayProducts(p, gamma,pip,pim,Kp,Km,nstable);
 53	  
 54	}
 55      }
 56    }
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      static const map<PdgId,unsigned int> & mode   = { { 553,1},{ 211,1}, {-211,1}};
 61      DecayedParticles UPS = apply<DecayedParticles>(event, "UPS");
 62      // loop over particles
 63      for(unsigned int ix=0;ix<UPS.decaying().size();++ix) {
 64	// check pi+pi- upslion(1S) decay mode
 65      	if (!UPS.modeMatches(ix,3,mode)) continue;
 66	const Particle  & ups1 = UPS.decayProducts()[ix].at( 553)[0];
 67	const Particle  & pips = UPS.decayProducts()[ix].at( 211)[0];
 68	const Particle  & pims = UPS.decayProducts()[ix].at(-211)[0];
 69	// boost to rest frame
 70	LorentzTransform boost;
 71	if (UPS.decaying()[ix].p3().mod() > 1*MeV)
 72	  boost = LorentzTransform::mkFrameTransformFromBeta(UPS.decaying()[ix].momentum().betaVec());
 73	FourMomentum ppipi = boost.transform(pips.momentum()+pims.momentum());
 74	Vector3 axis1 = ppipi.p3().unit();
 75	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(ppipi);
 76	FourMomentum ppi = boost1.transform(boost.transform(pips.momentum()));
 77	_h_pi->fill(abs(ppi.p3().unit().dot(axis1)));
 78	unsigned int nstable=0;
 79	Particles gamma,pip,pim,Kp,Km;
 80	findDecayProducts(ups1, gamma,pip,pim,Kp,Km,nstable);
 81	if(gamma.size()!=1 || nstable!=0) continue;
 82	// gamma pi+pi-
 83	if(Kp.empty()&&Km.empty()&&pip.size()==1&&pim.size()==1) {
 84	  ppipi = pip[0].momentum()+pim[0].momentum();
 85	  double mpipi = ppipi.mass();
 86	  _h_mass[0]->fill(mpipi);
 87	  axis1 *=-1;
 88	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(boost.transform(ups1.momentum()).betaVec());
 89	  FourMomentum pgamma = boost2.transform(boost.transform(gamma[0].momentum()));
 90	  Vector3 axis2 = pgamma.p3().unit();
 91	  double cGamma = axis2.dot(axis1);
 92	  ppipi = boost2.transform(boost1.transform(ppipi));
 93	  LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(ppipi.betaVec());
 94	  Vector3 axis3 = boost3.transform(boost2.transform(boost1.transform(pip[0].momentum()))).p3().unit();
 95	  double cH = axis3.dot(axis2);
 96	  int iloc=-1;
 97	  if(mpipi>0.6&&mpipi<1.)          iloc=0;
 98	  else if(mpipi>1.092&&mpipi<1.46) iloc=1;
 99	  if(iloc>=0) {
100	    _h_angle[iloc][0]->fill(cGamma);
101	    _h_angle[iloc][1]->fill(cH);
102	  }
103	}
104	// gamma K+K-
105	else if (pip.empty()&&pim.empty()&&Kp.size()==1&&Km.size()==1) {
106	  FourMomentum pKK = Kp[0].momentum()+Km[0].momentum(); 
107	  double mKK = pKK.mass();
108	  _h_mass[1]->fill(mKK);
109	  axis1 *=-1;
110	  LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(boost.transform(ups1.momentum()).betaVec());
111	  FourMomentum pgamma = boost2.transform(boost.transform(gamma[0].momentum()));
112	  Vector3 axis2 = pgamma.p3().unit();
113	  double cGamma = axis2.dot(axis1);
114	  pKK = boost2.transform(boost1.transform(pKK));
115	  LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pKK.betaVec());
116	  Vector3 axis3 = boost3.transform(boost2.transform(boost1.transform(Kp[0].momentum()))).p3().unit();
117	  double cH = axis3.dot(axis2);
118	  if(mKK>1.424 && mKK<1.62) {
119	    _h_angle[2][0]->fill(cGamma);
120	    _h_angle[2][1]->fill(cH);
121	  }
122	}
123      }
124    }
125
126
127    /// Normalise histograms etc., after the run
128    void finalize() {
129      for(unsigned int ix=0;ix<2;++ix)
130	normalize(_h_mass[ix],1.,false);
131      normalize(_h_pi,1.,false);
132      for(unsigned int ix=0;ix<3;++ix)
133	for(unsigned int iy=0;iy<2;++iy)
134	  normalize(_h_angle[ix][iy],1.,false);
135    }
136
137    /// @}
138
139
140    /// @name Histograms
141    /// @{
142    Histo1DPtr _h_mass[2];
143    Histo1DPtr _h_pi;
144    Histo1DPtr _h_angle[3][2];
145    /// @}
146
147  };
148
149
150  RIVET_DECLARE_PLUGIN(BABAR_2018_I1667191);
151
152}