Measurement of negatively charged pion spectra in inelastic p+p interactions at $p_{lab}$ = 20, 31, 40, 80 and 158 GeV/c

We present experimental results on inclusive spectra and mean multiplicities of negatively charged pions produced in inelastic p+p interactions at incident projectile momenta of 20, 31, 40, 80 and 158 GeV/c ($\sqrt{s} = $ 6.3, 7.7, 8.8, 12.3 and 17.3 GeV, respectively). The measurements were performed using the large acceptance NA61/SHINE hadron spectrometer at the CERN Super Proton Synchrotron. Two-dimensional spectra are determined in terms of rapidity and transverse momentum. Their properties such as the width of rapidity distributions and the inverse slope parameter of transverse mass spectra are extracted and their collision energy dependences are presented. The results on inelastic p+p interactions are compared with the corresponding data on central Pb+Pb collisions measured by the NA49 experiment at the CERN SPS. The results presented in this paper are part of the NA61/SHINE ion program devoted to the study of the properties of the onset of deconfinement and search for the critical point of strongly interacting matter. They are required for interpretation of results on nucleus-nucleus and proton-nucleus collisions.


1
Abstract We present experimental results on inclusive spectra and mean multiplicities of negatively charged pions produced in inelastic p+p interactions at incident projectile momenta of 20, 31, 40, 80 and 158 GeV/c ( √ s = 6.3, 7.7, 8.8, 12.3 and 17.3 GeV, respectively). The measurements were performed using the large acceptance NA61/SHINE hadron spectrometer at the CERN Super Proton Synchrotron.
Two-dimensional spectra are determined in terms of rapidity and transverse momentum. Their properties such as the width of rapidity distributions and the inverse slope parameter of transverse mass spectra are extracted and their collision energy dependences are presented. The results on inelastic p+p interactions are compared with the corresponding data on central Pb+Pb collisions measured by the NA49 experiment at the CERN SPS.
The results presented in this paper are part of the NA61/ SHINE ion program devoted to the study of the properties of the onset of deconfinement and search for the critical point of strongly interacting matter. They are required for interpretation of results on nucleus-nucleus and proton-nucleus collisions.

Introduction
This paper presents experimental results on inclusive spectra and mean multiplicities of negatively charged pions produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c. The measurements were performed by the multi-purpose NA61/SHINE -the SPS Heavy Ion and Neutrino Experiment [1] -at the CERN Super Proton Synchrotron (SPS). They are part of the NA61/SHINE ion program devoted to the study of the properties of the onset of deconfinement and search for the critical point of strongly interacting matter. Within this program data on p+p, Be+Be and p+Pb collisions were recorded and data on Ar+Ca and Xe+La collisions will be registered within the coming years. The started two dimensional scan in collision energy and size of colliding nuclei is mainly motivated by the observation of the onset of deconfinement [2,3] in central Pb+Pb collisions at about 30A GeV/c by the NA49 experiment at the CERN SPS. Recently the NA49 results were confirmed by the RHIC beam energy scan programme and their interpretation by the onset of deconfinement is supported by the LHC results (see Ref. [4] and references therein).
In addition to the ion programme, NA61/SHINE is conducting precise hadron production measurements for improving calculations of the initial neutrino flux in longbaseline neutrino oscillation experiments [5][6][7][8] as well as for more reliable simulations of cosmic-ray air showers [9,10].
An interpretation of the rich experimental results on nucleus-nucleus collisions relies to a large extent on a com-parison to the corresponding data on p+p and p+A interactions. However, the available data concern mainly basic features of unidentified charged hadron production and they are sparse. Many needed results on hadron spectra, fluctuations and correlations are missing. Detailed measurements of hadron spectra in a large acceptance in the beam momentum range covered by the data presented in this paper exist only for inelastic p+p interactions at 158 GeV/c [11][12][13]. Thus the new high precision measurements of hadron production properties in p+p and p+A interactions are necessary and they are performed in parallel with the corresponding measurements in nucleus-nucleus collisions. Among the many different hadrons produced in high energy collisions pions are the lightest and by far the most abundant ones. Thus, data on pion production properties are crucial for constraining basic properties of models of strong interactions. In particular, the most significant signals of the onset of deconfinement (the "kink" and "horn") [14] require precise measurements of the mean pion multiplicity at the same beam momenta per nucleon as the corresponding A+A data. Moreover, the NA61/SHINE data are taken with the same detector and the same acceptance.
In the CERN SPS beam momentum range of 10-450 GeV/c the mean multiplicity of negatively charged pions in inelastic p+p interactions increases from about 0.7 at 10 GeV/c to about 3.5 at 450 GeV/c [15]. Among three charged states of pions the most reliable measurements in the largest phase-space are usually possible for π − mesons. Neutral pions decay electromagnetically into two photons and thus measurements of their production properties require measurements of both photons and then extraction of the π 0 signal from the two-photon mass spectra. Charged pions are easy to detect by ionisation detectors as they decay weakly with a long lifetime (cτ = 7.8 m). A significant fraction of positively charged hadrons are protons (25%) and kaons (5%) [11][12][13]. Therefore measurements of π + mesons require their identification by measurements of the energy loss and/or time-of-flight. This identification is not as important for π − mesons because the contamination of negatively charged particles by K − mesons and anti-protons is below 10% [11][12][13] and can be subtracted reliably. The latter method is used in this paper and it allows to derive π − spectra in the broadest phase-space region in a uniform way. Results obtained using explicit pion identification are planned in future NA61/SHINE publications.
The paper is organised as follows. In Sec. 2 the NA61/ SHINE experimental set-up is described. Details on the beam, trigger and event selection are given in Sec. 3. Data reconstruction, simulation and detector performance are described in Sec. 4. Analysis techniques and final results are presented in Secs. 5 and 6, respectively. These results are compared with the corresponding data on central Pb+Pb collisions in Sec. 7. A summary in Sec. 8 closes the paper. The beam and trigger detector configuration used for data taking on p+p interactions in 2009 is shown. Alignment of the chosen coordinate system is shown on the plot; its origin lies in the middle of VTPC-2, on the beam axis. The nominal beam direction is along the z axis. The magnetic field bends charged particle trajectories in the x-z (horizontal) plane. The drift direction in the TPCs is along the y (vertical) axis.
The pion rapidity is calculated in the collision centre of mass system: y = atanh(β L ), where β L = p L /E is the longitudinal component of the velocity, p L and E are pion longitudinal momentum and energy given in the collision centre of mass system. The transverse component of the momentum is denoted as p T and the transverse mass m T is defined as m T = m 2 π + p 2 T , where m π is the charged pion mass. The nucleon mass and collision energy per nucleon pair in the centre of mass system are denoted as m N and √ s NN , respectively.

The NA61/SHINE facility
The NA61/SHINE experimental facility consists of a large acceptance hadron spectrometer located in the CERN North Area Hall 887 (EHN1) and the H2 beam-line to which beams accelerated in the CERN accelerator complex are delivered from the Super Proton Synchrotron. NA61/SHINE profits from the long development of the CERN proton and ion sources and the accelerator chain as well as the H2 beam line of the CERN North Area. The latter has recently been modified to also serve as a fragment separator as needed to produce the Be beams for NA61/SHINE. Numerous components of the NA61/SHINE set-up were inherited from its predecessor, the NA49 experiment [16].
The schematic layout of the NA61/SHINE detector is shown in Fig. 1.
A set of scintillation and Cherenkov counters as well as beam position detectors (BPDs) upstream of the spectrometer provide timing reference, identification and position measurements of incoming beam particles. The trigger scintillator counter S4 placed downstream of the target is used to select events with collisions in the target area. Details on this system are presented in Sec. 3.
The main tracking devices of the spectrometer are large volume Time Projection Chambers (TPCs). Two of them, the vertex TPCs (VTPC-1 and VTPC-2 in Fig. 1), are located in the magnetic fields of two super-conducting dipole magnets with a maximum combined bending power of 9 Tm which corresponds to about 1.5 T and 1.1 T fields in the upstream and downstream magnets, respectively. This field configuration was used for data taking on p+p interactions at 158 GeV/c. In order to optimise the acceptance of the detector at lower collision momenta, the field in both magnets was lowered in proportion to the beam momentum.
Two large TPCs (MTPC-L and MTPC-R) are positioned downstream of the magnets symmetrically to the beam line. The fifth small TPC (GAP-TPC) is placed between VTPC-1 and VTPC-2 directly on the beam line. It closes the gap between the beam axis and the sensitive volumes of the other TPCs.
The TPCs are filled with Ar:CO 2 gas mixtures in proportions 90:10 for the VTPCs and the GAP-TPC, and 95:5 for the MTPCs.
The particle identification capability of the TPCs based on measurements of the specific energy loss, dE/dx, is augmented by time-of-flight (tof ) measurements using Timeof-Flight (ToF) detectors. The high resolution forward calorimeter, the Projectile Spectator Detector (PSD), measures energy flow around the beam direction, which in nucleus-nucleus collisions is primarily given by the projectile spectators.
NA61/SHINE uses various solid nuclear targets and a liquid hydrogen target (see Sec. 3 for details). The targets are positioned about 80 cm upstream of the sensitive volume of VTPC-1.
The results presented in this paper were obtained using information from the Time Projection Chambers, the Beam Position Detectors as well as from the beam and trigger counters.

Beams, target, triggers and data samples
This section summarises basic information on the beams, target, triggers and recorded data samples which is relevant for the results presented in this paper.
Secondary beams of positively charged hadrons at 20, 31, 40, 80 and 158 GeV/c are produced from 400 GeV protons extracted from the SPS in a slow extraction mode with a flat-top of 10 seconds. The secondary beam momentum and intensity is adjusted by proper setting of the H2 beam-line magnet currents and collimators. The beam is transported along the H2 beam-line towards the experiment. The precision of the setting of the beam magnet currents is approximately 0.5%. This was verified by a direct measurement of the beam momentum at 31 GeV/c by bending the incoming beam particles into the TPCs with the maximum magnetic field [17]. The selected beam properties are given in Table 1.
The set-up of beam detectors is illustrated in the inset on Fig. 1. Protons from the secondary hadron beam are identified by two Cherenkov counters, a CEDAR [18] (either CEDAR-W or CEDAR-N) and a threshold counter (THC). The CEDAR counter, using a coincidence of six out of the eight photo-multipliers placed radially along the Cherenkov ring, provides positive identification of protons, while the THC, operated at pressure lower than the proton threshold, is used in anti-coincidence in the trigger logic. Due to their limited range of operation two different CEDAR counters were used, namely for beams at 20, 31, and 40 GeV/c the CEDAR-W counter and for beams at 80 and 158 GeV/c the CEDAR-N counter. The threshold counter was used for all beam energies. A selection based on signals from the Cherenkov counters allowed to identify beam protons with a purity of about 99%. A consistent value for the purity was found by bending the beam into the TPCs with the full magnetic field and using the dE/dx identification method [19]. The fraction of protons in the beams is given in Table 1. Two scintillation counters, S1 and S2, provide beam definition, together with the three veto counters V0, V1 and V1 p with a 1 cm diameter hole, which are defining the beam before the target. The S1 counter provides also the timing (start time for all counters). Beam protons are then selected by the coincidence S1 ∧ S2 ∧ V0 ∧ V1 ∧ V1 p ∧ CEDAR ∧ THC. Trajectories of individual beam particles were measured in a telescope of beam position detectors placed along the beam line (BPD-1/2/3 in Fig. 1). These counters are small (4.8 × 4.8 cm 2 ) proportional chambers with cathode strip readout, providing a resolution of about 100 µm in two orthogonal directions, see Ref. [20] for more details. The beam profile and divergence obtained from the BPD measurements are presented in Fig. 2. Due to properties of the H2 beam line both the beam width and divergence at the NA61/SHINE target increase with decreasing beam momentum.
For data taking on p+p interactions a liquid hydrogen target of 20.29 cm length (2.8% interaction length) and 3 cm diameter placed 88.4 cm upstream of VTPC-1 was used. The Liquid Hydrogen Target facility (LHT) filled the target cell with para-hydrogen obtained in a closed-loop liquefaction system which was operated at 75 mbar overpressure with respect to the atmosphere. At the atmospheric pressure of 965 mbar the liquid hydrogen density is ρ LH = 0.07 g/cm 3 . The boiling rate in the liquid hydrogen was not monitored during the data taking and thus the liquid hydrogen density is known only approximately. It has however no impact on the results presented in this paper as they are determined from particle yields per selected event and thus they are independent of the target density. Data taking with inserted and removed liquid hydrogen in the LHT was alternated in order to calculate a data-based correction for interactions with the material surrounding the liquid hydrogen.
Interactions in the target are selected by the trigger system by an anti-coincidence of the incoming beam protons with a small, 2 cm diameter, scintillation counter (S4) placed on the beam trajectory between the two vertex magnets (see    1). This minimum bias trigger is based on the disappearance of the incident proton. In addition, unbiased proton beam events were recorded with a frequency typically by a factor of 10 lower than the frequency of interaction events.

Data processing, simulation and detector performance
Detector parameters were optimised by a data-based calibration procedure which also took into account their time dependences. Small adjustments were determined in consecutive steps for: (i) detector geometry, TPC drift velocities and distortions due to the magnetic field inhomogeneities in the corners of VTPCs, (ii) magnetic field setting, (iii) specific energy loss measurements, (iv) time-of-flight measurements.
Each step involved reconstruction of the data required to optimise a given set of calibration constants and time dependent corrections followed by verification procedures. Details of the procedure and quality assessment are presented in Ref. [21]. The resulting performance in the measurements of quantities relevant for this paper is discussed below.
The main steps of the data reconstruction procedure are: (i) cluster finding in the TPC raw data, calculation of the cluster centre-of-gravity and total charge, (ii) reconstruction of local track segments in each TPC separately, (iii) matching of track segments into global tracks, (iv) track fitting through the magnetic field and determination of track parameters at the first measured TPC cluster, (v) determination of the interaction vertex using the beam trajectory (x and y coordinates) fitted in the BPDs and the trajectories of tracks reconstructed in the TPCs (z coordinate), (vi) refitting the particle trajectory using the interaction vertex as an additional point and determining the particle momentum at the interaction vertex, (vii) matching of ToF hits with the TPC tracks.
An example of a reconstructed p+p interaction at 40 GeV/c is shown in Fig. 3. Long tracks of one negatively charged and two positively charged particles are seen. All particles leave signals in the ToF detectors.
A simulation of the NA61/SHINE detector response is used to correct the reconstructed data. Several MC models were compared with the NA61/SHINE results on p+p, p+C and π+C interactions: FLUKA2008, URQMD1.3.1, VENUS4.12, EPOS1.99, GHEISHA2002, QGSJetII-3 and Sibyll2.1 [17,[22][23][24]. Based on these comparisons and taking into account continuous support and documentation from the developers the EPOS model [25] was selected for the MC simulation. The simulation consists of the following steps (see Ref. [26] for more details): (i) generation of inelastic p+p interactions using the EPOS model [25], (ii) propagation of outgoing particles through the detector material using the GEANT 3.21 package [27] which takes into account the magnetic field as well as relevant physics processes, such as particle interactions and decays, (iii) simulation of the detector response using dedicated NA61/SHINE packages which introduce distortions corresponding to all corrections applied to the real data, (iv) simulation of the interaction trigger selection by checking whether a charged particle hits the S4 counter, see Sec. 3, (v) storage of the simulated events in a file which has the same format as the raw data, (vi) reconstruction of the simulated events with the same reconstruction chain as used for the real data and (vii) matching of the reconstructed tracks to the simulated ones based on the cluster positions.
It should be underlined that only inelastic p+p interactions in the hydrogen in the target cell were simulated and reconstructed. Thus the Monte Carlo based corrections (see Sec. 5) can be applied only for inelastic events. The contribution of elastic events is removed by the event selection cuts (see Sec. 5.1), whereas the contribution of off-target interactions is subtracted based on the data (see Sec. 5.4).
Spectra of π − mesons presented in this paper were derived from spectra of all negatively charged hadrons corrected for a small contamination of mostly K − mesons and anti-protons. The typical acceptance in rapidity and transverse momentum is presented in Fig. 4 for p+p interactions at 20 and 158 GeV/c. This figure also shows acceptance regions for methods based on explicit pion identification using dE/dx and tof measurements. They are limited due to the geometrical acceptance of the ToF detectors, the finite resolution of the dE/dx measurements and limited data statistics.
The quality of measurements was studied by reconstructing masses of K 0 S particles from their V 0 decay topology. As an example the invariant mass distributions of K 0 S candidates found in p+p interactions at 20 and 158 GeV/c are plotted in Fig. 5. The differences between the measured peak positions and the literature value of the K 0 S mass [28] are smaller than 1 MeV/c 2 . The width of the distributions, related to the detector resolution, is about 25% smaller for the Monte Carlo than for the data. This implies that statistical and/or systematic uncertainties of track parameters reconstructed from the data are somewhat underestimated in the simulation. Systematic bias due to this imperfectness was estimated by varying the selection cuts and was found to be below 2% (see Sec. 5.7.2).
The track reconstruction efficiency and the resolution of kinematic quantities were calculated by matching reconstructed tracks to their generated partners. In only 0.1-0.2% of cases a single generated track is matched to more than one reconstructed partner, typically due to failure of matching reconstructed track segments. This effect is taken into account in the correction described in Sec. 5.6. As examples, the reconstruction efficiency as a function of rapidity and transverse momentum for negatively charged pions produced in p+p interactions at 20 and 158 GeV/c is shown in Fig. 6. The resolution of rapidity and transverse momentum measurements is illustrated in Fig. 7. The resolution was calculated as the FWHM of the distribution of the difference between the generated and reconstructed y and p T . These results were obtained for negatively charged pions passing the track selection criteria described in Sec. 5.2. Resolution of the transverse momentum is worse at low beam momenta. This is caused by the lower magnetic field and by the fact that the same rapidity region in the centre of mass frame corresponds to lower momenta in the laboratory frame. Figures 8 and 9 show further examples of the comparison between data and simulation. Distributions of the z coordinate of the fitted vertex are presented in Fig. 8. Distributions of the distance between the track trajectory extrapolated to the z coordinate of the vertex and the vertex in the x-y plane (b x and b y impact parameters) are given in Fig. 9. Differences visible in the tails of distributions are partially due to imperfect simulation of the detector response and, in case of the impact parameter, partially due to the contribution of background tracks from off-time beam particles which are not included in the simulation. The difference is smaller for events selected using more restrictive cuts on the off-time beam particles. A possible small bias due to these effects was estimated by varying the impact parameter cuts and was found to be below 1%.

Analysis technique
This section presents the procedures used for data analysis consisting of the following steps: The MC plot was normalised to the peak height of the data. The K 0 S candidates were selected within 0 < y < −1 and 0 < p T < 0.5 GeV/c for 20 GeV/c and −1 < y < 0 and 0 < p T < 0.5 GeV/c for 158 GeV/c. The distribution was fitted with the sum of a Lorentzian function (signal) and a second order polynomial (background).
(ii) determination of spectra of negatively charged hadrons using the selected events and tracks, (iii) evaluation of corrections to the spectra based on experimental data and simulations, (iv) calculation of the corrected spectra.
Corrections for the following biases were evaluated and applied: (i) geometrical acceptance, (ii) contribution of off-target interactions, (iii) contribution of particles other than negatively charged pions produced in inelastic p+p interactions, (iv) losses of inelastic p+p interactions as well as of negatively charged pions produced in accepted interactions due to the trigger and the event and track selection criteria employed in the analysis.
These steps are described in the successive subsections. The final results refer to π − mesons produced in inelastic p+p interactions by strong interaction processes and in electromagnetic decays of produced hadrons. Such pions are referred to as primary π − . The term primary will be used in the above meaning also for other particles.
The analysis was performed independently in (y, p T ) and (y, m T ) bins. The bin sizes were selected taking into account the statistical uncertainties as well as the resolution of the momentum reconstruction. Corrections as well as statistical and systematic uncertainties were calculated for each bin.

Event selection criteria
This section presents the event selection criteria. The number of events selected by the trigger (see Sec. 3) and used in the analysis is shown in Table 2. The fraction of events selected for the analysis increases with the interaction energy, mostly due to lower beam intensity (see Table 1) and resulting smaller off-time particle contamination, and smaller fraction of the low multiplicity events for which no tracks are found within the acceptance.
The following event selection criteria were applied to the events recorded with the interaction trigger: (i) no off-time beam particle is detected within ±2 µs around the trigger particle, (ii) the beam particle trajectory is measured in at least one of BPD-1 or BPD-2 and in the BPD-3 detector positioned just in front of the LHT, (iii) there is at least one track reconstructed in the TPCs and fitted to the interaction vertex, (iv) the vertex z position (fitted using the beam and TPC tracks) is not farther away than 40 cm from the centre of the LHT, (v) events with a single, well measured positively charged track with absolute momentum close to the beam momentum are rejected. The momentum thresholds are listed in Table 3.  The off-line (listed above) and on-line (the interaction trigger condition, see Sec. 3) event cuts select a large fraction of well measured (cuts (i) and (ii)) inelastic (cut (iii)) p+p interactions. The cut (iii) removes part of elastic interactions. However in some elastic events at beam momenta up to 80 GeV/c the beam particle is deflected enough to be measured in the detector. This is demonstrated in the momentum distributions shown in Fig. 10. Such events are removed by cut (v).
Moreover cut (iv) significantly suppresses interactions outside the hydrogen in the target cell. The corrections for the contribution of interactions outside the hydrogen in the target cell and the loss of inelastic events are presented in Secs. 5.4 and 5.6.

Track selection criteria
In order to select well-measured tracks of primary negatively charged hadrons as well as to reduce the contamination of tracks from secondary interactions, weak decays and off-time interactions the following track selection criteria were applied: acteristic for electrons are rejected. The electron contribution to particles with momenta above 20 GeV/c is corrected using the simulation. The electron selection procedure is visualised in Fig. 11.
The spectra of negatively charged particles after track and event selections were obtained in 2-dimensional bins of (y, p T ) and (y, m T ). The spectra were evaluated in the centreof-mass frame after rotation of the z axis into the proton beam direction measured event-by-event by the BPDs.

Correction for detector acceptance
The detection and reconstruction inefficiencies are corrected using the simulation described in Sec. 4. However, in order to limit the impact of possible inaccuracies of this simulation, only regions are accepted where the reconstruction efficiency (defined as the ratio of the number of reconstructed and matched MC π − tracks passing the track selection criteria to the number of generated π − ) equals at least . The acceptance is calculated as fraction of (y, p T , φ ) bins accepted for given y and p T for tracks with p x < 0 selected for this analysis.
90%. These regions were identified using a separate, statistically independent simulation in three-dimensional bins of y, p T or m T and the azimuthal angle φ (5 • bin width). The resulting acceptance maps are shown in Fig. 12. The acceptance calculated in the y < 0 region, not used for the final results, is shown also for comparison. We chose an upper limit of 1.5 GeV/c for the transverse momentum spectra, because beyond the admixture of background tracks reaches a level which cannot be handled by the correction procedures used in this paper. Future publications will be devoted to the high p T region.
Since neither target nor beam are polarized, we can assume a uniform distribution of particles in φ . The data falling into the accepted bins is summed over φ bins and the (y, p T /m T ) bin content is multiplied by a correction factor to compensate for the rejected φ ranges. The acceptance correction also compensates for the p x < 0 selection (see Sec. 5.2, point (iii)).
Even a small deviation of the beam direction from the nominal axis (z) results in a non-negligible bias in the reconstructed transverse momentum. The beam direction is measured in the BPDs, and the particle momenta are recalculated to the frame connected with the beam direction. However, the detector acceptance depends on the momentum in the detector frame. Therefore the acceptance selection is done in the detector frame, and the acceptance correction is applied as a weight to each track. The weights are used to obtain particle spectra corrected for the detector acceptance in the frame connected with the beam direction.

Correction for off-target interactions
The spectra were derived for events with liquid hydrogen in ( I ) and removed ( R ) from the LHT. The latter data set represents interactions outside the liquid hydrogen (interactions with materials downstream and upstream of the liquid hydrogen). The differential inclusive yield of negatively charged particles per event in interactions of beam protons with the liquid hydrogen inside the LHT (n T [h − ]) is calculated as: where: (i) n I [h − ] and n R [h − ] is the number of tracks in a given bin per event selected for the analysis (see Sec. 5.2) for the data with the liquid hydrogen inserted and removed, respectively, (ii) ε is the ratio of the interaction probabilities for the removed and inserted target operation.
ε was derived based on the distribution of the fitted z coordinate of the interaction vertex. All vertices far away from the target originate from interactions with the beamline and detector materials. Neglecting the beam attenuation in the target one gets: where N[z > −450 cm] is the number of events with fitted vertex z > −450 cm. Examples of distributions of z of the fitted vertex for events recorded with the liquid hydrogen inserted and removed are shown in Fig. 13. Values of ε are listed in Table 4.
The correction for the off-target interactions changes the yields obtained from the target inserted data by less than ±5%, except in the regions where the statistical uncertainty is high.

The correction for contamination of primary π − mesons
More than 90% of primary negatively charged particles produced in p+p interactions in the SPS energy range are π − mesons [11][12][13]. Thus π − meson spectra can be obtained by subtracting the estimated non-pion contribution from the spectra of negatively charged particles and additional particle identification is not required.
The simulation described in Sec. 4 was used to calculate corrections for the admixture of particles other than primary π − mesons to the reconstructed negatively charged particles. The dominating contributions are primary K − andp, and secondary π − from weak decays of Λ and K 0 S (feed-down) and from secondary interactions, incorrectly fitted to the primary vertex.
The EPOS spectra were adjusted based on the existing data [24,29]. Preliminary NA61/SHINE results were used to scale double differential spectra of K − , andp [30]. EPOS spectra of π − were replaced by the preliminary NA61/ SHINE results [31] normalised to the multiplicity from the world data compilation [32]. Spectra of Λ and K 0 S were scaled by a constant factor derived at each energy using the world data compilation [33] of total multiplicities. The impact of the adjustments on the final spectra ranges from −2% to +5% in most regions, except of the low p T region at the low beam momenta, where it reaches +20%.
As it was found in [24] the yields of K − andp are strongly correlated with the π − yield. Thus the correction for the contribution of primary hadrons is performed via the multiplicative factor c K . On the contrary the contribution due to weak decays and secondary interactions is mostly located in the low p T region, and it is weakly correlated with the primary pion yield in this region. Thus this feed-down contribution is corrected for using the additive correction c V . The total correction is applied in as: where The spectrum of a particle x is denoted as n[x] whereas n[other] stands for all primary and secondary particles other than K − ,p, π − and feed-down from Λ and K 0 S . The spectrum n[other] of all other particles originates mostly from secondary interactions with >90% occurring in the hydrogen target. This contribution was taken from the simulations without an additional adjustment. The superscript MC marks adjusted EPOS spectra. The subscript sel indicates that the event and track selection criteria were applied and then the correction for the detector acceptance was performed; the reconstructed tracks were identified by matching.

Correction for event as well as track losses and migration
The multiplicative correction c loss for losses of inelastic events as well as losses and bin-to-bin migration of primary π − mesons emitted within the acceptance is calculated using the Monte Carlo simulation as: where the subscript gen indicates the generated spectrum of primary π − mesons binned according to the generated momentum vector. Then the final, corrected π − meson spectrum in inelastic p+p interactions is calculated as The dominating effects contributing to the c loss correction are losses of inelastic events due to the trigger and off-line event selection, the pion migration between analysis bins, the pion reconstruction inefficiency.

Statistical uncertainties
Statistical errors receive contributions from the finite statistics of both the data as well as the simulated events used to obtain the correction factors. The dominating contribution is the uncertainty of the data which is calculated assuming a Poisson probability distribution for the number of entries in a bin. The Monte Carlo statistics was higher than the data statistics. Also the uncertainties of the Monte Carlo corrections are significantly smaller than the uncertainties of the number of entries in bins.

Systematic uncertainties
Systematic errors presented in this paper were calculated taking into account contributions from the following effects.  cut was changed by ±0.01 dE/dx units (where 1 corresponds to a minimum ionising particle, and 0.04 is a typical width of a single particle distribution), and the off-time interactions cut was varied from a ±1 µs to a ±3 µs time window. The assigned systematic uncertainty was calculated as the maximum of the absolute differences between the results obtained for both cut values. This contribution is drawn with a long-dashed red line (--) in Fig. 14. (ii) Uncertainty of the correction for contamination of the primary π − mesons. The systematic uncertainty σ ii of this correction was assumed as 20% (for 40, 80 and 158 GeV/c) and 40% (for 20 and 31 GeV/c) of the absolute value of the correction. At the low beam momenta there was less data available to adjust the simulated spectra, which was the reason to increase the uncertainty. This contribution is drawn with a dasheddotted blue line (-···-) in Fig. 14. The absolute correction is small thus the related systematic uncertainty is small also. (iii) Uncertainty of the correction for the event losses. The uncertainty was estimated using 20% of the correction magnitude and a comparison with the correction calculated using the VENUS [34] model: This contribution is drawn with a short-dashed green line (---) in Fig. 14. (iv) Uncertainty related to the track selection method. It was estimated by varying the track selection cuts: removing the impact parameter cut and decreasing the minimum number of required points to 25 (total) and 10 (in VTPCs) and by checking symmetries with respect to y = 0 and p T = 0. The potential bias is below 2% and the corresponding contribution was neglected.
The total systematic uncertainty was calculated by adding in quadrature the contributions σ sys = σ 2 i + σ 2 ii + σ 2 iii . This uncertainty is listed in the tables including numerical values and it is visualised by a shaded band around the data points in plots presenting the results. Statistical and systematic uncertainties in selected example regions are shown in Fig. 14. Systematic biases in different bins are correlated, whereas statistical fluctuations are almost independent. Figure 15 presents a comparison of the rapidity spectrum of π − mesons produced in inelastic p+p interactions at 158 GeV/c (for details see Sec. 6) from the present analysis with the corresponding spectrum measured by NA49 [11]. Statistical and systematic uncertainties of the NA49 spectrum are not explicitly given but the published information implies that the systematic uncertainty dominates and amounts to several %. The results agree within the systematic uncertainties of the NA61/SHINE spectra. The analysis method of p+p interactions at 158 GeV/c performed by NA49 [11] differed from the one used in this paper. In particular, pions were identified by dE/dx measurement and the NA49 event selection criteria did not include the selection according to the fitted z coordinate of the interaction vertex and the rejection of elastic interactions. Namely, all events passing the trigger selection and off-line quality cuts were used for the analysis. For comparison, this event selection procedure was applied to the NA61/SHINE data. As a result 20% more events were accepted. Approximately half of them were unwanted elastic and off-target interactions and half were wanted inelastic interactions. Then the corrections corresponding to the changed selection criteria were applied (the contribution of elastic events was subtracted using the estimate from Ref. [11]). The fully corrected rapidity spectrum obtained using this alternative analysis is also shown in Fig. 15. The differences between the results for the standard and alternative methods are below 0.5% at y < 2 and below 2% at higher y. Figure 16 shows a comparison of the rapidity distribution at 31 GeV/c with the MIRABELLE results at 32 GeV/c [35]. A parametrisation of the distribution and the total π − multiplicity are provided. As the parametrisation appears to be incorrectly normalised, we normalised it to the total multiplicity. The results agree within the NA61/SHINE systematic uncertainties. The spectra measured in p+p interactions should obey reflection symmetry with respect to mid-rapidity. As the NA61/SHINE acceptance extends somewhat below mid-rapidity a check of the reflection symmetry can be performed and used to validate the measurements. It was verified that the yields measured for y < 0 agree with those measured for y > 0 in the reflected acceptance within 1.5%. A similar agreement was also found at lower beam momenta. The measurements above mid-rapidity are taken as the final results. Nevertheless, for comparison the points at y < 0 were added in Figs. 15 and 20 in the regions where the p T acceptance extends to zero.

Results
This section presents results on inclusive π − meson spectra in inelastic p+p interactions at beam momenta of 20, 31, 40, 80 and 158 GeV/c. The spectra refer to pions produced by strong interaction processes and in electromagnetic decays of produced hadrons.
Numerical results corresponding to the plotted spectra as well as their statistical and systematic uncertainties are given in Ref. [36].
6.1 Double differential spectra The double differential inclusive spectra of π − mesons in rapidity and transverse momentum produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c are shown in Fig. 17. The transverse momentum distributions at 20 and 158 GeV/c are plotted in Fig. 18. Here d 2 n dy dp T or d 2 n dy dm T were calculated by dividing the fully corrected bin contents n[π − ] (see Sec. 5) by the bin size. The spectra in (y, m T ) are not shown here but they are given in the compilation of the numerical values [36].

Transverse mass spectra
The transverse mass spectra at mid-rapidity (0 < y < 0.2) are shown in Fig. 19 (left). A function was fitted in the range 0.2 < m T − m π < 0.7 GeV/c 2 and is indicated by lines in Fig. 19 (left). The fitted parameters were the normalisation A and the inverse slope T . They minimise the χ 2 function which was calculated using statistical errors only. In the χ 2 calculation a measured bin content (dn/dm T ) was compared with the integral of the fitted function in a bin divided by the bin width. Similar fits were performed to spectra in other rapidity bins containing data in the fit range. The rapidity dependence of the fitted inverse slope parameter T is presented in Fig. 19 (right). The T parameter decreases significantly when going from mid-rapidity to the projectile rapidity (y beam = 1.877, 2.094, 2.223, 2.569 and 2.909 at 20, 31, 40, 80 and 158 GeV/c, respectively).

Rapidity spectra
The rapidity spectra are shown in Fig. 20 (left). They were obtained by summing the measured m T spectra and using the exponential function Eq. (9). The function was fitted in the range ending at the maximum measured m T , and starting 0.9 GeV/c 2 below (note this is a different fit from the one shown in Fig. 19). The correction is typically below 0.2% and becomes significant (several %) only at y > 2.4. Half of the correction is added in quadrature to the systematic uncertainty in order to take into account a potential imperfectness of the exponential extrapolation. The pion yield increases with increasing collision energy at all measured rapidities.
The rapidity spectra are parametrised by the sum of two Gaussian functions symmetrically displaced with respect to  158 GeV/c (right) in various rapidity ranges. The legend provides the centres of the rapidity bins, y bin and the scaling factor c used to separate the spectra visually. mid-rapidity: where y 0 and σ 0 are fit parameters, and the total multiplicity π − (y 0 , σ 0 ) is calculated from the requirement that the integral over the measured spectrum equals the integral of the fitted function Eq. (10) in the range covered by the measurements. The χ 2 function was minimised in a similar way as in case of the m T spectra, namely using the integral of  Table 5.

Mean multiplicities
Mean multiplicities of π − mesons, π − , produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c were calculated as the integral of the fitted function Eq. (10). The extrapolation into the unmeasured region at large y contributes about 1%. Half of it is added in quadrature to the systematic uncertainty.
The dependence of the produced average π − multiplicity per inelastic p+p collision on the Fermi's energy measure [37], is plotted in Fig. 20 (right). The results of NA61/SHINE are in agreement with a compilation of the world data [11,32].

Comparison with central Pb+Pb collisions
In this section the NA61/SHINE results on inelastic p+p interactions are compared with the corresponding data on central Pb+Pb collisions published previously by NA49 [2,3]. Pion production properties which are different and similar in p+p interactions and central Pb+Pb collisions are identified. For completeness selected plots include the compilation of the world data on inelastic p+p interactions [11,32], as well as results on central Au+Au collisions from AGS [38,39] and RHIC [40][41][42][43][44], as processed in Ref. [3]. Figure 21 shows the ratio of transverse mass spectra of π − mesons produced at mid-rapidity (0 < y < 0.2) in central Pb+Pb collisions and p+p interactions at the same collision energy per nucleon. The spectra were normalised to unity before dividing. First, one observes that the ratio is not constant implying that the spectral shapes are different in p+p interactions and central Pb+Pb collisions. Second, it is seen that the ratio depends weakly, if at all, on collision energy. The ratio is higher than unity in the left (m T − m π < 0.1 GeV/c 2 ) and right (m T − m π > 0.5 GeV/c 2 ) parts of the m T range. It is below unity in the central region 0.1 < m T − m π < 0.5 GeV/c 2 .
The inverse slope parameter T of transverse mass spectra fitted in the range 0.2 < m T − m π < 0.7 GeV/c 2 is plotted versus the collision energy in Fig. 22 (left). The T parameter is larger by about 10-20 MeV/c 2 in central Pb+Pb collisions than in p+p interactions.
The transverse mass spectra measured by NA61/SHINE and NA49 allow a reliable calculation of mean transverse mass. A small correction to the measured value for the high m T region not covered by the measurements was applied based on the exponential extrapolation of the tail of the dis-  (11)). The results of NA61/SHINE are indicated by filled circles and the compilation of the world data [11,32] by open circles. The plotted statistical errors are smaller than the symbol size. The systematic uncertainties are indicated by the coloured band. Table 5: Numerical values of the parameters fitted to rapidity (see Eq. (10)) and transverse mass (see Eq. (9)) spectra of π − mesons produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c. In case of the rapidity fit parameters π − , σ , σ 0 and y 0 , the systematic uncertainty dominates. The uncertainties written in the table are the quadrature sum of the statistical and systematic uncertainties. All uncertainties are given numerically in [36]. For T and m T the statistical uncertainty is written first and the systematic one second. tributions. Half of the correction was added to the systematic uncertainty on m T . In spite of the different shapes of the m T spectra the mean transverse mass calculated for p+p interactions and central Pb+Pb collisions is similar, see Fig. 22 (right). This is because the differences shift the mean m T in opposite directions for different regions of m T and as a result leave it almost unchanged. Thus the mean transverse mass appears to be insensitive to the apparent changes of the pion production properties observed between p+p interactions and central Pb+Pb collisions. Figure 23 (left) presents the ratio of the normalised π − rapidity spectra produced in central Pb+Pb and inelastic p+p interactions at the same collision energy per nucleon. The spectra are plotted versus versus the rapidity scaled by the beam rapidity. Only weak, if any, energy dependence of the ratio is observed. Moreover, the ratio is close to unity in the central rapidity region (y/y beam < 0.6), whereas it is higher closer to beam rapidity (y/y beam > 0.6).
Consequently the r.m.s. width σ of rapidity distributions of π − mesons produced in p+p interactions is smaller than  the width in central Pb+Pb collisions. This is seen in Fig. 23 (right) where the energy dependence of σ is plotted. Additionally, p+p data from [11,35,45] are shown; they agree with the NA61/SHINE results.
Note, that when interpreting differences between results obtained for inelastic p+p interactions and central Pb+Pb collisions the isospin effects should be taken into account. This concerns both the spectra as well as the total multiplicities [11].
In order to reduce their influence the mean multiplicity of pions is obtained from a sum of mean multiplicities of negatively and positively charged pions using the phenomenological formula [15]: The results divided by the mean number of wounded nucleons (N W = 2 for p+p) are shown in Fig. 24 as a function of the Fermi energy measure F. The value of π + for the NA61/SHINE results on inelastic p+p interactions was estimated from the measured π − multiplicity assuming π + = π − + 2/3. This assumption is based on the compilation of the world data presented in Ref. [15] and the model presented therein. At beam momenta lower than 40A GeV/c the π / N W ratio is higher in p+p interactions than in central Pb+Pb collisions. The opposite relation holds for beam momenta higher than 40A GeV/c. The energy dependence for inelastic p+p interactions crosses the one for central Pb+Pb (Au+Au) collisions at about 40A GeV/c.

Summary
We presented experimental results on inclusive spectra and mean multiplicities of negatively charged pions produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c. Two dimensional spectra in transverse momentum and rapidity and parameters characterizing them were given. The results agree with existing sparse measurements, extend their range, accuracy and depth of detail.
The results on inelastic p+p interactions were compared with the corresponding data on central Pb+Pb collisions obtained by NA49. The spectra in p+p interactions are narrower both in rapidity and in transverse mass, which might be attributed to isospin effects. The mean pion multiplicity per wounded nucleon in p+p interactions increases more slowly with energy in the SPS range and crosses the corresponding dependence measured in the Pb+Pb collisions at about 40A GeV/c. Inverse slope parameter T of the transverse mass spectra at mid-rapidity (0 < y < 0.2) plotted against the collision energy per nucleon. The parameter T was fitted in the range 0.2 < m T − m π < 0.7 GeV/c 2 . The systematic uncertainty for the two lowest energy points for Pb+Pb, not given in [3] is assumed to be the same as for the higher energies [2]. Right: Mean transverse mass m T at mid-rapidity (0 < y < 0.2) versus the collision energy. The results on inelastic p+p interactions are compared with the corresponding data on central Pb+Pb (Au+Au) collisions.