1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Mass distribution in $e^+e^-\to e^+e^-D\bar{D}$ via $\gamma\gamma\to D\bar{D}$
 11  class BELLE_2006_I700451 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2006_I700451);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(FinalState(),"FS");
 26      declare(UnstableParticles(Cuts::abspid==411 || Cuts::abspid==421), "UFS");
 27      // histos
 28      for (unsigned int ix=0;ix<3;++ix) {
 29        book(_h_mass[ix],1,1,1+ix);
 30        if (ix<2) book(_h_mass_cut[ix],2,1,1+ix);
 31      }
 32      book(_h_angle,3,1,1);
 33    }
 34
 35    void findChildren(const Particle& p, map<long,int>& nRes, int &ncount) {
 36      for (const Particle& child : p.children()) {
 37        if (child.children().empty()) {
 38          --nRes[child.pid()];
 39          --ncount;
 40        }
 41        else {
 42          findChildren(child,nRes,ncount);
 43        }
 44      }
 45    }
 46
 47    bool findScattered(Particle beam, double& q2) {
 48      bool found = false;
 49      Particle scat = beam;
 50      while (!scat.children().empty()) {
 51        found = false;
 52        for (const Particle & p : scat.children()) {
 53          if (p.pid()==scat.pid()) {
 54            scat=p;
 55            found=true;
 56            break;
 57          }
 58        }
 59        if (!found) break;
 60      }
 61      if (!found) return false;
 62      q2 = -(beam.mom() - scat.mom()).mass2();
 63      return true;
 64    }
 65
 66    /// Perform the per-event analysis
 67    void analyze(const Event& event) {
 68      // find scattered leptons and calc Q2
 69      const Beam& beams = apply<Beam>(event, "Beams");
 70      double q12 = -1, q22 = -1;
 71      if (!findScattered(beams.beams().first,  q12)) vetoEvent;
 72      if (!findScattered(beams.beams().second, q22)) vetoEvent;
 73      // check the final state
 74      const FinalState & fs = apply<FinalState>(event, "FS");
 75      map<long,int> nCount;
 76      int ntotal(0);
 77      for (const Particle& p : fs.particles()) {
 78        nCount[p.pid()] += 1;
 79        ++ntotal;
 80      }
 81      // find the meson
 82      const FinalState& ufs = apply<FinalState>(event, "UFS");
 83      for (const Particle& p1 : ufs.particles()) {
 84        if (p1.children().empty() || p1.pid()<0) continue;
 85        bool matched=false;
 86        map<long,int> nRes = nCount;
 87        int ncount = ntotal;
 88        findChildren(p1,nRes,ncount);
 89        for (const Particle& p2 : ufs.particles(Cuts::pid==-p1.pid())) {
 90          map<long,int> nRes2 = nRes;
 91          int ncount2 = ncount;
 92          findChildren(p2,nRes2,ncount2);
 93          matched = true;
 94          for (const auto& val : nRes2) {
 95            if (abs(val.first)==11) {
 96              if (val.second!=1) {
 97                matched = false;
 98                break;
 99              }
100            }
101            else if (val.second!=0) {
102              matched = false;
103              break;
104            }
105          }
106          if (matched) {
107            FourMomentum pDD = p1.mom()+p2.mom();
108            double mDD = pDD.mass();
109            if (p1.abspid()==421) _h_mass[0]->fill(mDD);
110            else                  _h_mass[1]->fill(mDD);
111            _h_mass[2]->fill(mDD);
112            LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(pDD.betaVec());
113            FourMomentum pD = boost.transform(p1.mom());
114            double cTheta = abs(pD.z()/pD.p3().mod());
115            if (abs(cTheta)<0.5) _h_mass_cut[0]->fill(mDD);
116            else                 _h_mass_cut[1]->fill(mDD);
117            if (mDD>3.91 && mDD<3.95) _h_angle->fill(cTheta);
118            break;
119          }
120        }
121        if (matched) break;
122      }
123    }
124
125
126    /// Normalise histograms etc., after the run
127    void finalize() {
128	    normalize(_h_mass, 1.0, false);
129	    normalize(_h_mass_cut, 1.0, false);
130      normalize(_h_angle, 1.0, false);
131    }
132
133    /// @}
134
135
136    /// @name Histograms
137    /// @{
138    Histo1DPtr  _h_mass[3],_h_mass_cut[2],_h_angle;
139    /// @}
140
141
142  };
143
144
145  RIVET_DECLARE_PLUGIN(BELLE_2006_I700451);
146
147}