1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Bbar0 -> D*0 2pi+2pi-
 10  class CLEOII_2002_I557084 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_2002_I557084);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511);
 23      declare(ufs, "UFS");
 24      DecayedParticles B0(ufs);
 25      B0.addStable(PID::PI0);
 26      B0.addStable( 413);
 27      B0.addStable(-413);
 28      B0.addStable( 423);
 29      B0.addStable(-423);
 30      B0.addStable( 411);
 31      B0.addStable(-411);
 32      B0.addStable( 421);
 33      B0.addStable(-421);
 34      B0.addStable( 431);
 35      B0.addStable(-431);
 36      declare(B0, "B0");
 37      // histos
 38      for (unsigned int ix=0; ix<2; ++ix) {
 39       	book(_h_Dpi[ix],1,1,1+ix);
 40      	book(_h_4pi[ix],3,1,1+ix);
 41      }
 42      book(_h_3pi,2,1,1);
 43      book(_c,"TMP/nB");
 44    }
 45
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      // loop over particles
 50      DecayedParticles B0 = apply<DecayedParticles>(event, "B0");
 51      int imode = -1;
 52      for (unsigned int ix=0; ix<B0.decaying().size(); ++ix) {
 53        int sign = B0.decaying()[ix].pid()/B0.decaying()[ix].abspid();
 54        _c->fill();
 55        if ( (sign== 1 && B0.modeMatches(ix,5,mode1) ) ||
 56             (sign==-1 && B0.modeMatches(ix,5,mode1CC) ) ) {
 57          imode=0;
 58        }
 59        else if ( (sign== 1 && B0.modeMatches(ix,5,mode2) ) ||
 60                  (sign==-1 && B0.modeMatches(ix,5,mode2CC) ) ) {
 61          imode=1;
 62        }
 63        else {
 64          continue;
 65        }
 66        const Particles& pip = B0.decayProducts()[ix].at( sign*211);
 67        const Particles& pim = B0.decayProducts()[ix].at(-sign*211);
 68        FourMomentum pTotal;
 69        for (const Particle& p : pip) pTotal+=p.mom();
 70        for (const Particle& p : pim) pTotal+=p.mom();
 71        if (imode==0) {
 72          const Particle& pi0 = B0.decayProducts()[ix].at(111)[0];
 73          pTotal+= pi0.mom();
 74          _h_4pi[1]->fill(pTotal.mass2());
 75          continue;
 76        }
 77        _h_4pi[0]->fill(pTotal.mass2());
 78        // D* pi mases
 79        const Particle& Dstar = B0.decayProducts()[ix].at(-sign*423)[0];
 80        for (const Particle& p : pim) {
 81          FourMomentum pDpi = Dstar.mom()+p.mom();
 82          const double mDpi = pDpi.mass();
 83          _h_Dpi[0]->fill(mDpi);
 84          if (mDpi>2.3 and mDpi<2.6) {
 85            FourMomentum p3pi = B0.decaying()[ix].mom()-pDpi;
 86            _h_3pi->fill(p3pi.mass());
 87          }
 88        }
 89        for (const Particle& p : pip) {
 90          _h_Dpi[1]->fill((Dstar.mom()+p.mom()).mass());
 91        }
 92      }
 93    }
 94
 95
 96    /// Normalise histograms etc., after the run
 97    void finalize() {
 98      // first hist is differential BR/ BR for leptons
 99      const double br = .1033;
100      normalize(_h_Dpi, 1.0, false);
101      scale(_h_4pi, 1.0/br/ *_c);
102      normalize(_h_3pi,1.,false);
103    }
104
105    /// @}
106
107
108    /// @name Histograms
109    /// @{
110    Histo1DPtr _h_Dpi[2], _h_3pi, _h_4pi[2];
111    CounterPtr _c;
112    const map<PdgId,unsigned int> mode1   = { { -413,1}, { 211,2}, {-211,1}, {111,1}};
113    const map<PdgId,unsigned int> mode1CC = { {  413,1}, {-211,2}, { 211,1}, {111,1}};
114    const map<PdgId,unsigned int> mode2   = { { -423,1}, { 211,2}, {-211,2}};
115    const map<PdgId,unsigned int> mode2CC = { {  423,1}, { 211,2}, {-211,2}};
116    /// @}
117
118
119  };
120
121
122  RIVET_DECLARE_PLUGIN(CLEOII_2002_I557084);
123
124}