Rivet analyses

#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DirectFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/InvisibleFinalState.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {


  /// @brief Lepton differential ttbar analysis at 13 TeV
  class ATLAS_2025_I2971071 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2025_I2971071);

    void init() {

      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0 * MeV;

      // Get photons to dress leptons
      DirectFinalState photons(Cuts::pid == PID::PHOTON);

      // Projection to find the electrons
      DirectFinalState prompt_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
      LeptonFinder elecs(prompt_el, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 20 * GeV);
      LeptonFinder veto_elecs(prompt_el, photons, 0.1, eta_full);
      declare(elecs, "elecs");

      // Projection to find the muons
      DirectFinalState prompt_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
      LeptonFinder muons(prompt_mu, photons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 20 * GeV);
      LeptonFinder veto_muons(prompt_mu, photons, 0.1, eta_full);
      declare(muons, "muons");

      // needed to form jets
      const InvisibleFinalState neutrinos(OnlyPrompt::YES, TauDecaysAs::PROMPT);

      VetoedFinalState vfs;
      vfs.addVetoOnThisFinalState(veto_elecs);
      vfs.addVetoOnThisFinalState(veto_muons);
      vfs.addVetoOnThisFinalState(neutrinos);

      FastJets jet(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
      declare(jet, "Jets");

      // Book histograms
      bookHistos("lep_pt", 1);
      bookHistos("lep_eta", 2);
      bookHistos("dilep_pt", 3);
      bookHistos("dilep_mass", 4);
      bookHistos("dilep_rap", 5);
      bookHistos("dilep_dphi", 6);
      bookHistos("dilep_sumpt", 7);
      bookHistos("dilep_sumE", 8);
      bookHistos("lep_ptmax", 9);
      bookHistos("lep_ptmin", 10);

      // unrolled 2D distributions - 2nd-dim bin edges must be specified
      std::vector<double> massbins = {0., 80., 120., 200., 500.};

      bookHisto2D("lep_eta_mass", 21, massbins);
      bookHisto2D("dilep_rap_mass", 22, massbins);
      bookHisto2D("dilep_dphi_mass", 23, massbins);
    }

    void analyze(const Event& event) {
      DressedLeptons elecs = apply<LeptonFinder>(event, "elecs").dressedLeptons();
      DressedLeptons muons = apply<LeptonFinder>(event, "muons").dressedLeptons();

      if (elecs.empty() || muons.empty()) vetoEvent;
      if (elecs[0].charge() == muons[0].charge()) vetoEvent;

      // construct jets for embb final state
      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25 * GeV && Cuts::abseta < 2.5);
      // count b-tagged jets
      int nbjet = 0;
      for (const Jet& jet : jets) {
        if (jet.bTagged(Cuts::pT > 5.0 * GeV)) nbjet += 1;
      }
      // embb requires at least two b-jets
      bool embb = (nbjet >= 2);

      FourMomentum el = elecs[0].momentum();
      FourMomentum mu = muons[0].momentum();
      FourMomentum ll = elecs[0].momentum() + muons[0].momentum();

      // Fill histograms
      // include explicit overflow protection as last bins are inclusive
      fillHistos("lep_pt", min(el.pT() / GeV, 299.), embb);
      fillHistos("lep_pt", min(mu.pT() / GeV, 299.), embb);
      fillHistos("lep_eta", el.abseta(), embb);
      fillHistos("lep_eta", mu.abseta(), embb);
      fillHistos("dilep_pt", min(ll.pT() / GeV, 299.), embb);
      fillHistos("dilep_mass", min(ll.mass() / GeV, 599.), embb);
      fillHistos("dilep_rap", ll.absrap(), embb);
      fillHistos("dilep_dphi", deltaPhi(el, mu), embb);
      fillHistos("dilep_sumpt", min((el.pT() + mu.pT()) / GeV, 499.), embb);
      fillHistos("dilep_sumE", min((el.E() + mu.E()) / GeV, 799.), embb);
      fillHistos("lep_ptmax", min(max(el.pT(), mu.pT()), 299.), embb);
      fillHistos("lep_ptmin", min(min(el.pT(), mu.pT()), 199.), embb);

      // find mass bin variable, with overflow protection
      float massv = ll.mass() / GeV;
      if (massv > 499.) massv = 499.;
      // Fill unrolled 2D histograms vs mass
      fillHisto2D("lep_eta_mass", el.abseta(), massv, embb);
      fillHisto2D("lep_eta_mass", mu.abseta(), massv, embb);
      fillHisto2D("dilep_rap_mass", ll.absrap(), massv, embb);
      fillHisto2D("dilep_dphi_mass", deltaPhi(el, mu), massv, embb);
    }

    void finalize() {
      // Normalize to cross-section
      const double sf = crossSection() / femtobarn / sumOfWeights();

      // finalisation of 1D histograms
      for (auto& hist : _h) {
        const double norm = 1.0 / hist.second->integral();
        // histogram normalisation
        if (hist.first.find("norm") != string::npos)
          scale(hist.second, norm);
        else
          scale(hist.second, sf);
      }

      // finalisation of 2D histograms
      for (auto& hist : _h_multi) {
        if (hist.first.find("_norm") != std::string::npos) {
          // scaling for normalised distribution according integral of whole set
          hist.second->normalizeGroup(1.0, false);
        }
        else {
          // scaling for non-normalised distribution
          scale(hist.second, sf);
        }
      }
      divByGroupWidth(_h_multi);
    }


  private:

    /// @name Histogram helper functions
    /// @{
    void bookHistos(const std::string name, unsigned int index) {
      book(_h[name], index, 1, 1);
      book(_h["norm_" + name], index + 10, 1, 1);
      book(_h["embb_" + name], index, 1, 3);
      book(_h["norm_embb_" + name], index + 10, 1, 3);
    }

    void fillHistos(const std::string name, double value, bool embb) {
      _h[name]->fill(value);
      _h["norm_" + name]->fill(value);
      if (embb) {
        _h["embb_" + name]->fill(value);
        _h["norm_embb_" + name]->fill(value);
      }
    }

    void bookHisto2D(const std::string& name, unsigned int index, const std::vector<double>& massbins) {
      book(_h_multi[name], massbins);
      book(_h_multi[name + "_norm"], massbins);
      book(_h_multi[name + "_embb"], massbins);
      book(_h_multi[name + "_norm_embb"], massbins);
      for (size_t i = 1; i < _h_multi[name]->numBins() + 1; ++i) {
        book(_h_multi[name]->bin(i), index, 1, i);
        book(_h_multi[name + "_norm"]->bin(i), index + 3, 1, i);
        book(_h_multi[name + "_embb"]->bin(i), index, 1, 8 + i);
        book(_h_multi[name + "_norm_embb"]->bin(i), index + 3, 1, 8 + i);
      }
    }


    void fillHisto2D(const std::string& name, double val, double massval, bool embb) {
      _h_multi[name]->fill(massval, val);
      _h_multi[name + "_norm"]->fill(massval, val);
      if (embb) {
        _h_multi[name + "_embb"]->fill(massval, val);
        _h_multi[name + "_norm_embb"]->fill(massval, val);
      }
    }


    // pointers to 1D and 2D histograms
    map<string, Histo1DPtr> _h;
    map<string, Histo1DGroupPtr> _h_multi;
    /// @}
    // acceptance counter
  };

  RIVET_DECLARE_PLUGIN(ATLAS_2025_I2971071);
}