1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief phi -> 3 pion decay
 9  class SND_2001_I558279 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(SND_2001_I558279);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // Initialise and register projections
22      declare(UnstableParticles(), "UFS");
23      // Book histograms
24      book(_h_pm, 1, 1, 1);
25      book(_h_p0, 2, 1, 1);
26
27    }
28
29    void findDecayProducts(const Particle & mother, unsigned int & nstable, Particles &pip,
30			   Particles &pim, Particles &pi0) {
31      for(const Particle & p : mother.children()) {
32	int id = p.pid();
33	if (id == PID::PIMINUS ) {
34	  pim.push_back(p);
35	  ++nstable;
36	}
37        else if (id == PID::PIPLUS) {
38	  pip.push_back(p);
39	  ++nstable;
40	}
41	else if (id == PID::PI0) {
42	  pi0.push_back(p);
43	  ++nstable;
44	}
45	else if ( !p.children().empty() ) {
46          findDecayProducts(p, nstable, pip,pim,pi0);
47        }
48        else
49          ++nstable;
50      }
51    }
52
53    /// Perform the per-event analysis
54    void analyze(const Event& event) {
55      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==333)) {
56	Particles pip,pim,pi0;
57	unsigned int nstable(0);
58	findDecayProducts(p, nstable, pip,pim,pi0);
59	if(nstable==3 && pip.size()==1 && pim.size()==1&&pi0.size()==1) {
60	  _h_pm->fill((pip[0].momentum()+pim[0].momentum()).mass()/MeV);
61	  _h_p0->fill((pip[0].momentum()+pi0[0].momentum()).mass()/MeV);
62	  _h_p0->fill((pim[0].momentum()+pi0[0].momentum()).mass()/MeV);
63	}
64      }
65    }
66
67
68    /// Normalise histograms etc., after the run
69    void finalize() {
70
71      normalize(_h_pm);
72      normalize(_h_p0);
73
74    }
75
76    /// @}
77
78
79    /// @name Histograms
80    /// @{
81    Histo1DPtr _h_pm,_h_p0;
82    /// @}
83
84
85  };
86
87
88  RIVET_DECLARE_PLUGIN(SND_2001_I558279);
89
90
91}