1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief Photon spectrum in Upsilon(1,2S) decays
 9  class CUSB_1984_I204307 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(CUSB_1984_I204307);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // projections
22      declare(UnstableParticles(Cuts::pid==553 || Cuts::pid==100553), "UFS");
23      // histos
24      for (unsigned int ix=0; ix<2; ++ix) {
25	      book(_h[ix],1+ix,1,1);
26      }
27    }
28
29
30    /// Perform the per-event analysis
31    void analyze(const Event& event) {
32      // Find the Upsilons among the unstables
33      for (const Particle& ups : apply<UnstableParticles>(event, "UFS").particles()) {
34        unsigned int iups = ups.pid()/100000;
35      	unsigned int nhadron(0);
36      	Particles photons;
37      	for (const Particle& child : ups.children()) {
38      	  if (PID::isHadron(child.pid())) {
39      	    ++nhadron;
40          }
41      	  else if(child.pid()==PID::PHOTON) {
42      	    photons.push_back(child);
43          }
44      	}
45        if (nhadron==0 || photons.empty()) continue;
46        LorentzTransform boost;
47        if (ups.p3().mod() > 1*MeV) {
48          boost = LorentzTransform::mkFrameTransformFromBeta(ups.mom().betaVec());
49        }
50        for (const Particle& gamma : photons) {
51          const double E = boost.transform(gamma.mom()).E();
52          // E = randnorm(E,0.21*sqrt(E));
53          _h[iups]->fill(2.*E/ups.mass());
54        }
55      }
56    }
57
58
59    /// Normalise histograms etc., after the run
60    void finalize() {
61      normalize(_h, 1.0, false);
62    }
63
64    /// @}
65
66
67    /// @name Histograms
68    /// @{
69    Histo1DPtr _h[2];
70    /// @}
71
72
73  };
74
75
76  RIVET_DECLARE_PLUGIN(CUSB_1984_I204307);
77
78}