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 Add a short analysis description here
 10  class TPC_1988_I262143 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1988_I262143);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(Beam(), "Beams");
 23      declare(ChargedFinalState(), "FS");
 24
 25      // Book histograms
 26      book(_h["z_pi"],  1, 1, 1);
 27      book(_h["z_K"],   1, 1, 2);
 28      book(_h["z_p"],   1, 1, 3);
 29      book(_h["z_all"], 1, 1, 4);
 30
 31      book(_h["z2_pi"], 5, 1, 1);
 32      book(_h["z2_K"],  5, 1, 2);
 33      book(_h["z2_p"],  5, 1, 3);
 34
 35      tribook("pi", 6, 1, 1);
 36      tribook("K",  6, 1, 2);
 37      tribook("p",  6, 1, 3);
 38      tribook("2_K", 7, 1, 1);
 39      tribook("2_p", 7, 1, 2);
 40      tribook("3_p", 7, 1, 3);
 41    }
 42
 43    void tribook(const string& label, unsigned int d, unsigned int x, unsigned int y) {
 44      book(_nd["n_"+label], "TMP/n" + label, refData<YODA::BinnedEstimate<string>>(d, x, y));
 45      book(_nd["d_"+label], "TMP/d" + label, refData<YODA::BinnedEstimate<string>>(d, x, y));
 46      book(_ratio[label], d, x, y);
 47    }
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51
 52      if (_edges.empty()) {
 53        for (const auto& item : _ratio) {
 54          _edges[item.first] = _nd["n_" + item.first]->xEdges();
 55        }
 56      }
 57      // First, veto on leptonic events by requiring at least 4 charged FS particles
 58      const FinalState& fs = apply<FinalState>(event, "FS");
 59      const size_t numParticles = fs.particles().size();
 60
 61      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 62      if (numParticles < 2) {
 63        MSG_DEBUG("Failed leptonic event cut");
 64        vetoEvent;
 65      }
 66      MSG_DEBUG("Passed leptonic event cut");
 67
 68      // Get beams and average beam momentum
 69      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 70      const double meanBeamMom = (beams.first.p3().mod() + beams.second.p3().mod()) / 2.0;
 71      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 72      for (const Particle& p : fs.particles()) {
 73        double xP = p.p3().mod()/meanBeamMom;
 74        _h["z_all"]->fill(xP);
 75        histfill("d_pi", xP);
 76        histfill("d_K",  xP);
 77        histfill("d_p",  xP);
 78        int id = abs(p.pid());
 79        if (id==211) {
 80          _h["z_pi"]->fill(xP);
 81          _h["z2_pi"]->fill(xP, xP);
 82          histfill("n_pi", xP, 100.);
 83          histfill("d_2_K", xP);
 84          histfill("d_2_p", xP);
 85          histfill("d_3_p", xP);
 86        }
 87        else if (id==321) {
 88          _h["z_K"]->fill(xP);
 89          _h["z2_K"]->fill(xP, xP);
 90          histfill("n_K",  xP, 100.);
 91          histfill("n_2_K", xP);
 92          histfill("d_3_p", xP);
 93        }
 94        else if (id==2212) {
 95          _h["z_p"]->fill(xP);
 96          _h["z2_p"]->fill(xP, xP);
 97          histfill("n_p",  xP, 100.);
 98          histfill("n_2_p", xP);
 99          histfill("n_3_p", xP);
100        }
101      }
102    }
103
104    /// Normalise histograms etc., after the run
105    void finalize() {
106      scale(_h, 1./sumOfWeights());
107      for (auto& item : _ratio) {
108        divide(_nd["n_"+item.first], _nd["d_"+item.first], item.second);
109      }
110    }
111    
112    /// @}
113
114    void histfill(const string& label, const double value, const double weight = 1.0) {
115      string edge = "OTHER";
116      const string tag = label.substr(2);
117      size_t idx = _axis.index(value);
118      if (tag=="pi") {
119        if      (0.035 <= value && value <= 0.07)  idx -=  2;
120        else if (0.09  <= value && value <= 0.11)  idx -=  6;
121        else if (0.16  <= value && value <= 0.7 )  idx -= 10;
122        else idx=0;
123      }
124      if (tag == "K") {
125        if      (0.035 <= value && value <= 0.07)  idx -=  2;
126        else if (0.09  <= value && value <= 0.14)  idx -=  6;
127        else if (0.25  <= value && value <= 0.7 )  idx -= 11;
128        else idx=0;
129      }
130      if (tag == "p" || tag == "2_p") {
131        if      (0.035 <= value && value <= 0.11)  idx -=  2;
132        else if (0.25  <= value && value <= 0.7 )  idx -= 10;
133        else idx=0;
134      }
135      if (tag == "2_K") {
136        if      (value <= 0.07                  )  ;                 
137        else if (0.09 <= value && value <= 0.11 )  idx -=  4;           
138        else if (0.25  <= value && value <= 0.7 )  idx -= 12;
139        else idx=0;
140      }
141      if (tag == "3_p") {
142        if      (0.035 <= value && value <= 0.07)  idx -=  2;
143        else if (0.09  <= value && value <= 0.11)  idx -=  6;
144        else if (0.25  <= value && value <= 0.7 )  idx -= 14;
145        else idx=0;
146      }
147      if(idx && idx<= _edges[tag].size()) edge = _edges[tag][idx-1];
148      _nd[label]->fill(edge, weight);
149    }
150
151
152    /// @name Histograms
153    /// @{
154    map<string, Histo1DPtr> _h;
155    map<string, BinnedHistoPtr<string>> _nd;
156    map<string, BinnedEstimatePtr<string>> _ratio;
157    map<string, vector<string>> _edges;
158    YODA::Axis<double> _axis{0.025, 0.03, 0.035, 0.04, 0.045, 0.05, 0.055, 0.06, 0.06, 0.065,
159                             0.07, 0.075, 0.08, 0.085, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14,
160                             0.16, 0.18, 0.2, 0.22, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7};
161
162
163    /// }
164
165
166  };
167
168
169  RIVET_DECLARE_PLUGIN(TPC_1988_I262143);
170
171
172}