Rivet analyses referenceLHCF_2016_I1385877Longitudinal and transverse momentum distributions for neutral pions in the forward-rapidity regionExperiment: LHCF (LHC) Inspire ID: 1385877 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0); (1380.0, 1380.0); (4000.0, 328000.0) GeV Run details:
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. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: LHCF_2016_I1385877.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/UnstableParticles.hh"
5
6namespace Rivet {
7
8
9 /// @brief Longitudinal and transverse momenta of neutral pions in the forward-rapidity region
10 class LHCF_2016_I1385877 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(LHCF_2016_I1385877);
15
16
17 // In some models there can be very small-value pT but greater than 0.
18 // In order to avoid unphysical behavior in the first bin a cutoff is needed
19 // If you are sure the model does not have this problem you can set pt_cutoff to 0.
20 const double pt_cutoff = 0.01;
21
22
23 /// @name Analysis methods
24 /// @{
25 void bookHistosPP() {
26 if (isCompatibleWithSqrtS(7000*GeV, 1e-3)) {
27 book(_h_pi0_rap_pT, {8.8, 9., 9.2, 9.4, 9.6, 9.8, 10., 10.2, 10.4, 10.6, 10.8});
28 for (auto& b : _h_pi0_rap_pT->bins()) {
29 book(b, 1+b.index(), 1, 2);
30 }
31
32 book(_h_pi0_pT_pZ, {0., 0.2, 0.4, 0.6, 0.8, 1.});
33 for (auto& b : _h_pi0_pT_pZ->bins()) {
34 book(b, 11+b.index(), 1, 2);
35 }
36
37 book(_p_pi0_rap_apT, 1, 1, 2);
38 book(_h_pi0_rap, 21, 1, 2);
39 book(_p_pi0_raploss_apT, 22, 1, 2);
40 book(_h_pi0_raploss, 23, 1, 2);
41 }
42 else if (isCompatibleWithSqrtS(2760*GeV, 1e-3)) {
43 book(_p_pi0_rap_apT, 1, 1, 1);
44
45 book(_h_pi0_rap_pT, {8.8, 9., 9.2, 9.4, 9.6, 9.8});
46 for (auto& b : _h_pi0_rap_pT->bins()) {
47 book(b, 1+b.index(), 1, 1);
48 }
49
50 book(_h_pi0_pT_pZ, {0., 0.2, 0.4}, {"d12-x01-y01", "d13-x01-y01"});
51
52 book(_h_pi0_rap, 21, 1, 1);
53
54 book(_p_pi0_raploss_apT, 22, 1, 1);
55 book(_h_pi0_raploss, 23, 1, 1);
56 }
57 else {
58 MSG_WARNING("p-p collisions : energy out of range!");
59 }
60 }
61
62
63 void bookHistosPPb() {
64 if (isCompatibleWithSqrtS(sqrt(208)*5020*GeV, 1e-3)) {
65
66 book(_h_pi0_rap_pT, {8.8, 9., 9.2, 9.4, 9.6, 9.8, 10., 10.2, 10.4, 10.6, 10.8});
67 for (auto& b : _h_pi0_rap_pT->bins()) {
68 book(b, 1+b.index(), 1, 3);
69 }
70
71 book(_h_pi0_pT_pZ, {0., 0.2, 0.4, 0.6, 0.8, 1.0});
72 for (auto& b : _h_pi0_pT_pZ->bins()) {
73 book(b, 11+b.index(), 1, 3);
74 }
75
76 book(_p_pi0_rap_apT, 1, 1, 3);
77 book(_p_pi0_raploss_apT, 22, 1, 3);
78 }
79 else {
80 MSG_WARNING("p-Pb collisions : energy out of range!");
81 }
82 }
83
84
85 /// Book histograms and initialise projections before the run
86 void init() {
87
88 // Initialise and register projections
89 declare(UnstableParticles(), "UFS");
90 declare(Beam(), "Beam");
91
92 // Calculate beam rapidity
93 const Particle bm1 = beams().first;
94 const Particle bm2 = beams().second;
95 MSG_DEBUG("Beam 1 : momentum=" << bm1.pz() << " PID=" << bm1.pid() << " rapidity=" << bm1.rap());
96 MSG_DEBUG("Beam 2 : momentum=" << bm2.pz() << " PID=" << bm2.pid() << " rapidity=" << bm2.rap());
97 MSG_DEBUG("CM energy: " << sqrtS() );
98 _beam_rap = bm1.rap();
99
100 // Book histos for p-p or p-Pb mode
101 if (bm1.pid() == PID::PROTON && bm2.pid() == PID::PROTON) {
102 _isPP = true;
103 bookHistosPP();
104 } else if (bm1.pid() == PID::PROTON && bm2.pid() == PID::LEAD) {
105 _isPP = false;
106 bookHistosPPb();
107 } else MSG_WARNING("Beam PDG ID out of range --- should be pp or p-Pb");
108
109 }
110
111
112 /// Perform the per-event analysis
113 void analyze(const Event& event) {
114
115 // Select neutral pions
116 const UnstableParticles& ufs = apply<UnstableParticles> (event, "UFS");
117 const Particles pions = ufs.particles(Cuts::pz > 0 && Cuts::abspid == PID::PI0 && Cuts::pT > pt_cutoff*GeV);
118 for (const Particle& p : pions) {
119 const double pT = p.pT()/GeV;
120 const double rap = p.rap();
121 const double raploss = _beam_rap - p.rap();
122 _p_pi0_rap_apT->fill(rap, p.pT()/MeV);
123 _p_pi0_raploss_apT->fill(raploss, p.pT()/MeV);
124 _h_pi0_rap_pT->fill(rap, pT, 1.0/pT);
125 _h_pi0_pT_pZ->fill(pT, p.pz()/GeV, p.E()/GeV/pT);
126 if (_isPP) {
127 _h_pi0_rap->fill(rap);
128 _h_pi0_raploss->fill(raploss);
129 }
130 }
131
132 }
133
134
135 /// Normalise histograms etc., after the run
136 void finalize() {
137
138 const double inv_scale_factor = 1. / sumOfWeights() / (2.*PI);
139 const double pt_bin_width = 0.2;
140 for (auto& h : _h_pi0_pT_pZ->bins()) {
141 if (h->path() == "/LHCF_2016_I1385877/d12-x01-y01" ||
142 h->path() == "/LHCF_2016_I1385877/d12-x01-y02" ||
143 h->path() == "/LHCF_2016_I1385877/d12-x01-y03") {
144 h->scaleW( inv_scale_factor / (pt_bin_width-pt_cutoff) );
145 }
146 else {
147 h->scaleW( inv_scale_factor / pt_bin_width );
148 }
149 }
150
151 const double scale_factor = 1. / sumOfWeights() / (2.*PI);
152 const double rap_bin_width = 0.2;
153 for (auto& h : _h_pi0_rap_pT->bins()) {
154 const int cutoff_bin = h->indexAt(pt_cutoff);
155 if (cutoff_bin >= 0) {
156 const double cutoff_wdt = h->bin(cutoff_bin).xMax()-h->bin(cutoff_bin).xMin();
157 h->bin(cutoff_bin).scaleW((cutoff_wdt)/(cutoff_wdt-pt_cutoff));
158 }
159 h->scaleW( scale_factor / rap_bin_width );
160 }
161
162 if (_isPP) {
163 scale( _h_pi0_rap , 1/sumOfWeights() );
164 scale( _h_pi0_raploss , 1/sumOfWeights() );
165 }
166
167 }
168
169 /// @}
170
171
172
173 /// Flag for handling extra histograms in p-p runs
174 bool _isPP;
175 // Store the beam rapidity for rap-loss calculation (could just re-access this in analyze())
176 double _beam_rap;
177
178 /// @name Histograms
179 /// @{
180 Histo1DGroupPtr _h_pi0_pT_pZ, _h_pi0_rap_pT;
181 Profile1DPtr _p_pi0_rap_apT;
182 Histo1DPtr _h_pi0_rap;
183 Profile1DPtr _p_pi0_raploss_apT;
184 Histo1DPtr _h_pi0_raploss;
185 /// @}
186
187 };
188
189
190
191 RIVET_DECLARE_PLUGIN(LHCF_2016_I1385877);
192
193}
|