1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/ChargedFinalState.hh"
 4#include "Rivet/Projections/FastJets.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief jet rates at 91 GeV
10  class DELPHI_1990_I297698 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_1990_I297698);
15
16
17    /// @name Analysis methods
18    ///@{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      const ChargedFinalState cfs;
23      declare(cfs, "FS");
24      declare(FastJets(cfs, JetAlg::JADE), "JadeJets");
25      // histos
26      book(_h_2, 1, 1, 1);
27      book(_h_3, 1, 1, 2);
28      book(_h_4, 1, 1, 3);
29      book(_h_5, 1, 1, 4);
30    }
31
32
33    /// Perform the per-event analysis
34    void analyze(const Event& event) {
35      const FinalState& fs = apply<FinalState>(event, "FS");
36      const size_t numParticles = fs.particles().size();
37      if (numParticles < 2) {
38        MSG_DEBUG("Failed leptonic event cut");
39        vetoEvent;
40      }
41      MSG_DEBUG("Passed leptonic event cut");
42      const FastJets& jets = apply<FastJets>(event, "JadeJets");
43      if (jets.clusterSeq()) {
44        const double y_23 = jets.clusterSeq()->exclusive_ymerge_max(2);
45        const double y_34 = jets.clusterSeq()->exclusive_ymerge_max(3);
46        const double y_45 = jets.clusterSeq()->exclusive_ymerge_max(4);
47        const double y_56 = jets.clusterSeq()->exclusive_ymerge_max(5);
48        for (auto& b : _h_2->bins()) {
49          const double ycut = std::stod(b.xEdge());
50          if (y_23 < ycut) {
51            _h_2->fill(b.xEdge());
52          }
53        }
54        for (auto& b : _h_2->bins()) {
55          const double ycut = std::stod(b.xEdge());
56          if (y_34 < ycut && y_23 > ycut) {
57            _h_3->fill(b.xEdge());
58          }
59        }
60        for (auto& b : _h_4->bins()) {
61          const double ycut = std::stod(b.xEdge());
62          if (y_45 < ycut && y_34 > ycut) {
63            _h_4->fill(b.xEdge());
64          }
65        }
66        for (auto& b : _h_5->bins()) {
67          const double ycut = std::stod(b.xEdge());
68          if (y_56 < ycut && y_45 > ycut) {
69            _h_5->fill(b.xEdge());
70          }
71        }
72      }
73    }
74
75
76    /// Normalise histograms etc., after the run
77    void finalize() {
78      scale(_h_2, 100/sumOfWeights());
79      scale(_h_3, 100/sumOfWeights());
80      scale(_h_4, 100/sumOfWeights());
81      scale(_h_5, 100/sumOfWeights());
82    }
83
84    ///@}
85
86
87    /// @name Histograms
88    ///@{
89    BinnedHistoPtr<string> _h_2, _h_3, _h_4, _h_5;
90    ///@}
91
92
93  };
94
95
96  RIVET_DECLARE_PLUGIN(DELPHI_1990_I297698);
97
98}