1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FastJets.hh"
 4#include "Rivet/Projections/FinalState.hh"
 5#include "Rivet/Tools/HistoGroup.hh"
 6
 7namespace Rivet {
 8
 9  /// @brief Energy correlators inside jets at 13 TeV
10  class CMS_2024_I2760466 : public Analysis {
11    public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2024_I2760466);
15
16    /// Book histograms and initialise projections before the run
17    void init() {
18      // The final-state particles declared above are clustered using FastJet with
19      // the anti-kT algorithm and a jet-radius parameter 0.4
20      declare(FastJets(FinalState(), JetAlg::ANTIKT, 0.4), "jets");
21
22      // Book histograms
23      const vector<double> ptbins{97., 220., 330., 468., 638., 846., 1101., 1410., 1784.};
24      book(_e2c, ptbins);
25      book(_e3c, ptbins);
26      for (unsigned int i = 0; i < ptbins.size() - 1; ++i) { // i=[0, 7]
27        book(_e2c->bin(i+1), i+1, 1, 1);
28        book(_e3c->bin(i+1), i+9, 1, 1);
29        book(_ratios["ratio_"s+to_string(i+1)], 26 + i, 1, 1); // correspond to hepdata t[26, 33]
30      }
31    }
32
33    /// Perform the per-event analysis
34    void analyze(const Event& event) {
35      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
36      const Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 30 * GeV);
37
38      if (jets.size() < 2)  vetoEvent;  // at least 2 jets
39      // pT and eta requirement
40      for (size_t i = 0; i < 2; ++i) {
41        if (jets[i].abseta() > 2.1) vetoEvent;
42      }
43      // back to back, diffs between leading and subleading
44      const double dphi = deltaPhi(jets[0], jets[1]);
45      const double dpt = jets[0].pT() - jets[1].pT();
46      const double spt = jets[0].pT() + jets[1].pT();
47      if (abs(dphi) <= 2 || abs(dpt) / (spt) >= 0.3) vetoEvent;
48
49      // Calculate first-two-leading-jets e2c and e3c
50      for (size_t iJet = 0; iJet < 2; ++iJet) {
51        const Jet& jet = jets[iJet];
52        // 1GeV cut to sub-particles inside jet
53        const Particles particles = select(jet.particles(), Cuts::pT >= 1.0*GeV);
54        const double E = sum(particles, Kin::E, 0.0);
55
56        // Fill e2c and e3c histograms
57        for (const Particle& p1 : particles) {
58          for (const Particle& p2 : particles) {
59            const double dr1 = deltaR(p1, p2);
60            _e2c->fill(jet.pT(), dr1, p1.E() * p2.E() / sqr(E));
61            for (const Particle& p3 : particles) {
62              const double dr2 = deltaR(p3, p2);
63              const double dr3 = deltaR(p3, p1);
64              const double rl = max(max(dr1, dr2), dr3);
65              _e3c->fill(jet.pT(), rl, p1.E() * p2.E() * p3.E() / intpow(E, 3));
66            }
67          }
68        }
69      }
70    }
71
72    /// Normalise histograms etc., after the run
73    void finalize() {
74      normalize(_e2c);
75      normalize(_e3c);
76      for (const auto& b2 : _e2c->bins()) {
77        if (b2->sumW() <= 0.)  continue;
78        divide(_e3c->bin(b2.index()), b2, _ratios["ratio_"s+to_string(b2.index())]);
79      }
80    }
81
82    private:
83
84    Histo1DGroupPtr _e2c, _e3c;
85    map<string, Estimate1DPtr> _ratios;
86
87  };
88
89  RIVET_DECLARE_PLUGIN(CMS_2024_I2760466);
90
91}  // namespace Rivet