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 gamma gamma -> eta' pi+pi-
 10  class BELLE_2018_I1672149 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2018_I1672149);
 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      declare(FinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25      // counters
 26      for (unsigned int ix=0; ix<2; ++ix) {
 27        book(_sigma[ix],"TMP/sigma_"+toString(ix+1),refData(1+ix,1,1));
 28      }
 29    }
 30
 31    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 32      for (const Particle &child : p.children()) {
 33        if (child.children().empty()) {
 34          nRes[child.pid()]-=1;
 35          --ncount;
 36        }
 37        else {
 38          findChildren(child,nRes,ncount);
 39        }
 40      }
 41    }
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      const FinalState& fs = apply<FinalState>(event, "FS");
 46      // find the final-state particles
 47      map<long,int> nCount;
 48      int ntotal(0);
 49      for (const Particle& p : fs.particles()) {
 50        nCount[p.pid()] += 1;
 51        ++ntotal;
 52      }
 53      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 54      // first check for f_2 eta'
 55      bool foundf2(false);
 56      Particles f2s = ufs.particles(Cuts::pid==225);
 57      Particles etaps = ufs.particles(Cuts::pid==331);
 58      for (const Particle & f2 : f2s) {
 59        bool matched=false;
 60       	map<long,int> nRes=nCount;
 61       	int ncount = ntotal;
 62       	findChildren(f2,nRes,ncount);
 63        for (const Particle & etap : etaps) {
 64          map<long,int> nRes2=nRes;
 65          int ncount2 = ncount;
 66          findChildren(etap,nRes2,ncount2);
 67          if (ncount2 !=0 ) continue;
 68          matched = true;
 69          for (const auto& val : nRes2) {
 70            if (val.second!=0) {
 71              matched = false;
 72              break;
 73            }
 74          }
 75          if (matched) {
 76            foundf2 = true;
 77            break;
 78          }
 79        }
 80        if (foundf2) break;
 81      }
 82      // if we have the f_2 eta' state
 83      if (foundf2) {
 84        _sigma[0]->fill(sqrtS());
 85        if (sqrtS()>2.26*GeV) vetoEvent;
 86      }
 87      // see if we have eta' pi+pi-
 88      bool foundetap(false);
 89      for (const Particle & etap : etaps) {
 90        bool matched=false;
 91       	map<long,int> nRes=nCount;
 92       	int ncount = ntotal;
 93       	findChildren(etap,nRes,ncount);
 94        if (ncount !=2 ) continue;
 95        matched=true;
 96        for (const auto& val : nRes) {
 97          if (abs(val.first)==211) {
 98            if (val.second!=1) {
 99              matched = false;
100              break;
101            }
102          }
103          else if (val.second!=0) {
104            matched = false;
105            break;
106          }
107        }
108        if (matched) {
109          foundetap = true;
110          break;
111        }
112      }
113      // check there's no eta_c
114      if (foundetap && sqrtS()>2.62 && sqrtS()<3.06) {
115        for (const Particle& etac : ufs.particles(Cuts::pid==441)) {
116          bool matched=false;
117          map<long,int> nRes=nCount;
118          int ncount = ntotal;
119          findChildren(etac,nRes,ncount);
120          for(const auto& val : nRes) {
121            if (val.second!=0) {
122              matched = false;
123              break;
124            }
125          }
126          if (matched) {
127            foundetap=false;
128            break;
129          }
130        }
131      }
132      if (foundetap) _sigma[1]->fill(sqrtS());
133    }
134
135    /// Normalise histograms etc., after the run
136    void finalize() {
137      const double fact = crossSection()/nanobarn/sumOfWeights();
138      // loop over tables in paper
139      for(unsigned int ix=0;ix<2;++ix) {
140        scale(_sigma[ix],fact);
141        Estimate1DPtr tmp;
142        book(tmp, ix+1, 1, 1);
143        barchart(_sigma[ix],tmp);
144      }
145    }
146
147    /// @}
148
149
150    /// @name Histograms
151    /// @{
152    Histo1DPtr _sigma[2];
153    /// @}
154
155
156  };
157
158
159  RIVET_DECLARE_PLUGIN(BELLE_2018_I1672149);
160
161}