1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Delta++ rate and spectrum
  9  class ARGUS_1989_I278932 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1989_I278932);
 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      // histograms
 24      book(_r_ups, 1,1,2);
 25      book(_r_cont,2,1,2);
 26      book(_h_p,3,1,1);
 27      book(_w_ups ,"TMP/w_ups" );
 28      book(_w_cont,"TMP/w_cont");
 29    }
 30
 31
 32    /// Recursively walk the decay tree to find decay products of @a p
 33    void findDecayProducts(Particle mother, Particles& delta) {
 34      // dleta pdg code
 35      static const long id = 2224;
 36      for(const Particle & p: mother.children()) {
 37        if(p.abspid() == id) {
 38          delta.push_back(p);
 39        }
 40        else if(!p.children().empty()) {
 41          findDecayProducts(p, delta);
 42        }
 43      }
 44    }
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      UnstableParticles ufs = apply<UnstableParticles>(event, "UFS");
 49      Particles ups = ufs.particles(Cuts::pid==553);
 50      // continuum
 51      if (ups.empty()) {
 52        _w_cont->fill();
 53        _r_cont->fill(Ecm2, ufs.particles(Cuts::abspid==2224).size());
 54      }
 55      // upsilon decays
 56      else {
 57        for (const Particle &p : ups) {
 58          _w_ups->fill();
 59          Particles delta;
 60          findDecayProducts(p,delta);
 61          if (delta.empty())  continue;
 62          LorentzTransform boost;
 63          if (p.p3().mod() > 1*MeV) {
 64            boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 65          }
 66          for (const Particle& del : delta) {
 67            _r_ups->fill(Ecm1);
 68            double mom = boost.transform(del.momentum()).p3().mod();
 69            _h_p->fill(mom);
 70          }
 71        }
 72      }
 73    }
 74
 75
 76    /// Normalise histograms etc., after the run
 77    void finalize() {
 78      if (_w_cont->effNumEntries()>0) {
 79        scale(_r_cont, 1./ *_w_cont);
 80      }
 81      // direct upsilon decays, i.e. to gg gamma and ggg so need to renormalize
 82      // factor and values from arxiv:0612019
 83      if (_w_ups->effNumEntries()>0) {
 84        const double Bmumu=0.0249, rhad = 3.56;
 85        const double fact = 1.-(3.+rhad)*Bmumu;
 86        scale(_r_ups, 1./fact / *_w_ups);
 87        scale(_h_p  , 100./fact / *_w_ups);
 88      }
 89    }
 90
 91    /// @}
 92
 93
 94    /// @name Histograms
 95    /// @{
 96    BinnedHistoPtr<string> _r_ups, _r_cont;
 97    Histo1DPtr _h_p;
 98    CounterPtr _w_ups,_w_cont;
 99    const string Ecm1 = "9.46", Ecm2 = "10.2";
100    /// @}
101
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(ARGUS_1989_I278932);
107
108}