1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief e+ e- > pi0 pi0 hc
 10  class BESIII_2014_I1318650 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2014_I1318650);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // set the PDG code
 23      _pid = getOption<double>("PID", 9030443);
 24      // projections
 25      declare(FinalState(), "FS");
 26      declare(UnstableParticles(), "UFS");
 27      // histograms
 28      for (unsigned int ix=0; ix<2; ++ix) {
 29        book(_sigma[ix], 1, 1, 1+ix);
 30      }
 31
 32      for (const string& en : _sigma[0].binning().edges<0>()) {
 33        const double end = std::stod(en)*GeV;
 34        if (isCompatibleWithSqrtS(end)) {
 35          _ecms = en;
 36          break;
 37        }
 38      }
 39      if (_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
 40
 41    }
 42
 43    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 44      for (const Particle &child : p.children()) {
 45        if (child.children().empty()) {
 46          --nRes[child.pid()];
 47          --ncount;
 48        }
 49        else {
 50          findChildren(child,nRes,ncount);
 51        }
 52      }
 53    }
 54
 55    /// Perform the per-event analysis
 56    void analyze(const Event& event) {
 57      const FinalState& fs = apply<FinalState>(event, "FS");
 58      map<long,int> nCount;
 59      int ntotal(0);
 60      for (const Particle& p : fs.particles()) {
 61        nCount[p.pid()] += 1;
 62        ++ntotal;
 63      }
 64      const UnstableParticles & ufs = apply<UnstableParticles>(event, "UFS");
 65      // loop over any h_c
 66      for (const Particle& hc : ufs.particles(Cuts::pid==10443)) {
 67        if (hc.children().empty()) continue;
 68        // find the h_c
 69        map<long,int> nRes = nCount;
 70        int ncount = ntotal;
 71        findChildren(hc,nRes,ncount);
 72        // h_c pi0 pi0
 73        if(ncount!=2) continue;
 74        bool matched = true;
 75        for (const auto& val : nRes) {
 76          if (abs(val.first)==111) {
 77            if (val.second !=2) {
 78              matched = false;
 79              break;
 80            }
 81          }
 82          else if(val.second!=0) {
 83            matched = false;
 84            break;
 85          }
 86        }
 87        if (matched) {
 88          _sigma[0]->fill(_ecms);
 89          if (hc.parents().empty()) break;
 90          Particle Zc = hc.parents()[0];
 91          if (Zc.pid()==_pid && Zc.children().size()==2 &&
 92             (Zc.children()[0].pid()==PID::PI0 ||
 93              Zc.children()[1].pid()==PID::PI0)) _sigma[1]->fill(_ecms);
 94          break;
 95        }
 96      }
 97    }
 98
 99
100    /// Normalise histograms etc., after the run
101    void finalize() {
102      scale(_sigma, crossSection()/ sumOfWeights() /picobarn);
103    }
104
105    /// @}
106
107
108    /// @name Histograms
109    /// @{
110    BinnedHistoPtr<string> _sigma[2];
111    string _ecms;
112    int _pid;
113    /// @}
114
115
116  };
117
118
119  RIVET_DECLARE_PLUGIN(BESIII_2014_I1318650);
120
121}