1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief Add a short analysis description here
 9  class PLUTO_1979_I142517 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1979_I142517);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // Initialise and register projections
22      declare(FinalState(), "FS");
23      book(_c_hadrons, "sigma_hadrons", refData<YODA::BinnedEstimate<string>>(1,1,1));
24      book(_c_muons, "sigma_muons", refData<YODA::BinnedEstimate<string>>(1,1,1));
25      book(_mult, 1, 1, 1);
26      for (const string label : {"22.0", "27.6", "30.0", "31.6"}) {
27        const double E = std::stod(label);
28        if (isCompatibleWithSqrtS(E*GeV))  Ecm = label;
29      }
30    }
31
32
33    /// Perform the per-event analysis
34    void analyze(const Event& event) {
35      const FinalState& fs = apply<FinalState>(event, "FS");
36
37      map<long,int> nCount;
38      int ntotal(0);
39      for (const Particle& p : fs.particles()) {
40        nCount[p.pid()] += 1;
41        ++ntotal;
42      }
43      if (nCount[-13]==1 and nCount[13]==1 && ntotal==2+nCount[22]) {
44        _c_muons->fill(Ecm); // mu+mu- + photons
45      }
46      else {
47        _c_hadrons->fill(Ecm); // everything else
48      }
49    }
50
51
52    /// Normalise histograms etc., after the run
53    void finalize() {
54      const double fact = crossSection()/ sumOfWeights() /nanobarn;
55      scale({_c_hadrons, _c_muons}, fact);
56      divide(_c_hadrons, _c_muons, _mult);
57    }
58
59    /// @}
60
61
62    /// @name Histograms
63    /// @{
64    BinnedHistoPtr<string> _c_hadrons, _c_muons;
65    BinnedEstimatePtr<string> _mult;
66    string Ecm = "";
67    /// @}
68
69
70  };
71
72
73  RIVET_DECLARE_PLUGIN(PLUTO_1979_I142517);
74
75
76}