Rivet analyses

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

namespace Rivet {


  /// @brief Measurement of photonuclear jet production in ultra-peripheral Pb+Pb at 5.02 TeV
  class ATLAS_2024_I2829427 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2024_I2829427);

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

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

      // Jet reconstruction: anti-kT R=0.4
      declare(FastJets(FinalState(Cuts::abseta < 4.9), JetAlg::ANTIKT, 0.4,
                       JetMuons::NONE, JetInvisibles::NONE), "Jets");

      groupBook(2);
      groupBook(3);
    }


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

      // Retrieve clustered jets
      Jets js = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 15*GeV && Cuts::abseta < 4.4);

      // Require at least two jets in the event
      if (js.size() < 2) {
        MSG_DEBUG("Vetoing event: Did not find at least two jets passing cuts.");
        vetoEvent;
      }

      // Sum the 4-momenta of all selected jets
      FourMomentum p_sys;
      double hT = 0.0;
      for (const Jet& j : js) {
        p_sys += j.mom();
        hT += j.pT();
      }

      // Calculate the invariant mass of the jet system
      const double m_jets = p_sys.mass();

      // ASSUMPTION: The photon-going direction is +z
      // If the generator convention is different, this sign may need to be flipped
      const double y_jets = p_sys.rapidity();

      // Apply correlation cut between m_jets and HT
      if (m_jets <= 0.9 * hT || m_jets >= 4.0 * hT) {
        MSG_DEBUG("Vetoing event: Failed m_jets vs HT correlation cut.");
        vetoEvent;
      }

      // Apply H_T cut
      if (hT < 35*GeV || 275*GeV <= hT) {
        MSG_DEBUG("Vetoing event: HT cut.");
        vetoEvent;
      }

      // Calculate further kinematics variables
      const double sqrt_s_NN = 5020.0*GeV;
      const double xA = (m_jets / sqrt_s_NN) * exp(-y_jets);
      const double z_gamma = (m_jets / sqrt_s_NN) * exp(y_jets);

      // Fill the histograms
      groupFill(hT/GeV, y_jets, m_jets/GeV, 0);
      groupFill(hT/GeV, xA, z_gamma, 1);
    }


    /// Normalise histograms to cross-section
    void finalize() {
      scale(_h, crossSection() / microbarn / sumW());

      const double hTwidth = _xaxis[0].max(_xaxis[0].numBins()) - _xaxis[0].min(1);

      for (size_t ih=0; ih<_yaxis[0].numBins(); ++ih) {
        scale(_h[ih], 1.0/(hTwidth*_yaxis[0].width(ih+1)));
      }
      for (size_t ih=0; ih<_zaxis[0].numBins(); ++ih) {
        scale(_h[_yaxis[0].numBins()+ih], 1.0/(hTwidth*_zaxis[0].width(ih+1)));
      }

      for (size_t g : {0,1}) { // which observable group
        int ih = g? 100 : 400;
        vector<YODA::Axis<double>> axes{ _xaxis[g], _yaxis[g], _zaxis[g] };
        vector<size_t> N = {axes[0].numBins(), axes[1].numBins(), axes[2].numBins()};
        for (size_t idx : {0,1,2}) {
          for (size_t i0=0; i0<N[idx]; ++i0) {
            for (size_t i1=0; i1<N[(idx+1)%3]; ++i1) {
              if (_h.count(ih)) {
                scale(_h[ih], 1.0/(axes[idx].width(i0+1)*axes[(idx+1)%3].width(i1+1)));
              }
              ++ih;
            }
          }
        }
      }
    }

    /// @}


    /// @name Histo helpers
    /// @{

    void groupBook(unsigned int id) {

      bool isXAZgamma(id == 3);

      auto ref = refData<YODA::Estimate3D>(id, 1, 1);
      _xaxis[isXAZgamma] = YODA::Axis<double>(ref.xEdges());
      _yaxis[isXAZgamma] = YODA::Axis<double>(ref.yEdges());
      _zaxis[isXAZgamma] = YODA::Axis<double>(ref.zEdges());

      size_t Nx = _xaxis[isXAZgamma].numBins();
      size_t Ny = _yaxis[isXAZgamma].numBins();
      size_t Nz = _zaxis[isXAZgamma].numBins();

      size_t N = Nx*Ny + Nx*Nz + Ny*Nz;
      size_t offset = isXAZgamma? 100 : 400;
      for (size_t i = 0; i < N; ++i) {
        if (!hasRefData(100+offset+i, 1, 1))  continue;
        book(_h[offset+i], 100+offset+i, 1, 1);
      }

      if (!isXAZgamma) {
        // The paper also shows HT-integrated versions
        for (size_t i = 0; i < Ny*Nz; ++i) {
          if (!hasRefData(100+i, 1, 1))  continue;
          book(_h[i], 100+i, 1, 1);
        }
        // Keep track of which HT bins were actually measured
        const YODA::Estimate3D& matrix = refData<YODA::Estimate3D>(2,1,1);
        for (size_t ix = 0; ix < matrix.numBinsX(); ++ix) {
          for (size_t iy = 0; iy < matrix.numBinsY(); ++iy) {
            for (size_t iz = 0; iz < matrix.numBinsZ(); ++iz) {
              isMeasured[ix][iy][iz] = std::isfinite(matrix.bin(ix+1,iy+1,iz+1).val());
            }
          }
        }
      }
    }


    /// Fill a histogram group en mass
    void groupFill(double posX, double posY, double posZ, bool isXAZgamma) {
      const size_t Nx = _xaxis[isXAZgamma].numBins();
      const size_t Ny = _yaxis[isXAZgamma].numBins();
      const size_t Nz = _zaxis[isXAZgamma].numBins();

      const size_t ix = _xaxis[isXAZgamma].index(posX);
      const size_t iy = _yaxis[isXAZgamma].index(posY);
      const size_t iz = _zaxis[isXAZgamma].index(posZ);

      if (ix == 0 || iy == 0 || iz == 0)  return; // skip underflow
      if (ix > Nx || iy > Ny || iz > Nz)  return; // skip overflow

      size_t offset = isXAZgamma? 100 : 400;

      size_t ih = offset + (ix-1)*Ny + (iy-1);
      if (_h.count(ih)) _h[ih]->fill(posZ);

      ih = offset + Nx*Ny + (iy-1)*Nz + (iz-1);
      if (_h.count(ih)) _h[ih]->fill(posX);

      ih = offset + Nx*Ny + Ny*Nz + (iz-1)*Nx + (ix-1);
      if (_h.count(ih)) _h[ih]->fill(posY);

      if (!isXAZgamma) {
        // Fill m_jets and y_jets histograms within given H_T range
        ih = (iy-1);
        if (isMeasured[ix-1][iy-1][iz-1]) _h[ih]->fill(posZ);

        ih = Ny + (iz-1);
        if (isMeasured[ix-1][iy-1][iz-1]) _h[ih]->fill(posY);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    map<int,Histo1DPtr> _h;
    YODA::Axis<double> _xaxis[2], _yaxis[2], _zaxis[2];

    // Keep track of which 3D bins were actually measured
    bool isMeasured[8][10][9];
    /// @}

  };


  RIVET_DECLARE_PLUGIN(ATLAS_2024_I2829427);

}