1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Lambda(1520) production
  9  class ARGUS_1989_I262415 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1989_I262415);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(), "UFS");
 24
 25      // Book histograms
 26      book(_h_ups1_obs, 1, 1, 1);
 27      book(_h_cont_obs, 1, 1, 2);
 28      book(_h_ups1_all, 1, 2, 1);
 29      book(_h_cont_all, 1, 2, 2);
 30      book(_h_ups1    , 3, 1, 1);
 31      book(_h_cont    , 4, 1, 1);
 32      book(_weightSum_cont,"TMP/weightSum_cont");
 33      book(_weightSum_Ups1,"TMP/weightSum_Ups1");
 34
 35    }
 36
 37    /// Recursively walk the decay tree to find decay products of @a p
 38    void findDecayProducts(Particle mother, Particles& unstable) {
 39      for(const Particle & p: mother.children()) {
 40        const int id = p.pid();
 41	if(id == 3124) {
 42	  unstable.push_back(p);
 43	}
 44	if(!p.children().empty())
 45	  findDecayProducts(p, unstable);
 46      }
 47    }
 48
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      // Find the Upsilons among the unstables
 53      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 54      Particles upsilons = ufs.particles(Cuts::pid==553);
 55      // Continuum
 56      if (upsilons.empty()) {
 57        MSG_DEBUG("No Upsilons found => continuum event");
 58        _weightSum_cont->fill();
 59        for (const Particle& p : ufs.particles(Cuts::abspid==3124)) {
 60          const double xp = 2.*p.E()/sqrtS();
 61          const double beta = p.p3().mod() / p.E();
 62	  _h_cont->fill(xp,1./beta);
 63	  _h_cont_obs->fill(xp);
 64	  _h_cont_all->fill(xp);
 65	}
 66      }
 67      // Upsilon(s) found
 68      else {
 69        MSG_DEBUG("Upsilons found => resonance event");
 70        for (const Particle& ups : upsilons) {
 71	  _weightSum_Ups1->fill();
 72          Particles unstable;
 73          // Find the decay products we want
 74          findDecayProducts(ups, unstable);
 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 : unstable) {
 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	    _h_ups1->fill(xp,1./beta);
 85	    _h_ups1_obs->fill(xp);
 86	    _h_ups1_all->fill(xp);
 87	  }
 88	}
 89      }
 90    }
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      if (_weightSum_cont->val() > 0.) {
 95	scale(_h_cont    , 1. / *_weightSum_cont);
 96	scale(_h_cont_obs, 0.3/ *_weightSum_cont);
 97	scale(_h_cont_all, 1. / *_weightSum_cont);
 98      }
 99      if (_weightSum_Ups1->val() > 0.) {
100	scale(_h_ups1    , 1. / *_weightSum_Ups1);
101	scale(_h_ups1_obs, 0.3/ *_weightSum_Ups1);
102	scale(_h_ups1_all, 1. / *_weightSum_Ups1);
103      }
104
105    }
106
107    /// @}
108
109
110    /// @name Histograms
111    /// @{
112    Histo1DPtr _h_ups1, _h_cont;
113    Histo1DPtr _h_ups1_obs, _h_cont_obs;
114    Histo1DPtr _h_ups1_all, _h_cont_all;
115    CounterPtr _weightSum_cont,_weightSum_Ups1;
116    /// @}
117
118
119  };
120
121
122  RIVET_DECLARE_PLUGIN(ARGUS_1989_I262415);
123
124
125}