1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Lambda in Upsilon(1S) decay and nearby continuum
  9  class MD1_1994_I362599 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(MD1_1994_I362599);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(), "UFS");
 23      // histos
 24      book(_weightSum_cont,"TMP/weightSumcont");
 25      book(_weightSum_Ups1,"TMP/weightSumUps1");
 26      book(_h_spect,1,1,1);
 27      for (unsigned int ix=0; ix<2; ++ix) {
 28        for (unsigned int iy=0; iy<2; ++iy) {
 29          book(_mult[ix][iy],"/TMP/MULT_" +toString(ix) + "_" +toString(iy));
 30        }
 31      }
 32    }
 33
 34    /// Recursively walk the decay tree to find decay products of @a p
 35    void findDecayProducts(Particle mother, Particles& unstable) {
 36      for (const Particle & p: mother.children()) {
 37        const int id = p.abspid();
 38        if (id==PID::LAMBDA   || id==PID::XIMINUS) unstable.push_back(p);
 39        if (!p.children().empty()) findDecayProducts(p, unstable);
 40      }
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      // Find the upsilons
 47      // First in unstable final state
 48      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 49      Particles upsilons = ufs.particles(Cuts::pid==553);
 50      // continuum
 51      if (upsilons.empty()) {
 52        _weightSum_cont->fill();
 53        // Unstable particles
 54        size_t nLam = ufs.particles(Cuts::abspid==PID::LAMBDA).size();
 55        _mult[1][0]->fill(nLam);
 56      }
 57      else {
 58        for (const Particle& ups : upsilons) {
 59          _weightSum_Ups1->fill();
 60          Particles unstable;
 61          LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 62          // Find the decay products we want
 63          findDecayProducts(ups,unstable);
 64          for(const Particle & p : unstable)  {
 65            int id = p.abspid();
 66            if(id==PID::LAMBDA) {
 67              double xp = 2.*boost.transform(p.momentum()).p3().mod()/ups.mass();
 68              _h_spect->fill(xp);
 69              _mult[0][0]->fill();
 70            }
 71            else if(id==PID::XIMINUS) {
 72              _mult[0][1]->fill();
 73            }
 74          }
 75        }
 76      }
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82      if (_weightSum_Ups1->val() > 0.) {
 83        scale(_h_spect, 1./ *_weightSum_Ups1);
 84        for (size_t iy=0; iy<2; ++iy) {
 85          BinnedEstimatePtr<string> est;
 86          book(est, 2+iy*2, 1, 1);
 87          if (_weightSum_Ups1->val() > 0.) {
 88            scale(_mult[0][iy],1./ *_weightSum_Ups1);
 89            est->bin(1).set(_mult[0][iy]->val(), _mult[0][iy]->err());
 90          }
 91        }
 92      }
 93      if (_weightSum_cont->val() > 0.) {
 94        scale(_mult[1][0], 1./ *_weightSum_cont);
 95        Estimate1DPtr est;
 96        book(est, 3, 1, 1);
 97        est->bin(1).set(_mult[1][0]->val(), _mult[1][0]->err());
 98      }
 99    }
100
101    /// @}
102
103
104    /// @name Histograms
105    /// @{
106    Histo1DPtr _h_spect;
107    CounterPtr _weightSum_cont,_weightSum_Ups1;
108    CounterPtr _mult[2][2];
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(MD1_1994_I362599);
116
117}