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 B -> D* rho
 10  class CLEOII_2003_I611745 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_2003_I611745);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511 || Cuts::abspid==521);
 24      declare(ufs, "UFS");
 25      DecayedParticles BB(ufs);
 26      BB.addStable( 423);
 27      BB.addStable(-423);
 28      BB.addStable( 413);
 29      BB.addStable(-413);
 30      BB.addStable( 213);
 31      BB.addStable(-213);
 32      BB.addStable( 113);
 33      declare(BB, "BB");
 34      // histos
 35      for (unsigned int ix=0; ix<2; ++ix) {
 36        for (unsigned int iy=0; iy<3; ++iy) {
 37          book(_h[ix][iy], 2, 1+iy, 1+ix);
 38          if(iy<2) book(_p[ix][iy], "TMP/TMP_"+toString(ix)+"_"+toString(iy));
 39        }
 40      }
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      DecayedParticles BB = apply<DecayedParticles>(event, "BB");
 47      // loop over particles
 48      for (unsigned int ix=0; ix<BB.decaying().size(); ++ix) {
 49        int sign = BB.decaying()[ix].pid()/BB.decaying()[ix].abspid();
 50        unsigned int imode=0;
 51       	if ( ( BB.decaying()[ix].pid() == 521 && BB.modeMatches(ix,2,mode1  ) ) ||
 52             ( BB.decaying()[ix].pid() ==-521 && BB.modeMatches(ix,2,mode1CC) ) ) imode=0;
 53       	else if ( ( BB.decaying()[ix].pid() == 511 && BB.modeMatches(ix,2,mode2  ) ) ||
 54                  ( BB.decaying()[ix].pid() ==-511 && BB.modeMatches(ix,2,mode2CC) ) ) imode=1;
 55        else continue;
 56        // find D* decay products
 57        const Particle & Dstar = BB.decayProducts()[ix].at(-sign*(423-10*imode))[0];
 58        if (Dstar.children().size()!=2) continue;
 59        Particle D0;
 60        int ipi = imode==0 ?  111 : -sign*211;
 61        if (Dstar.children()[0].pid()   ==  ipi &&
 62            Dstar.children()[1].abspid()==-sign*421) {
 63          D0 = Dstar.children()[1];
 64        }
 65        else if (Dstar.children()[1].pid() == ipi &&
 66                 Dstar.children()[0].abspid()==-sign*421) {
 67          D0 = Dstar.children()[0];
 68        }
 69        else {
 70          continue;
 71        }
 72        // rho decay products
 73        const Particle& rho = BB.decayProducts()[ix].at( sign*213)[0];
 74        if (rho.children().size()!=2) continue;
 75        Particle pip;
 76        if (rho.children()[1].pid()==111 &&
 77            rho.children()[0].pid()==sign*211)
 78          pip = rho.children()[0];
 79        else if (rho.children()[0].pid()==111 &&
 80                 rho.children()[1].pid()==sign*211)
 81          pip = rho.children()[1];
 82        else {
 83          continue;
 84        }
 85        // boost to B rest frame
 86        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(BB.decaying()[ix]. mom().betaVec());
 87        FourMomentum prho   = boost.transform(rho.mom());
 88        FourMomentum pDstar = boost.transform(Dstar.mom());
 89        FourMomentum pPi = boost.transform(pip.mom());
 90        FourMomentum pD  = boost.transform(D0   .mom());
 91        const LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
 92        pD = boost2.transform(pD);
 93        double cD = pD.p3().unit().dot(pDstar.p3().unit());
 94        Vector3 trans1 = pD.p3().unit() - cD*pDstar.p3().unit();
 95        _h[imode][0]->fill(cD);
 96        const LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(prho.betaVec());
 97        pPi = boost3.transform(pPi);
 98        double cPi = pPi.p3().unit().dot(prho.p3().unit());
 99        Vector3 trans2 = pPi.p3().unit()-cPi*prho.p3().unit();
100        _h[imode][1]->fill(cPi);
101        double chi = atan2(trans1.cross(trans2).dot(prho.p3().unit()),trans1.dot(trans2));
102        _h[imode][2]->fill(chi);
103        _p[imode][0]->fill(-0.5*(1-5*sqr(cPi)));
104        _p[imode][1]->fill();
105      }
106    }
107
108
109    /// Normalise histograms etc., after the run
110    void finalize() {
111      for (unsigned int ix=0; ix<2; ++ix) {
112        normalize(_h[ix], 1.0,false);
113        Estimate0DPtr tmp;
114        book(tmp,1,1,1+ix);
115        divide(*_p[ix][0],*_p[ix][1],tmp);
116      }
117    }
118
119    /// @}
120
121
122    /// @name Histograms
123    /// @{
124    Histo1DPtr _h[2][3];
125    CounterPtr _p[2][2];
126    const map<PdgId,unsigned int> mode1   = { {-423,1}, { 213,1}};
127    const map<PdgId,unsigned int> mode1CC = { { 423,1}, {-213,1}};
128    const map<PdgId,unsigned int> mode2   = { {-413,1}, { 213,1}};
129    const map<PdgId,unsigned int> mode2CC = { { 413,1}, {-213,1}};
130    /// @}
131
132
133  };
134
135
136  RIVET_DECLARE_PLUGIN(CLEOII_2003_I611745);
137
138}