Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.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);

}