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- > eta omega and omega pi0
 10  class BESIII_2022_I2047667 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2047667);
 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      for(unsigned int ix=0;ix<2;++ix)
 26        book(_nMeson[ix], 1, 1, 1+ix);
 27      for (const string& en : _nMeson[0].binning().edges<0>()) {
 28        const double end = std::stod(en)*GeV;
 29        if (isCompatibleWithSqrtS(end)) {
 30          _ecms = en;
 31          break;
 32        }
 33      }
 34    }
 35
 36    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 37      for( const Particle &child : p.children()) {
 38	if(child.children().empty()) {
 39	  --nRes[child.pid()];
 40	  --ncount;
 41	}
 42	else
 43	  findChildren(child,nRes,ncount);
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      const FinalState& fs = apply<FinalState>(event, "FS");
 50
 51      map<long,int> nCount;
 52      int ntotal(0);
 53      for (const Particle& p : fs.particles()) {
 54      	nCount[p.pid()] += 1;
 55      	++ntotal;
 56      }
 57      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 58      for (const Particle& p : ufs.particles(Cuts::pid==223)) {
 59	if(p.children().empty()) continue;
 60	map<long,int> nRes = nCount;
 61	int ncount = ntotal;
 62	findChildren(p,nRes,ncount);
 63	if(ncount==1) {
 64	  bool matched=true;
 65	  for(auto const & val : nRes) {
 66	    if( val.first==111 ) {
 67	      if(val.second !=1) {
 68		matched = false;
 69		break;
 70	      }
 71	    }
 72	    else if(val.second!=0) {
 73	      matched = false;
 74	      break;
 75	    }
 76	  }
 77	  if(matched) {
 78	    _nMeson[0]->fill(_ecms);
 79	    continue;
 80	  }
 81	}
 82	else {
 83	  // check for eta
 84	  for (const Particle& p2 : ufs.particles(Cuts::pid==221)) {
 85	    map<long,int> nResB = nRes;
 86	    int ncountB = ncount;
 87	    findChildren(p2,nResB,ncountB);
 88	    if(ncountB!=0) continue;
 89	    bool matched = true;
 90	    for(auto const & val : nResB) {
 91	      if(val.second!=0) {
 92		matched = false;
 93		break;
 94	      }
 95	    }
 96	    if(matched) {
 97	      _nMeson[1]->fill(_ecms);
 98	      break;
 99	    }
100	  }
101	}
102      }
103    }
104
105
106    /// Normalise histograms etc., after the run
107    void finalize() {
108      double fact = crossSection()/picobarn/sumOfWeights();
109      for(unsigned int ix=0;ix<2; ++ix)
110        scale(_nMeson[ix],fact);
111    }
112
113    /// @}
114
115
116    /// @name Histograms
117    /// @{
118    BinnedHistoPtr<string> _nMeson[2];
119    string _ecms;
120    /// @}
121
122
123  };
124
125
126  RIVET_DECLARE_PLUGIN(BESIII_2022_I2047667);
127
128}