rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

LHCF_2015_I1351909

Measurement of very forward neutron energy spectra for 7 TeV proton-proton collisions at the Large Hadron Collider
Experiment: LHCF (LHC)
Inspire ID: 1351909
Status: VALIDATED
Authors:
  • Eugenio Berti
  • LHCf collaboration
References:
  • Phys. Lett. B 750 (2015) 360-366
Beams: p+ p+
Beam energies: (3500.0, 3500.0) GeV
Run details:
  • Differential production cross section of neutrons in p-p collisions in the very forward region expressed as a function of energy. Note that, being LHCf detector locate 140m away from IP, two additional effects must be taken into account. They are particles decay in the transport though the beam pipe and trajectories bending due to the dipole magnet. Because of them, the final energy spectra included about $0-6 \%$ of other hadrons, mainly $\Lambda^{0}$ and $K^{0}$. These effects are considered in the Rivet code making use of some approximations that are able to reproduce the model distributions shown in the paper within 15% in most cases, depending on the model, the pseudorapidity region and the energy bin.

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}