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_1979_I140294 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1979_I140294);
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
24      // Book histograms
25      book(_c_hadrons1, "/TMP/sigma_hadrons1", refData<YODA::BinnedEstimate<string>>(1,1,1));
26      book(_c_muons1,   "/TMP/sigma_muons1"  , refData<YODA::BinnedEstimate<string>>(1,1,1));
27      book(_c_hadrons2, "/TMP/sigma_hadrons2", refData<YODA::BinnedEstimate<int>>(2,1,1));
28      book(_c_muons2,   "/TMP/sigma_muons2"  , refData<YODA::BinnedEstimate<int>>(2,1,1));
29      for (const string& en : _c_hadrons1.binning().edges<0>()) {
30        double end = std::stod(en)*GeV;
31        if(isCompatibleWithSqrtS(end)) {
32          _ecms = en;
33          break;
34        }
35      }
36      if(_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
37    }
38
39
40    /// Perform the per-event analysis
41    void analyze(const Event& event) {
42      const FinalState& fs = apply<FinalState>(event, "FS");
43
44      map<long,int> nCount;
45      int ntotal(0);
46      for (const Particle& p : fs.particles()) {
47	nCount[p.pid()] += 1;
48	++ntotal;
49      }
50      // mu+mu- + photons
51      if(nCount[-13]==1 and nCount[13]==1 &&
52	 ntotal==2+nCount[22]) {
53	_c_muons1->fill(_ecms);
54	_c_muons2->fill(round(sqrtS()/GeV));
55      }
56      // everything else
57      else {
58	_c_hadrons1->fill(_ecms);
59	_c_hadrons2->fill(round(sqrtS()/GeV));
60      }
61    }
62
63
64    /// Normalise histograms etc., after the run
65    void finalize() {
66      BinnedEstimatePtr<string> mult1;
67      book(mult1, 1, 1, 1);
68      divide(_c_hadrons1,_c_muons1,mult1);
69      BinnedEstimatePtr<int> mult2;
70      book(mult2, 2, 1, 1);
71      divide(_c_hadrons2,_c_muons2,mult2);
72    }
73
74    /// @}
75
76
77    /// @name Histograms
78    /// @{
79    BinnedHistoPtr<string> _c_hadrons1, _c_muons1;
80    BinnedHistoPtr<int>    _c_hadrons2, _c_muons2;
81    string _ecms;
82    /// @}
83
84
85  };
86
87
88  RIVET_DECLARE_PLUGIN(PLUTO_1979_I140294);
89
90
91}