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 collisions
  class ALICE_2024_I2803487 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2024_I2803487);


    /// @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;
      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["p_t"], 1, 1, 1);
      book(_h["x_t"], 5, 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
      FourMomentum mom_in_EtCone;
      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();
      }

      // 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;
      }

      if (mom_in_EtCone.pT() > _iso_et_cut)  vetoEvent;

      _h["p_t"]->fill(photon.pT()/GeV);
      // scale event weight by 1/(2 pi * pT ) as in ALICE_2017_I1512110 (with d eta applied later)
      _h["x_t"]->fill(2* photon.pT() / sqrtS(), 1.0 / (2. * M_PI* photon.pT()));

    }

    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h["p_t"], crossSection()/nanobarn/sumOfWeights() / (2* _abseta)); // also divide by |eta| range
      scale(_h["x_t"], crossSection()/picobarn/sumOfWeights()* pow(sqrtS() / GeV, 4.5)/ (2* _abseta) / (sqrtS() / 2.) ); // scaled by sqrt(s)^4.5 and also divide by |eta| range. Further division by (sqrt(s)/2) to convert from xt to pT
    }

    /// @}

    // Analysis parameters
    bool _iso_charged = true;
    const double _abseta = 0.67;
    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_2024_I2803487);

}