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 Particle& p, 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 isParticleFromDecay(const Particle p, const Particles& parents) const {
38 return (parents.size() == ndecay);
39 }
40
41 bool isDeviated(Particle p, Particle parent) { //Select/Remove particles decayed between IP and LHCf
42 ConstGenVertexPtr pv = p.genParticle()->production_vertex();
43 assert(pv != nullptr);
44
45 const double decay_vertex = pv->position().z()/mm;
46
47 const double parent_charge = PID::charge(parent.pid());
48 const double descendant_charge = PID::charge(p.pid());
49
50 if(parent_charge == 0) { //Particles produced by neutral parent decay
51 if(descendant_charge == 0) {
52 return false;
53 } else {
54 if(decay_vertex >= D1_end)
55 return false;
56 else
57 return true; //Remove charged descendants produced from decay before end of D1
58 }
59 } else { //Particles produced by charged parent decay
60 if(decay_vertex <= D1_begin) {
61 if(descendant_charge == 0)
62 return false;
63 else
64 return true; //Remove charged descendants produced from decay before end of D1
65 } else {
66 return true; //Remove particles produced by charged parent decay after begin of D1
67 }
68 }
69
70 return false;
71 }
72
73 bool isSameParticle(Particle p1, Particle p2) {
74 if(p1.pid() == p2.pid() &&
75 mom(p1).t() == mom(p2).t() &&
76 mom(p1).x() == mom(p2).x() &&
77 mom(p1).y() == mom(p2).y() &&
78 mom(p1).z() == mom(p2).z())
79 return true;
80 else
81 return false;
82 }
83
84 bool isAlreadyProcessed(Particle p, vector<Particle> list) {
85 for(unsigned int ipar=0; ipar<list.size(); ++ipar)
86 if(isSameParticle(p, list[ipar]))
87 return true;
88
89 return false;
90 }
91
92 /// This method return a fake pseudorapidity to check id decayed particle is in LHCf acceptance
93 double RecomputeEta(Particle p) {
94 ConstGenVertexPtr pv = p.genParticle()->production_vertex();
95
96 const double x0 = pv->position().x()/mm;
97 const double y0 = pv->position().y()/mm;
98 const double z0 = pv->position().z()/mm;
99
100 const double px = p.px()/MeV;
101 const double py = p.py()/MeV;
102 const double pz = abs(p.pz()/MeV);
103
104 const double dist_to_lhcf = IPtoLHCf - z0;
105 const double x1 = x0 + (dist_to_lhcf * px/pz);
106 const double y1 = y0 + (dist_to_lhcf * py/pz);
107
108 const double r = sqrt(pow(x1, 2.)+pow(y1, 2.));
109 const double theta = atan(abs(r / IPtoLHCf));
110 const double pseudorapidity = - log (tan (theta/2.) );
111
112 return pseudorapidity;
113 }
114
115 /// Book histograms and initialise projections before the run
116 void init() {
117
118 // Initialise and register projections
119 // declare(FinalState("FS");
120 declare(FinalState(), "FS");
121 declare(Beam(), "Beams");
122
123 // Book histograms
124 book(_h_n_en_eta1, 1, 1, 1);
125 book(_h_n_en_eta2, 1, 1, 2);
126 book(_h_n_en_eta3, 1, 1, 3);
127
128 }
129
130 /// Perform the per-event analysis
131 void analyze(const Event& event) {
132
133 const FinalState &fs = applyProjection<FinalState> (event, "FS");
134 Particles fs_particles = fs.particles();
135 const Beam & beams = applyProjection<Beam> (event, "Beams");
136
137 vector<Particle> processed_parents;
138 processed_parents.clear();
139
140 for (Particle& p: fs_particles ) {
141
142 if(p.pz()/GeV<0.) continue;
143
144 double eta = 0.;
145 double en = 0.;
146
147 if(lhcf_like) {
148 //======================================================================
149 //========== LHCf-like analysis ========================================
150 //======================================================================
151
152 vector<Particle> parents = p.parents();
153
154 if(isParticleFromCollision(p, parents, beams)) { //Particles directly produced in collisions
155 if(!PID::isHadron(p.pid())) continue; //Remove non-hadron particles
156 if(PID::charge(p.pid()) != 0) continue; //Remove charged particles
157
158 eta = p.eta();
159 en = p.E()/GeV;
160 } else if(isParticleFromDecay(p, parents)) { //Particles produced from decay
161 ConstGenVertexPtr pv = p.genParticle()->production_vertex();
162 assert(pv != nullptr);
163
164 const double decay_vertex = pv->position().z()/mm;
165 Particle parent = parents[0];
166
167 if(decay_vertex < IPtoLHCf) { //If decay happens before LHCf we consider descendants
168 if(!PID::isHadron(p.pid())) continue; //Remove non-hadron descendants
169 if(isDeviated(p, parent)) continue; //Remove descendants deviated by D1
170
171 eta = RecomputeEta(p);
172 en = p.E()/GeV;
173 } else {//If decay happens after LHCf we consider parents
174 vector<Particle> ancestors;
175 ancestors.clear();
176
177 int ngeneration=0;
178 bool isValid=true;
179 bool isEnded=false;
180 while(!isEnded) //Loop over all generations in the decay
181 {
182 vector<Particle> temp_part;
183 temp_part.clear();
184 if(ngeneration==0) {
185 parent = parents[0];
186 temp_part = parent.parents();
187 }
188 else {
189 parent = ancestors[0];
190 temp_part = parent.parents();
191 }
192 ancestors.clear();
193 ancestors = temp_part;
194
195 Particle ancestor = ancestors[0];
196
197 if(isParticleFromCollision(parent, ancestors, beams)) { //if we found first particles produced in collisions we consider them
198 isEnded=true;
199
200 if(!PID::isHadron(parent.pid())) isValid=false; //Remove non-hadron ancestors/parents
201 if(PID::charge(parent.pid()) != 0) isValid=false; //Remove charged ancestors/parents
202 if(isAlreadyProcessed(parent, processed_parents))
203 isValid=false; //Remove already processed ancestors/parents when looping other descendants
204 else
205 processed_parents.push_back(parent); //Fill ancestors/parents in the list
206
207 eta = parent.eta();
208 en = parent.E()/GeV;
209 } else if (isParticleFromDecay(parent, ancestors)) { //if we found first particles produced entering LHCf we consider them
210 ConstGenVertexPtr pv_prev = parent.genParticle()->production_vertex();
211 assert(pv_prev != NULL);
212
213 const double previous_decay_vertex = pv_prev->position().z()/mm;
214
215 if(previous_decay_vertex < IPtoLHCf) {
216 isEnded=true;
217
218 if(!PID::isHadron(parent.pid())) isValid=false; //Remove non-hadron ancestors/parents
219 if(isDeviated(parent, ancestor)) isValid=false; //Remove ancestors/parents deviated by D1
220 if(isAlreadyProcessed(parent, processed_parents))
221 isValid=false; //Remove already processed ancestors/parents when looping other descendants
222 else
223 processed_parents.push_back(parent); //Fill ancestors/parents in the list
224
225 eta = RecomputeEta(parent);
226 en = parent.E()/GeV;
227 }
228 } else { //This condition should never happen
229 cout << "Looping over particles generation ended without match : Exit..." << endl;
230 vetoEvent;
231 }
232
233 ++ngeneration;
234 }
235
236 if(!isValid) continue;
237 }
238 } else { //This condition should never happen
239 cout << "Particle seems not to be produced in collision or decay : Exit..." << endl;
240 vetoEvent;
241 }
242
243 } else {
244 //======================================================================
245 //========== Only neutrons at IP =======================================
246 //======================================================================
247
248 vector<Particle> parents = p.parents();
249
250 //if(isParticleFromCollision(p, parents)) { //Particles directly produced in collisions
251 if(p.pid() != 2112 ) continue;
252
253 eta = p.eta();
254 en = p.E()/GeV;
255 //}
256 }
257
258 // Fill histograms
259 if( eta > 10.76 ){
260 _h_n_en_eta1->fill( en );
261
262 }else if(eta > 8.99 && eta < 9.22){
263 _h_n_en_eta2->fill( en );
264
265 }else if(eta > 8.81 && eta < 8.99){
266 _h_n_en_eta3->fill( en );
267
268 }
269 }
270 }
271
272
273 /// Normalise histograms etc., after the run
274 void finalize() {
275
276 scale(_h_n_en_eta1, crossSection()/millibarn/sumOfWeights()); // norm to cross section
277 scale(_h_n_en_eta2, crossSection()/millibarn/sumOfWeights()); // norm to cross section
278 scale(_h_n_en_eta3, crossSection()/millibarn/sumOfWeights()); // norm to cross section
279
280 }
281
282 //@}
283
284
285private:
286
287
288 /// @name Histograms
289 //@{
290 Histo1DPtr _h_n_en_eta1;
291 Histo1DPtr _h_n_en_eta2;
292 Histo1DPtr _h_n_en_eta3;
293 //@}
294
295
296};
297
298
299// The hook for the plugin system
300RIVET_DECLARE_PLUGIN(LHCF_2015_I1351909);
301
302
303}
|