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], "/TMP/NUps1");
 28      book(_n_Phi[1], "/TMP/NUps2");
 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    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      // Find the Upsilons among the unstables
 49      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 50      Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
 51      // Continuum
 52      if (upsilons.empty()) { 
 53        MSG_DEBUG("No Upsilons found => continuum event");
 54        _weightSum_cont->fill();
 55        for (const Particle& p : ufs.particles(Cuts::pid==333)) {
 56          const double xp = 2.*p.E()/sqrtS();
 57          const double beta = p.p3().mod() / p.E();
 58	  _h_cont->fill(xp,1./beta);
 59	}
 60      }
 61      // Upsilon(s) found
 62      else { 
 63        MSG_DEBUG("Upsilons found => resonance event");
 64        for (const Particle& ups : upsilons) {
 65          const int parentId = ups.pid();
 66	  if(parentId==553) {
 67	    _weightSum_Ups1->fill();
 68	  }
 69	  else {
 70	    _weightSum_Ups2->fill();
 71	  }
 72          Particles phis;
 73          // Find the decay products we want
 74          findDecayProducts(ups, phis);
 75          LorentzTransform cms_boost;
 76          if (ups.p3().mod() > 1*MeV)
 77            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 78          const double mass = ups.mass();
 79	  // loop over decay products
 80          for(const Particle& p : phis) {
 81            const FourMomentum p2 = cms_boost.transform(p.momentum());
 82            const double xp = 2.*p2.E()/mass;
 83            const double beta = p2.p3().mod()/p2.E();
 84	    if(parentId==553) {
 85	      _n_Phi[0]->fill();
 86	      _h_ups1->fill(xp,1./beta);
 87	    }
 88	    else {
 89	      _n_Phi[1]->fill();
 90	      _h_ups2->fill(xp,1./beta);
 91	    }
 92	  }
 93	}
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      if (_weightSum_cont->effNumEntries() > 0.)
101	scale(_h_cont, sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
102      if (_weightSum_Ups1->effNumEntries() > 0.) {
103	scale(_h_ups1, 1./ *_weightSum_Ups1);
104      }
105      if (_weightSum_Ups2->effNumEntries() > 0.) {
106	scale(_h_ups2, 1./ *_weightSum_Ups2);
107      }
108      // Counters
109      vector<CounterPtr> scales = {_weightSum_Ups1,_weightSum_Ups2};
110      for(unsigned int ix=0;ix<2;++ix) {
111	Scatter2DPtr scatter;
112	book(scatter, 3+ix, 1, 1, true);
113	scale(_n_Phi[ix], 1./ *scales[ix]);
114        scatter->point(0).setY(_n_Phi[ix]->val(),
115			       _n_Phi[ix]->err());
116      }
117    }
118
119    //@}
120
121
122    /// @name Histograms
123    //@{
124    Histo1DPtr _h_cont, _h_ups1, _h_ups2;
125    CounterPtr _n_Phi[2];
126    CounterPtr _weightSum_cont,_weightSum_Ups1,_weightSum_Ups2;
127    //@}
128
129
130  };
131
132
133  // The hook for the plugin system
134  RIVET_DECLARE_PLUGIN(ARGUS_1989_I262551);
135
136
137}