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 pi+pi-
 10  class DM2_1988_I264144 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(DM2_1988_I264144);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(FinalState(), "FS");
 25      declare(UnstableParticles(), "UFS");
 26      book(_numEtaPiPi, 1,1,1);
 27      for (const string& en : _numEtaPiPi.binning().edges<0>()) {
 28        const size_t idx = en.find("-");
 29        if(idx!=std::string::npos) {
 30          const double emin = std::stod(en.substr(0,idx));
 31          const double emax = std::stod(en.substr(idx+1,string::npos));
 32          if(inRange(sqrtS()/MeV, emin, emax)) {
 33            _ecms = en;
 34            break;
 35          }
 36        }
 37        else {
 38          const double end = std::stod(en)*MeV;
 39          if (isCompatibleWithSqrtS(end)) {
 40            _ecms = en;
 41            break;
 42          }
 43        }
 44      }
 45      if (_ecms.empty())
 46        MSG_ERROR("Beam energy incompatible with analysis.");
 47    }
 48
 49    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 50      for (const Particle &child : p.children()) {
 51	if(child.children().empty()) {
 52	  --nRes[child.pid()];
 53	  --ncount;
 54	}
 55	else
 56	  findChildren(child,nRes,ncount);
 57      }
 58    }
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      const FinalState& fs = apply<FinalState>(event, "FS");
 63
 64      map<long,int> nCount;
 65      int ntotal(0);
 66      for (const Particle& p : fs.particles()) {
 67	nCount[p.pid()] += 1;
 68	++ntotal;
 69      }
 70      const FinalState& ufs = apply<FinalState>(event, "UFS");
 71      for (const Particle& p : ufs.particles()) {
 72	if(p.children().empty()) continue;
 73	// find the omega
 74	if(p.pid()==221) {
 75	  map<long,int> nRes = nCount;
 76	  int ncount = ntotal;
 77	  findChildren(p,nRes,ncount);
 78	  // eta pi+pi-
 79	  if(ncount!=2) continue;
 80	  bool matched = true;
 81	  for(auto const & val : nRes) {
 82	    if(abs(val.first)==211) {
 83	      if(val.second !=1) {
 84		matched = false;
 85		break;
 86	      }
 87	    }
 88	    else if(val.second!=0) {
 89	      matched = false;
 90	      break;
 91	    }
 92	  }
 93	  if(matched)
 94	    _numEtaPiPi->fill(_ecms);
 95	}
 96      }
 97    }
 98
 99
100    /// Normalise histograms etc., after the run
101    void finalize() {
102      scale(_numEtaPiPi, crossSection()/ sumOfWeights() /nanobarn);
103    }
104
105    /// @}
106
107
108    /// @name Histograms
109    /// @{
110    BinnedHistoPtr<string> _numEtaPiPi;
111    string _ecms;
112    /// @}
113
114
115  };
116
117
118  RIVET_DECLARE_PLUGIN(DM2_1988_I264144);
119
120
121}