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) _b_p->fill();
 93	  else    _r_p->fill("10.58"s);
 94	}
 95	else {
 96	  if(ups) _b_0->fill();
 97	  else    _r_0->fill("10.58"s);
 98	  idXic=4132;
 99	}
100	// boost to Xi'_c rest frame
101	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
102	FourMomentum pXic = boost1.transform(Xi_c.momentum());
103	FourMomentum pXi  = boost1.transform(Xi[0].momentum());
104	// then to Xi_c rest frame
105	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pXic.betaVec());
106	Vector3 axis = pXic.p3().unit();
107	FourMomentum pp = boost2.transform(pXi);
108	// calculate angle
109	double cTheta = pp.p3().unit().dot(axis);
110	_h_ctheta->fill(cTheta);
111      }
112    }
113
114
115    /// Normalise histograms etc., after the run
116    void finalize() {
117      normalize(_h_p_0   );
118      normalize(_h_p_p   );
119      normalize(_h_ctheta);
120      if(_ups->effNumEntries()!=0) {
121	scale(_b_p,0.5/ *_ups);
122	scale(_b_0,0.5/ *_ups);
123      }
124      scale(_r_p,crossSection()/sumOfWeights()/femtobarn);
125      scale(_r_0,crossSection()/sumOfWeights()/femtobarn);
126    }
127
128    /// @}
129
130
131    /// @name Histograms
132    /// @{
133    Histo1DPtr _h_p_0,_h_p_p,_h_ctheta;
134    CounterPtr _b_p,_b_0;
135    BinnedHistoPtr<string> _r_p,_r_0;
136    CounterPtr _ups;
137    /// @}
138
139
140  };
141
142
143  RIVET_DECLARE_PLUGIN(BABAR_2007_I722622);
144
145}