Rivet analyses referenceLHCF_2015_I1351909Measurement of very forward neutron energy spectra for 7 TeV proton-proton collisions at the Large Hadron ColliderExperiment: LHCF (LHC) Inspire ID: 1351909 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
The Large Hadron Collider forward (LHCf) experiment is designed to use the LHC to verify the hadronic-interaction models used in cosmic-ray physics. Forward baryon production is one of the crucial points to understand the development of cosmic-ray showers. We report the neutron-energy spectra for LHC $\sqrt{s}$ = 7 TeV proton--proton collisions with the pseudo-rapidity $\eta$ ranging from 8.81 to 8.99, from 8.99 to 9.22, and from 10.76 to infinity. The measured energy spectra obtained from the two independent calorimeters of Arm1 and Arm2 show the same characteristic feature before unfolding the difference in the detector responses. We unfolded the measured spectra by using the multidimensional unfolding method based on Bayesian theory, and the unfolded spectra were compared with current hadronic-interaction models. The QGSJET II-03 model predicts a high neutron production rate at the highest pseudo-rapidity range similar to our results and the DPMJET 3.04 model describes our results well at the lower pseudo-rapidity ranges. However no model perfectly explains the experimental results in the whole pseudo-rapidity range. The experimental data indicate the most abundant neutron production rate relative to the photon production, which does not agree with predictions of the models. Source code: LHCF_2015_I1351909.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/Beam.hh"
5
6namespace Rivet {
7
8/// @brief Add a short analysis description here
9class LHCF_2015_I1351909 : public Analysis {
10public:
11
12 /// Constructor
13 RIVET_DEFAULT_ANALYSIS_CTOR(LHCF_2015_I1351909);
14
15 static constexpr bool lhcf_like = true;
16 static constexpr int ndecay = 1;
17 static constexpr int nbeam = 2;
18 static constexpr double D1_begin = 82000.; //mm 60000.; //mm
19 static constexpr double D1_end = 82000; //mm 90000.; //mm
20 static constexpr double IPtoLHCf = 141050.; //mm
21
22 /// @name Analysis methods
23
24 bool isParticleFromCollision(const Particles& parents, const Beam& beams) const {
25 bool beam[nbeam]={false};
26 if (parents.size()==nbeam) {
27 for ( int ipar=0; ipar < nbeam; ++ipar ) {
28 if ( parents[ipar].genParticle() == beams.beams().first.genParticle() ||
29 parents[ipar].genParticle() == beams.beams().second.genParticle() )
30 beam[ipar] = true;
31 }
32 if(beam[0] && beam[1]) return true;
33 }
34 return false;
35 }
36
37 bool isDeviated(Particle p, Particle parent) { //Select/Remove particles decayed between IP and LHCf
38 ConstGenVertexPtr pv = p.genParticle()->production_vertex();
39 assert(pv != nullptr);
40
41 const double decay_vertex = pv->position().z()/mm;
42
43 const double parent_charge = PID::charge(parent.pid());
44 const double descendant_charge = PID::charge(p.pid());
45
46 if(parent_charge == 0) { //Particles produced by neutral parent decay
47 if(descendant_charge == 0) {
48 return false;
49 } else {
50 if(decay_vertex >= D1_end)
51 return false;
52 else
53 return true; //Remove charged descendants produced from decay before end of D1
54 }
55 } else { //Particles produced by charged parent decay
56 if(decay_vertex <= D1_begin) {
57 if(descendant_charge == 0)
58 return false;
59 else
60 return true; //Remove charged descendants produced from decay before end of D1
61 } else {
62 return true; //Remove particles produced by charged parent decay after begin of D1
63 }
64 }
65
66 return false;
67 }
68
69 bool isSameParticle(Particle p1, Particle p2) {
70 if(p1.pid() == p2.pid() &&
71 mom(p1).t() == mom(p2).t() &&
72 mom(p1).x() == mom(p2).x() &&
73 mom(p1).y() == mom(p2).y() &&
74 mom(p1).z() == mom(p2).z())
75 return true;
76 else
77 return false;
78 }
79
80 bool isAlreadyProcessed(Particle p, vector<Particle> list) {
81 for(unsigned int ipar=0; ipar<list.size(); ++ipar)
82 if(isSameParticle(p, list[ipar]))
83 return true;
84
85 return false;
86 }
87
88 /// This method return a fake pseudorapidity to check id decayed particle is in LHCf acceptance
89 double RecomputeEta(Particle p) {
90 ConstGenVertexPtr pv = p.genParticle()->production_vertex();
91
92 const double x0 = pv->position().x()/mm;
93 const double y0 = pv->position().y()/mm;
94 const double z0 = pv->position().z()/mm;
95
96 const double px = p.px()/MeV;
97 const double py = p.py()/MeV;
98 const double pz = abs(p.pz()/MeV);
99
100 const double dist_to_lhcf = IPtoLHCf - z0;
101 const double x1 = x0 + (dist_to_lhcf * px/pz);
102 const double y1 = y0 + (dist_to_lhcf * py/pz);
103
104 const double r = sqrt(pow(x1, 2.)+pow(y1, 2.));
105 const double theta = atan(abs(r / IPtoLHCf));
106 const double pseudorapidity = - log (tan (theta/2.) );
107
108 return pseudorapidity;
109 }
110
111 /// Book histograms and initialise projections before the run
112 void init() {
113
114 // Initialise and register projections
115 // declare(FinalState("FS");
116 declare(FinalState(), "FS");
117 declare(Beam(), "Beams");
118
119 // Book histograms
120 book(_h_n_en_eta1, 1, 1, 1);
121 book(_h_n_en_eta2, 1, 1, 2);
122 book(_h_n_en_eta3, 1, 1, 3);
123
124 }
125
126 /// Perform the per-event analysis
127 void analyze(const Event& event) {
128
129 const FinalState &fs = apply<FinalState> (event, "FS");
130 Particles fs_particles = fs.particles();
131 const Beam & beams = apply<Beam> (event, "Beams");
132
133 vector<Particle> processed_parents;
134 processed_parents.clear();
135
136 for (Particle& p: fs_particles ) {
137
138 if(p.pz()/GeV<0.) continue;
139
140 double eta = 0.;
141 double en = 0.;
142
143 if(lhcf_like) {
144 //======================================================================
145 //========== LHCf-like analysis ========================================
146 //======================================================================
147
148 vector<Particle> parents = p.parents();
149
150 if(isParticleFromCollision(parents, beams)) { //Particles directly produced in collisions
151 if(!PID::isHadron(p.pid())) continue; //Remove non-hadron particles
152 if(PID::charge(p.pid()) != 0) continue; //Remove charged particles
153
154 eta = p.eta();
155 en = p.E()/GeV;
156 } else if(parents.size() == ndecay) { //Particles produced from decay
157 ConstGenVertexPtr pv = p.genParticle()->production_vertex();
158 assert(pv != nullptr);
159
160 const double decay_vertex = pv->position().z()/mm;
161 Particle parent = parents[0];
162
163 if(decay_vertex < IPtoLHCf) { //If decay happens before LHCf we consider descendants
164 if(!PID::isHadron(p.pid())) continue; //Remove non-hadron descendants
165 if(isDeviated(p, parent)) continue; //Remove descendants deviated by D1
166
167 eta = RecomputeEta(p);
168 en = p.E()/GeV;
169 } else {//If decay happens after LHCf we consider parents
170 vector<Particle> ancestors;
171 ancestors.clear();
172
173 int ngeneration=0;
174 bool isValid=true;
175 bool isEnded=false;
176 while(!isEnded) //Loop over all generations in the decay
177 {
178 vector<Particle> temp_part;
179 temp_part.clear();
180 if(ngeneration==0) {
181 parent = parents[0];
182 temp_part = parent.parents();
183 }
184 else {
185 parent = ancestors[0];
186 temp_part = parent.parents();
187 }
188 ancestors.clear();
189 ancestors = temp_part;
190
191 Particle ancestor = ancestors[0];
192
193 if(isParticleFromCollision(ancestors, beams)) { //if we found first particles produced in collisions we consider them
194 isEnded=true;
195
196 if(!PID::isHadron(parent.pid())) isValid=false; //Remove non-hadron ancestors/parents
197 if(PID::charge(parent.pid()) != 0) isValid=false; //Remove charged ancestors/parents
198 if(isAlreadyProcessed(parent, processed_parents))
199 isValid=false; //Remove already processed ancestors/parents when looping other descendants
200 else
201 processed_parents.push_back(parent); //Fill ancestors/parents in the list
202
203 eta = parent.eta();
204 en = parent.E()/GeV;
205 } else if (ancestors.size() == ndecay) { //if we found first particles produced entering LHCf we consider them
206 ConstGenVertexPtr pv_prev = parent.genParticle()->production_vertex();
207 assert(pv_prev != NULL);
208
209 const double previous_decay_vertex = pv_prev->position().z()/mm;
210
211 if(previous_decay_vertex < IPtoLHCf) {
212 isEnded=true;
213
214 if(!PID::isHadron(parent.pid())) isValid=false; //Remove non-hadron ancestors/parents
215 if(isDeviated(parent, ancestor)) isValid=false; //Remove ancestors/parents deviated by D1
216 if(isAlreadyProcessed(parent, processed_parents))
217 isValid=false; //Remove already processed ancestors/parents when looping other descendants
218 else
219 processed_parents.push_back(parent); //Fill ancestors/parents in the list
220
221 eta = RecomputeEta(parent);
222 en = parent.E()/GeV;
223 }
224 } else { //This condition should never happen
225 cout << "Looping over particles generation ended without match : Exit..." << endl;
226 vetoEvent;
227 }
228
229 ++ngeneration;
230 }
231
232 if(!isValid) continue;
233 }
234 } else { //This condition should never happen
235 cout << "Particle seems not to be produced in collision or decay : Exit..." << endl;
236 vetoEvent;
237 }
238
239 } else {
240 //======================================================================
241 //========== Only neutrons at IP =======================================
242 //======================================================================
243
244 vector<Particle> parents = p.parents();
245
246 if(p.pid() != 2112 ) continue;
247
248 eta = p.eta();
249 en = p.E()/GeV;
250 //}
251 }
252
253 // Fill histograms
254 if( eta > 10.76 ){
255 _h_n_en_eta1->fill( en );
256
257 }else if(eta > 8.99 && eta < 9.22){
258 _h_n_en_eta2->fill( en );
259
260 }else if(eta > 8.81 && eta < 8.99){
261 _h_n_en_eta3->fill( en );
262
263 }
264 }
265 }
266
267
268 /// Normalise histograms etc., after the run
269 void finalize() {
270
271 scale(_h_n_en_eta1, crossSection()/millibarn/sumOfWeights()); // norm to cross section
272 scale(_h_n_en_eta2, crossSection()/millibarn/sumOfWeights()); // norm to cross section
273 scale(_h_n_en_eta3, crossSection()/millibarn/sumOfWeights()); // norm to cross section
274
275 }
276
277 // @}
278
279
280private:
281
282
283 /// @name Histograms
284 // @{
285 Histo1DPtr _h_n_en_eta1;
286 Histo1DPtr _h_n_en_eta2;
287 Histo1DPtr _h_n_en_eta3;
288 // @}
289
290
291};
292
293
294RIVET_DECLARE_PLUGIN(LHCF_2015_I1351909);
295
296
297}
|