1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief R measurement
 9  class PLUTO_1982_I166799 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1982_I166799);
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      // counters for R
24      book(_c_hadrons[0], "/TMP/sigma_hadrons", refData<YODA::BinnedEstimate<string>>(1,1,1));
25      book(_c_hadrons[1], 2,1,1);
26      book(_c_muons, "/TMP/sigma_muons",refData<YODA::BinnedEstimate<string>>(1,1,1));
27      for(unsigned int ix=0;ix<2;++ix) {
28        for (const string& en : _c_hadrons[ix].binning().edges<0>()) {
29          const size_t idx = en.find("-");
30          if(idx!=std::string::npos) {
31            const double emin = std::stod(en.substr(0,idx));
32            const double emax = std::stod(en.substr(idx+1,string::npos));
33            if(inRange(sqrtS()/GeV, emin, emax)) {
34              _ecms[ix] = en;
35              break;
36            }
37          }
38          else {
39            const double end = std::stod(en)*GeV;
40            if (isCompatibleWithSqrtS(end)) {
41              _ecms[ix] = en;
42              break;
43            }
44          }
45        }
46      }
47      if(_ecms[0].empty() && _ecms[1].empty())
48        MSG_ERROR("Beam energy incompatible with analysis.");
49    }
50
51
52    /// Perform the per-event analysis
53    void analyze(const Event& event) {
54      const FinalState& fs = apply<FinalState>(event, "FS");
55
56      map<long,int> nCount;
57      int ntotal(0);
58      for (const Particle& p : fs.particles()) {
59        nCount[p.pid()] += 1;
60        ++ntotal;
61      }
62      if(nCount[-13]==1 and nCount[13]==1 && ntotal==2+nCount[22]) {
63        // mu+mu- + photons
64        _c_muons->fill(_ecms[0]);
65      }
66      else {
67        // everything else
68        _c_hadrons[0]->fill(_ecms[0]);
69        _c_hadrons[1]->fill(_ecms[1]);
70      }
71    }
72
73
74    /// Normalise histograms etc., after the run
75    void finalize() {
76      double fact = crossSection()/ sumOfWeights() /nanobarn;
77      scale(_c_hadrons[1],fact);
78      BinnedEstimatePtr<string> mult;
79      book(mult, 1, 1, 1);
80      divide(_c_hadrons[0],_c_muons,mult);
81    }
82
83    /// @}
84
85
86    /// @name Histograms
87    /// @{
88    BinnedHistoPtr<string> _c_hadrons[2], _c_muons;
89    string _ecms[2];
90    /// @}
91
92
93  };
94
95
96  RIVET_DECLARE_PLUGIN(PLUTO_1982_I166799);
97
98
99}