1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief Add a short analysis description here
 11  class TOPAZ_1995_I381900 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(TOPAZ_1995_I381900);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(Beam(), "Beams");
 26      declare(UnstableParticles(), "UFS");
 27      const ChargedFinalState cfs;
 28      declare(cfs, "CFS");
 29
 30      // Book histograms
 31      book(_h["charged"], 1, 1, 1);
 32      book(_h["pi"]     , 2, 1, 1);
 33      book(_h["Kp"]     , 2, 1, 2);
 34      book(_h["proton"] , 2, 1, 3);
 35      book(_h["K0"]     , 3, 1, 1);
 36      book(_wSum,"TMP/wSum");
 37
 38      _axes["charged"] = YODA::Axis<double>{22, 0.6, 5.0};
 39      _axes["pi"] = YODA::Axis<double>{1.17, 1.47, 1.77, 2.07, 2.345, 2.545, 2.695, 2.875, 3.195,
 40                                       3.46, 3.56, 3.69, 3.865, 4.05, 4.28, 4.565, 4.77, 4.89};
 41      _axes["Kp"] = YODA::Axis<double>{1.17, 1.47, 1.77, 2.27, 2.695, 3.09, 3.46, 3.56, 3.69, 3.865, 4.05, 4.28, 4.56};
 42      _axes["proton"] = YODA::Axis<double>{1.17, 1.47, 1.77, 2.27, 2.695, 3.09, 3.46, 3.56, 3.69, 3.85};
 43      _axes["K0"] = YODA::Axis<double>{1.575, 2.025, 2.4, 2.7, 3.0, 3.3, 3.75, 4.35};
 44    }
 45
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      if (_edges.empty()) {
 50        for (const auto& item : _h) {
 51          _edges[item.first] = item.second->xEdges();
 52        }
 53      }
 54      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 55      int nch = cfs.particles().size();
 56      if (nch<5)  vetoEvent;
 57      // Get beams and average beam momentum
 58      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 59      const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 60      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 61      _wSum->fill();
 62
 63      // neutral kaons
 64      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 65      for (const Particle & p : ufs.particles(Cuts::pid==130 || Cuts::pid==310)) {
 66        const double xi = -log(p.p3().mod()/meanBeamMom);
 67        fillhist("K0", xi);
 68      }
 69      // charged particles
 70      for (const Particle& p : cfs.particles()) {
 71        double xi = -log(p.p3().mod()/meanBeamMom);
 72        fillhist("charged", xi);
 73        int id = abs(p.pid());
 74        if (id==211) {
 75          fillhist("pi", xi);
 76        }
 77        else if (id==321) {
 78          fillhist("Kp", xi);
 79        }
 80        else if (id==2212) {
 81          fillhist("proton", xi);
 82        }
 83      }
 84    }
 85
 86    void fillhist(const string& label, const double value) {
 87      string edge = "OTHER";
 88      const size_t idx = _axes[label].index(value);
 89      if (idx && idx <= _edges[label].size())  edge = _edges[label][idx-1];
 90      _h[label]->fill(edge);
 91    }
 92
 93
 94    /// Normalise histograms etc., after the run
 95    void finalize() {
 96      scale(_h, 1./ *_wSum);
 97    }
 98
 99    /// @}
100
101
102    /// @name Histograms
103    /// @{
104    CounterPtr _wSum;
105    map<string, BinnedHistoPtr<string>> _h;
106    map<string, YODA::Axis<double>> _axes;
107    map<string, vector<string>> _edges;
108    /// @}
109
110
111  };
112
113
114  RIVET_DECLARE_PLUGIN(TOPAZ_1995_I381900);
115
116}