1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> Xc ell - nu_ell
  9  class BABAR_2009_I827787 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2009_I827787);
 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::abspid==511 ||
 23				Cuts::abspid==521),"UFS");
 24      // histograms
 25      for(unsigned int ix=0;ix<6;++ix) {
 26	book(_p_X[ix],1,1,1+ix);
 27	if(ix<3) book(_p_n[ix],2,1,1+ix);
 28      }
 29    }
 30    void findDecayProducts(Particle parent, Particles & em, Particles & ep,
 31			   Particles & nue, Particles & nueBar, bool & charm) {
 32      for(const Particle & p : parent.children()) {
 33	if(PID::isCharmHadron(p.pid())) {
 34	  charm=true;
 35	}
 36	else if(p.pid() == PID::EMINUS || p.pid()==PID::MUON) {
 37	  em.push_back(p);
 38	}
 39	else if(p.pid() == PID::EPLUS || p.pid()==PID::ANTIMUON) {
 40	  ep.push_back(p);
 41	}
 42	else if(p.pid() == PID::NU_E  || p.pid()==PID::NU_MU) {
 43	  nue.push_back(p);
 44	}
 45	else if(p.pid() == PID::NU_EBAR || p.pid()==PID::NU_MUBAR) {
 46	  nueBar.push_back(p);
 47	}
 48	else if(PID::isBottomHadron(p.pid())) {
 49	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 50	}
 51	else if(!PID::isHadron(p.pid())) {
 52	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 53	}
 54      }
 55    }
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59      if(_edges.empty()) _edges = _p_X[0]->xEdges();
 60      static const double Lambda=0.65*GeV;
 61      // find and loop over Upslion(4S)
 62      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 63	if(p.children().empty() ||
 64	   (p.children().size()==1 && p.children()[1].abspid()==p.abspid()))
 65	  continue;
 66	// find decay products
 67	bool charm = false;
 68	Particles em,ep,nue,nueBar;
 69	findDecayProducts(p,em,ep,nue,nueBar,charm);
 70	if(!charm) continue;
 71	FourMomentum pl,pnu;
 72	if(em.size()==1 && nueBar.size()==1 && em[0].pid()+1==-nueBar[0].pid()) {
 73	  pl  = em[0].momentum();
 74	  pnu = nueBar[0].momentum();
 75	}
 76	else if(ep.size()==1 && nue.size()==1 && nue[0].pid()==-ep[0].pid()+1) {
 77	  pl  = ep[0].momentum();
 78	  pnu = nue[0].momentum();
 79       	}
 80	else
 81	  continue;
 82	// boost to rest frame
 83	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 84	FourMomentum pX = boost.transform(p.momentum()-pl-pnu);
 85	pl = boost.transform(pl);
 86	double modp = pl.p3().mod();
 87	double mX  = pX.mass();
 88	double nX2 = sqr(mX)-2.*Lambda*pX.t()+sqr(Lambda);
 89	double m1=mX,m2=nX2;
 90	for(unsigned int ix=0;ix<6;++ix) {
 91          double Emin = 0.8;
 92          for(unsigned int iy=0;iy<_edges.size();++iy) {
 93            if(modp>Emin) _p_X[ix]->fill(_edges[iy],m1);
 94            Emin+=0.1;
 95          }
 96          m1 *=mX;
 97          if(ix>=3) continue;
 98          Emin = 0.8;
 99          for(unsigned int iy=0;iy<_edges.size();++iy) {
100	    if(modp>Emin) _p_n[ix]->fill(_edges[iy],m2);
101            Emin+=0.1;
102          }
103	  m2 *=nX2;
104	}
105      }
106    }
107
108
109    /// Normalise histograms etc., after the run
110    void finalize() {
111    }
112
113    /// @}
114
115
116    /// @name Histograms
117    /// @{
118    BinnedProfilePtr<string> _p_X[6],_p_n[3];
119    vector<string> _edges;
120    /// @}
121
122
123  };
124
125
126  RIVET_DECLARE_PLUGIN(BABAR_2009_I827787);
127
128}