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