1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/TriggerUA5.hh"
 4#include "Rivet/Projections/ChargedFinalState.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// UA5 charged multiplicity measurements at 546 GeV
10  class UA5_1987_I244829 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(UA5_1987_I244829);
15
16
17    /// Book histograms and initialise projections before the run
18    void init() {
19      declare(TriggerUA5(), "Trigger");
20      declare(ChargedFinalState(Cuts::abseta < 5.0), "CFS");
21
22      book(_hist_mean_nch, 1, 1, 1);
23      book(_hist_nch,      3, 1, 1);
24      book(_sumWPassed, "SumW");
25
26      vector<double> tmp; tmp.reserve(41);
27      for (int i = 1; i < 41; ++i)  tmp.push_back(2*i-1);
28      tmp.push_back(82);
29      _axis = YODA::Axis<double>(tmp);
30    }
31
32
33    /// Perform the per-event analysis
34    void analyze(const Event& event) {
35      if (_edges.empty())  _edges = _hist_nch->xEdges();
36      // Trigger
37      const TriggerUA5& trigger = apply<TriggerUA5>(event, "Trigger");
38      if (!trigger.nsdDecision()) vetoEvent;
39
40      _sumWPassed->fill();
41
42      // Count final state particles in several eta regions
43      const int Nch = apply<ChargedFinalState>(event, "CFS").size();
44
45      // Fill histograms
46      _hist_nch->fill(map2string(Nch));
47      _hist_mean_nch->fill(546, Nch);
48
49    }
50
51    string map2string(const size_t nch) const {
52      if ( 82 <= nch && nch <=  84)  return   "82.0 - 84.0";
53      if ( 86 <= nch && nch <=  88)  return   "86.0 - 88.0";
54      if ( 90 <= nch && nch <=  92)  return   "90.0 - 92.0";
55      if ( 94 <= nch && nch <=  96)  return   "94.0 - 96.0";
56      if ( 98 <= nch && nch <= 100)  return  "98.0 - 100.0";
57      if (102 <= nch && nch <= 110)  return "102.0 - 110.0";
58      if (112 <= nch && nch <= 120)  return "112.0 - 120.0";
59      const size_t idx = _axis.index(nch);
60      if (idx && idx <= _edges.size())  return _edges[idx-1];
61      return "OTHER";
62    }
63
64
65    /// Normalise histograms etc., after the run
66    void finalize() {
67
68      scale(_hist_nch, 1.0 / *_sumWPassed);
69      scale(_hist_mean_nch, 1.0 / *_sumWPassed);
70
71      for (auto& b : _hist_nch->bins()) {
72        const size_t idx = b.index();
73        double sf = idx < 46? 2.0 : 8.0;
74        if (idx == 40)  sf = 3.0;
75        b.scaleW(1.0/sf);
76      }
77
78    }
79
80
81  private:
82
83    CounterPtr _sumWPassed;
84
85    BinnedHistoPtr<int> _hist_mean_nch;
86    BinnedHistoPtr<string> _hist_nch;
87    vector<string> _edges;
88    YODA::Axis<double> _axis;
89
90  };
91
92
93
94  RIVET_DECLARE_ALIASED_PLUGIN(UA5_1987_I244829, UA5_1987_S1640666);
95
96}