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 TASSO_1984_I199468 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1984_I199468);
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(2,1,1));
28      book(_c_muons2,   "/TMP/sigma_muons2"  , refData(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() &&
37         !inRange(sqrtS()/GeV,39.8, 43.1) &&
38         !inRange(sqrtS()/GeV,43.2, 45.2))
39        MSG_ERROR("Beam energy incompatible with analysis.");
40    }
41
42
43    /// Perform the per-event analysis
44    void analyze(const Event& event) {
45
46      const FinalState& fs = apply<FinalState>(event, "FS");
47
48      map<long,int> nCount;
49      int ntotal(0);
50      for (const Particle& p : fs.particles()) {
51	nCount[p.pid()] += 1;
52	++ntotal;
53      }
54      // mu+mu- + photons
55      if(nCount[-13]==1 and nCount[13]==1 &&
56	 ntotal==2+nCount[22]) {
57	_c_muons1->fill(_ecms);
58	_c_muons2->fill(sqrtS()/GeV);
59      }
60      // everything else
61      else {
62	_c_hadrons1->fill(_ecms);
63	_c_hadrons2->fill(sqrtS()/GeV);
64      }
65    }
66
67    /// Normalise histograms etc., after the run
68    void finalize() {
69      BinnedEstimatePtr<string> mult1;
70      book(mult1, 1, 1, 1);
71      divide(_c_hadrons1,_c_muons1,mult1);
72      Estimate1DPtr mult2;
73      book(mult2, 2, 1, 1);
74      divide(_c_hadrons2,_c_muons2,mult2);
75    }
76    /// @}
77
78    /// @name Histograms
79    /// @{
80    BinnedHistoPtr<string> _c_hadrons1, _c_muons1;
81    Histo1DPtr _c_hadrons2, _c_muons2;
82    string _ecms;
83    /// @}
84
85  };
86
87
88  RIVET_DECLARE_PLUGIN(TASSO_1984_I199468);
89
90
91}