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