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