1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief pi, K and proton spectra at 14,22 and 34 GeV
 10  class TASSO_1983_I181470 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1983_I181470);
 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(FinalState(), "FS");
 26
 27      vector<int> hist1,hist2;
 28      sqs = 1.;
 29      if(isCompatibleWithSqrtS(14*GeV)) {
 30	hist1 = {19,21,23};
 31	hist2 = {20,22,24};
 32	sqs = 14.;
 33      }
 34      else if (isCompatibleWithSqrtS(22*GeV)) {
 35	hist1 = {25,27,11};
 36	hist2 = {26,10,12};
 37	sqs = 22.;
 38      }
 39      else if (isCompatibleWithSqrtS(34*GeV)) {
 40	hist1 = {13,15,17};
 41	hist2 = {14,16,18};
 42	sqs = 34.;
 43      }
 44      else
 45        MSG_WARNING("CoM energy of events sqrt(s) = " << sqrtS()/GeV
 46          << " doesn't match any available analysis energy .");
 47
 48      book(_h_p_pi , hist1[0],1,1);
 49      book(_h_p_K  , hist1[1],1,1);
 50      book(_h_p_p  , hist1[2],1,1);
 51      book(_h_x_pi , hist2[0],1,1);
 52      book(_h_x_K  , hist2[1],1,1);
 53      book(_h_x_p  , hist2[2],1,1);
 54
 55    }
 56
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60      // First, veto on leptonic events by requiring at least 4 charged FS particles
 61      const FinalState& fs = apply<FinalState>(event, "FS");
 62      const size_t numParticles = fs.particles().size();
 63
 64      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 65      if (numParticles < 2) {
 66        MSG_DEBUG("Failed leptonic event cut");
 67        vetoEvent;
 68      }
 69      MSG_DEBUG("Passed leptonic event cut");
 70
 71      // Get beams and average beam momentum
 72      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 73      const double meanBeamMom = ( beams.first.p3().mod() +
 74                                   beams.second.p3().mod() ) / 2.0;
 75      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 76
 77      for (const Particle& p : fs.particles()) {
 78	double xE = p.E()/meanBeamMom;
 79	if(abs(p.pid())==211) {
 80	  _h_p_pi->fill(p.p3().mod());
 81	  _h_x_pi->fill(xE          );
 82	}
 83	else if(abs(p.pid())==321) {
 84	  _h_p_K->fill(p.p3().mod());
 85	  _h_x_K->fill(xE          );
 86	}
 87	else if(abs(p.pid())==2212) {
 88	  _h_p_p->fill(p.p3().mod());
 89	  _h_x_p->fill(xE          );
 90	}
 91      }
 92    }
 93
 94
 95    /// Normalise histograms etc., after the run
 96    void finalize() {
 97
 98      double fact1 = crossSection()/nanobarn/sumOfWeights();
 99      double fact2 = sqr(sqs)/GeV2*crossSection()/microbarn/sumOfWeights();
100
101      scale(_h_p_pi, fact1);
102      scale(_h_p_K , fact1);
103      scale(_h_p_p , fact1);
104
105      scale(_h_x_pi, fact2);
106      scale(_h_x_K , fact2);
107      scale(_h_x_p , fact2);
108    }
109
110    /// @}
111
112
113    /// @name Histograms
114    /// @{
115    Histo1DPtr   _h_p_pi,_h_p_K,_h_p_p;
116    Histo1DPtr   _h_x_pi,_h_x_K,_h_x_p;
117    double sqs;
118    /// @}
119
120
121  };
122
123
124  RIVET_DECLARE_PLUGIN(TASSO_1983_I181470);
125
126}