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_2003_I616446 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(CMD2_2003_I616446);
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(_numOmegaPi, "TMP/OmegaPi");
27    }
28
29
30    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
31      for (const Particle &child : p.children()) {
32	if(child.children().empty()) {
33	  nRes[child.pid()]+=1;
34	  ++ncount;
35	}
36	else
37	  findChildren(child,nRes,ncount);
38      }
39    }
40
41    /// Perform the per-event analysis
42    void analyze(const Event& event) {
43      const FinalState& fs = apply<FinalState>(event, "FS");
44
45      map<long,int> nCount;
46      int ntotal(0);
47      for (const Particle& p : fs.particles()) {
48	nCount[p.pid()] += 1;
49	++ntotal;
50      }
51      // three particles (pi0 pi0 gamma)
52      if(ntotal!=3) vetoEvent;
53      const FinalState& ufs = apply<FinalState>(event, "UFS");
54      for (const Particle& p : ufs.particles()) {
55	if(p.children().empty()) continue;
56	// find the omega
57	if(p.pid()==223) {
58	  map<long,int> nRes;
59	  int ncount(0);
60	  findChildren(p,nRes,ncount);
61	  // only omega to pi0 gamma mode
62	  if(ncount!=2) continue;
63	  if(nRes[111]!=1 || nRes[22]!=1) continue;
64	  // omega pi0
65	  if(nCount[111]-nRes[111]==1)
66	    _numOmegaPi->fill();
67	}
68      }
69    }
70
71
72    /// Normalise histograms etc., after the run
73    void finalize() {
74
75      double sigma = _numOmegaPi->val();
76      double error = _numOmegaPi->err();
77      sigma *= crossSection()/ sumOfWeights() /nanobarn;
78      error *= crossSection()/ sumOfWeights() /nanobarn;
79      Estimate1DPtr mult;
80      book(mult, 1, 1, 1);
81      for (auto& b : mult->bins()) {
82        if (inRange(sqrtS()/MeV, b.xMin(), b.xMax())) {
83          b.set(sigma, error);
84        }
85      }
86    }
87
88    /// @name Histograms
89    /// @{
90    CounterPtr _numOmegaPi;
91    /// @}
92
93  };
94
95
96  RIVET_DECLARE_PLUGIN(CMD2_2003_I616446);
97
98
99}