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 BELLE_2007_I732595 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2007_I732595);
 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 || Cuts::abspid==521),"UFS");
 23      // histograms
 24      for (unsigned int ix=0; ix<2; ++ix) {
 25        book(_p[ix], 1, 1, 1+2*ix);
 26      }
 27    }
 28
 29    void findDecayProducts(const 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 || p.pid()==PID::MUON) {
 36          em.push_back(p);
 37        }
 38        else if (p.pid() == PID::EPLUS || p.pid()==PID::ANTIMUON) {
 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()) {
 59        _edges = _p[0]->xEdges();
 60        for(const string & en : _edges)
 61          _eCut.push_back(std::stod(en)*GeV);
 62      }
 63      // find and loop over B mesons
 64      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 65        if (p.children().empty() ||
 66           (p.children().size()==1 && p.children()[1].abspid()==p.abspid())) continue;
 67        // find decay products
 68        bool charm = false;
 69        Particles em,ep,nue,nueBar;
 70        findDecayProducts(p,em,ep,nue,nueBar,charm);
 71        if (!charm) continue;
 72        FourMomentum pl,pnu;
 73        if (em.size()==1 && nueBar.size()==1 && em[0].pid()+1==-nueBar[0].pid()) {
 74          pl  = em[0].mom();
 75          pnu = nueBar[0].mom();
 76        }
 77        else if (ep.size()==1 && nue.size()==1 && nue[0].pid()==-ep[0].pid()+1) {
 78          pl  = ep[0].mom();
 79          pnu = nue[0].mom();
 80        }
 81        else {
 82          continue;
 83        }
 84        // boost to rest frame
 85        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.mom().betaVec());
 86        double mX2 = (p.mom()-pl-pnu).mass2();
 87        pl = boost.transform(pl);
 88        double Estar = pl.E();
 89        for (unsigned int ix=0;ix<2;++ix) {
 90          for(unsigned int ix=0;ix<_eCut.size();++ix) {
 91            if(Estar>_eCut[ix]) {
 92              _p[0]->fill(_edges[ix],mX2);
 93              _p[1]->fill(_edges[ix],sqr(mX2));
 94            }
 95          }
 96        }
 97      }
 98    }
 99
100
101    /// Normalise histograms etc., after the run
102    void finalize() {
103      // compute <(mx2-<mx2>)^2> = <mx4>-<mx2>^2
104      BinnedEstimatePtr<string> tmp;
105      book(tmp,1,1,2);
106      for (unsigned int ix=0; ix<_p[0]->numBins(); ++ix) {
107        const double value = _p[1]->bin(ix+1).mean(2)-sqr(_p[0]->bin(ix+1).mean(2));
108        const double error = value*sqrt(sqr(_p[1]->bin(ix+1).relStdErr(2))+4.*sqr(_p[0]->bin(ix+1).relStdErr(2)));
109        tmp->bin(ix+1).set(value,error);
110      }
111    }
112
113    /// @}
114
115
116    /// @name Histograms
117    /// @{
118    BinnedProfilePtr<string> _p[2];
119    vector<string> _edges;
120    vector<double> _eCut;
121    /// @}
122
123
124  };
125
126
127  RIVET_DECLARE_PLUGIN(BELLE_2007_I732595);
128
129}