Rivet analyses
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 sNN = 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++ -*-
#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);
}