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