1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> e spectrum
  9  class CLEOII_2004_I647288 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_2004_I647288);
 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(),"UFS");
 23      // histograms
 24      for (unsigned int ix=0; ix<2; ++ix) {
 25        book(_h[ix],2,1,1+ix);
 26        book(_p[ix],1,1,1+ix);
 27      }
 28      book(_h[2],3,1,1);
 29      book(_c,"TMP/nb");
 30    }
 31
 32    void findDecayProducts(const Particle& parent,
 33			   vector<Particles>&  em, vector<Particles>& ep,
 34			   vector<Particles>& nue, vector<Particles>& nueBar,
 35			   bool& charm, bool secondary) {
 36      for (const Particle& p : parent.children()) {
 37      	if (PID::isCharmHadron(p.pid())) {
 38      	  charm=true;
 39          findDecayProducts(p,em,ep,nue,nueBar,charm,true);
 40      	}
 41      	else if (p.pid() == PID::EMINUS || p.pid()==PID::MUON) {
 42      	  em[secondary].push_back(p);
 43      	}
 44      	else if (p.pid() == PID::EPLUS || p.pid()==PID::ANTIMUON) {
 45      	  ep[secondary].push_back(p);
 46      	}
 47      	else if (p.pid() == PID::NU_E  || p.pid()==PID::NU_MU) {
 48      	  nue[secondary].push_back(p);
 49      	}
 50      	else if (p.pid() == PID::NU_EBAR || p.pid()==PID::NU_MUBAR) {
 51      	  nueBar[secondary].push_back(p);
 52      	}
 53      	else if (PID::isBottomHadron(p.pid())) {
 54      	  findDecayProducts(p,em,ep,nue,nueBar,charm,false);
 55      	}
 56      	else if (!PID::isHadron(p.pid())) {
 57      	  findDecayProducts(p,em,ep,nue,nueBar,charm,secondary);
 58      	}
 59      }
 60    }
 61
 62    /// Perform the per-event analysis
 63    void analyze(const Event& event) {
 64      if(_edges.empty()) _edges = _p[0]->xEdges();
 65      // find and loop over Upslion(4S)
 66      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 67      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
 68        // boost to rest frame
 69        LorentzTransform cms_boost;
 70        if (p.p3().mod() > 1*MeV) {
 71          cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 72        }
 73        // loop over decay products
 74        for (const Particle& p2 : p.children()) {
 75          if (p2.abspid()==511 || p2.abspid()==521) {
 76            _c->fill();
 77            // find decay products
 78            bool charm = false;
 79            vector<Particles> em(2),ep(2),nue(2),nueBar(2);
 80            findDecayProducts(p2,em,ep,nue,nueBar,charm,false);
 81            for (unsigned int isec=0;isec<2;++isec) {
 82              FourMomentum pl,pnu;
 83              if (em[isec].size()==1 && nueBar[isec].size()==1 && em[isec][0].pid()+1==-nueBar[isec][0].pid()) {
 84                pl  = cms_boost.transform(em[isec][0].momentum());
 85                pnu = cms_boost.transform(nueBar[isec][0].momentum());
 86              }
 87              else if (ep[isec].size()==1 && nue[isec].size()==1 && nue[isec][0].pid()==-ep[isec][0].pid()+1) {
 88                pl  = cms_boost.transform(ep[isec][0].momentum());
 89                pnu = cms_boost.transform(nue[isec][0].momentum());
 90              }
 91              else {
 92                continue;
 93              }
 94              _h[isec]->fill(pl.p3().mod());
 95              if (charm && isec==0) {
 96                _h[2]->fill(pl.p3().mod());
 97                double Emin=0.6;
 98                for(size_t i=0; i<_p[0]->numBins();++i) {
 99                  if(pl.E()>Emin) {
100                    _p[0]->fill(_edges[i],pl.E());
101                    _p[1]->fill(_edges[i],sqr(pl.E()));
102                  }
103                  Emin+=0.1;
104                }
105              }
106            }
107          }
108        }
109      }
110    }
111
112
113    /// Normalise histograms etc., after the run
114    void finalize() {
115      // 0.5 from counting e and mu
116      scale(_h, 0.5/ *_c);
117      // compute <(el^2-<el>^2>
118      BinnedEstimatePtr<string> tmp;
119      book(tmp,1,1,3);
120      for (const auto& b0 : _p[0]->bins()) {
121        const auto& b1 = _p[1]->bin(b0.index());
122        const double val = b1.yMean()-sqr(b0.yMean());
123        const double err = val*sqrt(sqr(b1.relErrW())+4.*sqr(b0.relErrW()));
124        tmp->bin(b0.index()).set(val, err);
125      }
126    }
127
128    /// @}
129
130
131    /// @name Histograms
132    /// @{
133    Histo1DPtr _h[3];
134    BinnedProfilePtr<string> _p[2];
135    CounterPtr _c;
136    vector<string> _edges;
137    /// @}
138
139
140  };
141
142
143  RIVET_DECLARE_PLUGIN(CLEOII_2004_I647288);
144
145}