Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"

namespace Rivet {


  /// @brief Measurement of the inclusive isolated photon production cross section in pp and pPb collisions
  class ALICE_2025_I2895564: public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2025_I2895564);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // This option exists because JetPhox can not do isolation on charged particles only
      _iso_charged = getOption("ISO_CHARGED", true);
      _iso_et_cut = _iso_charged ? _iso_et_cut_charged : _iso_et_cut_all;

      if (!_iso_charged) MSG_WARNING("LIMITED VALIDITY - measured data was obtained for ISO_CHARGED=true case only.");

      const FinalState fs;
      // calorimeter particles
      VisibleFinalState visFS(fs);
      VetoedFinalState calo_fs(visFS);
      declare(calo_fs, "calo");

      // photon
      PromptFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < _abseta);
      declare(photonfs, "photons");

      // Book histograms
      // specify custom binning
      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
      book(_h["pp_8_p_t"], 1, 1, 1);
      book(_h["pPb_8_16_p_t"], 2, 1, 1);
      book(_h["pPb_5_02_p_t"], 3, 1, 1);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Get the photon
      const Particles& photons = apply<PromptFinalState>(event, "photons").particlesByPt();
      if (photons.empty())  vetoEvent;
      // We only look at the hardest/leading photon
      const FourMomentum photon = photons[0].momentum();


      // Compute photon isolation with a standard ET cone
      double mom_in_EtCone = 0.0;
      const Particles calo_fs = apply<VetoedFinalState>(event, "calo").particles();
      for (const Particle& p : calo_fs) {
        // Possibly skip non-charged particles
        if (_iso_charged && !p.isCharged()) continue;
        // Check if it's in the cone of .4
        if (deltaR(p, photon) >= _iso_dr) continue;
        // Increment sum
        mom_in_EtCone += p.momentum().pT();
      }

      // photon itself is not charged, so only remove it if we consider all particles
      if (!_iso_charged) {
        // Remove the photon energy from the isolation
        mom_in_EtCone -= photon.pT();
      }

      if (mom_in_EtCone > _iso_et_cut)  vetoEvent;

      if (beamsMatch(PID::PROTON, PID::PROTON,4000*GeV,4000*GeV)) {
        _h["pp_8_p_t"]->fill(photon.pT()/GeV);
      } else if (beamsMatch(PID::LEAD, PID::PROTON,533000*GeV,6500*GeV)) {
        _h["pPb_8_16_p_t"]->fill(photon.pT()/GeV);
      } else if (beamsMatch(PID::LEAD, PID::PROTON,328000*GeV,4000*GeV)) {
        _h["pPb_5_02_p_t"]->fill(photon.pT()/GeV);
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h, crossSection()/nanobarn/sumOfWeights() / (2* _abseta)); // also divide by |eta| range
    }

    /// @}

    // Analysis parameters

    bool _iso_charged = true;
    const double _abseta = 0.7;
    const double _iso_dr = 0.4;
    double _iso_et_cut;
    const double _iso_et_cut_charged = 1.5*GeV;
    const double _iso_et_cut_all = 2.0*GeV;

    /// @name Histograms
    /// @{
    map<string, Histo1DPtr> _h;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(ALICE_2025_I2895564);

}