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.37, 2.545, 2.695, 2.85,3.1,3.315,
 40                                       3.46, 3.56, 3.69, 3.865, 4.05, 4.28, 4.565, 4.77, 4.89};
 41      _axes["K0"] = YODA::Axis<double>{1.5, 2.1, 2.4, 2.7, 3.0, 3.3, 3.6, 4.5};
 42    }
 43
 44
 45    /// Perform the per-event analysis
 46    void analyze(const Event& event) {
 47      if (_edges.empty()) {
 48        for (const auto& item : _h) {
 49          _edges[item.first] = item.second->xEdges();
 50        }
 51      }
 52      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 53      int nch = cfs.particles().size();
 54      if (nch<5)  vetoEvent;
 55      // Get beams and average beam momentum
 56      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 57      const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 58      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 59      _wSum->fill();
 60
 61      // neutral kaons
 62      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 63      for (const Particle & p : ufs.particles(Cuts::pid==130 || Cuts::pid==310)) {
 64        const double xi = -log(p.p3().mod()/meanBeamMom);
 65        fillhist("K0", "K0", xi);
 66      }
 67      // charged particles
 68      for (const Particle& p : cfs.particles()) {
 69        double xi = -log(p.p3().mod()/meanBeamMom);
 70        fillhist("charged", "charged", xi);
 71        int id = abs(p.pid());
 72        if (id==211) {
 73          fillhist("pi", "pi", xi);
 74        }
 75        else if (id==321) {
 76          fillhist("Kp", "pi", xi);
 77        }
 78        else if (id==2212) {
 79          fillhist("proton", "pi", xi);
 80        }
 81      }
 82    }
 83
 84    void fillhist(const string& label_hist, const string& label_axis, const double value) {
 85      string edge = "OTHER";
 86      size_t idx = _axes[label_axis].index(value);
 87      // binning for pi, Kp and proton same but different masked bins (pion only 1 masked bin 9)
 88      if ( label_axis=="pi" ) {
 89        if (idx==9) idx=0;
 90        else if (idx>9) --idx;
 91      }
 92      if (idx && idx <= _edges[label_axis].size())  edge = _edges[label_axis][idx-1];
 93      _h[label_hist]->fill(edge);
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      scale(_h, 1./ *_wSum);
100      for( auto & hist : _h) {
101        for(auto & b: hist.second->bins()) {
102          size_t idx = b.index();
103          if (hist.first=="pi" || _axes.find(hist.first)==_axes.end()) {
104            auto it = std::find(_edges["pi"].begin(), _edges["pi"].end(), _edges[hist.first][idx-1]);
105            idx = std::distance(_edges["pi"].begin(), it)+1;
106            if(idx>=9) ++idx;
107            b.scaleW(1./_axes["pi"].width(idx));
108          }
109          else
110            b.scaleW(1./_axes[hist.first].width(idx));
111        }
112      }
113    }
114
115    /// @}
116
117
118    /// @name Histograms
119    /// @{
120    CounterPtr _wSum;
121    map<string, BinnedHistoPtr<string>> _h;
122    map<string, YODA::Axis<double>> _axes;
123    map<string, vector<string>> _edges;
124    /// @}
125
126
127  };
128
129
130  RIVET_DECLARE_PLUGIN(TOPAZ_1995_I381900);
131
132}