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 Strange hadrons at 35GeV
 12  class CELLO_1990_I283026 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CELLO_1990_I283026);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      declare(Beam(), "Beams");
 26      declare(ChargedFinalState(), "FS");
 27      declare(UnstableParticles(), "UFS");
 28
 29      // Book histograms
 30      book(_h["K0"],     1, 1, 1);
 31      book(_h["Kstar"],  2, 1, 1);
 32      book(_h["Lambda"], 3, 1, 1);
 33
 34      _axes["K0"] = YODA::Axis<double>({0.0315, 0.0365, 0.0535, 0.0745, 0.0955, 0.1165,
 35                                        0.1415, 0.175, 0.2375, 0.282, 0.3775, 0.652});
 36      _axes["Kstar"] = YODA::Axis<double>({0.055, 0.195, 0.335, 0.537, 0.883});
 37      _axes["Lambda"] = YODA::Axis<double>({0.067, 0.085, 0.103, 0.131, 0.175, 0.269, 0.447, 0.723});
 38
 39    }
 40
 41
 42    /// Perform the per-event analysis
 43    void analyze(const Event& event) {
 44
 45      if (_edges.empty()) {
 46        for (const auto& item : _h) {
 47          _edges[item.first] = item.second->xEdges();
 48        }
 49      }
 50
 51      // First, veto on leptonic events by requiring at least 4 charged FS particles
 52      const FinalState& fs = apply<FinalState>(event, "FS");
 53      const size_t numParticles = fs.particles().size();
 54
 55      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 56      if (numParticles < 2) {
 57        MSG_DEBUG("Failed leptonic event cut");
 58        vetoEvent;
 59      }
 60      MSG_DEBUG("Passed leptonic event cut");
 61
 62      // Get beams and average beam momentum
 63      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 64      const double meanBeamMom = 0.5*(beams.first.p3().mod() + beams.second.p3().mod());
 65      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 66
 67      // Final state of unstable particles to get particle spectra
 68      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 69
 70      for (const Particle& p : ufs.particles()) {
 71        const int id = p.abspid();
 72        if (id == PID::K0S || id == PID::K0L) {
 73          discfill("K0", p.E()/meanBeamMom);
 74        }
 75        else if(abs(id)==323) {
 76          discfill("Kstar", p.E()/meanBeamMom);
 77        }
 78        else if(abs(id)==3122) {
 79          discfill("Lambda", p.E()/meanBeamMom);
 80        }
 81      }
 82    }
 83
 84    void discfill(const string& name, const double value) {
 85      string edge = "OTHER";
 86      const size_t idx = _axes[name].index(value);
 87      if (idx && idx <= _edges[name].size())  edge = _edges[name][idx-1];
 88      _h[name]->fill(edge);
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      scale(_h, 1./sumOfWeights());
 95      for( auto & hist : _h) {
 96        for(auto & b: hist.second->bins()) {
 97          const size_t idx = b.index();
 98          b.scaleW(1./_axes[hist.first].width(idx));
 99        }
100      }
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    map<string, BinnedHistoPtr<string>> _h;
109    map<string, YODA::Axis<double>> _axes;
110    map<string, vector<string>> _edges;
111    /// @}
112
113
114  };
115
116
117  RIVET_DECLARE_PLUGIN(CELLO_1990_I283026);
118
119
120}