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 B0 -> D*0 omega
 10  class BABAR_2011_I920989 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2011_I920989);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projection
 23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511);
 24      declare(ufs, "UFS");
 25      // histos
 26      for(unsigned int ix=0;ix<4;++ix)
 27	for(unsigned int iy=0;iy<2;++iy)
 28	  book(_h[ix][iy],1,1+ix,1+iy);
 29    }
 30    
 31    void findChildren(const Particle & p, Particles & pim, Particles & pip,
 32		      Particles & pi0, unsigned int &ncount) {
 33      for( const Particle &child : p.children()) {
 34	if(child.pid()==PID::PIPLUS) {
 35	  pip.push_back(child);
 36	  ncount+=1;
 37	}
 38	else if(child.pid()==PID::PIMINUS) {
 39	  pim.push_back(child);
 40	  ncount+=1;
 41	}
 42	else if(child.pid()==PID::PI0) {
 43	  pi0.push_back(child);
 44	  ncount+=1;
 45	}
 46	else if(child.children().empty()) {
 47	  ncount+=1;
 48	}
 49    	else
 50    	  findChildren(child,pim,pip,pi0,ncount);
 51      }
 52    }
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56      for(const Particle & B0 : apply<UnstableParticles>(event,"UFS").particles()) {
 57	if(B0.children().size()!=2) continue;
 58	Particle Dstar,omega;
 59	if(B0.children()[0].abspid()==423 &&
 60	   B0.children()[1].pid()==223) {
 61	  Dstar = B0.children()[0];
 62	  omega = B0.children()[1];
 63	}
 64	else if (B0.children()[1].abspid()==423 &&
 65		 B0.children()[0].pid()==223) {
 66	  Dstar = B0.children()[1];
 67	  omega = B0.children()[0];
 68	}
 69	else
 70	  continue;
 71	// check the no of decay products
 72	if(Dstar.children().size()!=2 || omega.children().size()!=3)
 73	  continue;
 74	// find the children of the D* meson
 75	Particle D0;
 76	if(Dstar.children()[0].pid()==111 &&
 77	   Dstar.children()[1].abspid()==421)
 78	  D0 = Dstar.children()[1];
 79	else if(Dstar.children()[1].pid()==111 &&
 80		Dstar.children()[0].abspid()==421)
 81	  D0 = Dstar.children()[0];
 82	else
 83	  continue;
 84	// children of the omega
 85	unsigned int ncount=0;
 86	Particles pip,pim,pi0;
 87	findChildren(omega,pim,pip,pi0,ncount);
 88	if( ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) continue;
 89	// boost to B rest frame
 90	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(B0.momentum().betaVec());
 91	FourMomentum pDstar = boost1.transform(Dstar.momentum());
 92	FourMomentum pD0    = boost1.transform(D0   .momentum());
 93	FourMomentum pomega = boost1.transform(omega.momentum());
 94	FourMomentum pPip   = boost1.transform(pip[0].momentum());
 95	FourMomentum pPim   = boost1.transform(pim[0].momentum());
 96	// boost to D* frame
 97	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
 98	pD0 = boost2.transform(pD0);
 99	double c1 = pD0.p3().unit().dot(pDstar.p3().unit());
100	// boost to omega frame
101	LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pomega.betaVec());
102	pPip   = boost3.transform(pPip);
103	pPim   = boost3.transform(pPim);
104	Vector3 axisOmega = pPip.p3().cross(pPim.p3()).unit();
105	double c2 = pomega.p3().unit().dot(axisOmega);
106	for(unsigned int ix=0;ix<4;++ix) {
107	  _h[ix][0]->fill(c1);
108	  _h[ix][1]->fill(c2);
109	}
110      }
111    }
112
113
114    /// Normalise histograms etc., after the run
115    void finalize() {
116      for(unsigned int ix=0;ix<4;++ix)
117	for(unsigned int iy=0;iy<2;++iy)
118	  normalize(_h[ix][iy],1.,false);
119    }
120
121    /// @}
122
123
124    /// @name Histograms
125    /// @{
126    Histo1DPtr _h[4][2];
127    /// @}
128
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(BABAR_2011_I920989);
134
135}