1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief phi spectrum at Upsilon 4 and 5s
 9  class CLEO_2007_I728872 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2007_I728872);
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
25      // Book histograms
26      book(_h_4S, 1, 1, 5);
27      book(_h_5S, 2, 1, 5);
28      book(_c_4S, "/TMP/N4S");
29      book(_c_5S, "/TMP/N5S");
30    }
31
32    void findDecayProducts(const Particle & p, Particles & phi) {
33      for(const Particle & child : p.children()) {
34	if(child.pid()==333)
35	  phi.push_back(child);
36	else if(!child.children().empty())
37	  findDecayProducts(child,phi);
38      }
39    }
40
41    /// Perform the per-event analysis
42    void analyze(const Event& event) {
43      // Find the Upsilons among the unstables
44      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
45      Particles upsilons = ufs.particles(Cuts::pid==300553 || Cuts::pid==400553|| Cuts::pid==9000553);
46      for (const Particle& ups : upsilons) {
47        LorentzTransform cms_boost;
48        if (ups.p3().mod() > 1*MeV)
49          cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
50	Particles phis;
51        findDecayProducts(ups, phis);
52
53	if(ups.pid()==300553)
54	  _c_4S->fill();
55	else
56	  _c_5S->fill();
57
58	for(const Particle & phi : phis) {
59          FourMomentum p2 = cms_boost.transform(phi.momentum());
60	  double x = 2.*p2.p3().mod()/ups.mass();
61	if(ups.pid()==300553)
62	  _h_4S->fill(x);
63	else
64	  _h_5S->fill(x);
65	}
66      }
67    }
68
69
70    /// Normalise histograms etc., after the run
71    void finalize() {
72      if(_h_4S->effNumEntries()!=0)
73	scale(_h_4S   ,100./ *_c_4S);
74      if(_h_5S->effNumEntries()!=0)
75	scale(_h_5S   ,100./ *_c_5S);
76    }
77
78    /// @}
79
80
81    /// @name Histograms
82    /// @{
83    Histo1DPtr _h_4S,_h_5S;
84    CounterPtr _c_4S,_c_5S;
85    /// @}
86
87
88  };
89
90
91  RIVET_DECLARE_PLUGIN(CLEO_2007_I728872);
92
93
94}