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- > pi pi psi_2
10  class BESIII_2023_I2158340 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2023_I2158340);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      //check CMS energy
23      if (!inRange(sqrtS(),4.23,4.70)) {
24        throw Error("Unexpected sqrtS,must be between 4.23 and 4.70 GeV");
25      }
26      // Initialise and register projections
27      declare(FinalState(), "FS");
28      declare(UnstableParticles(), "UFS");
29      // histos
30      for (unsigned int ix=0; ix<2; ++ix) {
31        book(_n[ix], "TMP/n_"+toString(ix+1),refData(1,1,1));
32      }
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()]-=1;
39          --ncount;
40        }
41        else {
42          findChildren(child,nRes,ncount);
43        }
44      }
45    }
46
47    /// Perform the per-event analysis
48    void analyze(const Event& event) {
49      const FinalState& fs = apply<FinalState>(event, "FS");
50      map<long,int> nCount;
51      int ntotal(0);
52      for (const Particle& p : fs.particles()) {
53      	nCount[p.pid()] += 1;
54      	++ntotal;
55      }
56      const FinalState& ufs = apply<FinalState>(event, "UFS");
57      for (const Particle& chi : ufs.particles(Cuts::pid==20445)) {
58        if (chi.children().empty()) continue;
59        map<long,int> nRes=nCount;
60        int ncount = ntotal;
61        findChildren(chi,nRes,ncount);
62        if (ncount!=2) continue;
63        if (nRes.find(211)!=nRes.end() && nRes.find(-211)!=nRes.end() && nRes[211]==1 && nRes[-211]==1) {
64          _n[0]->fill(sqrtS());
65          break;
66        }
67        else if (nRes.find(111)!=nRes.end() && nRes[111]==2) {
68          _n[1]->fill(sqrtS());
69          break;
70        }
71      }
72    }
73
74
75    /// Normalise histograms etc., after the run
76    void finalize() {
77      Estimate1DPtr mult;
78      book(mult, 1, 1, 1);
79      divide(_n[1],_n[0],mult);
80    }
81
82    /// @}
83
84
85    /// @name Histograms
86    /// @{
87    Histo1DPtr _n[2];
88    /// @}
89
90
91  };
92
93
94  RIVET_DECLARE_PLUGIN(BESIII_2023_I2158340);
95
96}