Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/LeptonFinder.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/PromptFinalState.hh"

namespace Rivet {


  /// @brief WZgamma production in 3-lepton channels at 13 TeV
  class CMS_2025_I2905870 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2025_I2905870);


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

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

      // Initialise and register projections

      // The basic final-state projection:
      // all final-state particles within
      // the given eta acceptance
      const FinalState fs(Cuts::abseta < 4.9);
      declare(fs, "finalstate");

      // 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 jetfs(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
      declare(jetfs, "jets");

      // FinalState of prompt photons and bare muons and electrons in the event
      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
      PromptFinalState bare_leps(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);

      // Dress the bare prompt leptons with prompt photons within dR < 0.1,
      // and apply some fiducial cuts on the dressed leptons
      Cut lepton_cuts = Cuts::abseta < 2.5 && Cuts::pT > 15*GeV;
      // Prompt photons used for 'dressed leptons'
      PromptFinalState prompt_photons(Cuts::abspid == PID::PHOTON);
      prompt_photons.acceptMuonDecays(true);
      prompt_photons.acceptTauDecays(true);
      LeptonFinder dressed_leps(bare_leps, prompt_photons, 0.1, lepton_cuts);
      declare(prompt_photons, "prompt_photons");
      declare(photons, "photons");
      declare(dressed_leps, "leptons");

      // Book histograms
      // specify custom binning
      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
      book(_h_xsec_exp, 7, 1, 1);
    }


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

      // Retrieve dressed leptons, sorted by pT
      Particles leptons = apply<LeptonFinder>(event, "leptons").particlesByPt();
      Particles muons = select(leptons, Cuts::abspid == PID::MUON);
      Particles electrons = select(leptons, Cuts::abspid == PID::ELECTRON);
      Particles photons = apply<PromptFinalState>(event,"photons").particlesByPt(Cuts::abseta < 2.5 && Cuts::pT > 15*GeV);
      const Particles& fsparticles = apply<FinalState>(event, "finalstate").particles();
      Particles iso_photons;

      if (leptons.size() != 3) vetoEvent;

      if (photons.size() == 0) vetoEvent;

      // Overlap removal: discard photons too close to any lepton
      double _lepton_photon_dr_cut = 0.3;
      idiscardIfAnyDeltaRLess(photons, leptons, _lepton_photon_dr_cut);

      if (photons.size() == 0) vetoEvent;

      // photon et and dr cut
      double _phoIsoConeSize = 0.4;
      double _phoMaxRelIso = 0.5;
      for (const auto& photon : photons) {

        double photonptsum = 0;
        for (const auto& fsparticle : fsparticles) {
          if (isSame(fsparticle, photon))  continue;

          if (deltaR(fsparticle, photon) > _phoIsoConeSize)  continue;

          photonptsum += fsparticle.pt();
        }

        if (photonptsum/photon.pt() > _phoMaxRelIso)  continue;

        iso_photons.push_back(photon);
      }

      if (iso_photons.size() == 0) vetoEvent;

      // extract mZ
      double mZ = -1;
      double temp_mZ1 = -1;
      double temp_mZ2 = -1;
      // emumu and muee final states
      if ((muons.size() == 2) && (electrons.size() == 1)) {
        if (muons[0].pid() * muons[1].pid() > 0) vetoEvent;
        mZ = (muons[0].mom() + muons[1].mom()).mass();
      }
      if ((muons.size() == 1) && (electrons.size() == 2)) {
        if (electrons[0].pid() * electrons[1].pid() > 0) vetoEvent;
        mZ = (electrons[0].mom() + electrons[1].mom()).mass();
      }

      // eee and mumumu final states
      if (muons.size() == 3) {
        if ((muons[0].pid() * muons[1].pid() > 0) && (muons[0].pid() * muons[2].pid() > 0)) vetoEvent;
        if (muons[0].pid() * muons[1].pid() < 0) {
          temp_mZ1 = (muons[0].mom() + muons[1].mom()).mass();
          if (muons[0].pid() * muons[2].pid() < 0) {
            temp_mZ2 = (muons[0].mom() + muons[2].mom()).mass();
          } else {
            temp_mZ2 = (muons[1].mom() + muons[2].mom()).mass();
          }
        } else {
          temp_mZ1 = (muons[0].mom() + muons[2].mom()).mass();
          temp_mZ2 = (muons[1].mom() + muons[2].mom()).mass();
        }
        if (fabs(temp_mZ1 - 91.188*GeV) > fabs(temp_mZ2 - 91.188*GeV)){
          mZ = temp_mZ2;
        } else {
          mZ = temp_mZ1;
        }
      }

      if (electrons.size() == 3) {
        if ((electrons[0].pid() * electrons[1].pid() > 0) && (electrons[0].pid() * electrons[2].pid() > 0)) vetoEvent;
        if (electrons[0].pid() * electrons[1].pid() < 0) {
          temp_mZ1 = (electrons[0].mom() + electrons[1].mom()).mass();
          if (electrons[0].pid() * electrons[2].pid() < 0) {
            temp_mZ2 = (electrons[0].mom() + electrons[2].mom()).mass();
          } else {
            temp_mZ2 = (electrons[1].mom() + electrons[2].mom()).mass();
          }
        } else {
          temp_mZ1 = (electrons[0].mom() + electrons[2].mom()).mass();
          temp_mZ2 = (electrons[1].mom() + electrons[2].mom()).mass();
        }
        if (fabs(temp_mZ1 - 91.188*GeV) > fabs(temp_mZ2 - 91.188*GeV)){
          mZ = temp_mZ2;
        } else {
          mZ = temp_mZ1;
        }
      }

      if ((mZ > 120*GeV) || (mZ < 60*GeV)) vetoEvent;

      _h_xsec_exp->fill();

    }


    /// Normalise histograms etc., after the run
    void finalize() {

      scale(_h_xsec_exp, crossSection() / femtobarn / sumOfWeights()); // norm to generated cross-section (after cuts)

    }

    /// @}


    /// @name Histograms
    /// @{
    CounterPtr _h_xsec_exp;
    /// @}


  };


  RIVET_DECLARE_PLUGIN(CMS_2025_I2905870);

}