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 Cross section for  $e^+e^-\to\omega \chi_{c(1,2)}$ and $X(3872)$ for $\sqrt{s}=4.66$ to 4.95 GeV
 10  class BESIII_2024_I2779452 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2024_I2779452);
 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      // histos
 28      for(unsigned int ix=0;ix<3;++ix)
 29        book(_h[ix],1+ix,1,1);
 30      for(const string& en : _h[0].binning().edges<0>()) {
 31        const double end = std::stod(en)*GeV;
 32        if (isCompatibleWithSqrtS(end)) {
 33          _ecms = en;
 34          break;
 35        }
 36      }
 37      if (_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
 38    }
 39
 40    void findChildren2(const Particle& p, map<long,int> & nRes, int &ncount) {
 41      for (const Particle &child : p.children()) {
 42        if(child.pid()==443 ) {
 43          nRes[child.pid()] += 1;
 44          ++ncount;
 45        }
 46        else if (child.children().empty()) {
 47          ++nRes[child.pid()];
 48          ++ncount;
 49        }
 50        else {
 51          findChildren(child,nRes,ncount);
 52        }
 53      }
 54    }
 55
 56    void findChildren(const Particle& p, map<long,int> & nRes, int &ncount) {
 57      for (const Particle &child : p.children()) {
 58        if (child.children().empty()) {
 59          --nRes[child.pid()];
 60          --ncount;
 61        }
 62        else {
 63          findChildren(child,nRes,ncount);
 64        }
 65      }
 66    }
 67
 68    /// Perform the per-event analysis
 69    void analyze(const Event& event) {
 70      const FinalState& fs = apply<FinalState>(event, "FS");
 71      map<long,int> nCount;
 72      int ntotal(0);
 73      for (const Particle& p: fs.particles()) {
 74      	nCount[p.pid()] += 1;
 75      	++ntotal;
 76      }
 77      const FinalState& ufs = apply<FinalState>(event, "UFS");
 78      for (const Particle& XX : ufs.particles(Cuts::pid==_pid || Cuts::pid==20443 || Cuts::pid==445)) {
 79      	if (XX.children().empty()) continue;
 80        bool matched = false;
 81       	map<long,int> nRes = nCount;
 82       	int ncount = ntotal;
 83       	findChildren(XX,nRes,ncount);
 84        for (const Particle & omega : ufs.particles(Cuts::pid==223)) {
 85          if (omega.parents()[0].pid()==_pid || omega.children().empty()) continue;
 86          map<long,int> nRes2 = nRes;
 87          int ncount2 = ncount;
 88          findChildren(omega,nRes2,ncount2);
 89          matched = true;
 90          for (const auto& val : nRes2) {
 91            if (val.second!=0) {
 92              matched = false;
 93              break;
 94            }
 95          }
 96          if(matched) {
 97            if     (XX.pid()==20443) _h[1]->fill(_ecms);
 98            else if(XX.pid()==  445) _h[2]->fill(_ecms);
 99            // check decay mode
100            else {
101              int nOut=0;
102              map<long,int> nDecay;
103              findChildren2(XX, nDecay,nOut);
104              if(nOut==3 && nDecay[443]==1 && nDecay[211]==1 && nDecay[-211]==1)
105                _h[0]->fill(_ecms);
106            }
107            break;
108          }
109        }
110        if(matched) break;
111      }
112    }
113
114
115    /// Normalise histograms etc., after the run
116    void finalize() {
117      scale(_h, crossSection()/ sumOfWeights() /picobarn);
118    }
119
120    /// @}
121
122
123    /// @name Histograms
124    /// @{
125    int _pid;
126    BinnedHistoPtr<string> _h[3];
127    string _ecms;
128    /// @}
129
130
131  };
132
133
134  RIVET_DECLARE_PLUGIN(BESIII_2024_I2779452);
135
136}