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- -> hadrons near phi
 10  class CMD2_1995_I406880 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CMD2_1995_I406880);
 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(_nKpKm      , "TMP/KpKm"    , refData(1, 1, 1));
 27      book(_nK0K0      , "TMP/K0K0"    , refData(1, 1, 2));
 28      book(_n3pi       , "TMP/3pi"     , refData(1, 1, 3));
 29      book(_numEtaGamma, "TMP/EtaGamma", refData(1, 1, 4));
 30    }
 31
 32    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 33      for (const Particle &child : p.children()) {
 34	if(child.children().empty()) {
 35	  --nRes[child.pid()];
 36	  --ncount;
 37	}
 38	else
 39	  findChildren(child,nRes,ncount);
 40      }
 41    }
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      const FinalState& fs = apply<FinalState>(event, "FS");
 46
 47      map<long,int> nCount;
 48      int ntotal(0);
 49      for (const Particle& p : fs.particles()) {
 50	nCount[p.pid()] += 1;
 51	++ntotal;
 52      }
 53      if(ntotal==2) {
 54	if(nCount[321]==1 && nCount[-321]==1)
 55	  _nKpKm->fill(sqrtS()/MeV);
 56	else if(nCount[130]==1 && nCount[310]==1)
 57	  _nK0K0->fill(sqrtS()/MeV);
 58      }
 59      else if(ntotal==3 && nCount[211] == 1 && nCount[-211] == 1 && nCount[111] == 1)
 60	_n3pi->fill(sqrtS()/MeV);
 61
 62      const FinalState& ufs = apply<FinalState>(event, "UFS");
 63      for (const Particle& p : ufs.particles()) {
 64	if(p.children().empty()) continue;
 65	// find the omega
 66	if(p.pid()==221) {
 67	  map<long,int> nRes = nCount;
 68	  int ncount = ntotal;
 69	  findChildren(p,nRes,ncount);
 70	  // eta pi+pi-
 71	  if(ncount!=1) continue;
 72	  bool matched = true;
 73          for(auto const & val : nRes) {
 74	    if(val.first==22) {
 75	      if(val.second !=1) {
 76		matched = false;
 77		break;
 78	      }
 79	    }
 80	    else if(val.second!=0) {
 81	      matched = false;
 82	      break;
 83	    }
 84	  }
 85	  if(matched)
 86	    _numEtaGamma->fill(sqrtS()/MeV);
 87	}
 88      }
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      double fact = crossSection()/ sumOfWeights() /nanobarn;
 95      Estimate1DPtr mult;
 96      scale(_nKpKm,fact);
 97      book(mult, 1, 1, 1);
 98      barchart(_nKpKm,mult);
 99      scale(_nK0K0,fact);
100      book(mult, 1, 1, 2);
101      barchart(_nK0K0,mult);
102      scale(_n3pi,fact);
103      book(mult, 1, 1, 3);
104      barchart(_n3pi,mult);
105      scale(_numEtaGamma,fact);
106      book(mult, 1, 1, 4);
107      barchart(_numEtaGamma,mult);
108    }
109    /// @}
110
111
112    /// @name Histograms
113    /// @{
114    Histo1DPtr _nKpKm,_nK0K0,_n3pi,_numEtaGamma;
115    /// @}
116
117
118  };
119
120
121  RIVET_DECLARE_PLUGIN(CMD2_1995_I406880);
122
123
124}