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      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 59      for(const Particle& p : ufs.particles()) {
 60	if(p.children().empty()) continue;
 61	if(p.children()[0].abspid()==p.abspid()) continue;
 62	_nB->fill();
 63	bool charm = false;
 64	Particles em,ep,nue,nueBar;
 65	findDecayProducts(p,em,ep,nue,nueBar,charm);
 66	if(charm) continue;
 67	FourMomentum pl,pnu;
 68	if(em.size()==1 && nueBar.size()==1) {
 69	  pl  = em[0].momentum();
 70	}
 71	else if(ep.size()==1 && nue.size()==1) {
 72	  pl  = ep[0].momentum();
 73	}
 74	else
 75	  continue;
 76	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 77	pl  = boost.transform(pl );
 78	double pp = pl.p3().mod();
 79	for(const auto & bin : _h->bins())
 80	  if(bin.xMin()<pp) _h->fill(bin.xMid());
 81      }
 82    }
 83
 84
 85    /// Normalise histograms etc., after the run
 86    void finalize() {
 87      // 1e3 due br normalization and 0.1 to remove bin width
 88      scale(_h, 0.1*1e3/ *_nB);
 89    }
 90
 91    /// @}
 92
 93
 94    /// @name Histograms
 95    /// @{
 96    Histo1DPtr _h;
 97    CounterPtr _nB;
 98    /// @}
 99
100
101  };
102
103
104  RIVET_DECLARE_PLUGIN(BABAR_2012_I1079912);
105
106}