1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6#include "Rivet/Projections/Beam.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief kaon spectra at LEP1
 12  class DELPHI_1995_I377487 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1995_I377487);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24      declare(Beam(), "Beams");
 25      declare(ChargedFinalState(), "FS");
 26      declare(UnstableParticles(), "UFS");
 27
 28      // Book histograms
 29      book(_h_K0_x , 8, 1, 1);
 30      book(_h_K0_xi, 9, 1, 1);
 31      book(_h_Ks_x ,10, 1, 1);
 32
 33    }
 34
 35
 36    /// Perform the per-event analysis
 37    void analyze(const Event& event) {
 38
 39      if (_edges.empty())  _edges = _h_K0_x->xEdges();
 40
 41      // First, veto on leptonic events by requiring at least 4 charged FS particles
 42      const FinalState& fs = apply<FinalState>(event, "FS");
 43      const size_t numParticles = fs.particles().size();
 44
 45      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 46      if (numParticles < 2) {
 47        MSG_DEBUG("Failed leptonic event cut");
 48        vetoEvent;
 49      }
 50      MSG_DEBUG("Passed leptonic event cut");
 51
 52      // Get beams and average beam momentum
 53      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 54      const double meanBeamMom = ( beams.first.p3().mod() +
 55                                   beams.second.p3().mod() ) / 2.0;
 56      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 57
 58      // Final state of unstable particles to get particle spectra
 59      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 60
 61      for (const Particle& p : ufs.particles(Cuts::pid==130 or Cuts::pid==310 or Cuts::abspid==323)) {
 62        const double xp = p.p3().mod()/meanBeamMom;
 63        if (abs(p.pid())==323) {
 64          _h_Ks_x->fill(xp);
 65        }
 66        else {
 67          _h_K0_x->fill(map2string(xp));
 68          _h_K0_xi->fill(-log(xp));
 69        }
 70      }
 71    }
 72
 73    string map2string(const double val) const {
 74      const size_t idx = _axis.index(val);
 75      if (idx || idx <= _edges.size())  return _edges[idx-1];
 76      return "OTHER";
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82
 83      scale(_h_K0_x,  1./sumOfWeights());
 84      for(auto & b: _h_K0_x->bins()) {
 85        const size_t idx = b.index();
 86        b.scaleW(1./_axis.width(idx));
 87      }
 88      scale(_h_K0_xi, 1./sumOfWeights());
 89      scale(_h_Ks_x,  1./sumOfWeights());
 90    }
 91
 92    /// @}
 93
 94
 95    /// @name Histograms
 96    /// @{
 97    BinnedHistoPtr<string> _h_K0_x;
 98    Histo1DPtr _h_K0_xi,_h_Ks_x;
 99    vector<string> _edges;
100    YODA::Axis<double> _axis{ 0.00030213628351706805, 0.00048341905738733416, 0.00048343557923922857,
101                              0.00026205300091453314, 0.0007702529509799692, 0.0009481642553664134,
102                              0.0017226600969315643, 0.002004423179522298, 0.0016843611112658217,
103                              0.0026292281438344, 0.00267627755270744, 0.0036267300396739185,
104                              0.004237796021633787, 0.005212931632136056, 0.007189937374782032,
105                              0.007726421785666127, 0.010282046710587495, 0.012196841637666128,
106                              0.014664716763387292, 0.018701111778413465, 0.02210348200534462,
107                              0.027403036058393532, 0.03380578808779788, 0.04112055882855764,
108                              0.049671035882778436, 0.22518836390597363, 1.0 };
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(DELPHI_1995_I377487);
116
117
118}