1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief kaon spetra at Upsilon(4S) and nearby continuum
  9  class CLEO_1982_I180612 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1982_I180612);
 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      for (unsigned int ix=0; ix<2; ++ix) {
 25        book(_c[ix],"TMP/c_"+toString(ix+1));
 26        book(_h_mult[ix], 1, 1, ix+1);
 27        for (unsigned int iy=0; iy<2; ++iy) {
 28          book(_h[ix][iy], 2, 1+iy, 1+ix);
 29        }
 30      }
 31    }
 32
 33    /// Recursively walk the decay tree to find decay products of @a p
 34    void findDecayProducts(Particle mother, Particles& unstable) {
 35      for (const Particle & p: mother.children()) {
 36        const int id = abs(p.pid());
 37        if (id == 310 || id == 130 || id == 321) {
 38          unstable.push_back(p);
 39        }
 40        else if(!p.children().empty()) {
 41          findDecayProducts(p, unstable);
 42        }
 43      }
 44    }
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 49      Particles upsilons = ufs.particles(Cuts::pid==300553);
 50      // continuum
 51      if (upsilons.empty()) {
 52        _c[0]->fill();
 53        for (const Particle& p : ufs.particles(Cuts::pid==130 ||
 54                                               Cuts::pid==310 ||
 55                                               Cuts::abspid==321)) {
 56          int id = p.abspid();
 57          double modp = p.p3().mod();
 58          if (id==321) {
 59            _h[0][0]->fill(modp);
 60            _h_mult[0]->fill(_edges[0]);
 61          }
 62          else {
 63            _h[0][1]->fill(modp);
 64            _h_mult[0]->fill(_edges[1]);
 65          }
 66        }
 67      }
 68      /// found an upsilon
 69      else {
 70        for (const Particle& ups : upsilons) {
 71          _c[1]->fill();
 72          Particles unstable;
 73          // Find the decay products we want
 74          findDecayProducts(ups,unstable);
 75          LorentzTransform cms_boost;
 76          if (ups.p3().mod() > 0.001) {
 77            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 78          }
 79          for( const Particle & p : unstable) {
 80            int id = p.abspid();
 81            double modp = cms_boost.transform(p.momentum()).p3().mod();
 82            if (id==321) {
 83              _h[1][0]->fill(modp);
 84              _h_mult[1]->fill(_edges[0]);
 85            }
 86            else {
 87              _h[1][1]->fill(modp);
 88              _h_mult[1]->fill(_edges[1]);
 89            }
 90          }
 91        }
 92      }
 93    }
 94
 95
 96    /// Normalise histograms etc., after the run
 97    void finalize() {
 98      for (unsigned int ix=0;ix<2;++ix) {
 99        scale(_h_mult[ix], 1./ *_c[ix]);
100        for (unsigned int iy=0;iy<2;++iy) {
101          scale(_h[ix][iy],crossSection()/nanobarn/sumOfWeights());
102        }
103      }
104    }
105
106    /// @}
107
108
109    /// @name Histograms
110    /// @{
111    Histo1DPtr _h[2][2];
112    BinnedHistoPtr<string> _h_mult[2];
113    CounterPtr _c[2];
114    vector<string> _edges = { "E+ E- --> (K+ + K-) X", "E+ E- --> K0 X" };
115    /// @}
116
117
118  };
119
120
121  RIVET_DECLARE_PLUGIN(CLEO_1982_I180612);
122
123}