1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief chi_c1 -> gamma +V
  9  class BESIII_2011_I894356 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2011_I894356);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(Cuts::pid==20443), "UFS");
 23      // histograms
 24      for(unsigned int ix=0;ix<3;++ix)
 25	book(_h[ix],1,1,1+ix);
 26    }
 27    
 28    void findChildren(const Particle & p, Particles & pim, Particles & pip,
 29		      Particles & pi0, unsigned int &ncount) {
 30      for( const Particle &child : p.children()) {
 31	if(child.pid()==PID::PIPLUS) {
 32	  pip.push_back(child);
 33	  ncount+=1;
 34	}
 35	else if(child.pid()==PID::PIMINUS) {
 36	  pim.push_back(child);
 37	  ncount+=1;
 38	}
 39	else if(child.pid()==PID::PI0) {
 40	  pi0.push_back(child);
 41	  ncount+=1;
 42	}
 43	else if(child.children().empty()) {
 44	  ncount+=1;
 45	}
 46    	else
 47    	  findChildren(child,pim,pip,pi0,ncount);
 48      }
 49    }
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53      for (const Particle& p :  apply<UnstableParticles>(event, "UFS").particles()) {
 54	if(p.children().size()!=2) continue;
 55	Particle gamma,meson;
 56	// get the photon and the meson
 57	if(p.children()[0].pid()==PID::PHOTON) {
 58	  gamma = p.children()[0];
 59	  meson = p.children()[1];
 60	}
 61	else if(p.children()[1].pid()==PID::PHOTON) {
 62	  gamma = p.children()[1];
 63	  meson = p.children()[0];
 64	}
 65	else
 66	  continue;
 67	// check the meson
 68	int imode=-1;
 69	if(meson.pid()==333)      imode = 0;
 70	else if(meson.pid()==113) imode = 1;
 71	else if(meson.pid()==223) imode = 2;
 72	else continue;
 73	Particle d1,d2,d3;
 74	if(imode==0) {
 75	  if(meson.children().size()!=2) continue;
 76	  if(meson.children()[0].pid()==-meson.children()[1].pid() &&
 77	     meson.children()[0].abspid()==321) {
 78	    d1 = meson.children()[0];
 79	    d2 = meson.children()[0];
 80	  }
 81	  else
 82	    continue;
 83	}
 84	else if(imode==1) {
 85	  if(meson.children().size()!=2) continue;
 86	  if(meson.children()[0].pid()==-meson.children()[1].pid() &&
 87	     meson.children()[0].abspid()==211) {
 88	    d1 = meson.children()[0];
 89	    d2 = meson.children()[0];
 90	  }
 91	  else
 92	    continue;
 93	}
 94	else if(imode==2) {
 95	  unsigned int ncount=0;
 96	  Particles pip,pim,pi0;
 97	  findChildren(meson,pim,pip,pi0,ncount);
 98	  if(ncount==3 && pim.size()==1 && pip.size()==1 && pi0.size()==1) {
 99	    d1=pip[0];
100	    d2=pim[0];
101	    d3=pi0[0];
102	  }
103	  else continue;
104	}
105	if(d1.pid()<0) swap(d1,d2);
106	// boost in chi_c1 frame
107	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
108	FourMomentum pMeson = boost1.transform(meson.momentum());
109	FourMomentum p1=boost1.transform(d1.momentum());
110	FourMomentum p2=boost1.transform(d2.momentum());
111	Vector3 axis1 = pMeson.p3().unit();
112	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pMeson.betaVec());
113	p1=boost2.transform(p1);
114	p2=boost2.transform(p2);
115	// axis in meson rest frame
116	Vector3 axis2;
117	if(imode<2) {
118	  axis2=p1.p3().unit();
119	}
120	else {
121	  FourMomentum p3 = boost1.transform(d3.momentum());
122	  p3 = boost2.transform(p3);
123	  axis2=p1.p3().cross(p2.p3()).unit();
124	}
125	_h[imode]->fill(axis1.dot(axis2));
126      }
127    }
128
129
130    /// Normalise histograms etc., after the run
131    void finalize() {
132      for(unsigned int ix=0;ix<3;++ix)
133	normalize(_h[ix],1.,false);
134    }
135
136    /// @}
137
138
139    /// @name Histograms
140    /// @{
141    Histo1DPtr _h[3];
142    /// @}
143
144
145  };
146
147
148  RIVET_DECLARE_PLUGIN(BESIII_2011_I894356);
149
150}