1 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 30\(A\) 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 long-baseline neutrino oscillation experiments [58] 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 comparison 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  [1113]. 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 \(\pi ^{-}\) 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 \(\pi ^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\tau = 7.8\) m). A significant fraction of positively charged hadrons are protons (25 %) and kaons (5 %) [1113]. Therefore measurements of \(\pi ^{+}\) mesons require their identification by measurements of the energy loss and/or time-of-flight (tof). This identification is not as important for \(\pi ^{-}\) mesons because the contamination of negatively charged particles by \({\mathrm {K}^-}\) mesons and anti-protons is below 10 % [1113] and can be subtracted reliably. The latter method is used in this paper and it allows to derive \(\pi ^{-}\) 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 Sect. 2 the NA61/SHINE experimental set-up is described. Details on the beam, trigger and event selection are given in Sect. 3. Data reconstruction, simulation and detector performance are described in Sect. 4. Analysis techniques and final results are presented in Sects. 5 and 6, respectively. These results are compared with the corresponding data on central Pb+Pb collisions in Sect. 7. A summary in Sect. 8 closes the paper.

The pion rapidity is calculated in the collision centre of mass system: \(y = \mathrm{atanh}(\beta _\mathrm {L})\), where \(\beta _\mathrm {L} = p_\mathrm {L}/E\) is the longitudinal component of the velocity, \(p_\mathrm {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_\mathrm {T}}}\) and the transverse mass \({{m_\mathrm {T}}}\) is defined as \({{m_\mathrm {T}}} = \sqrt{m_{\pi }^2 + {{p_\mathrm {T}}}^2}\), where \(m_{\pi }\) is the charged pion mass. The nucleon mass and collision energy per nucleon pair in the centre of mass system are denoted as \(m_\mathrm {N}\) and \(\sqrt{s_\mathrm {NN}}\), respectively.

2 The NA61/SHINE facility

The NA61/SHINE experimental facility [20] 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.

Fig. 1
figure 1

The schematic layout of the NA61/SHINE experiment at the CERN SPS (horizontal cut, not to scale). 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 in the figure; 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

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 Sect. 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 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, \(\mathrm {d}E/\mathrm {d}x \), is augmented by tof measurements using Time-of-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 LHT (see Sect. 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 TPCs, the BPDs as well as from the beam and trigger counters.

3 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.

Table 1 Basic beam properties and numbers of events recorded for p+p interactions at 20, 31, 40, 80 and 158 GeV/c. The first column gives the beam momentum. The second and third columns list typical numbers of beam particles at NA61/SHINE per spill (about 10 s) and the fraction of protons in the beam, respectively

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 \(\mathrm {d}E/\mathrm {d}x \) 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\(^\mathrm {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 \(\mathrm{S1} \wedge \mathrm{S2} \wedge \overline{\mathrm{V0}} \wedge \overline{\mathrm{V1}} \wedge \overline{\mathrm{V1}^\mathrm{p}} \wedge \mathrm{CEDAR} \wedge \overline{\mathrm{THC}}\). Trajectories of individual beam particles were measured in a telescope of BPDs placed along the beam line (BPD-1/2/3 in Fig. 1). These counters are small (\(4.8\times 4.8\) cm\(^2\)) proportional chambers with cathode strip readout, providing a resolution of about 100 \(\upmu \)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.

Fig. 2
figure 2

Top The beam spot as measured by BPD-3 after the \(\overline{\mathrm{V1 }}\) cut described in the text for 20 GeV/c (left) and 158 GeV/c (right) beams. Bottom The beam divergence in \(x\) and \(y\) for 20 GeV/c (left) and 158 GeV/c (right) beams. All distributions were arbitrarily scaled to the full colour scale. Widths of the distributions are given in the legend

For data taking on p+p interactions a LHT 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 \(\rho _\mathrm {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 Fig. 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.

4 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:

  1. (i)

    detector geometry, TPC drift velocities and distortions due to the magnetic field inhomogeneities in the corners of VTPCs,

  2. (ii)

    magnetic field setting,

  3. (iii)

    specific energy loss measurements,

  4. (iv)

    tof 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:

  1. (i)

    cluster finding in the TPC raw data, calculation of the cluster centre-of-gravity and total charge,

  2. (ii)

    reconstruction of local track segments in each TPC separately,

  3. (iii)

    matching of track segments into global tracks,

  4. (iv)

    track fitting through the magnetic field and determination of track parameters at the first measured TPC cluster,

  5. (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),

  6. (vi)

    refitting the particle trajectory using the interaction vertex as an additional point and determining the particle momentum at the interaction vertex,

  7. (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.

Fig. 3
figure 3

An example of a p+p interaction at 40 GeV/c measured in the NA61/SHINE detector. The measured points (green) are used to fit tracks (red lines) to the interaction point. The grey dots show the noise clusters. Due to the central gap of the VTPCs only a small part of the trajectory of the negatively charged particle is seen in VTPC-1

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 \(\pi \)+C interactions: FLUKA2008, URQMD1.3.1, VENUS4.12, EPOS1.99, GHEISHA2002, QGSJetII-3 and Sibyll2.1 [17, 2224]. 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):

  1. (i)

    generation of inelastic p+p interactions using the EPOS model [25],

  2. (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,

  3. (iii)

    simulation of the detector response using dedicated NA61/SHINE packages which introduce distortions corresponding to all corrections applied to the real data,

  4. (iv)

    simulation of the interaction trigger selection by checking whether a charged particle hits the S4 counter, see Sect. 3,

  5. (v)

    storage of the simulated events in a file which has the same format as the raw data,

  6. (vi)

    reconstruction of the simulated events with the same reconstruction chain as used for the real data and

  7. (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 Sect. 5) can be applied only for inelastic events. The contribution of elastic events is removed by the event selection cuts (see Sect. 5.1), whereas the contribution of off-target interactions is subtracted based on the data (see Sect. 5.4).

Spectra of \(\pi ^{-}\) mesons presented in this paper were derived from spectra of all negatively charged hadrons corrected for a small contamination of mostly \({\mathrm {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 \(\mathrm {d}E/\mathrm {d}x \) and tof measurements. They are limited due to the geometrical acceptance of the ToF detectors, the finite resolution of the \(\mathrm {d}E/\mathrm {d}x \) measurements and limited data statistics.

Fig. 4
figure 4

Typical acceptance regions for \(\pi ^{-}\) meson spectra in p+p interactions at 20 GeV/c (left) and 158 GeV/c (right) for different analysis methods: the method used in this paper which does not require an explicit pion identification, the method which identifies pions via their energy loss (\(\mathrm {d}E/\mathrm {d}x \)) and, in addition, their time-of-flight (tof)

The quality of measurements was studied by reconstructing masses of \({\mathrm {K}^0_\mathrm {S}}\) particles from their V\(^0\) decay topology. As an example the invariant mass distributions of \({\mathrm {K}^0_\mathrm {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 \({\mathrm {K}^0_\mathrm {S}}\) mass [28] are smaller than 1 \(\text {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 Sect. 5.7.2).

Fig. 5
figure 5

Invariant mass distribution of reconstructed \({\mathrm {K}^0_\mathrm {S}}\) candidates in p+p interactions at 20 (top) and 158 GeV/c (bottom) for the measured data and EPOS model based Monte Carlo simulations. The MC plot was normalised to the peak height of the data. The \({\mathrm {K}^0_\mathrm {S}}\) candidates were selected within \(0<y<-1\) and \(0<{{p_\mathrm{T}}}<0.5\,\text {GeV/}{ c} \) for 20 GeV/c and \(-1<y<0\) and \(0<{{p_\mathrm{T}}}<0.5\,\text {GeV/}{ c} \) for 158 GeV/c. The distribution was fitted with the sum of a Lorentzian function (signal) and a second order polynomial (background)

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 Sect. 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_\mathrm {T}}}\). These results were obtained for negatively charged pions passing the track selection criteria described in Sect. 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.

Fig. 6
figure 6

Reconstruction efficiency of negatively charged pions produced in p+p interactions at 20 (left) and 158 GeV/c (right) as a function of rapidity and transverse momentum. It was calculated by dividing the number of tracks passing the track selection cuts specified in Sect. 5.2 by the number of the generated tracks. The selection criteria include the requirement of at least 90 % reconstruction efficiency

Fig. 7
figure 7

Resolution of rapidity (top, scaled by \(10^3\)) and transverse momentum (bottom) measurements for negatively charged pions produced in p+p interactions at 20 (left) and 158 GeV/c (right) as a function of pion rapidity and transverse momentum. The results are obtained using the track selection cuts specified in Sect. 5.2

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 %.

Fig. 8
figure 8

Distribution of fitted vertex \(z\) coordinate for p+p interactions at 20 (top) and 158 GeV/c (bottom). The black line shows the data after target removed subtraction (see Sect. 5.4). The filled area shows the distribution for the reconstructed Monte Carlo simulation. This distribution was normalised to the total integral of the data plot. The dashed vertical lines show the \(z\) vertex selection range

Fig. 9
figure 9

Distribution of the impact parameter in the \(x\) (left) and \(y\) (right) coordinate for p+p interactions at 20 (top) and 158 GeV/c (bottom). The black line shows the data after target removed subtraction (see Sect. 5.4). The filled area shows the reconstructed Monte Carlo simulation. The dashed vertical lines show the accepted range (see Sect. 5.2). The dotted blue line in the 20 GeV/c plots show the distribution obtained using the alternative event selection. Namely only events with no off-time beam particles within the time window of \(\pm 6\ \upmu \)s around the trigger particle time were accepted

5 Analysis technique

This section presents the procedures used for data analysis consisting of the following steps:

  1. (i)

    applying event and track selection criteria,

  2. (ii)

    determination of spectra of negatively charged hadrons using the selected events and tracks,

  3. (iii)

    evaluation of corrections to the spectra based on experimental data and simulations,

  4. (iv)

    calculation of the corrected spectra.

Corrections for the following biases were evaluated and applied:

  1. (i)

    geometrical acceptance,

  2. (ii)

    contribution of off-target interactions,

  3. (iii)

    contribution of particles other than negatively charged pions produced in inelastic p+p interactions,

  4. (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 \(\pi ^{-}\) 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 \(\pi ^{-}\). The term primary will be used in the above meaning also for other particles.

The analysis was performed independently in (\(y\), \({{p_\mathrm {T}}}\)) and (\(y\), \({{m_\mathrm {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.

5.1 Event selection criteria

This section presents the event selection criteria. The number of events selected by the trigger (see Sect. 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.

Table 2 Number of events recorded with the interaction trigger (all) and selected for the analysis (selected)

The following event selection criteria were applied to the events recorded with the interaction trigger:

  1. (i)

    no off-time beam particle is detected within \(\pm 2\ \upmu \)s around the trigger particle,

  2. (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,

  3. (iii)

    there is at least one track reconstructed in the TPCs and fitted to the interaction vertex,

  4. (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,

  5. (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.

Table 3 Momentum thresholds used to reject elastic interactions (cut (v))

The off-line (listed above) and on-line (the interaction trigger condition, see Sect. 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).

Fig. 10
figure 10

Momentum distributions at 20 (left), 40 (middle) and 158 GeV/c (right) of the positively charged tracks in events passing selection cuts (i)–(iv), containing a single track, which is positively charged and measured in the GAP TPC and MTPC. The distributions were normalised to all events. The vertical dashed lines at 20 and 40 GeV/c show the momentum threshold used to remove elastic events (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 Sects. 5.4 and 5.6.

5.2 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:

  1. (i)

    the track momentum fit at the interaction vertex should have converged,

  2. (ii)

    the fitted track charge is negative,

  3. (iii)

    the fitted track momentum component \(p_x\) is negative. This selection minimises the angle between the track trajectory and the TPC pad direction for the chosen magnetic field direction. This reduces statistical and systematic uncertainties of the cluster position, energy deposit and track parameters,

  4. (iv)

    the total number of reconstructed points on the track should be greater than 30,

  5. (v)

    the sum of the number of reconstructed points in VTPC-1 and VTPC-2 should be greater than 15 or the number of reconstructed points in the GAP-TPC should be greater than 4,

  6. (vi)

    the distance between the track extrapolated to the interaction plane and the interaction point (impact parameter) should be smaller than 4 cm in the horizontal (bending) plane and 2 cm in the vertical (drift) plane,

  7. (vii)

    the track should be measured in a high (\(\ge \)90 %) TPC acceptance region (see Sect. 5.3),

  8. (viii)

    tracks with \(\mathrm {d}E/\mathrm {d}x \) and total momentum values characteristic 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.

    Fig. 11
    figure 11

    Left Distribution of particle energy loss as a function of the logarithm of total momentum, for negatively charged particles produced in p+p interactions at 40 GeV/c. The black contour shows the electron selection region. Right \(\mathrm {d}E/\mathrm {d}x \) distribution in the momentum range indicated in the figure and marked with vertical dashed lines in the left panel. Electrons and negatively charged hadrons show separated peaks

The spectra of negatively charged particles after track and event selections were obtained in 2-dimensional bins of (\(y\), \({{p_\mathrm {T}}}\)) and (\(y\), \({{m_\mathrm {T}}}\)). The spectra were evaluated in the centre-of-mass frame after rotation of the \(z\) axis into the proton beam direction measured event-by-event by the BPDs.

5.3 Correction for detector acceptance

The detection and reconstruction inefficiencies are corrected using the simulation described in Sect. 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 \(\pi ^{-}\) tracks passing the track selection criteria to the number of generated \(\pi ^{-}\)) equals at least 90 %. These regions were identified using a separate, statistically independent simulation in three-dimensional bins of \(y\), \({{p_\mathrm {T}}}\) or \({{m_\mathrm {T}}}\) and the azimuthal angle \(\phi \) (\(5^\circ \) 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_\mathrm {T}}}\) region.

Fig. 12
figure 12

Detector acceptance at 20 (left) and 158 GeV/c (right). The acceptance is calculated as fraction of (\(y\), \({{p_\mathrm {T}}}\), \(\phi \)) bins accepted for given \(y\) and \({{p_\mathrm {T}}}\) for tracks with \(p_x<0\) selected for this analysis

Since neither target nor beam are polarized, we can assume a uniform distribution of particles in \(\phi \). The data falling into the accepted bins is summed over \(\phi \) bins and the (\(y\), \({{p_\mathrm {T}}}\)/\({{m_\mathrm {T}}}\)) bin content is multiplied by a correction factor to compensate for the rejected \(\phi \) ranges. The acceptance correction also compensates for the \(p_x < 0\) selection (see Sect. 5.2, point (5.2)).

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.

5.4 Correction for off-target interactions

The spectra were derived for events with liquid hydrogen in (\(_\mathrm {I}\)) and removed (\(_\mathrm {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 _\mathrm {T}[{\mathrm {h}^-}]\)) is calculated as:

$$\begin{aligned} n _\mathrm {T}[{\mathrm {h}^-}] = \frac{1}{1-\varepsilon } \cdot \left( n _\mathrm {I}[{\mathrm {h}^-}] - \varepsilon \cdot n _\mathrm {R}[{\mathrm {h}^-}] \right) , \end{aligned}$$
(1)

where:

  1. (i)

    \(n _\mathrm {I}[{\mathrm {h}^-}]\) and \(n _\mathrm {R}[{\mathrm {h}^-}]\) is the number of tracks in a given bin per event selected for the analysis (see Sect. 5.2) for the data with the liquid hydrogen inserted and removed, respectively,

  2. (ii)

    \(\varepsilon \) is the ratio of the interaction probabilities for the removed and inserted target operation.

\(\varepsilon \) 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 beam-line and detector materials. Neglecting the beam attenuation in the target one gets:

$$\begin{aligned} \epsilon = \frac{N_\mathrm {R}}{N_\mathrm {I}} \cdot \frac{N_\mathrm {I}[z>-450\mathrm {\ cm}]}{N_\mathrm {R}[z>-450\mathrm {\ cm}]}\ , \end{aligned}$$
(2)

where \(N[z>-450\mathrm {\ 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 \(\epsilon \) are listed in Table 4.

Table 4 The ratio of the interaction probabilities, \(\varepsilon \), for the removed and inserted target operation for data taking on p+p interactions at 20, 31, 40, 80 and 158 GeV/c
Fig. 13
figure 13

Distribution of fitted vertex \(z\) coordinate for data on p+p interactions at 40 GeV/c. The distribution for the data recorded with the removed liquid hydrogen was multiplied by a factor of \(N_\mathrm {I}[z>-450\mathrm {\ cm}]/ N_\mathrm {R} [z>-450\mathrm {\ cm}]\)

The correction for the off-target interactions changes the yields obtained from the target inserted data by less than \(\pm 5~\%\), except in the regions where the statistical uncertainty is high.

5.5 The correction for contamination of primary \(\pi ^{-}\) mesons

More than 90 % of primary negatively charged particles produced in p+p interactions in the SPS energy range are \(\pi ^{-}\) mesons [1113]. Thus \(\pi ^{-}\) 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 Sect. 4 was used to calculate corrections for the admixture of particles other than primary \(\pi ^{-}\) mesons to the reconstructed negatively charged particles. The dominating contributions are primary \({\mathrm {K}^-}\) and \({\bar{p }} \), and secondary \(\pi ^{-}\) from weak decays of \({\mathrm {\Lambda }}\) and \({\mathrm {K}^0_\mathrm {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 \({\mathrm {K}^-}\), and \({\bar{p }} \) [30]. EPOS spectra of \(\pi ^{-}\) were replaced by the preliminary NA61/SHINE results [31] normalised to the multiplicity from the world data compilation [32]. Spectra of \({\mathrm {\Lambda }}\) and \({\mathrm {K}^0_\mathrm {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_\mathrm {T}}}\) region at the low beam momenta, where it reaches \(+20~\%\).

As it was found in [24] the yields of \({\mathrm {K}^-}\) and \({\bar{p }} \) are strongly correlated with the \(\pi ^{-}\) yield. Thus the correction for the contribution of primary hadrons is performed via the multiplicative factor \(c_\mathrm {K}\). On the contrary the contribution due to weak decays and secondary interactions is mostly located in the low \({{p_\mathrm {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_\mathrm {V}\). The total correction is applied in as:

$$\begin{aligned} n _\mathrm {prim}[\pi ^{-}] = \left( n _\mathrm {T}[{\mathrm {h}^-}] - c_\mathrm {V} \right) \cdot c_\mathrm {K}, \end{aligned}$$
(3)

where

$$\begin{aligned} c_\mathrm {V} = \left( n [\pi ^{-}_{\mathrm {\Lambda }}]+n [\pi ^{-}_{{\mathrm {K}^0_\mathrm {S}}}]+n [\mathrm {other}] \right) ^\mathrm {MC}_\mathrm {sel}, \end{aligned}$$
(4)
$$\begin{aligned} c_\mathrm {K} = \left( \frac{n [\pi ^{-}]}{n [{\mathrm {K}^-}]+n [{\bar{p }} ]+n [\pi ^{-}]} \right) ^\mathrm {MC}_\mathrm {sel}. \end{aligned}$$
(5)

The spectrum of a particle \(x\) is denoted as \(n [x]\) whereas \(n [\mathrm {other}]\) stands for all primary and secondary particles other than \({\mathrm {K}^-}\), \({\bar{p }} \), \(\pi ^{-}\) and feed-down from \({\mathrm {\Lambda }}\) and \({\mathrm {K}^0_\mathrm {S}}\). The spectrum \(n [\mathrm {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 \(^{\mathrm {MC}}\) marks adjusted EPOS spectra. The subscript \(_{\mathrm {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.

5.6 Correction for event as well as track losses and migration

The multiplicative correction \(c_\mathrm {loss}\) for losses of inelastic events as well as losses and bin-to-bin migration of primary \(\pi ^{-}\) mesons emitted within the acceptance is calculated using the Monte Carlo simulation as:

$$\begin{aligned} c_\mathrm {loss} = n [\pi ^{-}]_\mathrm {gen}^\mathrm {MC}\ /\ n [\pi ^{-}]_\mathrm {sel}^\mathrm {MC}, \end{aligned}$$
(6)

where the subscript \(_{\mathrm {gen}}\) indicates the generated spectrum of primary \(\pi ^{-}\) mesons binned according to the generated momentum vector. Then the final, corrected \(\pi ^{-}\) meson spectrum in inelastic p+p interactions is calculated as

$$\begin{aligned} n [\pi ^{-}] = c_\mathrm {loss}\cdot n _\mathrm {prim}[\pi ^{-}]\ . \end{aligned}$$
(7)

The dominating effects contributing to the \(c_\mathrm {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.

5.7 Statistical and systematic uncertainties

5.7.1 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.

5.7.2 Systematic uncertainties

Systematic errors presented in this paper were calculated taking into account contributions from the following effects.

  1. (i)

    Possible biases due to event and track cuts which are not corrected for. These are:

    • a possible bias due to the \(\mathrm {d}E/\mathrm {d}x \) cut applied to remove electrons,

    • a possible bias related to the removal of events with off-time beam particles close in time to the trigger particle.

    The magnitude \(\sigma _\mathrm {i}\) of possible biases was estimated by varying values of the corresponding cuts. The \(\mathrm {d}E/\mathrm {d}x \) cut was changed by \(\pm 0.01\) \(\mathrm {d}E/\mathrm {d}x \) 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 \(\pm 1~\upmu \)s to a \(\pm 3~\upmu \)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.

    Fig. 14
    figure 14

    Statistical and systematic uncertainties in selected bins of \({{p_\mathrm {T}}}\) for 20 (left) and 158 GeV/c (right) p+p data. The shaded band shows the statistical uncertainty. The coloured thin lines show the contributions to the systematic uncertainty listed in Sect. 5.7.2. The thick black lines show the total systematic uncertainty, which was calculated by adding the contributions in quadrature (continuous line) or linearly (dashed/dotted line, shown for comparison)

  2. (ii)

    Uncertainty of the correction for contamination of the primary \(\pi ^{-}\) mesons. The systematic uncertainty \(\sigma _\mathrm {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 dashed-dotted blue line in Fig. 14. The absolute correction is small thus the related systematic uncertainty is small also.

  3. (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, 35] model:

    $$\begin{aligned}&\sigma _\mathrm {iii} = 0.2\cdot \left| 1 - c_\mathrm {loss}^\mathrm {EPOS}\right| + \left| c_\mathrm {loss}^\mathrm {EPOS}-c_\mathrm {loss}^\mathrm {VENUS}\right| . \end{aligned}$$
    (8)

    This contribution is drawn with a short-dashed green line in Fig. 14.

  4. (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_\mathrm {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 \(\sigma _\mathrm {sys} = \sqrt{\sigma _\mathrm {i}^2+\sigma _\mathrm {ii}^2+\sigma _\mathrm {iii}^2}\). 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 \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 158 GeV/c (for details see Sect. 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.

Fig. 15
figure 15

Rapidity distribution of \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 158 GeV/c. The big blue points show the results obtained with an alternative method: without vertex fit requirement and rejection of events with a single very high momentum positively charged track. The results of NA61/SHINE (this paper, red dots) are compared with the NA49 measurements [11] (black squares). The open symbols show points reflected with respect to mid-rapidity. A single NA61/SHINE point measured at \(y<0\) is also shown for comparison. The shaded band shows the NA61/SHINE systematic uncertainty

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 \(\mathrm {d}E/\mathrm {d}x \) 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  [36, 37]. A parametrisation of the distribution and the total \(\pi ^{-}\) 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.

Fig. 16
figure 16

Rapidity distribution of \(\pi ^{-}\) mesons produced in inelastic p+p interactions. The NA61/SHINE results at 31 GeV/c (blue points) are compared with the MIRABELLE measurement (parametrised by the black line) at 32 GeV/c. The shaded band shows the NA61/SHINE systematic uncertainty

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_\mathrm {T}}}\) acceptance extends to zero.

6 Results

This section presents results on inclusive \(\pi ^{-}\) 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. [38].

6.1 Double differential spectra

The double differential inclusive spectra of \(\pi ^{-}\) 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 \(\frac{\mathrm {d}^2n}{\mathrm {d}y\,\mathrm {d}{{p_\mathrm {T}}}}\) or \(\frac{\mathrm {d}^2n}{\mathrm {d}y\,\mathrm {d}{{m_\mathrm {T}}}}\) were calculated by dividing the fully corrected bin contents \(n [\pi ^{-}]\) (see Sect. 5) by the bin size. The spectra in (\(y\), \({{m_\mathrm {T}}}\)) are not shown here but they are given in the compilation of the numerical values [38].

Fig. 17
figure 17

Double differential spectra \(\mathrm {d}^2n/(\mathrm {d}y\,\mathrm {d}{{p_\mathrm{T}}})[(\text {GeV}/{c})^{-1}]\) of \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c

Fig. 18
figure 18

Transverse momentum spectra of \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 20 (left) and 158 GeV/c (right) in various rapidity ranges. The legend provides the centres of the rapidity bins, \(y_\mathrm {bin}\) and the scaling factor \(c\) used to separate the spectra visually

6.2 Transverse mass spectra

The transverse mass spectra at mid-rapidity (\(0 < y < 0.2\)) are shown in Fig. 19 (left). A function

$$\begin{aligned} \frac{\mathrm {d}n}{\mathrm {d}{{m_\mathrm {T}}}} = A\cdot {{m_\mathrm {T}}}\cdot \exp \left( -\frac{{{m_\mathrm {T}}}}{T}\right) \end{aligned}$$
(9)

was fitted in the range \(0.2 < {{m_\mathrm {T}}} - m_\pi < 0.7\,\text {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 \(\chi ^2\) function which was calculated using statistical errors only. In the \(\chi ^2\) calculation a measured bin content (\(\mathrm {d}n/\mathrm {d}{{m_\mathrm {T}}}\)) was compared with the integral of the fitted function in a bin divided by the bin width.

Fig. 19
figure 19

Left Transverse mass spectra at mid-rapidity (\(0 < y < 0.2\)). The fitted exponential function Eq. (9) is indicated by solid lines in the fit range \(0.2 < {{m_\mathrm {T}}} - m_\pi < 0.7\,\text {GeV/}c^2 \) and dashed lines outside the fit range. The data points for different beam momenta were scaled for better readability. Right The inverse slope parameter \(T\) of the transverse mass spectra as a function of rapidity divided by the projectile rapidity. The fit range is \(0.2 <{{m_\mathrm {T}}} - m_\pi < 0.7\,\text {GeV/}c^2 \). The results refer to \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c

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_\mathrm {beam}\) = 1.877, 2.094, 2.223, 2.569 and 2.909 at 20, 31, 40, 80 and 158 GeV/c, respectively).

6.3 Rapidity spectra

The rapidity spectra are shown in Fig. 20 (left). They were obtained by summing the measured \({{m_\mathrm {T}}}\) spectra and using the exponential function Eq. (9). The function was fitted in the range ending at the maximum measured \({{m_\mathrm {T}}}\), and starting 0.9 \(\text {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.

Fig. 20
figure 20

Left Rapidity spectra obtained from sums of the measured and extrapolated \({{m_\mathrm {T}}}\) spectra. Closed symbols indicate measured points, open points are reflected with respect to mid-rapidity. The measured points at \(y<0\) are shown for systematic comparison only. The plotted statistical errors are smaller than the symbol size. The systematic uncertainties are indicated by the coloured bands. The lines indicate fits of the sum of two symmetrically displaced Gaussian functions (see Eq. (10)) to the spectra. The results refer to \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c. Right Dependence of the mean total multiplicity of \(\pi ^{-}\) mesons produced in inelastic p+p interactions on Fermi’s energy measure \(F\) (see Eq. (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

The rapidity spectra are parametrised by the sum of two Gaussian functions symmetrically displaced with respect to mid-rapidity:

$$\begin{aligned} \frac{\mathrm {d}n}{\mathrm {d}y}&= \frac{\langle {\pi ^{-}}\rangle (y_0,\sigma _0)}{2\sigma _0\sqrt{2\pi }}\cdot \nonumber \\&\cdot \left[ \exp \left( -\frac{(y-y_0)^2}{2\sigma _0^2}\right) +\exp \left( -\frac{(y+y_0)^2}{2\sigma _0^2}\right) \right] ,\nonumber \\ \end{aligned}$$
(10)

where \(y_0\) and \(\sigma _0\) are fit parameters, and the total multiplicity \(\langle {\pi ^{-}}\rangle (y_0,\sigma _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 \(\chi ^2\) function was minimised in a similar way as in case of the \({{m_\mathrm {T}}}\) spectra, namely using the integral of the function in a given bin. The numerical values of the fitted parameters as well as the r.m.s. width \(\upsigma = \sqrt{y_0^2 + \upsigma _0^2}\) are given in Table 5.

Table 5 Numerical values of the parameters fitted to rapidity (see Eq. (10)) and transverse mass (see Eq. (9)) spectra of \(\pi ^{-}\) mesons produced in inelastic p+p interactions at 20, 31, 40, 80 and 158 GeV/c. In case of the rapidity fit parameters \(\langle {\pi ^{-}}\rangle , \, \sigma , \, \sigma _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 [38]. For \(T\) and \(\langle {{m_\mathrm {T}}}\rangle \) the statistical uncertainty is written first and the systematic one second

6.4 Mean multiplicities

Mean multiplicities of \(\pi ^{-}\) mesons, \(\langle {\pi ^{-}}\rangle \), 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 \(\langle {\pi ^{-}}\rangle \) multiplicity per inelastic p+p collision on the Fermi’s energy measure [39],

$$\begin{aligned} F \equiv \left[ \frac{(\sqrt{s_\mathrm {NN}} - 2 m_\mathrm {N})^3}{\sqrt{s_\mathrm {NN}}} \right] ^{1/4} \, \end{aligned}$$
(11)

is plotted in Fig. 20 (right). The results of NA61/SHINE are in agreement with a compilation of the world data [11, 32].

7 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 [40, 41] and RHIC [4246], as processed in Ref. [3].

Figure 21 shows the ratio of transverse mass spectra of \(\pi ^{-}\) 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_\mathrm {T}}} - m_{\pi } < 0.1\,\text {GeV/}c^2 \)) and right (\({{m_\mathrm {T}}} - m_{\pi } > 0.5\,\text {GeV/}c^2 \)) parts of the \({{m_\mathrm {T}}}\) range. It is below unity in the central region \(0.1 < {{m_\mathrm {T}}} - m_{\pi } < 0.5\,\text {GeV/}c^2 \).

Fig. 21
figure 21

The ratio of the normalised transverse mass spectra of \(\pi ^{-}\) mesons at mid-rapidity produced in central Pb+Pb collisions and inelastic p+p interactions at the same collision energy per nucleon. The coloured bands represent the systematic uncertainty of the p+p data

The inverse slope parameter \(T\) of transverse mass spectra fitted in the range \(0.2 < {{m_\mathrm {T}}} - m_\pi < 0.7\,\text {GeV/}c^2 \) is plotted versus the collision energy in Fig. 22 (left). The \(T\) parameter is larger by about 10–20\(\text {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_\mathrm {T}}}\) region not covered by the measurements was applied based on the exponential extrapolation of the tail of the distributions. Half of the correction was added to the systematic uncertainty on \(\langle {{{m_\mathrm {T}}}}\rangle \). In spite of the different shapes of the \({{m_\mathrm {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_\mathrm {T}}}\) in opposite directions for different regions of \({{m_\mathrm {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.

Fig. 22
figure 22

Left 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_\mathrm {T}}}-m_\pi < 0.7\,\text {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 \(\langle {{{m_\mathrm {T}}}}\rangle \) 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

Figure 23 (left) presents the ratio of the normalised \(\pi ^{-}\) rapidity spectra produced in central Pb+Pb and inelastic p+p interactions at the same collision energy per nucleon. The spectra are plotted 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_\mathrm {beam} < 0.6\)), whereas it is higher closer to beam rapidity (\(y/y_\mathrm {beam} > 0.6\)).

Consequently the r.m.s. width \(\sigma \) of rapidity distributions of \(\pi ^{-}\) 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 \(\sigma \) is plotted. Additionally, p+p data from [11, 36, 37, 47] are shown; they agree with the NA61/SHINE results.

Fig. 23
figure 23

Left The ratio of normalised rapidity spectra of \(\pi ^{-}\) mesons produced in central Pb+Pb collisions and inelastic p+p interactions at the same collision energy per nucleon plotted versus the rapidity scaled by the beam rapidity. The coloured bands represent the NA61/SHINE systematic uncertainty. Right Energy dependence of the width of the rapidity distribution of \(\pi ^{-}\) mesons produced in p+p interactions and central Pb+Pb collisions. The systematic uncertainty for the Pb+Pb points is not given

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]:

$$\begin{aligned} \langle {\pi }\rangle = \frac{3}{2}\left( \langle {\pi ^{+}}\rangle + \langle {\pi ^{-}}\rangle \right) . \end{aligned}$$
(12)

The results divided by the mean number of wounded nucleons (\(N_\mathrm {W} = 2\) for p+p) are shown in Fig. 24 as a function of the Fermi energy measure \(F\). The value of \(\langle {\pi ^{+}}\rangle \) for the NA61/SHINE results on inelastic p+p interactions was estimated from the measured \(\langle {\pi ^{-}}\rangle \) multiplicity assuming \(\langle {\pi ^{+}}\rangle = \langle {\pi ^{-}}\rangle + 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 40\(A\) GeV/c the \(\langle {\pi }\rangle /\langle {N_\mathrm {W}}\rangle \) ratio is higher in p+p interactions than in central Pb+Pb collisions. The opposite relation holds for beam momenta higher than 40\(A\) GeV/c. The energy dependence for inelastic p+p interactions crosses the one for central Pb+Pb (Au+Au) collisions at about 40\(A\) GeV/c.

Fig. 24
figure 24

Mean multiplicity of all pions per wounded nucleon produced in inelastic p+p interactions and central Pb+Pb (Au+Au) collisions. The vertical lines show the total uncertainty

8 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 40\(A\) GeV/c.