1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief pi, K and proton spectra at 12 and 30 GeV
 10  class TASSO_1980_I153656 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1980_I153656);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(ChargedFinalState(), "FS");
 26
 27      // Book histograms
 28      _iHist=-1;
 29      sqs = 1.;
 30      if (isCompatibleWithSqrtS(12*GeV)) {
 31        _iHist = 0;
 32        sqs = 12.;
 33      }
 34      else if (isCompatibleWithSqrtS(30*GeV)) {
 35        _iHist = 1;
 36        sqs = 30.;
 37      }
 38      else
 39        MSG_ERROR("Beam energy " << sqrtS() << " GeV not supported!");
 40
 41      book(_h["p_pi"],3*_iHist+2,1,1);
 42      book(_h["p_K"] ,3*_iHist+3,1,1);
 43      book(_h["p_p"] ,3*_iHist+4,1,1);
 44
 45      tribook("pi", 3*_iHist+ 8, 1, 1);
 46      tribook("K",  3*_iHist+ 9, 1, 1);
 47      tribook("p",  3*_iHist+10, 1, 1);
 48
 49      if (_iHist) {
 50        _axes["pi"] = YODA::Axis<double>({0.325, 0.375, 0.425, 0.5, 0.575, 0.7, 0.9, 1.1, 1.3, 1.5});
 51        _axes["K"]  = YODA::Axis<double>({0.4, 0.5, 0.575, 0.7, 0.9, 1.1});
 52        _axes["p"]  = YODA::Axis<double>({0.5, 0.7, 1.2325, 2.0975});
 53        _axes["r"]  = YODA::Axis<double>({0.4, 0.5, 0.575, 0.7, 0.9, 1.1});
 54      }
 55      else {
 56        _axes["pi"] = YODA::Axis<double>({0.3, 0.4, 0.5, 0.6, 1.0, 1.6});
 57        _axes["K"]  = YODA::Axis<double>({0.4, 0.5, 0.6, 1.0});
 58        _axes["p"]  = YODA::Axis<double>({0.5, 0.7, 0.9, 2.43});
 59        _axes["r"]  = YODA::Axis<double>({0.4, 0.5, 0.6, 1.0});
 60      }
 61      _axes["rp"]  = YODA::Axis<double>({0.475, 0.725, 1.0, 2.2});
 62    }
 63
 64    void tribook(const string& label, unsigned int d, unsigned int x, unsigned int y) {
 65      book(_h["n_"+label], "TMP/n_"+label, refData<YODA::BinnedEstimate<string>>(d, x, y));
 66      book(_h["d_"+label], "TMP/d_"+label, refData<YODA::BinnedEstimate<string>>(d, x, y));
 67      book(_r[label], d, x, y);
 68    }
 69
 70
 71    /// Perform the per-event analysis
 72    void analyze(const Event& event) {
 73      if (_edges.empty()) {
 74        for (const auto& item : _h) {
 75          _edges[item.first] = item.second->xEdges();
 76        }
 77        for (const auto& item : _r) {
 78          _edges["r"+item.first] = _h["n_"+item.first]->xEdges();
 79        }
 80      }
 81      // First, veto on leptonic events by requiring at least 4 charged FS particles
 82      const ChargedFinalState& fs = apply<ChargedFinalState>(event, "FS");
 83      const size_t numParticles = fs.particles().size();
 84
 85      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 86      if (numParticles < 2) {
 87        MSG_DEBUG("Failed leptonic event cut");
 88        vetoEvent;
 89      }
 90      MSG_DEBUG("Passed leptonic event cut");
 91
 92      // Get beams and average beam momentum
 93      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 94      const double meanBeamMom = ( beams.first.p3().mod() +
 95                                   beams.second.p3().mod() ) / 2.0;
 96      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 97
 98      for (const Particle& p : fs.particles()) {
 99        double modp = p.p3().mod();
100        fillND("d_pi", modp);
101        fillND("d_K", modp);
102        fillND("d_p", modp);
103        if (abs(p.pid())==211) {
104          fillhist("p_pi", modp);
105          fillND("n_pi", modp);
106        }
107        else if (abs(p.pid())==321) {
108          fillhist("p_K", modp);
109          fillND("n_K", modp);
110        }
111        else if (abs(p.pid())==2212) {
112          fillhist("p_p", modp);
113          fillND("n_p", modp);
114        }
115      }
116    }
117
118    void fillhist(const string& label, const double value, const double weight = 1.0) {
119      string edge = "OTHER";
120      const string tag = label.substr(2);
121      const size_t idx = _axes[tag].index(value);
122      if (idx && idx <= _edges[label].size())  edge = _edges[label][idx-1];
123      _h[label]->fill(edge, weight);
124    }
125
126    void fillND(const string& label, const double value) {
127      string edge = "OTHER";
128      const string tag = label.substr(2);
129      const size_t idx = _axes[(tag == "p")? "rp" : "r"].index(value);
130      if (idx && idx <= _edges["r"+tag].size())  edge = _edges["r"+tag][idx-1];
131      _h[label]->fill(edge);
132    }
133
134    /// Normalise histograms etc., after the run
135    void finalize() {
136
137      scale(_h["p_pi"], crossSection()/nanobarn/sumOfWeights());
138      scale(_h["p_K"],  crossSection()/nanobarn/sumOfWeights());
139      scale(_h["p_p"],  crossSection()/nanobarn/sumOfWeights());
140      vector<string> s1={"pi","K","p"};
141      for ( auto & key : s1 ) {
142        for(auto & b: _h["p_"+key]->bins()) {
143          const size_t idx = b.index();
144          b.scaleW(1./_axes[key].width(idx));
145        }
146      }
147      for (auto& item : _r) {
148        divide(_h["n_"+item.first], _h["d_"+item.first], item.second);
149      }
150
151    }
152
153    /// @}
154
155
156    /// @name Histograms
157    /// @{
158    map<string,BinnedHistoPtr<string>> _h;
159    map<string,BinnedEstimatePtr<string>> _r;
160    map<string, YODA::Axis<double>> _axes;
161    map<string, vector<string>> _edges;
162    int _iHist;
163    double sqs;
164    /// @}
165
166
167  };
168
169
170  RIVET_DECLARE_PLUGIN(TASSO_1980_I153656);
171
172
173}