1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Electron spectrum in B decays
  9  class BABAR_2017_I1498564 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2017_I1498564);
 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::pid==300553),"UFS");
 23      // Book histograms
 24      // specify custom binning
 25      book(_h_light, 1,1,1);
 26      book(_h_all  , 1,1,2);
 27      book(_nB,"/TMP/nB");
 28    }
 29
 30
 31    void findDecayProducts(Particle parent, Particles & em, Particles & ep,
 32			   Particles & nue, Particles & nueBar, bool & charm) {
 33      for(const Particle & p : parent.children()) {
 34	if(PID::isCharmHadron(p.pid())) {
 35	  charm=true;
 36	}
 37	else if(p.pid() == PID::EMINUS) {
 38	  em.push_back(p);
 39	}
 40	else if(p.pid() == PID::EPLUS) {
 41	  ep.push_back(p);
 42	}
 43	else if(p.pid() == PID::NU_E) {
 44	  nue.push_back(p);
 45	}
 46	else if(p.pid() == PID::NU_EBAR) {
 47	  nueBar.push_back(p);
 48	}
 49	else if(PID::isBottomHadron(p.pid())) {
 50	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 51	}
 52	else if(!PID::isHadron(p.pid())) {
 53	  findDecayProducts(p,em,ep,nue,nueBar,charm);
 54	}
 55      }
 56    }
 57
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& event) {
 61      // find and loop over Upslion(4S)
 62      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 63      	// boost to rest frame
 64      	LorentzTransform cms_boost;
 65      	if (p.p3().mod() > 1*MeV)
 66      	  cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 67      	// loop over decay products
 68      	for(const Particle & p2 : p.children()) {
 69      	  if(p2.abspid()==511 || p2.abspid()==521) {
 70	    _nB->fill();
 71	    bool charm = false;
 72      	    Particles em,ep,nue,nueBar;
 73	    findDecayProducts(p2,em,ep,nue,nueBar,charm);
 74	    if(em.size()==1 && nueBar.size()==1) {
 75	      double pmod = cms_boost.transform(em[0].momentum()).p3().mod();
 76	      _h_all->fill(pmod);
 77	      if(!charm) _h_light->fill(pmod);
 78	    }
 79	    else if(ep.size()==1 && nue.size()==1) {
 80	      double pmod = cms_boost.transform(ep[0].momentum()).p3().mod();
 81	      _h_all->fill(pmod);
 82	      if(!charm) _h_light->fill(pmod);
 83	    }
 84      	  }
 85      	}
 86      }
 87    }
 88
 89
 90    /// Normalise histograms etc., after the run
 91    void finalize() {
 92
 93      // normalize to number of B decays
 94      scale(_h_all, 1./ *_nB);
 95      scale(_h_light, 1./ *_nB);
 96    }
 97
 98    //@}
 99
100
101    /// @name Histograms
102    //@{
103    Histo1DPtr _h_light,_h_all;
104    CounterPtr _nB;
105    //@}
106
107
108  };
109
110
111  RIVET_DECLARE_PLUGIN(BABAR_2017_I1498564);
112
113}