1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/ChargedFinalState.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief Charged particle spectra between 2.2 and 4.8 GeV
 9  class BESII_2004_I622224 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(BESII_2004_I622224);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      const ChargedFinalState fs;
22      declare(fs, "FS");
23      unsigned int iloc(0);
24      if(isCompatibleWithSqrtS(2.2 , 1E-3)) iloc = 1;
25      else if(isCompatibleWithSqrtS(2.6 , 1E-3)) iloc = 2;
26      else if(isCompatibleWithSqrtS(3.0 , 1E-3)) iloc = 3;
27      else if(isCompatibleWithSqrtS(3.2 , 1E-3)) iloc = 4;
28      else if(isCompatibleWithSqrtS(4.6 , 1E-3)) iloc = 5;
29      else if(isCompatibleWithSqrtS(4.8 , 1E-3)) iloc = 6;
30      else MSG_ERROR("Beam energy not supported!");
31      assert(iloc != 0);
32      book(_h_ln, iloc, 1, 1);
33      book(_h_weight, "TMP/Weight");
34    }
35
36    string map2string(const double xi) const {
37      const size_t idx = axis.index(xi) - 1;
38      if (idx < edges.size())  return edges[idx];
39      return "OTHER";
40    }
41
42    /// Perform the per-event analysis
43    void analyze(const Event& event) {
44      if (edges.empty())  edges = _h_ln->xEdges();
45      const ChargedFinalState& fs = apply<ChargedFinalState>(event, "FS");
46      if (fs.particles().size()==2 &&
47          abs(fs.particles()[0].pid())==13 &&
48          abs(fs.particles()[1].pid())==13) vetoEvent;
49      for (const Particle& p : fs.particles()) {
50        const Vector3 mom3 = p.p3();
51        double pp = mom3.mod();
52        double xi = -log(2.*pp/sqrtS());
53        _h_ln->fill(map2string(xi));
54      }
55      _h_weight->fill();
56    }
57
58
59    /// Normalise histograms etc., after the run
60    void finalize() {
61      scale(_h_ln, 1./dbl(*_h_weight));
62      vector<string> edges;
63    }
64
65    /// @}
66
67
68    /// @name Histograms
69    /// @{
70    BinnedHistoPtr<string> _h_ln;
71    CounterPtr _h_weight;
72    vector<string> edges;
73    YODA::Axis<double> axis{0.0, 0.3875, 0.5875, 0.7, 0.8, 0.9, 1.0,
74                            1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
75                            2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.1, 3.4};
76    /// @}
77
78
79  };
80
81
82  RIVET_DECLARE_PLUGIN(BESII_2004_I622224);
83
84}