1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/Beam.hh"
 4#include "Rivet/Projections/UnstableParticles.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief Kaon spectra at 3.63, 4.03 and 4.5 GeV
10  class PLUTO_1977_I118873 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1977_I118873);
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(UnstableParticles(Cuts::pid==PID::K0S), "UFS");
24      book(_h_sigma,1,1,1);
25      for (const string& en :_h_sigma.binning().edges<0>()) {
26        const size_t idx = en.find("-");
27        if(idx!=std::string::npos) {
28          const double emin = std::stod(en.substr(0,idx));
29          const double emax = std::stod(en.substr(idx+1,string::npos));
30          if(inRange(sqrtS()/GeV, emin, emax)) {
31            _ecms = en;
32            break;
33          }
34        }
35        else {
36          const double end = std::stod(en)*GeV;
37          if (isCompatibleWithSqrtS(end)) {
38            _ecms = en;
39            break;
40          }
41        }
42      }
43      if (isCompatibleWithSqrtS(3.63*GeV)) {
44        book(_h_spectrum, 2, 1, 1);
45      }
46      else if (isCompatibleWithSqrtS(4.03*GeV)) {
47        book(_h_spectrum, 3, 1, 1);
48      }
49      else if (isCompatibleWithSqrtS(4.5*GeV)) {
50        book(_h_spectrum, 4, 1, 1);
51      }
52      else if(_ecms.empty())
53        MSG_ERROR("Beam energy not supported!");
54    }
55
56
57    /// Perform the per-event analysis
58    void analyze(const Event& event) {
59      // Get beams and average beam momentum
60      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
61      const double meanBeamMom = 0.5*(beams.first.p3().mod() + beams.second.p3().mod());
62      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
63      // unstable particles
64      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
65        _h_sigma->fill(_ecms);
66        const double xp = p.E()/meanBeamMom;
67        const double beta = p.p3().mod()/p.E();
68        if(_h_spectrum) _h_spectrum->fill(xp,1./beta);
69      }
70    }
71
72
73    /// Normalise histograms etc., after the run
74    void finalize() {
75      if(_h_spectrum)
76        scale(_h_spectrum, sqr(sqrtS())*crossSection()/nanobarn/sumOfWeights());
77      scale(_h_sigma, crossSection()/nanobarn/sumOfWeights());
78    }
79
80    /// @}
81
82
83    /// @name Histograms
84    /// @{
85    Histo1DPtr _h_spectrum;
86    BinnedHistoPtr<string> _h_sigma;
87    string _ecms;
88    /// @}
89
90
91  };
92
93
94  RIVET_DECLARE_PLUGIN(PLUTO_1977_I118873);
95
96
97}