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	Scatter2D temphisto(refData(1, 1, 1+ix));
106	Scatter2DPtr  mult;
107	book(mult, 1, 1, 1+ix);
108	for (size_t b = 0; b < temphisto.numPoints(); b++) {
109	  const double x  = temphisto.point(b).x();
110	  pair<double,double> ex = temphisto.point(b).xErrs();
111	  pair<double,double> ex2 = ex;
112	  if(ex2.first ==0.) ex2. first=0.0001;
113	  if(ex2.second==0.) ex2.second=0.0001;
114	  if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
115	    mult->addPoint(x, sigma, ex, make_pair(error,error));
116	  }
117	  else {
118	    mult->addPoint(x, 0., ex, make_pair(0.,.0));
119	  }
120	}
121      }
122    }
123
124    /// @}
125
126
127    /// @name Histograms
128    /// @{
129    CounterPtr _nMeson[2];
130    /// @}
131
132
133  };
134
135
136  RIVET_DECLARE_PLUGIN(BESIII_2022_I2047667);
137
138}