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      static const double Lambda=0.65*GeV;
 60      // find and loop over Upslion(4S)
 61      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 62	if(p.children().empty() ||
 63	   (p.children().size()==1 && p.children()[1].abspid()==p.abspid()))
 64	  continue;
 65	// find decay products
 66	bool charm = false;
 67	Particles em,ep,nue,nueBar;
 68	findDecayProducts(p,em,ep,nue,nueBar,charm);
 69	if(!charm) continue;
 70	FourMomentum pl,pnu;
 71	if(em.size()==1 && nueBar.size()==1 && em[0].pid()+1==-nueBar[0].pid()) {
 72	  pl  = em[0].momentum();
 73	  pnu = nueBar[0].momentum();
 74	}
 75	else if(ep.size()==1 && nue.size()==1 && nue[0].pid()==-ep[0].pid()+1) {
 76	  pl  = ep[0].momentum();
 77	  pnu = nue[0].momentum();
 78       	}
 79	else
 80	  continue;
 81	// boost to rest frame
 82	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 83	FourMomentum pX = boost.transform(p.momentum()-pl-pnu);
 84	pl = boost.transform(pl);
 85	double modp = pl.p3().mod();
 86	double mX  = pX.mass();
 87	double nX2 = sqr(mX)-2.*Lambda*pX.t()+sqr(Lambda);
 88	double m1=mX,m2=nX2;
 89	for(unsigned int ix=0;ix<6;++ix) {
 90	  for(const auto & bin : _p_X[ix]->bins())
 91	    if(bin.xMin()<modp) _p_X[ix]->fill(bin.xMid(),m1);
 92	  m1 *=mX;
 93	  if(ix>=3) continue;
 94	  for(const auto & bin : _p_n[ix]->bins())
 95	    if(bin.xMin()<modp) _p_n[ix]->fill(bin.xMid(),m2);
 96	  m2 *=nX2;
 97	}
 98      }
 99    }
100
101
102    /// Normalise histograms etc., after the run
103    void finalize() {
104    }
105
106    /// @}
107
108
109    /// @name Histograms
110    /// @{
111    Profile1DPtr _p_X[6],_p_n[3];
112
113    /// @}
114
115
116  };
117
118
119  RIVET_DECLARE_PLUGIN(BABAR_2009_I827787);
120
121}