1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief e- spectrum in semileptonic B decays
  9  class BELLE_2005_I665011 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2005_I665011);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projection
 22      declare(UnstableParticles(Cuts::abspid==511 ||
 23				Cuts::abspid==521),"UFS");
 24      // histos
 25      for (unsigned int ix=0;ix<2;++ix) {
 26        book(_h[ix],1,1,1+ix);
 27        book(_c[ix],"TMP/c_"+toString(ix+1));
 28      }
 29    }
 30
 31    void findDecayProducts(const 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    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      // find and loop over B mesos
 61      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 62        unsigned int iB = (p.abspid()-511)/10;
 63        _c[iB]->fill();
 64        // find decay products
 65        bool charm = false;
 66        Particles em,ep,nue,nueBar;
 67        findDecayProducts(p,em,ep,nue,nueBar,charm);
 68        if (!charm) continue;
 69        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 70        double pmod(0.);
 71        if (em.size()==1 && nueBar.size()==1) {
 72          pmod = boost.transform(em[0].momentum()).p3().mod();
 73        }
 74        else if (ep.size()==1 && nue.size()==1) {
 75          pmod = boost.transform(ep[0].momentum()).p3().mod();
 76        }
 77        else {
 78          continue;
 79        }
 80        _h[iB]->fill(pmod);
 81      }
 82    }
 83
 84
 85    /// Normalise histograms etc., after the run
 86    void finalize() {
 87      // normalize hists
 88      for (unsigned int ix=0; ix<2; ++ix) {
 89	      scale(_h[ix],1./ *_c[ix]);
 90      }
 91      // ratio
 92      Estimate1DPtr ratio;
 93      book(ratio,1,1,3);
 94      divide(_h[1], _h[0], ratio);
 95    }
 96
 97    /// @}
 98
 99
100    /// @name Histograms
101    /// @{
102    Histo1DPtr _h[2];
103    CounterPtr _c[2];
104    /// @}
105
106
107  };
108
109
110  RIVET_DECLARE_PLUGIN(BELLE_2005_I665011);
111
112}