1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5
 6namespace Rivet {
 7
 8
 9  /// @brief e+e- > KS0 KL0
10  class SND_2006_I720035 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(SND_2006_I720035);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22
23      declare(FinalState(), "FS");
24      book(_nK0K0, 1, 1, 1);
25      for (const string& en : _nK0K0.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 (_ecms.empty())
44        MSG_ERROR("Beam energy incompatible with analysis.");
45    }
46
47
48    /// Perform the per-event analysis
49    void analyze(const Event& event) {
50
51      const FinalState& fs = apply<FinalState>(event, "FS");
52
53      map<long,int> nCount;
54      int ntotal(0);
55      for (const Particle& p : fs.particles()) {
56	nCount[p.pid()] += 1;
57	++ntotal;
58      }
59      if(ntotal==2 &&
60	 nCount[130]==1 && nCount[310]==1)
61	_nK0K0->fill(_ecms);
62
63    }
64
65
66    /// Normalise histograms etc., after the run
67    void finalize() {
68      scale(_nK0K0, crossSection()/ sumOfWeights() /nanobarn);
69    }
70
71    /// @}
72
73
74    /// @name Histograms
75    /// @{
76    BinnedHistoPtr<string> _nK0K0;
77    string _ecms;
78    /// @}
79
80
81  };
82
83
84  RIVET_DECLARE_PLUGIN(SND_2006_I720035);
85
86
87}