Rivet analyses

title: LHCF_2016_I1385877

Longitudinal and transverse momentum distributions for neutral pions in the forward-rapidity region

Experiment: LHCF (LHC)

Inspire ID: 1385877

Status: VALIDATED

Authors: - Eugenio Berti - LHCf collaboration

References: - Phys. Rev. D 94 (2016) 032007

Beams: p+ p+, p+ p+, p+ 1000822080

Beam energies: (3500.0, 3500.0); (1380.0, 1380.0); (4000.0, 328000.0)GeV

Run details: - Inclusive production cross section of neutral pion in p-p and p-Pb collisions in the very forward region expressed as a function of transverse and longitudinal momentum.

The differential cross sections for inclusive neutral pions as a function of transverse and longitudinal momentum in the very forward-rapidity region have been measured at the LHC with the LHC forward detector in proton-proton collisions at $s=2.76$ and $7\,$TeV and in proton-lead collisions at nucleon-nucleon center-of-mass energies of $s_\textrm{NN}=5.02\,$TeV. Such differential cross sections in proton-proton collisions are compatible with the hypotheses of limiting fragmentation and Feynman scaling. Comparing proton-proton with proton-lead collisions, we find a sizeable suppression of the production of neutral pions in the differential cross sections after subtraction of ultraperipheral proton-lead collisions. This suppression corresponds to the nuclear modification factor value of about $0.1-0.3$. The experimental measurements presented in this paper provide a benchmark for the hadronic interaction Monte Carlo simulation codes that are used for the simulation of cosmic-ray air showers.

Source code:LHCF_2016_I1385877.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief Longitudinal and transverse momenta of neutral pions in the forward-rapidity region class LHCF_2016_I1385877 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(LHCF_2016_I1385877);

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

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

  // Initialise and register projections
  declare(UnstableParticles(), "UFS");
  declare(Beam(), "Beam");

  // Calculate beam rapidity
  const Particle bm1 = beams().first;
  const Particle bm2 = beams().second;
  MSG_DEBUG("Beam 1 : momentum=" << bm1.pz() << " PID=" << bm1.pid() << " rapidity=" << bm1.rap());
  MSG_DEBUG("Beam 2 : momentum=" << bm2.pz() << " PID=" << bm2.pid() << " rapidity=" << bm2.rap());
  MSG_DEBUG("CM energy: " << sqrtS());
  _beam_rap = bm1.rap();

  // Book histos for p-p or p-Pb mode
  for (double eVal : allowedEnergies()) {
    const string en = toString(round(eVal));
    if (isCompatibleWithSqrtS(eVal)) _sqs = en;
    if (en == "7000"s) {
      book(_g[en + "rap_pT"], {8.8, 9., 9.2, 9.4, 9.6, 9.8, 10., 10.2, 10.4, 10.6, 10.8});
      for (auto& b : _g[en + "rap_pT"]->bins()) {
        book(b, 1 + b.index(), 1, 2);
      }
      book(_g[en + "pT_pZ"], {0., 0.2, 0.4, 0.6, 0.8, 1.});
      for (auto& b : _g[en + "pT_pZ"]->bins()) {
        book(b, 11 + b.index(), 1, 2);
      }
      book(_p[en + "rap_apT"], 1, 1, 2);
      book(_h[en + "rap"], 21, 1, 2);
      book(_p[en + "raploss_apT"], 22, 1, 2);
      book(_h[en + "raploss"], 23, 1, 2);
    }
    if (en == "2760"s) {
      book(_p[en + "rap_apT"], 1, 1, 1);
      book(_g[en + "rap_pT"], {8.8, 9., 9.2, 9.4, 9.6, 9.8});
      for (auto& b : _g[en + "rap_pT"]->bins()) {
        book(b, 1 + b.index(), 1, 1);
      }
      book(_g[en + "pT_pZ"], {0., 0.2, 0.4}, {"d12-x01-y01", "d13-x01-y01"});
      book(_h[en + "rap"], 21, 1, 1);
      book(_p[en + "raploss_apT"], 22, 1, 1);
      book(_h[en + "raploss"], 23, 1, 1);
    }
    else if (en == "1044160"s) {
      book(_g[en + "rap_pT"], {8.8, 9., 9.2, 9.4, 9.6, 9.8, 10., 10.2, 10.4, 10.6, 10.8});
      for (auto& b : _g[en + "rap_pT"]->bins()) {
        book(b, 1 + b.index(), 1, 3);
      }
      book(_g[en + "pT_pZ"], {0., 0.2, 0.4, 0.6, 0.8, 1.0});
      for (auto& b : _g[en + "pT_pZ"]->bins()) {
        book(b, 11 + b.index(), 1, 3);
      }
      book(_p[en + "rap_apT"], 1, 1, 3);
      book(_p[en + "raploss_apT"], 22, 1, 3);
    }
  }
  raiseBeamErrorIf(_sqs.empty());
}


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

  // Select neutral pions
  const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
  const Particles pions = ufs.particles(Cuts::pz > 0 && Cuts::abspid == PID::PI0
                                        && Cuts::pT > 0.01 * GeV);
  for (const Particle& p : pions) {
    const double pT = p.pT() / GeV;
    const double rap = p.rap();
    const double raploss = _beam_rap - p.rap();
    _p[_sqs + "rap_apT"]->fill(rap, p.pT() / MeV);
    _p[_sqs + "raploss_apT"]->fill(raploss, p.pT() / MeV);
    _g[_sqs + "rap_pT"]->fill(rap, pT, 1.0 / pT);
    _g[_sqs + "pT_pZ"]->fill(pT, p.pz() / GeV, p.E() / GeV / pT);
    if (_sqs != "1044160"s) {
      _h[_sqs + "rap"]->fill(rap);
      _h[_sqs + "raploss"]->fill(raploss);
    }
  }
}


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

  const double sf = 1.0 / sumOfWeights() / TWOPI;
  const double bw = 0.2;
  for (auto& grp : _g) {
    for (auto& hist : grp.second->bins()) {
      if (grp.first.find("pT_") != string::npos) {
        if (hist.index() == 1)
          hist->scaleW(sf / (bw - 0.01));
        else
          hist->scaleW(sf / bw);
      }
      else {
        size_t cutoff_bin = hist->indexAt(0.01);
        if (cutoff_bin > 0) {
          const double cutoff_wdt = hist->bin(cutoff_bin).xWidth();
          hist->bin(cutoff_bin).scaleW((cutoff_wdt) / (cutoff_wdt - 0.01));
        }
        hist->scaleW(sf / bw);
      }
    }
  }
  scale(_h, 1.0 / sumOfWeights());
}

/// @}


// Store the beam rapidity for rap-loss calculation (could just re-access this in analyze())
double _beam_rap;

string _sqs = "";

/// @name Histograms
/// @{
map<string, Histo1DPtr> _h;
map<string, Profile1DPtr> _p;
map<string, Histo1DGroupPtr> _g;
/// @}

};

RIVET_DECLARE_PLUGIN(LHCF_2016_I1385877); } ```