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 Add a short analysis description here
 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, "TMP/OmegaPiPi");
 27      book(_numEtaPiPi, "TMP/EtaPiPi");
 28      book(_num5Pi, "TMP/5Pi");
 29
 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
 46      const FinalState& fs = apply<FinalState>(event, "FS");
 47
 48      map<long,int> nCount;
 49      int ntotal(0);
 50      for (const Particle& p : fs.particles()) {
 51	nCount[p.pid()] += 1;
 52	++ntotal;
 53      }
 54      const FinalState& ufs = apply<FinalState>(event, "UFS");
 55      bool foundRes = false;
 56      for (const Particle& p : ufs.particles()) {
 57	if(p.children().empty()) continue;
 58	// find the eta
 59	if(p.pid()==221) {
 60	  map<long,int> nRes = nCount;
 61	  int ncount = ntotal;
 62	  findChildren(p,nRes,ncount);
 63	  // eta pi+pi-
 64	  if(ncount!=2) continue;
 65	  bool matched = true;
 66	  for(auto const & val : nRes) {
 67	    if(abs(val.first)==211) {
 68	      if(val.second !=1) {
 69		matched = false;
 70		break;
 71	      }
 72	    }
 73	    else if(val.second!=0) {
 74	      matched = false;
 75	      break;
 76	    }
 77	  }
 78	  if(matched) {
 79	    _numEtaPiPi->fill();
 80	    foundRes = true;
 81	  }
 82	}
 83	// find the omega
 84	else if(p.pid()==223) {
 85	  map<long,int> nRes = nCount;
 86	  int ncount = ntotal;
 87	  findChildren(p,nRes,ncount);
 88	  // eta pi+pi-
 89	  if(ncount!=2) continue;
 90	  bool matched = true;
 91	  for(auto const & val : nRes) {
 92	    if(abs(val.first)==211) {
 93	      if(val.second !=1) {
 94		matched = false;
 95		break;
 96	      }
 97	    }
 98	    else if(val.second!=0) {
 99	      matched = false;
100	      break;
101	    }
102	  }
103	  if(matched) {
104	    _numOmegaPiPi->fill();
105	    foundRes = true;
106	  }
107	}
108      }
109
110      if(foundRes) vetoEvent;
111      if(nCount[-211]==2&&nCount[211]==2&&nCount[111]==1)
112	_num5Pi->fill();
113    }
114
115    /// Normalise histograms etc., after the run
116    void finalize() {
117      for(unsigned int ix=1;ix<4;++ix) {
118	double sigma,error;
119	if(ix==1) {
120	  sigma = _numOmegaPiPi->val();
121	  error = _numOmegaPiPi->err();
122	}
123	else if(ix==2) {
124	  sigma = _numEtaPiPi->val();
125	  error = _numEtaPiPi->err();
126	}
127	else {
128	  sigma = _num5Pi->val();
129	  error = _num5Pi->err();
130	}
131	sigma *= crossSection()/ sumOfWeights() /nanobarn;
132	error *= crossSection()/ sumOfWeights() /nanobarn; 
133	Scatter2D temphisto(refData(ix, 1, 1));
134	Scatter2DPtr mult;
135	book(mult, ix, 1, 1);
136	for (size_t b = 0; b < temphisto.numPoints(); b++) {
137	  const double x  = temphisto.point(b).x();
138	  pair<double,double> ex = temphisto.point(b).xErrs();
139	  pair<double,double> ex2 = ex;
140	  if(ex2.first ==0.) ex2. first=0.0001;
141	  if(ex2.second==0.) ex2.second=0.0001;
142	  if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
143	    mult->addPoint(x, sigma, ex, make_pair(error,error));
144	  }
145	  else {
146	    mult->addPoint(x, 0., ex, make_pair(0.,.0));
147	  }
148	}
149      }
150    }
151
152    //@}
153
154    /// @name Histograms
155    //@{
156    CounterPtr _numEtaPiPi,_numOmegaPiPi,_num5Pi;
157    //@}
158
159  };
160
161
162  // The hook for the plugin system
163  RIVET_DECLARE_PLUGIN(CMD2_2000_I532970);
164
165
166}