Measurement of forward neutral pion transverse momentum spectra for $\sqrt{s}$ = 7TeV proton-proton collisions at LHC

The inclusive production rate of neutral pions in the rapidity range greater than $y=8.9$ has been measured by the Large Hadron Collider forward (LHCf) experiment during LHC $\sqrt{s}=7$\,TeV proton-proton collision operation in early 2010. This paper presents the transverse momentum spectra of the neutral pions. The spectra from two independent LHCf detectors are consistent with each other and serve as a cross check of the data. The transverse momentum spectra are also compared with the predictions of several hadronic interaction models that are often used for high energy particle physics and for modeling ultra-high-energy cosmic-ray showers.

The inclusive production rate of neutral pions in the rapidity range greater than y = 8.9 has been measured by the Large Hadron Collider forward (LHCf) experiment during √ s = 7 TeV protonproton collision operation in early 2010. This paper presents the transverse momentum spectra of the neutral pions. The spectra from two independent LHCf detectors are consistent with each other and serve as a cross-check of the data. The transverse momentum spectra are also compared with the predictions of several hadronic interaction models that are often used for high-energy particle physics and for modeling ultra-high-energy cosmic-ray showers.

I. INTRODUCTION
One of the important tasks of strong-interaction physics described by Quantum Chromodynamics (QCD) is to provide a detailed understanding of forward particle production in hadronic interactions. QCD involves two types of limiting processes: "hard" and "soft".
Hard processes occur in the range characterized by a large four-momentum transfer t, where |t| should be larger than 1 GeV 2 . Note that units used in this report are c = k (Boltzmann constant) = 1. Deep inelastic scattering that is accompanied by the exchange of virtual photons or vector bosons, or jets produced by large transverse momentum (p T ) partons are typical phenomena that are categorized as hard processes. The hard processes have been successfully described by perturbation theory, owing to the asymptotic freedom of QCD at high energy.
On the other hand, soft processes occur when the fourmomentum transfer |t| is smaller than 1 GeV 2 . These processes, which correspond to a large impact parameter, have a large QCD coupling constant and cannot be calculated by perturbative QCD. Gribov-Regge theory is applicable for describing soft processes [1,2] and the Pomeron contribution, as a component of the Gribov-Regge approach to high-energy hadronic interactions, increases with increasing energy [3] and should dominate at the TeV energy scale. However there still exists a problem for the theories that involve these virtual quasiparticles. Since the treatment of the Pomeron differs amongst the model theories they predict different results for particle production. Thus a deeper understanding of soft processes is needed and soft processes are mostly equivalent to forward or large rapidity particle production in hadronic interactions. However experimental data for large rapidity are meager. Moreover the experimental data that do exist have so far been carried out at relatively low energy, for example ISR [4] at √ s = 53 GeV and UA7 [5] at √ s = 630 GeV. The Large Hadron Collider forward (LHCf) experiment [6] has been designed to measure the hadronic production cross sections of neutral particles emitted in very forward angles in proton-proton collisions at the LHC, including zero degrees. The LHCf detectors have the capability for precise measurements of forward high-energy inclusive-particle-production cross sections of photons, neutrons, and possibly other neutral mesons and baryons. Among the many secondary neutral particles that LHCf can detect, the π 0 mesons are the most sensitive to the details of the proton-proton interactions. Thus a high priority has been given to analyzing forward π 0 production data in order to provide key information for an as yet un-established hadronic interaction theory at the TeV energy scale. The analysis in this paper concentrates on obtaining the inclusive production rate for π 0 s in the rapidity range larger than y = 8.9 as a function of the π 0 transverse momentum.
In addition to the aim described above, this work is also motivated by an application to the understanding of Ultra-High-Energy Cosmic-Ray (UHECR) phenomena, which are sensitive to the details of soft π 0 production at extreme energy. It is known that the lack of knowledge about forward particle production in hadronic collisions hinders the interpretation of observations of UHECR [7,8]. Although UHECR observations have made notable advances in the last few years [9][10][11][12][13][14][15], critical parts of the analysis depend on Monte Carlo (MC) simulations of air shower development that are sensitive to the choice of the hadronic interaction model. It should also be remarked that the LHC has reached 7 TeV collision energy, which in the laboratory frame of UHECR observations is equivalent to 2.6 × 10 16 eV, and this energy is above the "knee" region of the primary cosmic ray energy spectrum (∼ 4 × 10 15 eV) [16]. The data provided by LHCf should then provide a useful bench mark for the MC codes that are used for the simulation of UHECR atmospheric showers.
This paper is organized as follows. In Sec. II the LHCf detectors are described. Sec. III summarizes the conditions for taking data and the MC simulation methodology. In Sec. IV the analysis framework is described. The factors that contribute to the systematic uncertainty of the results are explained in Sec. V and the analysis results are then presented in Sec. VI. Sec. VII discusses the results that have been obtained and compare these with the predictions of several hadronic interaction models. Finally, concluding remarks are found in Sec. VIII.

II. THE LHCF DETECTORS
Two independent LHCf detectors, called Arm1 and Arm2, have been installed in the instrumentation slots of the target neutral absorbers (TANs) [17] located ±140 m from the ATLAS interaction point (IP1) and at zero degree collision angle. Fig. 1 shows schematic views of the Arm1 (left) and Arm2 (right) detectors. Inside a TAN the beam-vacuum-chamber makes a Y-shaped transition from a single common beam tube facing IP1 to two separate beam tubes joining to the arcs of the LHC. Charged particles produced at IP1 and directed towards the TAN are swept aside by the inner beam separation dipole magnet D1 before reaching the TAN. Consequently only neutral particles produced at IP1 enter the LHCf detector. At this location the LHCf detectors cover the pseudorapidity range from 8.7 to infinity for zero degree beam crossing angle. With a maximum beam crossing angle of 140 µrad, the pseudorapidity range can be extended to 8.4 to infinity.
Each LHCf detector has two sampling and imaging calorimeters composed of 44 radiation lengths (X 0 ) of tungsten and 16 sampling layers of 3 mm thick plastic scintillator. The transverse sizes of the calorimeters are 20 ×20 mm 2 and 40 ×40 mm 2 in Arm1, and 25 ×25 mm 2 and 32 ×32 mm 2 in Arm2. The smaller calorimeters cover zero degree collision angle. Four X-Y layers of position sensitive detectors are interleaved with the layers of tungsten and scintillator in order to provide the transverse positions of the showers. Scintillating fiber (SciFi) belts are used for the Arm1 position sensitive layers and silicon micro-strip sensors are used for Arm2. Readout pitches are 1 mm and 0.16 mm for Arm1 and Arm2, respectively.
More detail on the scientific goals and the construction and performance of the detectors can be found in previous reports [18][19][20][21][22].  The experimental data used for the analysis of this paper were obtained on May 15th and 16th 2010 during proton-proton collisions at √ s=7 TeV with zero degree beam crossing angle (LHC Fill 1104). Data taking was carried out in two different runs: the first run was on May 15th from 17:45 to 21:23, and the second run was on May 16th from 00:47 to 14:05. The events that were recorded during a luminosity optimization scan and a calibration run were removed from the data set for this analysis. The range of total luminosity of the three crossing bunch pairs was L = (6.3−6.5)×10 28 cm −2 s −1 for the first run and L = (4.8−5.9)×10 28 cm −2 s −1 for the second run. These ranges of luminosity were ideal for the LHCf data acquisition system. The integrated luminosities for the data analysis reported in this paper were derived from the counting rate of the LHCf Front Counters [23], and were 2.53 nb −1 (Arm1) and 1.90 nb −1 (Arm2) after taking the live time percentages into account. The average live time percentages for the first/second run were 85.7 %/81.1 % for Arm1 and 67.0 %/59.7 % for Arm2. The live time percentages for the second run were smaller than the first run owing to a difference in the trigger schemes. In both runs the trigger efficiency achieved was >99% for photons with energy E > 100 GeV [24].
The events containing more than one collision in a single bunch crossing (pile-up events) could potentially cause a bias in the p T spectra. For example combinatorial single-hits from different collisions within a single bunch crossing might be identified as multi-hit events from a single collision and removed from the analysis. (Multi-hit events have two showers in a single calorimeter and are eliminated from the data analysis. The production rates are later corrected for this cut. See Fig. 2 and related discussion.) However it can be shown that pile-up events are negligible for the LHCf data taking conditions of this report. Given that a collision has occurred, the probability of pile-up (P pileup ) is calculated from the Poisson probability distribution for n collisions P poi (n) according to P pileup = P poi (n ≥ 2)/P poi (n ≥ 1). With the highest bunch luminosity L = 2.3 × 10 28 cm −2 s −1 used in this analysis, an inelastic cross section σ inel = 73.6 mb and the revolution frequency of LHC f rev = 11.2 kHz, the pile-up probability is P pileup ∼ 0.07. However considering that the acceptance of the LHCf calorimeter for inelastic collisions is ∼0.03, only 0.2% of events have more than one shower event in a single calorimeter due to pile-up and this is negligible.
Detailed discussions of pile-up effects and background events from collisions between the beam and residual gas molecules in the beam tube can be found in previous reports [6,24].

B. Methodology for performing Monte Carlo simulations
MC simulation consists of three steps: (1) protonproton interaction event generation at IP1, (2) transport from IP1 to the LHCf detectors and (3) the response of the LHCf detectors.
Proton-proton interaction events at √ s = 7 TeV and the resulting flux of secondary particles and their kinematics are simulated with cosmos (version 8.81). cosmos acts as the front end for the external hadronic interaction models (qgsjet II-03 [25], dpmjet 3.04 [26], sibyll 2.1 [27] and epos 1.99 [28]) that describe the proton-proton interactions. While pythia 8.145 [29,30] serves as its own front end for the generation of protonproton interaction events.
Next, the generated secondary particles are transported in the beam pipe from IP1 to the TAN, taking account of the deflection of charged particles by the Q1 quadrupole and D1 beam separation dipole, particle decay, and particle interaction with the beam pipe and the Y-shaped beam-vacuum-chamber transition made of copper (1 X 0 projected thickness in front of the LHCf detectors). Charged particles are swept away by the D1 magnet before reaching the LHCf detectors. This simulation uses the epics library [31] (version 7.49) and a part of cosmos. epics deals with the transport of secondary particles. Particle interactions with the residual gas molecules inside the beam pipe are not simulated.
Contamination from beam-gas background events in the data set used for analysis is estimated to be only ∼ 0.1 % and has no significant impact on the p T spectra reported.
Finally the simulations of the showers produced in the LHCf detectors and their response are carried out for the particles arriving at the TAN using the cosmos and epics libraries. The survey data for detector position and random fluctuations equivalent to electrical noise are taken into account in this step. The Landau-Pomeranchuk-Migdal effect [32,33] that longitudinally lengthens an electromagnetic shower at high energy is also considered. A change of the p T spectra caused by LPM effects is only at the 1% level since the reconstruction of energy deposited in the calorimeters is carried out to a sufficiently deep layer whereby the energy of electromagnetic showers is almost perfectly deposited within the calorimeter.
The simulations of the LHCf detectors are tuned to test beam data taken at the CERN SPS in 2007 [20]. The validity of the detector simulation was checked by comparing the shower development and deposited energy for each calorimeter layer to the results obtained by the fluka library [34].
In order to validate the reconstruction algorithms and to estimate a possible reconstruction bias beyond the energy range of the SPS test beam results, the MC simulations are generated for 1.0 × 10 8 inelastic collisions, where the secondary particles are generated by the epos 1.99 [28] hadronic interaction model. This MC simulation is referred to as the "reference MC simulation" in the following text.
Similarly the "toy MC simulations" discussed below are performed in order to determine various correction factors to use in the event reconstruction processes. In the toy MC simulations, a single-photon with a given fixed energy is directed at the LHCf detectors.

A. Event reconstruction and selection
Observation of π 0 mesons by a LHCf detector is illustrated in Fig. 2. The π 0 s are identified by their decay into two photons. Since the π 0 s decay very close to their point of creation at IP1, the opening angle (θ) between the two photons is the transverse distance between photon impact points at the LHCf detectors divided by the distance from IP1 (z = ±141.05 m). Consequently the opening angle for the photons from π 0 decay that are detected by a LHCf detector is constrained by θ 0.4 mrad for Arm1 and θ 0.6 mrad for Arm2. Other kinematic variables of the π 0 s (energy, p T , and rapidity) are also reconstructed by using the photon energy and incident position measured by each calorimeter. Note that for the analysis of this paper events having two photons entering the same calorimeter (multi-hit events) are rejected (right panel of Fig. 2). The accuracy of energy recon-struction for such events is still under investigation. The final inclusive production rates reported in this paper are corrected for this cut. In order to ensure good event reconstruction efficiency, the range of the π 0 rapidity and p T are limited to 8.9 < y < 11.0 and p T < 0.6 GeV, respectively. All particles other than photons from π 0 decay are ignored in this analysis. Thus, also according to the multi-hit π 0 correction described in detail in Sec. IV F, the reported production rates are inclusive. The standard reconstruction algorithms are described in this section and systematic uncertainties will be discussed in Sec. V.

Hit position reconstruction
The transverse impact positions of particles entering the calorimeters are determined using the information provided by the position sensitive layers. In this analysis, the transverse impact position of the core of an electromagnetic shower is taken from the position of the peak signal on the position sensitive layer that has the largest energy deposited amongst all the position sensitive layers.
Hit positions that fall within 2 mm of the edges of the calorimeters are removed from analysis due to the large uncertainty in the energy determination of such events owing to shower leakage. For the toy MC simulations, the position reconstruction resolution is defined as the one standard deviation difference between the true primary photon position and the reconstructed position of the shower axis. The estimated resolution using the toy MC simulations and test beam data for a single photon with energy E > 100 GeV is better than 200 µm and 100 µm for Arm1 and Arm2, respectively [21,35].
Multi-hit events defined to have more than one photon registered in a single calorimeter are eliminated from the analysis in this paper. Multi-hit candidates that have two distinct peaks in the lateral shower impact distribution are searched for using the algorithm that has been implemented in the TSpectrum [36] class in root [37]. When the separation between peaks is greater than 1 mm and the lower energy photon has more than 5 % of the en-ergy of the nearby photon, the MC simulation estimated efficiencies for identifying multi-hit events are larger than 70 % and 90 % for Arm1 and Arm2, respectively [24]. The efficiency for Arm2 is better than that for Arm1 owing to the finer readout pitches of the silicon micro-strip sensors. The subtraction of the remaining contamination by multi-hit events is discussed in Sec. IV C.
On the other hand for single-hit events not having two identifiable peaks, the MC simulation estimated efficiency for correctly identifying true single photon events with energy E > 100 GeV is better than 98 % both for Arm1 and Arm2, although the precise percentage depends slightly on the photon energy.

Energy reconstruction
The charge information in each scintillation layer is converted to a deposited energy by using calibration factors obtained from the SPS electron test beam data taken below 200 GeV [19]. In this analysis the deposited energy is scaled to the number of minimum ionizing shower particles with a coefficient 1 MIP = 0.453 MeV that corresponds to the most probable deposited energy by a 150 GeV muon passing through a 3 mm thick plastic scintillator. The sum of the energy deposited in the 2 nd to 13 th scintillation layers (dE [MIP]) is then converted to the primary photon energy E[GeV] using a polynomial function The are determined from the response of the calorimeters to single photons by the toy MC simulations. The validity of this method has been confirmed with the SPS beam tests. The MC estimated energy resolution for single photons above 100 GeV considering the LHC data taking situation is given by the expression Corrections for shower leakage effects [18,19] are carried out during the energy reconstruction process. Corrections are applied for leakage of particles out of the calorimeters and for leakage of particles in that have escaped from the adjacent calorimeter. Both of the leakage effects depend on the transverse location of the shower axis in the calorimeters. The correction factors have been estimated from the toy MC simulations. The light-yield collection efficiency of the plastic scintillation layers [19] is also a function of the transverse location of the shower axis and corrected for in this step.
Events having a reconstructed energy below 100 GeV are eliminated from the analysis firstly to reject particles produced by interaction of collision products with the beam pipe, and secondly to avoid errors due to trigger inefficiency (see Sec. III A).

Particle identification
The particle identification (PID) process is applied in order to select pure electromagnetic showers, specifically photons from π 0 decay, and to reduce hadron contamination, specifically from neutrons. A parameter L 90% is defined for this purpose. L 90% is the longitudinal distance, in units of radiation length, measured from the 1st tungsten layer of a calorimeter to the position where the energy deposition integral reaches 90 % of the total shower energy deposition. Fig. 3 shows the distribution of L 90% for the 20 mm calorimeter of the Arm1 detector for events having a reconstructed energy in the range 500 GeV< E <1 TeV. Experimental data (black dots) and the MC simulations based on qgsjet II-03 (shaded areas) are shown. The normalization factors of pure photon and pure hadron incident events are modified to get the best agreement between the L 90% distributions of the experimental data and the MC simulations. The best agreement is obtained by a chi-square test of the L 90% distribution of the experimental data relative to the MC simulation. The two distinct peaks correspond to photon (L 90% 20 X 0 ) and hadron (L 90% 20 X 0 ) events.
PID criteria that depend on the energy of the individual photons are defined in terms of the L 90% distribution in order to keep the π 0 selection efficiency at approximately 90 % over the entire p T range. These criteria f L90% (E 1 , E 2 ) are expressed as a function of the photon energies measured by the small (E 1 ) and large (E 2 ) calorimeters and have been determined by the toy MC simulations for each Arm. The remaining hadron contamination is removed by background subtraction introduced in Sec. IV C. The unavoidable selection inefficiency of 10 % is corrected for in the unfolding process to be discussed later (Sec. IV D). Table I summarizes the π 0 event selection criteria that are applied prior to reconstruction of the π 0 kinematics.  Candidates for π 0 events are selected using the characteristic peak in the two-photon invariant mass distribution corresponding to the π 0 rest mass. Reconstruction of the invariant mass m γγ is done using the incident positions and energies information of the photon pair, where q i and E i are the energy-momentum 4-vectors and energies of the decay photons in the laboratory frame, respectively. θ is the opening angle between the two photons in the laboratory frame. The last approximation in Eq. (3) is valid since the π 0 s decay very close to IP1 (mean π 0 flight path 1 mm). This approximation and the reconstruction algorithm for π 0 events have been verified by analysis of the reference MC simulations of the energy, rapidity and p T of the π 0 s. The reconstructed invariant mass is concentrated near peaks at 135.2±0.2 MeV in Arm1 and 134.8±0.2 MeV in Arm2, thus reproducing the π 0 mass. The uncertainties given for the mass peaks are statistical only.
It should be noted however that in the π 0 analysis of the experimental LHCf data energy scale corrections are needed so the π 0 mass peaks for Arm1 and Arm2 occur at the proper value. With no energy scale corrections applied to the LHCf data, the reconstructed invariant mass peaks using gain calibration constants determined by test beam data occur at 145.8±0.1 MeV (Arm1) and 139.9±0.1 MeV (Arm2). Therefore energy scale corrections of −8.1 % (Arm1) and −3.8 % (Arm2) applied to the raw measured photon energies are needed to bring the reconstructed π 0 rest mass into agreement with the world averaged π 0 rest mass [38]. The cause of these energy scale corrections is probably due to a temperature dependent shift of PMT gain. However at this point the temperature dependent shift of PMT gain is only qualitatively understood. Note that the typical uncertainty in opening angle is estimated to be less than 1 % relative to the reconstructed invariant mass by the position determination resolution and the alignment of the position sensitive detectors.

C. Background subtraction
Background contamination of two-photon π 0 events by hadron events and the accidental coincidence of two pho-tons not coming from the decay of a single π 0 are subtracted using the so called "sideband" method. Fig. 4 shows an example of the reconstructed twophoton invariant mass distribution of the experimental data of Arm1 in the rapidity range from 9.0 to 9.2. The energy scale correction discussed in the previous section has been applied. The sharp peak around 135 MeV is due to π 0 events. The solid curve represents the best-fit of a composite physics model to the invariant mass distribution of the data. The model consists of an asymmetric Gaussian distribution (also known as a bifurcated Gaussian distribution) for the signal component and a 3rd order Chebyshev polynomial function for the background component. The dashed curve indicates the background component.
Using the expected mean (m) and 1 σ deviations (σ l for lower side and σ u for upper side) of the signal component, the signal window is defined as the invariant mass region within the two solid arrows shown in Fig. 4, where the lower and upper limits are given bym − 3σ l andm + 3σ u , respectively. The background window is constructed from the two sideband regions, [m − 6σ l ,m − 3σ l ] and [m+3σ u ,m+6σ u ], that are defined as the invariant mass regions within the dashed arrows in Fig. 4.
The rapidity and p T distributions of the signal (f (y, p T ) Sig ) are then obtained by subtracting the background distribution (f (y, p T ) BG ), estimated by the background window, from the signal-rich distribution (f (y, p T ) Sig+BG ) selected from the signal window. The fraction of the background component included in the signal window can be estimated using the likelihood function (L BG (y, p T , m γγ )) characterized by the best-fit 3rd order Chebyshev polynomial function. For simplicity, L BG (y, p T , m γγ ) is shortened as L BG in the following text. Thus the signal distribution with background subtracted is given by where R(y, p T ,m, σ l , σ u ) is the normalization for the background distribution and written as

D. Unfolding of spectra
The raw rapidity -p T distributions must be corrected for unavoidable reconstruction inefficiency and for the smearing caused by finite position and energy resolutions. An iterative Bayesian method [39,40] is used to simultaneously correct for both effects. The advantages of an iterative Bayesian method with respect to other unfolding algorithms are discussed in another report [39]. The unfolding procedure for the data is organized as follows. First, the response of the LHCf detectors to single π 0 events is simulated by toy MC calculations. In the toy MC simulations, two photons from the decay of π 0 s and low energy background particles such as those originating in a prompt photon event or a beam-pipe interaction are traced through the detector and then reconstructed with the event reconstruction algorithm introduced above. Note that the single π 0 kinematics that are simulated within the allowed phase space are independent of the particular interaction model that is being used. The background particles are simulated by a hadronic interaction model which is discussed later, since the amount of background particles is not directly measured by the LHCf detector.
The detector response to π 0 events depends on rapidity and p T , since the performance of the particle identification algorithm and the selection efficiency of events with a single photon hit in both calorimeters depend upon the energy and the incident position of a particle. The reconstructed rapidity -p T distributions for given true rapidity -p T distributions then lead to the calculation of the response function. Then the reconstructed rapidity and p T spectra are corrected with the response function which is equivalent to the likelihood function in Bayes' theorem. The corrections are carried out iteratively whereby the starting point of the current iteration is the ending point of the previous iteration. Statistical uncertainty is also propagated from the first iteration to the last. Iteration is stopped at or before the 4th iteration to obtain a regularization of the unfolded events.
Validation of the unfolding procedure is checked by applying the response function to the reference MC simulation samples. The default response function is determined with two photons from π 0 decay and the low en-ergy (E < 100 GeV) background particles generated by epos 1.99. Validity of the choice of epos 1.99 is tested by comparing two corrected spectra, one generated by epos 1.99 and another by pythia 8.145. No statistically significant difference between the corrected spectra is found. A chi-square test of the corrected spectra based on the default response function against the true spectra ensures the chi-square probability is greater than 60 %. Thus it is concluded that with the background subtraction and unfolding methods used in this analysis there is no significant bias and the statistical uncertainty is correctly quoted. Accordingly no systematic uncertainty related to the choice of the hadronic interaction models for the reference MC simulations is considered in the analysis that follows.

E. Acceptance and branching ratio correction
The apertures of the LHCf calorimeters do not cover the full 2π azimuthal angle over the entire rapidity range that is sampled. A correction for this is applied to the data before it is compared with theoretical expectations.
The correction is done using the rapidity -p T phase space. Correction coefficients are determined as follows. First, using a toy MC simulation, a single π 0 is generated at IP1 and the decay photons are propagated to the LHCf detectors. The energy-momentum 4-vectors of the π 0 s are randomly chosen so that they cover the rapidity range that the LHCf detectors are able to measure. The beam pipe shadow on the calorimeter and the actual detector positions are taken into account using survey data.
Next fiducial area cuts in the transverse X-Y plane are applied to eliminate particles that do not fall within the acceptance of the calorimeters. In the fiducial area cuts, a systematic shift of the proton beam axis is applied according to the reconstruction of the beam-axis during LHC operation. In addition a cut is applied to eliminate photons with energy less than 100 GeV. This corresponds to the treatment of the actual data for reducing the background contamination by particle interactions with the beam pipe.
Finally two phase space distributions of π 0 s are produced; one is for all π 0 s generated at IP1 and the other is for π 0 s accepted by the calorimeters. The ratio of the distribution of accepted π 0 s divided by the distribution of all π 0 s is then the geometrical acceptance efficiency. Fig. 5 shows the acceptance efficiency as a function of the π 0 rapidity and p T and dashed curves indicate lines of constant π 0 energy, E = 1 TeV, 2 TeV and 3 TeV. The left and right panels indicate the acceptance efficiency for Arm1 and Arm2, respectively. The final rapidity and p T spectra are obtained by applying the acceptance map shown in Fig. 5 to the acceptance uncorrected data. Note that the correction maps in Fig. 5 are purely kinematic and do not depend upon the particular hadronic interaction model that has been used. The uncertainty of the acceptance map caused by the finite statistics of the MC simulations is negligible.
The branching ratio of π 0 decay into two photons is 98.8 % and then inefficiency due to π 0 decay into channels other than two photons (1.2 %) is taken into account by increasing the acceptance efficiency in rapidity -p T phase space by 1.2 % everywhere and is independent of the particular hadronic interaction model.  The detected events have been classified into two types of events: single-hit π 0 and multi-hit π 0 events. The former class consists of two photons, one in each of the calorimeters of an Arm1 or Arm2 detector. A multihit π 0 event is defined as a single π 0 accompanied with at least one additional background particle (photon or neutron) in one of the calorimeters. In this analysis, only single-hit π 0 events are considered, and multi-hit π 0 events are rejected in the single-hit selection process (Sec. IV A 1) when the energy of the additional background particle is beyond the energy threshold of the cut.
The loss of multi-hit π 0 events is corrected for with the help of event generators. A range of ratios of multi-hit plus single-hit to single-hit π 0 events is estimated using several hadronic interaction models in each rapidity range. The observed p T spectra are then multiplied by the mean of these ratios and also contribute a systematic uncertainty corresponding to the variation among the interaction models. In this way the single-hit π 0 spectra are corrected so they represent inclusive π 0 production spectra. The p T dependent range of the flux of multi-hit π 0 events has been estimated using qgsjet II-03, dpmjet 3.04, sibyll 2.1, epos 1.99 and pythia 8.145, and resulted in a range of 0 %-10 % of the flux of single-hit π 0 events.

A. Energy scale
The known rest mass of the π 0 s is 134.9766 ± 0.0006 MeV [38] whereas the peak of the two-photon invariant mass measured by the two LHCf detectors occurs at 145.8±0.1 MeV (Arm1) and 139.9±0.1 MeV (Arm2) where the ± 0.1 MeV uncertainties are statistical. The mass excess error is +8.1 % for Arm1 and +3.8 % for Arm2. According to Eq. 3 there are two possible sources for mass excess error; (1) systematic over estimates of the energies E 1 and E 2 of the two decay photons and (2) systematic over estimate of the opening angle between the two photons. As discussed in Sec. IV B the typical uncertainty in opening angle is less than 1 %, too small to explain the observed mass excesses. This leaves measurement of the photon energies as the source of mass excess error.
The uncertainty in measurement of photon energy has also been investigated in a beam test at SPS and calibration with a radiation source. The estimated uncertainty of photon energy from these tests is 3.5 %. The 3.5 % uncertainty is dominated by the uncertainties in factors converting measured charge to deposited energy [20]. Note that the linearity of each PMT was carefully tested before detector assembly over a wide range of signal amplitude by exciting the scintillator with a 337 nm UV laser pulse [6,19]. The difference of reconstructed energy between the reconstruction algorithm with and without non-linearity correction of PMTs for 3 TeV photons is only 0.5 % at maximum, nevertheless the measured nonlinear response functions have been applied in the analysis.
The systematic uncertainties estimated by the beam test data at SPS (3.5 % for both Arms) are considered as uncorrelated among the p T bins, while the systematic uncertainties owing to the mass excess errors (8.1 % for Arm1 and 3.8 % for Arm2) are considered as correlated between each p T bin. The systematic shift of bin contents due to the energy scale uncertainties is estimated using two energy spectra by artificially scaling the energy with the two extremes. The ratios of the two extreme spectra to the non-scaled spectrum are assigned as systematic shifts in each bin.

B. Particle identification
The L 90% distribution described in Sec. IV A 3 is used to select LHCf π 0 events for the p T spectra presented in Sec. VI. Some disagreements in the L 90% distribution are found between the LHCf data and the MC simulations. This may be caused by residual errors of the channel-tochannel calibrations of the LHCf detector relative to the LHCf detector simulation.
The corresponding systematic uncertainty of the L 90% distribution is evaluated by comparing the L 90% distri-bution of the LHCf π 0 candidate events of the measured data with the MC simulation. The L 90% distribution for LHCf π 0 events is increased by at most one radiation length compared to the MC simulation. The systematic shifts of p T spectra bin contents are taken from the ratio of p T spectra with artificial shifts of the L 90% distribution to the p T spectra without any L 90% shift. This effect may distort the measured p T spectra by 0-20 % depending on p T .

C. Offset of beam axis
In the geometrical analysis of the data, the projected position of the zero degree collision angle at the LHCf detectors (beam center) can vary from fill to fill owing to slightly different beam transverse position and crossing angles at IP1. The beam center at the LHCf detectors can be determined by two methods; first by using the distribution of particle incident positions measured by the LHCf detectors and second by using the information from the Beam Position Monitors (BPMSW) installed ±21 m from IP1 [41]. Consistent results for the beam center are obtained by the two methods applied to LHC fills 1089-1134 within 1 mm accuracy. The systematic shifts to p T spectra bin contents are evaluated by taking the ratio of spectra with the beam-center displaced by 1 mm to spectra with no displacement as determined by the distribution of particle incident positions measured by the LHCf detectors. Owing to the fluctuations of the beam-center position, the p T spectra are modified by 5-20 % depending on the rapidity range.

D. Single-hit selection
Since energy reconstruction is degraded when more than one photon hits a given calorimeter, only single-hit events are used in the analysis. Owing to selection efficiency greater than 98 % for single-hit events and rejection of contamination by multi-hit events by the invariant mass cut, the systematic shift caused by the uncertainty in single-hit selection to bin contents is 3 %.

E. Position dependent correction
As described in Sec. IV A 2, energy reconstruction of the photons is sensitive to shower leakage effects which are a function of the photon incident position. Systematic uncertainties related to the leakage-out and leakagein effects arise from residual errors of calorimeter response when tuning of the LHCf detector simulation to the calibration data taken at SPS [20] that then lead to a mis-reconstruction of energy. Another source of uncertainties in energy reconstruction is an error in light-yield collection efficiency which is also dependent on the photon incident position. The systematic uncertainty due to position dependent effects is estimated by comparing two distributions of the energy deposited at each incident position bin. The first distribution is taken from the beam tests at SPS and the second distribution is generated by toy MC simulations that assume the upstream geometry of the test beam at SPS. Shifts of reconstructed p T attributed to the residual errors in calorimeter response between these two energy distributions are assigned as the systematic uncertainties. The typical systematic shifts of Arm1 (Arm2) are 5 % (5 %) for low p T and 40 % (30 %) for large p T . Owing to the light guide geometry, the systematic uncertainty of the Arm1 detector is larger than the Arm2 detector.

F. Luminosity
The instantaneous luminosity is derived from the counting rate of the Front Counters (FC). The calibration of the FC counting rates to the instantaneous luminosity was made during the Van der Meer scans on April 26th and May 9th 2010 [23]. The calibration factors obtained from two Van der Meer scans differ by 2.1 %. The estimated luminosities by the two FCs for the May 15th data differ by 2.7 %. Considering the uncertainty of ±5.0 % in the beam intensity measurement during the Van der Meer scans [42], we estimate an uncertainty of ±6.1 % in the luminosity determination.

VI. RESULTS OF ANALYSIS
The p T spectra derived from the independent analyses of the Arm1 and Arm2 detectors are presented in Fig. 6 for six ranges of rapidity y: 8.9 to 9.0, 9.0 to 9.2, 9.2 to 9.4, 9.4 to 9.6, 9.6 to 10.0 and 10.0 to 11.0. The spectra in Fig. 6 are after all corrections discussed in previous sections have been applied. The inclusive production rate of neutral pions is given by the expression σ inel is the inelastic cross section for proton-proton collisions at √ s = 7 TeV. Ed 3 σ/dp 3 is the inclusive cross section of π 0 production. The number of inelastic collisions, N inel , used for normalizing the production rates of Fig. 6 has been calculated from N inel = σ inel Ldt, assuming the inelastic cross section σ inel = 73.6 mb. This value for σ inel has been derived from the best COMPETE fits [38] and the TOTEM result for the elastic scattering cross section [43]. Using the integrated luminosities reported in Sec. III A, N inel is 1.85×10 8 for Arm1 and 1.40×10 8 for Arm2. d 2 N (p T , y) is the number of π 0 s detected in the transverse momentum interval (dp T ) and the rapidity interval (dy) with all corrections applied.
In Fig. 6, the red dots and blue triangles represent the results from Arm1 and Arm2, respectively. The error bars and shaded rectangles indicate the one standard deviation statistical and total systematic uncertainties, respectively. The total systematic uncertainties are given by adding all uncertainty terms except for the luminosity in quadrature. The vertical dashed lines shown in the rapidity range below 9.2 indicate the p T threshold of the Arm2 detector owing to the photon energy threshold and the geometrical acceptance. The p T threshold of the Arm1 detector occurs at a higher value of p T than Arm2 due to its smaller acceptance. A general agreement between the Arm1 and Arm2 p T spectra within statistical and systematic uncertainties is evident in Fig. 6. Fig. 7 presents the combined p T spectra of the Arm1 and Arm2 detectors (black dots). The 68 % confidence intervals incorporating the statistical and systematic uncertainties are indicated by the shaded green rectangles. The combined spectra below the p T threshold of Arm1 are taken from the Arm2 spectra alone. Above the p T threshold of Arm1, experimental p T spectra of the Arm1 and Arm2 detectors have been combined following the "pull method" [44] and the combined spectra have accordingly been obtained by minimizing the value of the chi-square function defined as where the index i represents the p T bin number running from 1 to n (the total number of p T bins), N obs a,i is the number of events and σ a,i is the uncertainty of the Arm-a analysis calculated by quadratically adding the statistical uncertainty and the energy scale uncertainty estimated by test beam data at SPS. The S a,i denotes the systematic correction to the number of events in the i-th bin of Arm-a: The coefficient f j a,i is the systematic shift of i-th bin content due to the j-th systematic uncertainty term. The systematic uncertainty is assumed fully uncorrelated between the Arm1 and Arm2 detectors, and consists of six uncertainties related to energy scale owing to the invariant mass shift, PID, beam center position, single-hit, position dependent correction, and contamination by mulithit π 0 events. Coefficients ε j a , which should follow a Gaussian distribution, can be varied to achieve the minimum χ 2 value in each chi-square test, while they are constrained by the penalty term The π 0 production rates for the combined data of LHCf are summarized in Tables IV-IX. Note that the uncertainty in the luminosity determination ±6.1%, that is not included in Fig. 7, can make a p T independent shift of all spectra.
For comparison, the p T spectra predicted by various hadronic interaction models are also shown in Fig. 7. The hadronic interaction models that have been used in Fig. 7 are dpmjet 3.04 (solid, red), qgsjet II-03 (dashed, blue), sibyll 2.1 (dotted, green), epos 1.99 (dashed dotted, magenta), and pythia 8.145 (default parameter set, dashed double-dotted, brown). In these MC simulations, π 0 s from short lived particles that decay within 1 m from IP1, for example η → 3π 0 , are also counted to be consistent with the treatment of the experimental data. Note that, since the experimental p T spectra have been corrected for the influences of the detector responses, event selection efficiencies and geometrical acceptance efficiencies, the p T spectra of the interaction models may be compared directly to the experimental spectra as presented in Fig. 7. Fig. 8 presents the ratios of p T spectra predicted by the various hadronic interaction models to the combined p T spectra. Error bars have been taken from the statistical and systematic uncertainties. A slight step found around p T = 0.3 GeV in 8.9 < y < 9.0 is due to low p T cutoff of the Arm1 data. The ratios are summarized in Tables X-XV.

A. Transverse momentum spectra
Several points can be made about Fig. 8. First, dpmjet 3.04 and pythia 8.145 show overall agreement with the LHCf data for 9.2 < y < 9.6 and p T < 0.2 GeV, while the expected π 0 production rates by both models exceed the LHCf data as p T becomes large. The latter observation can be explained by the baryon/meson production mechanism that has been employed in both models. More specifically, the "popcorn model" [45,46] is used to produce baryons and mesons through string breaking, and this mechanism tends to lead to hard pion spectra. sibyll 2.1, which is also based on the popcorn model, also predicts harder pion spectra than the experimental data, although the expected π 0 yield is generally small.
On the other hand, qgsjet II-03 predicts π 0 spectra that are softer than the LHCf data and the other models. This might be due to the fact that only one quark exchange is allowed in the qgsjet model. The remnants produced in a proton-proton collision are likewise baryons with relatively small mass, so fewer pions with large energy are produced.
Among hadronic interaction models tested in this anal- ysis, epos 1.99 shows the best overall agreement with the LHCf data. However epos 1.99 behaves softer than the data in the low p T region, p T 0.4 GeV in 9.0 < y < 9.4 and p T 0.3 GeV in 9.4 < y < 9.6, and behaves harder in the large p T region. Specifically a dip found in the ratio of epos 1.99 to the LHCf data for y > 9.0 can be attributed to the transition between two pion production mechanisms: string fragmentation via cut Pomeron process (low energy ∼ low p T for the fixed rapidity) and remnants of projectile/target (high energy ∼ large p T for the fixed rapidity) [47].

B. Average transverse momentum
According to the scaling law proposed by several authors [48][49][50], the average transverse momentum as a function of rapidity should be independent of the center of mass energy in the projectile fragmentation region. Average transverse momentum, p T , can be obtained by fitting an empirical function to the p T spectra in each rapidity range. In this analysis, among several ansatz proposed for fitting the p T spectra, an exponential distribution has been first chosen with the form This distribution is motivated by a thermodynamical model [51]. The parameter A [GeV −2 ] is a normalization factor and T [GeV] is the temperature of π 0 s with a given transverse momentum p T . Using Eq. 10, p T is derived as a function of T : where K α (m π 0 /T ) is the modified Bessel function.
Best-fit results for T and p T are summarized in Table II. The worse fit quality values are found for 9.2 < y < 9.4 (χ 2 /dof = 3.6) and 9.4 < y < 9.6 (χ 2 /dof = 11.1). These are caused by data points near p T = 0.25 GeV which exceed the best-fit exponential distribution and the experimental p T spectra decreasing more rapidly than Eq. (10) for p T > 0.3 GeV. The upper panels in Fig. 9 show the experimental p T spectra (black dots and green shaded rectangles) and the best-fit of Eq. (10) (dashed curve) in the rapidity range 9.2 < y < 9.4 and 9.4 < y < 9.6. The bottom panels in Fig. 9 show the ratio of the best-fit distribution to the experimental data (blue triangles). Shaded rectangles indicate the statistical and systematic uncertainties. Even though the minimum χ 2 /dof values are large, the bestfit T values are consistent with temperatures that are typical of soft QCD processes and the predictions of the thermodynamical model (T 180 MeV) [51] for y > 8.9. Another possibility is that the p T distributions in Fig. 7 can also be described by a Gaussian distribution: The Gaussian width σ Gauss determines the mean square p T of the p T spectra. p T is derived as a function of σ Gauss according to: where f (p T ) is given by Eq. (12). Best-fit results for σ Gauss and p T are summarized in Table II. In this case good fit quality values are found for all rapidity ranges. The best-fit of Eq. (12) (dotted curve) and the ratio of the best-fit Gaussian distribution to the experimental data (red open boxes) are found in Fig. 9. A third approach for estimating p T is simply numerically integrating the p T spectra. With this approach p T is given by where f (p T ) is the measured spectrum given in Fig. 7 for each of the six ranges of rapidity. In this analysis, p T is obtained over the rapidity range 9.2 < y < 11.0 where the p T spectra are available down to 0 GeV. Although the upper limits of numerical integration are actually finite, p upper T ≦ 0.6 GeV, the contribution of the high p T tail to p T is negligible. p upper T and the obtained p T are summarized in Table II.
The values of p T obtained by the three methods discussed above are in general agreement. When a specific values of p T are needed for this paper the values chosen ( p T LHCf ) are defined as follows. For the rapidity range 8.9 < y < 9.2, p T LHCf is taken from the weighted mean of p T obtained by the exponential fit of Eq. (11) and the Gaussian fit of Eq. (13). The systematic uncertainty related to a possible bias of the p T extraction methods is estimated by the difference of p T derived from these two different fitting functions. The estimated systematic uncertainty is ±6 % for both rapidity bins. For the rapidity range 9.2 < y < 11.0, the results obtained by the Gaussian fit and numerical integration are used to calculate the weighted mean of p T LHCf in order to avoid the poor quality of fit of the exponential function in this rapidity range. Systematic uncertainty is estimated to be ±3 % and ±2 % for 9.2 < y < 9.4 and 9.4 < y < 11.0, respectively. The values of p T LHCf obtained by the above calculation are summarized in Table III. Ratios of the best-fit exponential or Gaussian distribution to the experimental data (blue triangles or red open boxes) and the statistical and systematic uncertainties (green shaded areas). For both the upper and bottom panels, the rapidity ranges 9.2 < y < 9.4 and 9.4 < y < 9.6 are shown on the left and right panels, respectively.  The values of p T that have been obtained in this analysis, shown in Table III, are compared in Fig. 10 with the results from UA7 at SppS ( √ s = 630 GeV) [5] and the predictions of several hadronic interaction models. In Fig. 10 p T is presented as a function of rapidity loss ∆y ≡ y beam − y, where beam rapidity y beam is 8.92 for √ s = 7 TeV and 6.50 for √ s = 630 GeV. This shift of rapidity scales the results with beam energy and it allows a direct comparison between LHCf results and past experimental results at different collision energies. The black dots and the red diamonds indicate the LHCf data and the UA7 results, respectively. Although the LHCf and UA7 data in Fig. 10 have limited overlap and the systematic errors of the UA7 data are relatively large, the p T spectra for LHCf and UA7 in Fig. 10   II-03 gives softer π 0 spectra (smaller p T ) than the experimental data. For each prediction, solid and dashed lines indicate p T at the center of mass energy at SppS and the LHC, respectively. Of the three models the predictions by epos 1.99 show the smallest dependence of p T on the two center of mass energies, and this tendency is consistent with the LHCf and UA7 results except for the UA7 data at ∆y = −0.15 and 0.25. It is also evident in Fig. 10 that amongst the three models the best agreement with the LHCf data is obtained by epos 1.99.

VIII. CONCLUSIONS
The inclusive production of neutral pions in the rapidity range larger than y = 8.9 has been measured by the LHCf experiment in proton-proton collisions at the LHC in early 2010. Transverse momentum spectra of neutral pions have been measured by two independent LHCf detectors, Arm1 and Arm2, and give consistent results. The combined Arm1 and Arm2 spectra have been compared with the predictions of several hadronic interaction models. dpmjet 3.04, epos 1.99 and pythia 8.145 agree with the LHCf combined results in general for the rapidity range 9.0 < y < 9.6 and p T < 0.2 GeV. qgsjet II-03 has poor agreement with LHCf data for 8.9 < y < 9.4, while it agrees with LHCf data for y > 9.4. Among the hadronic interaction models tested in this paper, epos 1.99 shows the best overall agreement with the LHCf data even for y > 9.6.
The average transverse momentum, p T , of the combined p T spectra is consistent with typical values for soft QCD processes. The p T spectra for LHCf and UA7 in Fig. 10 mostly appear to lie along a common curve. The p T spectra derived by LHCf agrees with the expectation of epos 1.99. Additional experimental data are needed to establish the dependence, or independence, of p T on the center of mass collision energy.      XV: Ratio of π 0 production rate of MC simulation to data in the rapidity range 10.0<y<11.0.