Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Tools/HistoGroup.hh"

namespace Rivet {

  /// @brief Jet cross-section ratios at 13 TeV
  class ATLAS_2024_I2791854 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2024_I2791854);


    double maxDY(const Jets& jets) const {
      double biggest_y_gap = 0;
      for (size_t iJet = 0; iJet < jets.size(); ++iJet) {
        for (size_t jJet = iJet + 1; jJet < jets.size(); ++jJet) {
          double y_gap = std::abs(jets[iJet].rap() - jets[jJet].rap());
          if (y_gap > biggest_y_gap) {
            biggest_y_gap = y_gap;
          }
        }
      }
      return biggest_y_gap;
    }

    double maxMjj(const Jets& jets) const {
      double maxMass = 0;
      for (size_t iJet = 0; iJet < jets.size(); ++iJet) {
        for (size_t jJet = iJet + 1; jJet < jets.size(); ++jJet) {
          double mass = (jets[iJet].mom() + jets[jJet].mom()).mass();
          if (mass > maxMass) {
            maxMass = mass;
          }
        }
      }
      return maxMass;
    }


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      const FinalState fs(Cuts::abseta < 4.9);

      // The final-state particles declared above are clustered using FastJet with
      // the anti-kT algorithm and a jet-radius parameter 0.4
      // muons and neutrinos are excluded from the clustering
      FastJets jet4(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
      declare(jet4, "Jets");

      // booking pT3-binned HT2 results
      for (size_t pt3bin = 0; pt3bin < _pt3_cuts.size(); ++pt3bin) {
        HistoGroupPtr<int, double> tmp;
        _ht2.emplace_back(book(tmp, {2, 3, 4, 5}));
        for (auto& b : _ht2[pt3bin]->bins()) {
          book(b, 33 + b.index() + 4 * pt3bin, 1, 1);
        }
      }

      // booking all other results
      for (size_t njet = 0; njet < 4; ++njet) {
        Histo1DPtr tmp;
        size_t hid = 54 + njet;

        if (njet < 3) {
          _h_ptjet.emplace_back(book(tmp, hid, 1, 1));
        }

        hid += 3;
        _h_mjjmax.emplace_back(book(tmp, hid, 1, 1));

        hid += 4;
        _h_mjj.emplace_back(book(tmp, hid, 1, 1));

        hid += 4;
        _h_dy.emplace_back(book(tmp, hid, 1, 1));

        hid += 4;
        _h_dymax.emplace_back(book(tmp, hid, 1, 1));
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {


      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
      const Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 60 * GeV && Cuts::abseta < 4.5);

      if (jets.size() < 2) vetoEvent; // Dijet events

      const double mjj = (jets[0].mom() + jets[1].mom()).mass();
      const double mjjmax = maxMjj(jets);
      const double dy = std::abs(jets[0].rap() - jets[1].rap());
      const double dymax = maxDY(jets);
      const double ht2 = jets[0].pt() + jets[1].pt();

      if (ht2 < 250) vetoEvent;

      const double pt3 = jets.size() > 2 ? jets[2].pt() : 0.0;
      for (size_t nj = std::min(5ul, jets.size()); nj >= 2; --nj) {
        for (size_t pt3bin = 0; pt3bin < _pt3_cuts.size(); ++pt3bin) {
          if (nj > 2 && pt3 < _pt3_cuts[pt3bin] * ht2) continue;
          _ht2[pt3bin]->fill(nj, ht2);
        }
        _h_mjj[nj - 2]->fill(mjj / GeV);
        _h_mjjmax[nj - 2]->fill(mjjmax / GeV);
        _h_dy[nj - 2]->fill(dy);
        _h_dymax[nj - 2]->fill(dymax);
      }

      size_t jetSize = jets.size();
      for (size_t j = 0; j < 3; ++j) {
        for (int i = std::min(j + 1, jetSize - 1); i >= 0; --i) {
          _h_ptjet[j]->fill(jets[i].pt() / GeV);
        }
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      const double sf = crossSectionPerEvent();
      scale(_ht2, sf);
      scale(_h_dy, sf);
      scale(_h_dymax, sf);
      scale(_h_mjj, sf);
      scale(_h_mjjmax, sf);
      scale(_h_ptjet, sf);

      Estimate1DPtr ratio;
      for (size_t ir = 0; ir < 4; ++ir) {
        size_t hid = 10 * ir + 6;
        // Temporary hack for missing hepdata
        if (ir == 1) hid = 10 + 4;
        if (ir == 2) hid = 18 + 3;
        if (ir == 3) hid = 25 + 4;

        if (ir < 2) {
          book(ratio, hid, 1, 1);
          divide(_h_ptjet[ir + 1], _h_ptjet[ir], ratio);
        }

        book(ratio, ++hid, 1, 1);
        divide(_h_dy[ir == 3 ? 2 : ir + 1], _h_dy[ir == 3 ? 0 : ir], ratio);

        book(ratio, ++hid, 1, 1);
        divide(_h_dymax[ir == 3 ? 2 : ir + 1], _h_dymax[ir == 3 ? 0 : ir], ratio);

        book(ratio, ++hid, 1, 1);
        divide(_h_mjj[ir == 3 ? 2 : ir + 1], _h_mjj[ir == 3 ? 0 : ir], ratio);

        book(ratio, ++hid, 1, 1);
        divide(_h_mjjmax[ir == 3 ? 2 : ir + 1], _h_mjjmax[ir == 3 ? 0 : ir], ratio);
      }

      for (size_t pt3bin = 0; pt3bin < _pt3_cuts.size(); ++pt3bin) {
        for (size_t ir = 0; ir < 4; ++ir) {
          size_t hid = pt3bin + 10 * ir + 1;
          if (ir == 2) hid = pt3bin + 10 + 8 + 1;
          if (ir == 3) hid = pt3bin + 10 + 8 + 7 + 1;

          if ((pt3bin == 1 || pt3bin == 3) && ir > 0) continue;
          if (ir > 0 && pt3bin == 2) hid -= 1;
          if (ir > 0 && pt3bin == 4) hid -= 2;

          book(ratio, hid, 1, 1);
          MSG_DEBUG("ratio HT2: " << hid << " for pt3bin = " << pt3bin << " and ir = " << ir);
          divide(_ht2[pt3bin]->bin(ir == 3 ? 3 : ir + 2), _ht2[pt3bin]->bin(ir == 3 ? 1 : ir + 1), ratio);
        }
      }
    }

    /// @name Histograms
    /// @{

    // for pT-inclusive histograms
    vector<HistoGroupPtr<int, double>> _ht2;
    vector<Histo1DPtr> _h_dymax, _h_dy, _h_mjjmax, _h_mjj, _h_ptjet;

    const vector<double> _pt3_cuts{0.0, 0.05, 0.1, 0.2, 0.3};

    /// @}
  };


  RIVET_DECLARE_PLUGIN(ATLAS_2024_I2791854);

}