1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/DecayedParticles.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Cross sections for $e^+e^-\to J/\psi + \pi^+\pi^-$, $\pi^0\pi^0$ or $K^+K^-$
 11  class CLEO_2006_I710864 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2006_I710864);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      _mode = 0;
 24      if ( getOption("MODE") == "SIGMA" ) _mode = 0;
 25      else if ( getOption("MODE") == "DECAY" ) _mode = 1;
 26      // set the PDG code
 27      _pid = getOption<int>("PID", 9030443);
 28      // setup for cross sections
 29      if (_mode==0) {
 30        // projections
 31        declare(FinalState(), "FS");
 32        declare(UnstableParticles(Cuts::pid== 443), "UFS");
 33        // histograms
 34        for (unsigned int ix=0; ix<3; ++ix) {
 35          book(_sigma[ix],1,1,1+ix);
 36        }
 37      for (const string& en : _sigma[0].binning().edges<0>()) {
 38        const size_t idx = en.find("-");
 39        if(idx!=std::string::npos) {
 40          const double emin = std::stod(en.substr(0,idx));
 41          const double emax = std::stod(en.substr(idx+1,string::npos));
 42          if(inRange(sqrtS()/GeV, emin, emax)) {
 43            _ecms = en;
 44            break;
 45          }
 46        }
 47        else {
 48          const double end = std::stod(en)*GeV;
 49          if (isCompatibleWithSqrtS(end)) {
 50            _ecms = en;
 51            break;
 52          }
 53        }
 54      }
 55      if (_ecms.empty())
 56        MSG_ERROR("Beam energy incompatible with analysis.");
 57      }
 58      // setup for decays
 59      else {
 60        // projections
 61        UnstableParticles ufs = UnstableParticles(Cuts::abspid==_pid);
 62        declare(ufs, "UFS");
 63        DecayedParticles PSI(ufs);
 64        PSI.addStable(443);
 65        declare(PSI, "PSI");
 66        book(_h,2,1,1);
 67      }
 68    }
 69
 70    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 71      for (const Particle& child : p.children()) {
 72        if (child.children().empty()) {
 73          --nRes[child.pid()];
 74          --ncount;
 75        }
 76        else {
 77          findChildren(child,nRes,ncount);
 78        }
 79      }
 80    }
 81
 82    /// Perform the per-event analysis
 83    void analyze(const Event& event) {
 84      // cross section
 85      if (_mode==0) {
 86        const FinalState& fs = apply<FinalState>(event, "FS");
 87        map<long,int> nCount;
 88        int ntotal(0);
 89        for (const Particle& p : fs.particles()) {
 90          nCount[p.pid()] += 1;
 91          ++ntotal;
 92        }
 93        for (const Particle& psi : apply<UnstableParticles>(event, "UFS").particles()) {
 94          if (psi.children().empty()) continue;
 95          map<long,int> nRes = nCount;
 96          int ncount = ntotal;
 97          findChildren(psi,nRes,ncount);
 98          if (ncount!=2) continue;
 99          if (nRes.find(211)!=nRes.end() && nRes.find(-211)!=nRes.end() &&
100             nRes[211]==1 && nRes[-211]==1) {
101            _sigma[0]->fill(_ecms);
102            break;
103          }
104          else if (nRes.find(111)!=nRes.end() && nRes[111]==2) {
105            _sigma[1]->fill(_ecms);
106            break;
107          }
108          else if (nRes.find(321)!=nRes.end() && nRes.find(-321)!=nRes.end() &&
109                   nRes[321]==1 && nRes[-321]==1) {
110            _sigma[2]->fill(_ecms);
111            break;
112          }
113        }
114      }
115      // decays
116      else {
117        DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
118        for (unsigned int ix=0; ix<PSI.decaying().size(); ++ix) {
119          if (!PSI.modeMatches(ix,3,mode)) continue;
120          const Particle & pip  = PSI.decayProducts()[ix].at( 211)[0];
121          const Particle & pim  = PSI.decayProducts()[ix].at(-211)[0];
122          double mpipi = (pip.mom()+pim.mom()).mass();
123          _h->fill(mpipi);
124        }
125      }
126    }
127
128
129    /// Normalise histograms etc., after the run
130    void finalize() {
131      // sigma
132      if (_mode==0) {
133        scale(_sigma, crossSection()/ sumOfWeights() /picobarn);
134      }
135      // decay
136      else {
137        normalize(_h, 1.0, false);
138      }
139    }
140
141    /// @}
142
143
144    /// @name Histograms
145    /// @{
146    int _mode;
147    int _pid;
148    Histo1DPtr _h;
149    BinnedHistoPtr<string> _sigma[3];
150    string _ecms;
151    const map<PdgId,unsigned int> mode = { { 211,1}, {-211,1}, { 443,1}};
152    /// @}
153
154
155  };
156
157
158  RIVET_DECLARE_PLUGIN(CLEO_2006_I710864);
159
160}