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      for(auto & b: _h_ln->bins()) {
63        const size_t idx = b.index();
64        b.scaleW(1./axis.width(idx));
65      }
66    }
67
68    /// @}
69
70
71    /// @name Histograms
72    /// @{
73    BinnedHistoPtr<string> _h_ln;
74    CounterPtr _h_weight;
75    vector<string> edges;
76    YODA::Axis<double> axis{0.05, 0.45, 0.6, 0.7, 0.8, 0.9, 1.0,
77                            1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
78                            2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.5};
79    /// @}
80
81
82  };
83
84
85  RIVET_DECLARE_PLUGIN(BESII_2004_I622224);
86
87}