1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Xi_c' spectrum
  9  class BABAR_2007_I722622 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2007_I722622);
 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      book(_h_p_0,3,1,2);
 25      book(_h_p_p,3,1,1);
 26      book(_h_ctheta,4,1,1);
 27      book(_b_p,1,1,1);
 28      book(_b_0,1,1,2);
 29      book(_r_p,2,1,1);
 30      book(_r_0,2,1,2);
 31      book(_ups,"/TMP/ups");
 32    }
 33
 34    void findChildren(Particle parent, unsigned int & nStable,
 35		      Particles &Xi, unsigned int &nPi) {
 36      for(const Particle & p : parent.children()) {
 37	if(p.abspid()==PID::PIPLUS) {
 38	  ++nPi;
 39	  ++nStable;
 40	}
 41	else if(p.abspid()==PID::XIMINUS) {
 42	  Xi.push_back(p);
 43	  ++nStable;
 44	}
 45	else if(!p.children().empty()) {
 46	  findChildren(p,nStable,Xi,nPi);
 47	}
 48	else {
 49	  ++nStable;
 50	}
 51      }
 52    }
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56      static const int id0 = 4312, idp = 4322;
 57      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 58      bool ups = !ufs.particles(Cuts::pid==300553).empty();
 59      if(ups) _ups->fill();
 60      for (const Particle& p : ufs.particles(Cuts::abspid==idp or Cuts::abspid==id0)) {
 61	int idXic=4232;
 62	if(p.abspid()==idp) {
 63	  _h_p_p->fill(p.momentum().p3().mod());
 64	}
 65	else {
 66	  _h_p_0->fill(p.momentum().p3().mod());
 67	  idXic=4132;
 68	}
 69	// first find the Xi_c in the decay
 70	if(p.children().size()!=2) continue;
 71	int sign=p.pid()/p.abspid();
 72	Particle Xi_c;
 73	if(p.children()[0].pid()==sign*idXic &&
 74	   p.children()[1].pid()==22) {
 75	  Xi_c = p.children()[0];
 76	}
 77	else if(p.children()[1].pid()==sign*idXic &&
 78		p.children()[0].pid()==22) {
 79	  Xi_c = p.children()[1];
 80	}
 81	else
 82	  continue;
 83	// and the children of the Xi_c
 84	Particles Xi;
 85	unsigned int nStable(0),nPi(0);
 86	findChildren(Xi_c,nStable,Xi,nPi);
 87	if(Xi.size()!=1) continue;
 88	if( !(p.abspid()==idp && nPi==2 && nStable==3) &&
 89	    !(p.abspid()==id0 && nPi==1 && nStable==2) )
 90	  continue;
 91	if(p.abspid()==idp) {
 92	  if(ups)
 93	    _b_p->fill(0.5);
 94	  else
 95	    _r_p->fill(0.5);
 96	}
 97	else {
 98	  if(ups)
 99	    _b_0->fill(0.5);
100	  else
101	    _r_0->fill(0.5);
102	  idXic=4132;
103	}
104	// boost to Xi'_c rest frame
105	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
106	FourMomentum pXic = boost1.transform(Xi_c.momentum());
107	FourMomentum pXi  = boost1.transform(Xi[0].momentum());
108	// then to Xi_c rest frame
109	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pXic.betaVec());
110	Vector3 axis = pXic.p3().unit();
111	FourMomentum pp = boost2.transform(pXi);
112	// calculate angle
113	double cTheta = pp.p3().unit().dot(axis);
114	_h_ctheta->fill(cTheta);
115      }
116    }
117
118
119    /// Normalise histograms etc., after the run
120    void finalize() {
121      normalize(_h_p_0   );
122      normalize(_h_p_p   );
123      normalize(_h_ctheta);
124      if(_ups->effNumEntries()!=0) {
125	scale(_b_p,0.5/ *_ups);
126	scale(_b_0,0.5/ *_ups);
127      }
128      scale(_r_p,crossSection()/sumOfWeights()/femtobarn);
129      scale(_r_0,crossSection()/sumOfWeights()/femtobarn);
130    }
131
132    /// @}
133
134
135    /// @name Histograms
136    /// @{
137    Histo1DPtr _h_p_0,_h_p_p,_h_ctheta;
138    Histo1DPtr _b_p,_b_0;
139    Histo1DPtr _r_p,_r_0;
140    CounterPtr _ups;
141    /// @}
142
143
144  };
145
146
147  RIVET_DECLARE_PLUGIN(BABAR_2007_I722622);
148
149}