1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> X_u ell nu
  9  class BABAR_2012_I1079912 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2012_I1079912);
 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 or
 23				Cuts::abspid==521),"UFS");
 24      // histograms
 25      book(_h,1,1,1);
 26      book(_nB,"/TMP/nB");
 27    }
 28
 29    void findDecayProducts(Particle parent, Particles & em, Particles & ep,
 30			   Particles & nue, Particles & nueBar, bool & charm) {
 31      for(const Particle & p : parent.children()) {
 32	if(PID::isCharmHadron(p.pid())) {
 33	  charm=true;
 34	}
 35	else if(p.pid() == PID::EMINUS) {
 36	  em.push_back(p);
 37	}
 38	else if(p.pid() == PID::EPLUS) {
 39	  ep.push_back(p);
 40	}
 41	else if(p.pid() == PID::NU_E || p.pid()==PID::NU_MU) {
 42	  nue.push_back(p);
 43	}
 44	else if(p.pid() == PID::NU_EBAR || p.pid()==PID::NU_MUBAR) {
 45	  nueBar.push_back(p);
 46	}
 47	else if(PID::isBottomHadron(p.pid())) {
 48	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 49	}
 50	else if(!PID::isHadron(p.pid())) {
 51	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 52	}
 53      }
 54    }
 55
 56    /// Perform the per-event analysis
 57    void analyze(const Event& event) {
 58      if(_edges.empty()) _edges=_h->xEdges();
 59      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 60      for(const Particle& p : ufs.particles()) {
 61	if(p.children().empty()) continue;
 62	if(p.children()[0].abspid()==p.abspid()) continue;
 63	_nB->fill();
 64	bool charm = false;
 65	Particles em,ep,nue,nueBar;
 66	findDecayProducts(p,em,ep,nue,nueBar,charm);
 67	if(charm) continue;
 68	FourMomentum pl,pnu;
 69	if(em.size()==1 && nueBar.size()==1) {
 70	  pl  = em[0].momentum();
 71	}
 72	else if(ep.size()==1 && nue.size()==1) {
 73	  pl  = ep[0].momentum();
 74	}
 75	else
 76	  continue;
 77	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 78	pl  = boost.transform(pl );
 79	double pp = pl.p3().mod();
 80        double Emin=1.0;
 81        for(unsigned int ix=0;ix<_edges.size();++ix) {
 82	  if(pp>Emin) _h->fill(_edges[ix]);
 83          Emin+=0.1;
 84        }
 85      }
 86    }
 87
 88
 89    /// Normalise histograms etc., after the run
 90    void finalize() {
 91      // 1e3 due br normalization and 0.1 to remove bin width
 92      scale(_h, 1e3/ *_nB);
 93    }
 94
 95    /// @}
 96
 97
 98    /// @name Histograms
 99    /// @{
100    BinnedHistoPtr<string> _h;
101    CounterPtr _nB;
102    vector<string> _edges;
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(BABAR_2012_I1079912);
110
111}