1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief pi0 and eta at Upsilon 1 and continuum
  9  class CRYSTAL_BALL_1991_I297905 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CRYSTAL_BALL_1991_I297905);
 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_cont_pi[0] ,3, 1, 1);
 25      book(_h_cont_pi[1] ,3, 1, 2);
 26      book(_h_ups1_pi    ,4, 1, 1);
 27      book(_h_cont_eta[0],5, 1, 1);
 28      book(_h_cont_eta[1],5, 1, 2);
 29      book(_h_ups1_eta ,6, 1, 1);
 30      // counters
 31      book(_n_Pi[0] , 1, 1, 1);
 32      book(_n_Pi[1] , 1, 1, 2);
 33      book(_n_Eta[0], 2, 1, 1);
 34      book(_n_Eta[1], 2, 1, 2);
 35      book(_weightSum_cont,"/TMP/weightSum_cont");
 36      book(_weightSum_Ups1,"/TMP/weightSum_Ups1");
 37    }
 38
 39    /// Recursively walk the decay tree to find decay products of @a p
 40    void findDecayProducts(Particle mother, Particles& unstable) {
 41      for (const Particle & p: mother.children()) {
 42        const int id = p.pid();
 43        if (id == 111 or id == 221) {
 44          unstable.push_back(p);
 45        }
 46        if (!p.children().empty()) {
 47          findDecayProducts(p, unstable);
 48        }
 49      }
 50    }
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      // Find the Upsilons among the unstables
 55      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 56      Particles upsilons = ufs.particles(Cuts::pid==553);
 57      // Continuum
 58      if (upsilons.empty()) {
 59        MSG_DEBUG("No Upsilons found => continuum event");
 60        _weightSum_cont->fill();
 61        for (const Particle& p : ufs.particles(Cuts::pid==111 or Cuts::pid==221)) {
 62          const int id = p.pid();
 63          const double xp = 2.*p.E()/sqrtS();
 64          const double beta = p.p3().mod() / p.E();
 65          if (id==111) {
 66            _n_Pi[0]->fill(sqrtS()/GeV);
 67            for (unsigned int ix=0; ix<2; ++ix) {
 68              _h_cont_pi[ix]->fill(xp,1./beta);
 69            }
 70          }
 71          else {
 72            _n_Eta[0]->fill(Ecm);
 73            for (unsigned int ix=0; ix<2; ++ix) {
 74              _h_cont_eta[ix]->fill(xp,1./beta);
 75            }
 76          }
 77        }
 78      }
 79      // Upsilon(s) found
 80      else {
 81        MSG_DEBUG("Upsilons found => resonance event");
 82        for (const Particle& ups : upsilons) {
 83          _weightSum_Ups1->fill();
 84          Particles unstable;
 85          // Find the decay products we want
 86          findDecayProducts(ups, unstable);
 87          LorentzTransform cms_boost;
 88          if (ups.p3().mod() > 1*MeV)
 89            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 90            const double mass = ups.mass();
 91            // loop over decay products
 92            for (const Particle& p : unstable) {
 93              const int id = p.pid();
 94              const FourMomentum p2 = cms_boost.transform(p.momentum());
 95              const double xp = 2.*p2.E()/mass;
 96              const double beta = p2.p3().mod()/p2.E();
 97              if (id==111) {
 98                _n_Pi[1]->fill(sqrtS()/GeV);
 99                _h_ups1_pi->fill(xp,1./beta);
100              }
101              else if (id==221) {
102                _n_Eta[1]->fill(Ecm);
103                _h_ups1_eta->fill(xp,1./beta);
104            }
105          }
106        }
107      }
108    }
109
110
111    /// Normalise histograms etc., after the run
112    void finalize() {
113      // Scale histos
114      if (_weightSum_cont->val() > 0.) {
115        scale(_h_cont_pi[0] , 1./ *_weightSum_cont);
116        scale(_h_cont_pi[1] , sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
117        scale(_h_cont_eta[0], 1./ *_weightSum_cont);
118        scale(_h_cont_eta[1], sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
119      	scale(_n_Pi[0],  1./ *_weightSum_cont);
120      	scale(_n_Eta[0], 1./ *_weightSum_cont);
121      }
122      if (_weightSum_Ups1->val() > 0.) {
123        scale(_h_ups1_pi, 1./ *_weightSum_Ups1);
124        scale(_h_ups1_eta, 1./ *_weightSum_Ups1);
125      	scale(_n_Pi[1],  1./ *_weightSum_Ups1);
126      	scale(_n_Eta[1], 1./ *_weightSum_Ups1);
127      }
128    }
129
130    /// @}
131
132
133    /// @name Histograms
134    /// @{
135    Histo1DPtr _h_cont_pi[2] , _h_ups1_pi;
136    Histo1DPtr _h_cont_eta[2], _h_ups1_eta;
137    BinnedHistoPtr<string> _n_Eta[2];
138    CounterPtr _n_Pi[2];
139    CounterPtr _weightSum_cont,_weightSum_Ups1;
140    const string Ecm = "9.8";
141    /// @}
142
143  };
144
145
146  RIVET_DECLARE_PLUGIN(CRYSTAL_BALL_1991_I297905);
147
148}