1 Introduction

The study of hadron production has a long history in high-energy particle and nuclear physics, as well as in cosmic-ray physics. The absolute yields and the transverse momentum (\(p_{\mathrm {T}} \)) spectra of identified hadrons in high-energy hadron–hadron collisions are among the most basic physical observables. They can be used to test the predictions for non-perturbative quantum chromodynamics (QCD) processes like hadronization and soft-parton interactions, and the validity of their implementation in Monte Carlo (MC) event generators. Spectra of identified particles in proton–nucleus collisions also constitute an important reference for studies of high-energy heavy-ion collisions, where final-state effects are known to modify the spectral shape and yields of different hadron species [17].

The present analysis focuses on the measurement of the \(p_{\mathrm {T}} \) spectra of charged hadrons, identified mostly via their energy deposits in silicon detectors, in pPb collisions at \(\sqrt{s_{NN}} =\) 5.02\(\,\text {TeV}\). The analysis procedures are similar to those previously used in the measurement of pion, kaon, and proton production in pp collisions at several center-of-mass energies [8]. Results on \({\pi }\), \(\mathrm {K}\), and \(\mathrm {p}\) production in pPb collisions have been also reported by the ALICE Collaboration [9].

A detailed description of the CMS (Compact Muon Solenoid) detector can be found in Ref. [10]. The CMS experiment uses a right-handed coordinate system, with the origin at the nominal interaction point (IP) and the \(z\) axis along the counterclockwise-beam direction. The pseudorapidity \(\eta \) and rapidity \(y\) of a particle (in the laboratory frame) with energy \(E\), momentum \(p\), and momentum along the \(z\) axis \(p_z\) are defined as \(\eta = -\ln [\tan (\theta /2)]\), where \(\theta \) is the polar angle with respect to the \(z\) axis and \(y = \frac{1}{2}\ln [(E+p_z)/(E-p_z)]\), respectively. The central feature of the CMS apparatus is a superconducting solenoid of 6\(\,\text {m}\) internal diameter. Within the 3.8 T field volume are the silicon pixel and strip tracker, the crystal electromagnetic calorimeter, and the brass/scintillator hadron calorimeter. The tracker measures charged particles within the pseudorapidity range \(|\eta | < 2.4\). It has 1440 silicon pixel and 15 148 silicon strip detector modules, ordered in 13 tracking layers in the \(y\) region studied here. In addition to the barrel and endcap detectors, CMS has extensive forward calorimetry. Steel/quartz-fiber forward calorimeters (HF) cover \(3 < |\eta | < 5\). Beam Pick-up Timing for the eXperiments (BPTX) devices were used to trigger the detector readout. They are located around the beam pipe at a distance of 175\(\,\text {m}\) from the IP on either side, and are designed to provide precise information on the Large Hadron Collider (LHC) bunch structure and timing of the incoming beams.

The reconstruction of charged particles in CMS is bounded by the acceptance of the tracker (\(|\eta | < \) 2.4) and by the decreasing tracking efficiency at low momentum (greater than about 60 % for \(p > 0.05\), 0.10, 0.20, and 0.40\({\,\text {GeV/}c}\) for \(\mathrm {e}\), \({\pi }\), \(\mathrm {K}\), and \(\mathrm {p}\), respectively). Particle identification capabilities using specific ionization are restricted to \(p < 0.15{\,\text {GeV/}c} \) for electrons, \(p < 1.20{\,\text {GeV/}c} \) for pions, \(p < 1.05{\,\text {GeV/}c} \) for kaons, and \(p < 1.70{\,\text {GeV/}c} \) for protons. Pions are identified up to a higher momentum than kaons because of their high relative abundance. In view of the \((y,p_{\mathrm {T}})\) regions where pions, kaons, and protons can all be identified (\(p = p_{\mathrm {T}} \cosh y\)), the band \(-1 < y < 1\) (in the laboratory frame) was chosen for this measurement, since it is a good compromise between the \(p_{\mathrm {T}}\) range and \(y\) coverage.

In this paper, comparisons are made to predictions from three MC event generators. The Hijing [11] event generator is based on a two-component model for hadron production in high-energy nucleon and nuclear collisions. Hard parton scatterings are assumed to be described by perturbative QCD and soft interactions are approximated by string excitations with an effective cross section. In version 2.1 [12], in addition to modification of initial parton distributions, multiple scatterings inside a nucleus lead to transverse momentum broadening of both initial and final-state partons. This is responsible for the enhancement of intermediate-\(p_{\mathrm {T}} \) (2–6\({\,\text {GeV/}c}\)) hadron spectra in proton–nucleus collisions, with respect to the properly scaled spectra of proton–proton collisions (Cronin effect). The Ampt [13] event generator is a multi-phase transport model. It starts from the same initial conditions as Hijing, contains a partonic transport phase, the description of the bulk hadronization, and finally a hadronic rescattering phase. These processes lead to hydrodynamic-like effects in simulated nucleus–nucleus collisions, but not necessarily in proton–nucleus collisions. The latest available version (1.26/2.26) is used. The Epos [14] event generator uses a quantum mechanical multiple scattering approach based on partons and strings, where cross sections and particle production are calculated consistently, taking into account energy conservation in both cases. Nuclear effects related to transverse momentum broadening, parton saturation, and screening have been introduced. The model can be used both for extensive air shower simulations and accelerator physics. Epos Lhc [15] is an improvement of version 1.99 (v3400) and contains a three-dimensional viscous event-by-event hydrodynamic treatment. This is a major difference with respect to the Hijing and Ampt models for proton–nucleus collisions.

2 Data analysis

The data were taken in September 2012 during a 4-h-long pPb run with very low probability of multiple interactions (0.15 % “pileup”). A total of 2.0 million collisions were collected, corresponding to an integrated luminosity of approximately \(1\,\upmu \mathrm {b}^{-1} \). The dominant uncertainty for the reported measurements is systematic in nature. The beam energies were 4\(\,\text {TeV}\) for protons and 1.58\(\,\text {TeV}\) per nucleon for lead nuclei, resulting in a center-of-mass energy per nucleon pair of \(\sqrt{s_{NN}} =\) 5.02\(\,\text {TeV}\). Due to the asymmetric beam energies the nucleon-nucleon center-of-mass in the pPb collisions was not at rest with respect to the laboratory frame but was moving with a velocity \(\beta = -0.434\) or rapidity \(-0.465\). Since the higher-energy proton beam traveled in the clockwise direction, i.e. at \(\theta = \pi \), the rapidity of a particle emitted at \(y_\mathrm{cm}\) in the nucleon-nucleon center-of-mass frame is detected in the laboratory frame with a shift, \(y - y_\mathrm{cm} = -0.465\), i.e. a particle with \(y=0\) moves with rapidity 0.465 in the Pb-beam direction in the center-of-mass system. The particle yields reported in this paper have been measured for laboratory rapidity \(|y | < 1\) to match the experimentally accessible region.

The event selection consisted of the following requirements:

  • at the trigger level, the coincidence of signals from both BPTX devices, indicating the presence of both proton and lead bunches crossing the interaction point; in addition, at least one track with \(p_{\mathrm {T}} > 0.4{\,\text {GeV/}c} \) in the pixel tracker;

  • offline, the presence of at least one tower with energy above 3\(\,\text {GeV}\) in each of the HF calorimeters; at least one reconstructed interaction vertex; beam-halo and beam-induced background events, which usually produce an anomalously large number of pixel hits [16], are suppressed.

The efficiencies for event selection, tracking, and vertexing were evaluated using simulated event samples produced with the Hijing 2.1 MC event generator, where the CMS detector response simulation was based on Geant4 [17]. Simulated events were reconstructed in the same way as collision data events. The final results were corrected to a particle level selection applied to the direct MC output, which is very similar to the data selection described above: at least one particle (proper lifetime \(\tau > 10^{-18}\,\text {s} \)) with \(E > 3\,\text {GeV} \) in the range \(-5 < \eta < -3\) and at least one in the range \(3 < \eta <5\); this selection is referred to in the following as the “double-sided” (DS) selection. These requirements are expected to suppress single-diffractive collisions in both the data and MC samples. From the MC event generators studied in this paper, the DS selection efficiency for inelastic, hadronic collisions is found to be 94–97 %.

The simulated ratio of the data selection efficiency to the DS selection efficiency is shown as a function of the reconstructed track multiplicity in the top panel of Fig. 1. The ratio is used to correct the measured events. The results are also corrected for the fraction of DS events without a reconstructed track. This fraction, as given by the simulation, is about 0.1 %.

Fig. 1
figure 1

Top the ratio of selected events to double-sided (DS) events (ratio of the corresponding efficiencies in the inelastic sample), according to Epos Lhc and Hijing MC simulations, as a function of the reconstructed primary charged-particle multiplicity. Bottom acceptance, tracking efficiency (left scale), and misreconstructed-track rate (right scale) in the range \(|\eta | < 2.4\) as a function of \(p_{\mathrm {T}} \) for positively charged pions, kaons, and protons

The extrapolation of particle spectra into the unmeasured \((y,p_{\mathrm {T}})\) regions is model dependent, particularly at low \(p_{\mathrm {T}}\). A high-precision measurement therefore requires reliable track reconstruction down to the lowest possible \(p_{\mathrm {T}}\). The present analysis extends to \(p_{\mathrm {T}} \approx 0.1{\,\text {GeV/}c} \) by exploiting special tracking algorithms [18], used in previous studies [8, 16, 19], to provide high reconstruction efficiency and low background rate. The charged-pion mass was assumed when fitting particle momenta.

The acceptance of the tracker (\(C_\mathrm{a}\)) is defined as the fraction of primary charged particles leaving at least two hits in the pixel detector. It is flat in the region \(-2 < \eta < 2\) and \(p_{\mathrm {T}} > 0.4{\,\text {GeV/}c} \), and its value is 96–98 % (Fig. 1, bottom panel). The loss of acceptance at \(p_{\mathrm {T}} < 0.4{\,\text {GeV/}c} \) is caused by energy loss and multiple scattering of particles, both depending on the particle mass. Likewise, the reconstruction efficiency (\(C_\mathrm{e}\)) is about 75–85 %, degrading at low \(p_{\mathrm {T}}\), also in a mass-dependent way. The misreconstructed-track rate (\(C_f\)) is very small, reaching 1 % only for \(p_{\mathrm {T}} <\) 0.2\({\,\text {GeV/}c}\). The probability of reconstructing multiple tracks (\(C_\mathrm{m}\)) from a single true track is about 0.1 %, mostly due to particles spiralling in the strong magnetic field of the CMS solenoid. The efficiencies and background rates do not depend on the charged-multiplicity of the event. They largely factorize in \(\eta \) and \(p_{\mathrm {T}}\), but for the final corrections an \((\eta ,p_{\mathrm {T}})\) matrix is used.

The region where pPb collisions occur (beam spot) is measured by reconstructing vertices from many events. Since the bunches are very narrow in the transverse direction, the \(xy\) location of the interaction vertices is well constrained; conversely, their \(z\) coordinates are spread over a relatively long distance and must be determined on an event-by-event basis. The vertex position is determined using reconstructed tracks which have \(p_{\mathrm {T}} > 0.1{\,\text {GeV/}c} \) and originate from the vicinity of the beam spot, i.e. their transverse impact parameters \(d_\mathrm{T}\) satisfy the condition \(d_\mathrm{T} < 3\,\sigma _T\). Here \(\sigma _\mathrm{T}\) is the quadratic sum of the uncertainty in the value of \(d_\mathrm{T}\) and the root-mean-square of the beam spot distribution in the transverse plane. The agglomerative vertex-reconstruction algorithm [20] was used, with the \(z\) coordinates (and their uncertainties) of the tracks at the point of closest approach to the beam axis as input. For single-vertex events, there is no minimum requirement on the number of tracks associated with the vertex, even one-track vertices are allowed. Only tracks associated with a primary vertex are used in the analysis. If multiple vertices are present, the tracks from the highest multiplicity vertex are used. The resultant bias is negligible since the pileup rate is extremely small.

The vertex reconstruction resolution in the \(z\) direction is a strong function of the number of reconstructed tracks and it is always smaller than 0.1\(\,\text {cm}\). The distribution of the \(z\) coordinates of the reconstructed primary vertices is Gaussian, with a standard deviation of 7.1\(\,\text {cm}\). The simulated data were reweighted so as to have the same vertex \(z\) coordinate distribution as the data.

The hadron spectra were corrected for particles of non-primary origin (\(\tau > 10^{-12}\,\text {s} \)). The main sources of secondary particles are weakly decaying particles, mostly \(\mathrm {K^0_S}\), \(\mathrm {\Lambda }\)/\(\mathrm {\overline{\Lambda }}\), and \(\mathrm {\Sigma ^+}\)/\(\mathrm {\overline{\Sigma }^-}\). While the correction (\(C_\mathrm{s}\)) is around 1 % for pions, it rises up to 15 % for protons with \(p_{\mathrm {T}} \approx 0.2{\,\text {GeV/}c} \). As none of the mentioned weakly decaying particles decay into kaons, the correction for kaons is small. Based on studies comparing reconstructed \(\mathrm {K^0_S}\), \(\mathrm {\Lambda }\), and \(\mathrm {\overline{\Lambda }}\) spectra and predictions from the Hijing event generator, the corrections are reweighted by a \(p_{\mathrm {T}}\)-dependent factor.

For \(p < 0.15{\,\text {GeV/}c} \), electrons can be clearly identified. The overall \(\mathrm {e}^\pm \) contamination of the hadron yields is below 0.2 %. Although muons cannot be separated from pions, their fraction is very small, below 0.05 %. Since both contaminations are negligible, no corrections are applied for them.

3 Estimation of energy loss rate and yield extraction

In this paper an analytical parametrization [21] has been used to approximate the energy loss of charged particles in the silicon detectors. The method provides the probability density \(P(\Delta |\varepsilon , l)\) of energy deposit \(\Delta \), if the most probable energy loss rate \(\varepsilon \) at a reference path-length \(l_0 = 450\,\upmu \text {m} \) and the path-length \(l\) are known. It was used in conjunction with a maximum likelihood method, for the estimate of \(\varepsilon \).

For pixel clusters, the energy deposits were calculated as the sum of individual pixel deposits. In the case of strips, the energy deposits were corrected for capacitive coupling and cross-talk between neighboring strips. The readout threshold, the coupling parameter, and the standard deviation of the Gaussian noise for strips were determined from data, using tracks with close-to-normal incidence.

For an accurate determination of \(\varepsilon \), the response of all readout chips was calibrated with multiplicative gain correction factors. The measured energy deposit spectra were compared to the energy loss parametrization and hit-level corrections (affine transformation of energy deposits using scale factors and shifts) were introduced. The corrections were applied to individual hits during the determination of the \(\ln \varepsilon \) fit templates (described below).

The best value of \(\varepsilon \) for each track was calculated with the corrected energy deposits by minimizing the joint energy deposit negative log-likelihood of all hits on the trajectory (index \(i\)), \(\chi ^2 = -2 \sum _i \ln P(\Delta _i|\varepsilon ,l_i)\). Hits with incompatible energy deposits (contributing more than 12 to the joint \(\chi ^2\)) were excluded. At most one hit was removed; this affected about 1.5 % of the tracks.

Distributions of \(\ln \varepsilon \) as a function of total momentum \(p\) for positive particles are plotted in the top panel of Fig. 2 and compared to the predictions of the energy loss method [21] for electrons, pions, kaons, and protons. The remaining deviations were taken into account by means of track-level corrections mentioned above (affine transformation of templates using scale factors and shifts, \(\ln \varepsilon \rightarrow \alpha \ln \varepsilon + \delta \)).

Fig. 2
figure 2

Top distribution of \(\ln \varepsilon \) as a function of total momentum \(p\), for positively charged particles (\(\varepsilon \) is the most probable energy loss rate at a reference path length \(l_0 = 450\,\upmu \text {m} \)). The \(z\) scale is shown in arbitrary units and is linear. The curves show the expected \(\ln \varepsilon \) for electrons, pions, kaons, and protons (Eq. (30.11) in Ref. [22]). Bottom example \(\ln \varepsilon \) distribution at \(\eta = 0.35\) and \(p_{\mathrm {T}} = 0.775{\,\text {GeV/}c} \), with bin widths \(\Delta \eta = 0.1\) and \(\Delta p_{\mathrm {T}} = 0.05 {\,\text {GeV/}c} \). Scale factors (\(\alpha \)) and shifts (\(\delta \)) are indicated (see text). The inset shows the distribution with logarithmic vertical scale

Low-momentum particles can be identified unambiguously and can therefore be counted. Conversely, at high momentum, the \(\ln \varepsilon \) bands overlap (above about 0.5\({\,\text {GeV/}c}\) for pions and kaons and 1.2\({\,\text {GeV/}c}\) for protons); the particle yields therefore need to be determined by means of a series of template fits in \(\ln \varepsilon \), in bins of \(\eta \) and \(p_{\mathrm {T}}\) (Fig. 2, bottom panel). Finally, fit templates, giving the expected \(\ln \varepsilon \) distributions for all particle species (electrons, pions, kaons, and protons), were built from tracks. All kinematical parameters and hit-related observables were kept, but the energy deposits were regenerated by sampling from the analytical parametrization. For a less biased determination of track-level residual corrections, enhanced samples of each particle type were employed. These were used for setting starting values of the fits. For electrons and positrons, photon conversions in the beam-pipe and innermost first pixel layer were used. For high-purity \({\pi }\) and enhanced \(\mathrm {p}\) samples, weakly decaying hadrons were selected (\(\mathrm {K^0_S}\), \(\mathrm {\Lambda }\)/\(\mathrm {\overline{\Lambda }}\)). The relations and constraints described in Ref. [8] were also exploited, this way better constraining the parameters of the fits: fitting the \(\ln \varepsilon \) distributions in number of hits (\({n_\mathrm{hits}} \)) and track-fit \(\chi ^2/\hbox {ndf}\) slices simultaneously; fixing the distribution \({n_\mathrm{hits}} \) of particle species, relative to each other; using the expected continuity for refinement of track-level residual corrections, in neighboring \((\eta ,p_{\mathrm {T}})\) bins; using the expected convergence for track-level residual corrections, as the \(\ln \varepsilon \) values of two particle species approach each other.

The results of the (iterative) \(\ln \varepsilon \) fits are the yields for each particle species and charge in bins of \((\eta ,p_{\mathrm {T}})\) or \((y,p_{\mathrm {T}})\), both inclusive and divided into classes of reconstructed primary charged-track multiplicity. In the end, the histogram fit \(\chi ^2/\hbox {ndf}\) values were usually close to unity. Although pion and kaon yields could not be determined for \(p > 1.30{\,\text {GeV/}c} \), their sum was measured. This information is an important constraint when fitting the \(p_{\mathrm {T}}\) spectra.

The statistical uncertainties for the extracted yields are given by the fits. The observed local variations of parameters in the \((\eta ,p_{\mathrm {T}})\) plane for track-level corrections cannot be attributed to statistical fluctuations and indicate that the average systematic uncertainties in the scale factors and shifts are about \(10^{-2}\) and \(2 \times 10^{-3}\), respectively. These scale factors and shifts agree with those seen in the high-purity samples to well within a factor of two. The systematic uncertainties in the yields in each bin were obtained by refitting the histograms with the parameters changed by these amounts.

4 Corrections and systematic uncertainties

The measured yields in each \((\eta ,p_{\mathrm {T}})\) bin, \(\Delta N_\mathrm{measured}\), were first corrected for the misreconstructed-track rate (\(C_\mathrm{f}\)) and the fraction of secondary particles (\(C_\mathrm{s}\)):

$$\begin{aligned} \Delta N' = \Delta N_\mathrm{measured} \cdot (1 - C_\mathrm{f}) \cdot (1 - C_\mathrm{s}). \end{aligned}$$
(1)

The distributions were then unfolded to take into account the finite \(\eta \) and \(p_{\mathrm {T}}\) resolutions. The \(\eta \) distribution of the tracks is almost flat and the \(\eta \) resolution is very good. Conversely, the \(p_{\mathrm {T}}\) distribution is steep in the low-momentum region and separate corrections in each \(\eta \) bin were necessary. An unfolding procedure with linear regularization [23] was used, based on response matrices obtained from MC samples for each particle species.

The corrected yields were obtained by applying corrections for acceptance (\(C_\mathrm{a}\)), efficiency (\(C_\mathrm{e}\)), and multiple track reconstruction rate (\(C_\mathrm{m}\)):

$$\begin{aligned} \frac{1}{N_\mathrm{ev}} \frac{\text {d}^2 N}{\text {d}\eta \, \text {d}p_{\mathrm {T}}}_\mathrm{corrected} = \frac{1}{C_\mathrm{a} \cdot C_\mathrm{e} \cdot (1 + C_\mathrm{m})} \frac{\Delta N'}{N_\mathrm{ev} \Delta \eta \Delta p_{\mathrm {T}}}, \end{aligned}$$
(2)

where \(N_\mathrm{ev}\) is the corrected number of DS events (Fig. 1). Bins with acceptance smaller than 50 %, efficiency smaller than 50 %, multiple-track rate greater than 10 %, or containing less than 80 tracks were not used.

Finally, the differential yields \(\text {d}^2N/\text {d}\eta \,\text {d}p_{\mathrm {T}} \) were transformed to invariant yields \(\text {d}^2N/\text {d}y\,\text {d}p_{\mathrm {T}} \) by multiplying with the Jacobian \(E/p\) and the \((\eta ,p_{\mathrm {T}})\) bins were mapped into a \((y,p_{\mathrm {T}})\) grid. As expected, there is a small (5–10 %) \(y\) dependence in the narrow region considered (\(|y |<1\)), depending on event multiplicity. The yields as a function of \(p_{\mathrm {T}} \) were obtained by averaging over rapidity.

The systematic uncertainties are very similar to those in Ref. [8] and are summarized in Table 1. The uncertainties of the corrections related to the event selection and pileup are fully or mostly correlated and were treated as normalization uncertainties: 3.0 % uncertainty on the yields and 1.0 % on the average \(p_{\mathrm {T}}\). In order to study the influence of the high \(p_{\mathrm {T}}\) extrapolation on \(\langle \text {d}N/\text {d}y\rangle \) and \(\langle p_{\mathrm {T}} \rangle \), the \(1/n\) parameter of the fitted Tsallis–Pareto function (Sect. 5) was varied. While keeping the function in the measured range, \(1/n\) was increased and decreased by \(\pm 0.1\) above the highest \(p_{\mathrm {T}}\) measured point, ensuring that the two function pieces are continuous both in value and derivative. The choice of the magnitude for the variation was motivated by the fitted \(1/n\) values and their distance from a Boltzmann distribution. (The resulting functions are plotted in Fig. 3 as dotted lines.) The high \(p_{\mathrm {T}}\) extrapolation introduces sizeable systematic uncertainties, 4–6 % for \(\langle \text {d}N/\text {d}y\rangle \), and 9–15 % for \(\langle p_{\mathrm {T}} \rangle \) in case of the DS selection.

Fig. 3
figure 3

Transverse momentum distributions of identified charged hadrons (pions, kaons, protons, sum of pions and kaons) in the range \(|y |<1\), for positively (top) and negatively (bottom) charged particles. Kaon and proton distributions are scaled as shown in the legends. Fits to Eqs. (3) and (5) are superimposed. Error bars indicate the uncorrelated statistical uncertainties, while boxes show the uncorrelated systematic uncertainties. The fully correlated normalization uncertainty (not shown) is 3.0 %. Dotted lines illustrate the effect of varying the \(1/n\) value of the Tsallis–Pareto function by \(\pm 0.1\) above the highest measured \(p_{\mathrm {T}}\) point

Table 1 Summary of the systematic uncertainties affecting the \(p_{\mathrm {T}} \) spectra. Values in parentheses indicate uncertainties in the \(\langle p_{\mathrm {T}} \rangle \) measurement. The systematic uncertainty related to the low \(p_{\mathrm {T}}\) extrapolation is small compared to the contributions from other sources and therefore not included in the combined systematic uncertainty of the measurement. Representative, particle-specific uncertainties (\({\pi }\), \(\mathrm {K}\), \(\mathrm {p}\)) are given for \(p_{\mathrm {T}} =0.6{\,\text {GeV/}c} \) in the third group of systematic uncertainties

The tracker acceptance and the track reconstruction efficiency generally have small uncertainties (1 and 3 %, respectively), but change rapidly at very low \(p_{\mathrm {T}}\) (bottom panel of Fig. 1), leading to a 6 % uncertainty on the yields in that range. For the multiple-track and misreconstructed-track rate corrections, the uncertainty is assumed to be 50 % of the correction, while for the case of the correction for secondary particles it was estimated to be 20 %. These mostly uncorrelated uncertainties are due to the imperfect modeling of the detector: regions with mismodeled efficiency in the tracker, alignment uncertainties, and channel-by-channel varying hit efficiency. These circumstances can change frequently in momentum space, so can be treated as uncorrelated.

The systematic uncertainties originating from the unfolding procedure were studied. Since the \(p_{\mathrm {T}} \) response matrices are close to diagonal, the unfolding of \(p_{\mathrm {T}} \) distributions did not introduce substantial systematics. At the same time the inherited uncertainties were properly propagated. The introduced correlations between neighboring \(p_{\mathrm {T}} \) bins were neglected, hence statistical uncertainties were regarded as uncorrelated while systematic uncertainties were expected to be locally correlated in \(p_{\mathrm {T}} \). The systematic uncertainty of the fitted yields is in the range 1–10 % depending mostly on total momentum.

5 Results

In previously published measurements of unidentified and identified particle spectra [16, 24], the following form of the Tsallis–Pareto-type distribution [25, 26] was fitted to the data:

$$\begin{aligned} \frac{\text {d}^2 N}{\text {d}y \, \text {d}p_{\mathrm {T}}} = \frac{\text {d}N}{\text {d}y} \cdot C \cdot p_{\mathrm {T}} \left[ 1 + \frac{m_{\mathrm {T}}- m}{nT} \right] ^{-n}, \end{aligned}$$
(3)

where

$$\begin{aligned} C = \frac{(n-1)(n-2)}{nT[nT + (n-2) m]} \end{aligned}$$
(4)

and \(m_{\mathrm {T}} = \sqrt{m^2 + p_{\mathrm {T}} ^2}\) (factors of \(c\) are omitted from the preceding formulae). The free parameters are the integrated yield \(\text {d}N/\text {d}y\), the exponent \(n\), and parameter \(T\). The above formula is useful for extrapolating the spectra to zero \(p_{\mathrm {T}}\) and very high \(p_{\mathrm {T}} \) and for extracting \(\langle p_{\mathrm {T}} \rangle \) and \(\text {d}N/\text {d}y\). Its validity for different multiplicity bins was cross-checked by fitting MC spectra in the \(p_{\mathrm {T}} \) ranges where there are data points, and verifying that the fitted values of \(\langle p_{\mathrm {T}} \rangle \) and \(\text {d}N/\text {d}y\) were consistent with the generated values. Nevertheless, for a more robust estimation of both quantities (\(\langle p_{\mathrm {T}} \rangle \) and \(\langle \text {d}N/\text {d}y \rangle \)), the data points and their uncertainties were used in the measured range and the fitted functions only for the extrapolation in the unmeasured regions. According to some models of particle production based on non-extensive thermodynamics [26], the parameter \(T\) is connected with the average particle energy, while \(n\) characterizes the “non-extensivity” of the process, i.e. the departure of the spectra from a Boltzmann distribution (\(n = \infty \)).

As discussed earlier, pions and kaons cannot be unambiguously distinguished at higher momenta. Because of this, the pion-only, the kaon-only, and the joint pion and kaon \(\text {d}^2N/\text {d}y \, \text {d}p_{\mathrm {T}} \) distributions were fitted for \(|y | < 1\) and \(p < 1.20{\,\text {GeV/}c} \), \(|y | < 1\) and \(p < 1.05{\,\text {GeV/}c} \), and \(|\eta | <1\) and \(1.05 < p < 1.7{\,\text {GeV/}c} \), respectively. Since the ratio \(p/E\) for the pions (which are more abundant than kaons) at these momenta can be approximated by \(p_{\mathrm {T}}/m_{\mathrm {T}} \) at \(\eta \approx 0\), Eq. (3) becomes:

$$\begin{aligned} \frac{\text {d}^2 N}{\text {d}\eta \, \text {d}p_{\mathrm {T}}} \approx \frac{\text {d}N}{\text {d}y} \cdot C \cdot \frac{p_{\mathrm {T}} ^2}{m_{\mathrm {T}}} \left( 1 + \frac{m_{\mathrm {T}}-m}{nT} \right) ^{-n}. \end{aligned}$$
(5)

The approximate fractions of particles outside the measured \(p_{\mathrm {T}}\) range depend on track multiplicity; they are 15–30 % for pions, 40–50 % for kaons, and 20–35 % for protons. The average transverse momentum \(\langle p_{\mathrm {T}} \rangle \) and its uncertainty were obtained using data points in the measured range complemented by numerical integration of Eq. (3) with the fitted parameters in the unmeasured regions, under the assumption that the particle yield distributions follow the Tsallis–Pareto function in the low-\(p_{\mathrm {T}} \) and high-\(p_{\mathrm {T}} \) regions.

The results discussed in the following are for laboratory rapidity \(|y | < 1\). In all cases, error bars indicate the uncorrelated statistical uncertainties, while boxes show the uncorrelated systematic uncertainties. The fully correlated normalization uncertainty is not shown. For the \(p_{\mathrm {T}}\) spectra, the average transverse momentum, and the ratio of particle yields, the data are compared to Ampt 1.26/2.26 [13], Epos Lhc [14, 15], and Hijing 2.1 [11] MC event generators. Numerical results corresponding to the plotted spectra, fit results, as well as their statistical and systematic uncertainties are given in Ref. [27].

5.1 Inclusive measurements

The transverse momentum distributions of positively and negatively charged hadrons (pions, kaons, protons) are shown in Fig. 3, along with the results of the fits to the Tsallis–Pareto parametrization (Eqs. (3) and (5)). The fits are of good quality with \(\chi ^2/\hbox {ndf}\) values in the range 0.4–2.8 (Table 2). Figure 4 presents the data compared to the Ampt, Epos Lhc, and Hijing predictions. Epos Lhc gives a good description, while other generators predict steeper \(p_{\mathrm {T}}\) distributions than found in data.

Fig. 4
figure 4

Transverse momentum distributions of identified charged hadrons (pions, kaons, protons) in the range \(|y |<1\), for positively (top) and negatively (bottom) charged particles. Measured values (same as in Fig. 3) are plotted together with predictions from Ampt, Epos Lhc, and Hijing. Error bars indicate the uncorrelated statistical uncertainties, while boxes show the uncorrelated systematic uncertainties. The fully correlated normalization uncertainty (not shown) is 3.0 %

Table 2 Fit results (\(\text {d}N/\text {d}y\), \(1/n\), and \(T\)) and goodness-of-fit values for the DS selection shown together with calculated averages (\(\langle \text {d}N/\text {d}y\rangle \), \(\langle p_{\mathrm {T}} \rangle \)) for charged pions, kaons, and protons. The systematic uncertainty related to the low \(p_{\mathrm {T}}\) extrapolation is small compared to the contributions from other sources and therefore not included in the combined systematic uncertainty of the measurement. Combined uncertainties are given

Ratios of particle yields as a function of the transverse momentum are plotted in Fig. 5. While the \(\mathrm {K}/{\pi }\) ratios are well described by the Ampt simulation, only Epos Lhc is able to predict both \(\mathrm {K}/{\pi }\) and \(\mathrm {p}/{\pi }\) ratios. The ratios of the yields for oppositely charged particles are close to one, as expected for LHC energies at midrapidity.

Fig. 5
figure 5

Ratios of particle yields as a function of transverse momentum. \(\mathrm {K}\)/\({\pi }\) and \(\mathrm {p}\)/\({\pi }\) values are shown in the top panel, and opposite-charge ratios are plotted in the bottom panel. Error bars indicate the uncorrelated statistical uncertainties, while boxes show the uncorrelated systematic uncertainties. In the top panel, curves indicate predictions from Ampt, Epos Lhc, and Hijing

5.2 Multiplicity dependent measurements

A study of the dependence on track multiplicity is motivated partly by the intriguing hadron correlations measured in pp and pPb collisions at high track multiplicities [2831], suggesting possible collective effects in “central” pp and pPb collisions at the LHC. At the same time, it was seen that in pp collisions the characteristics of particle production (\(\langle p_{\mathrm {T}} \rangle \), ratios) at LHC energies are strongly correlated with event particle multiplicity rather than with the center-of-mass energy of the collision [8]. The strong dependence on multiplicity (or centrality) was also seen in dAu collisions at RHIC [6, 7]. In addition, the multiplicity dependence of particle yield ratios is sensitive to various final-state effects (hadronization, color reconnection, collective flow) implemented in MC models used in collider and cosmic-ray physics [32].

The event multiplicity \(N_\mathrm{rec}\) is obtained from the number of reconstructed tracks with \(|\eta |<2.4\), where the tracks are reconstructed using the same algorithm as for the identified charged hadrons [18]. (The multiplicity variable \(N_\mathrm{trk}^\mathrm{offline}\), used in Ref. [29], is obtained from a different track reconstruction configuration and a value of \(N_\mathrm{trk}^\mathrm{offline} = 110\) corresponds roughly to \(N_\mathrm{rec} = 170\).) The event multiplicity was divided into 19 classes, defined in Table 3. To facilitate comparisons with models, the corresponding corrected charged particle multiplicity in the same acceptance of \(|\eta |<2.4\) (\(N_\mathrm{tracks}\)) is also determined. For each multiplicity class, the correction from \(N_\mathrm{rec}\) to \(N_\mathrm{tracks}\) uses the efficiency estimated with the Hijing simulation in \((\eta ,p_{\mathrm {T}})\) bins. The corrected data are then integrated over \(p_{\mathrm {T}}\), down to zero yield at \(p_{\mathrm {T}} = 0\) (with a linear extrapolation below \(p_{\mathrm {T}} = 0.1 {\,\text {GeV/}c} \)). Finally, the integrals for each eta slice are summed. The average corrected charged-particle multiplicity \(\langle N_\mathrm{tracks} \rangle \), and also its values with the condition \(p_{\mathrm {T}} > 0.4{\,\text {GeV/}c} \), are shown in Table 3 for each event multiplicity class. The value of \(\langle N_\mathrm{tracks} \rangle \) is used to identify the multiplicity class in Figs. 6, 7, 8, 9, and 10.

Fig. 6
figure 6

Transverse momentum distributions of charged pions, kaons, and protons, normalized such that the fit integral is unity, in every second multiplicity class (\(\langle N_\mathrm{tracks} \rangle \) values are indicated) in the range \(|y |<1\), fitted with the Tsallis–Pareto parametrization (solid lines). For better visibility, the result for any given \(\langle N_\mathrm{tracks} \rangle \) bin is shifted by 0.3 units with respect to the adjacent bins. Error bars indicate the uncorrelated statistical uncertainties, while boxes show the uncorrelated systematic uncertainties. Dotted lines illustrate the effect of varying the \(1/n\) value of the Tsallis–Pareto function by \(\pm 0.1\) above the highest measured \(p_{\mathrm {T}}\) point

Fig. 7
figure 7

Ratios of particle yields in the range \(|y |<1\) as a function of the corrected track multiplicity for \(|\eta |<2.4\). \(\mathrm {K}\)/\({\pi }\) and \(\mathrm {p}\)/\({\pi }\) values are shown in the top panel, and opposite-charge ratios are plotted in the bottom panel. Error bars indicate the uncorrelated combined uncertainties, while boxes show the uncorrelated systematic uncertainties. In the top panel, curves indicate predictions from Ampt, Epos Lhc, and Hijing

Fig. 8
figure 8

Average transverse momentum of identified charged hadrons (pions, kaons, protons) in the range \(|y |<1\), as a function of the corrected track multiplicity for \(|\eta |<2.4\), computed assuming a Tsallis–Pareto distribution in the unmeasured range. Error bars indicate the uncorrelated combined uncertainties, while boxes show the uncorrelated systematic uncertainties. The fully correlated normalization uncertainty (not shown) is 1.0 %. Curves indicate predictions from Ampt, Epos Lhc, and Hijing

Fig. 9
figure 9

Average transverse momentum of identified charged hadrons (pions, kaons, protons; top panel) and ratios of particle yields (bottom panel) in the range \(|y |<1\) as a function of the corrected track multiplicity for \(|\eta |<2.4\), for pp collisions (open symbols) at several energies [8], and for pPb collisions (filled symbols) at \(\sqrt{s_{NN}} =\) 5.02\(\,\text {TeV}\). Both \(\langle p_{\mathrm {T}} \rangle \) and yield ratios were computed assuming a Tsallis–Pareto distribution in the unmeasured range. Error bars indicate the uncorrelated combined uncertainties, while boxes show the uncorrelated systematic uncertainties. For \(\langle p_{\mathrm {T}} \rangle \) the fully correlated normalization uncertainty (not shown) is 1.0 %. In both plots, lines are drawn to guide the eye (gray solid pp 0.9\(\,\text {TeV}\), gray dotted pp 2.76\(\,\text {TeV}\), black dash-dotted pp 7\(\,\text {TeV}\), colored solid pPb 5.02\(\,\text {TeV}\)). The ranges of \(\langle p_{\mathrm {T}} \rangle \), \(\mathrm {K}/{\pi }\) and \(\mathrm {p}/{\pi }\) values measured by ALICE in various centrality PbPb collisions (see text) at \(\sqrt{s_{NN}} = 2.76\,\text {TeV} \) [33] are indicated with horizontal bands

Table 3 Relationship between the number of reconstructed tracks (\(N_\mathrm{rec}\)) and the average number of corrected tracks (\(\langle N_\mathrm{tracks} \rangle \)) in the region \(|\eta | < 2.4\), and also with the condition \(p_{\mathrm {T}} > 0.4{\,\text {GeV/}c} \) (used in Ref. [29]), in the 19 multiplicity classes considered

Transverse-momentum distributions of identified charged hadrons, normalized such that the fit integral is unity, in selected multiplicity classes for \(|y | < 1\) are shown in Fig. 6 for pions, kaons, and protons. The distributions of negatively and positively charged particles have been summed. The distributions are fitted with the Tsallis–Pareto parametrization with \(\chi ^2/\hbox {ndf}\) values in the range 0.8–4.0 for pions, 0.1–1.1 for kaons, and 0.1–0.7 for protons. For kaons and protons, the parameter \(T\) increases with multiplicity, while for pions \(T\) slightly increases and the exponent \(n\) slightly decreases with multiplicity (not shown).

The ratios of particle yields are displayed as a function of track multiplicity in Fig. 7. The \(\mathrm {K}/{\pi }\) and \(\mathrm {p}/{\pi }\) ratios are flat, or slightly rising, as a function of \(\langle N_\mathrm{tracks} \rangle \). While none of the models is able to precisely reproduce the track multiplicity dependence, the best and worst matches to the overall scale are given by Epos Lhc and Hijing, respectively. The ratios of yields of oppositely charged particles are independent of \(\langle N_\mathrm{tracks} \rangle \) as shown in the bottom panel of Fig. 7. The average transverse momentum \(\langle p_{\mathrm {T}} \rangle \) is shown as a function of multiplicity in Fig. 8. As expected from the discrepancies between theory and data shown in Fig. 4, Epos Lhc again gives a reasonable description, while the other event generators presented here underpredict the measured values. For the dependence of \(T\) on multiplicity (not shown), the predictions match the pion data well; the kaon and proton values are much higher than in Ampt or Hijing.

5.3 Comparisons to pp and PbPb data

The comparison with pp data taken at various center-of-mass energies (0.9, 2.76, and 7\(\,\text {TeV}\)) [8] is shown in Fig. 9, where the dependence of \(\langle p_{\mathrm {T}} \rangle \) and the particle yield ratios (\(\mathrm {K}/{\pi }\) and \(\mathrm {p}/{\pi }\)) on the track multiplicity is shown. The plots also display the ranges of these values measured by ALICE in PbPb collisions at \(\sqrt{s_{NN}} =\) 2.76\(\,\text {TeV}\) for centralities from peripheral (80–90 % of the inelastic cross-section) to central (0–5 %) [33]. These ALICE PbPb data cover a much wider range of \(N_\mathrm{tracks}\) than is shown in the plot. Although PbPb data are not available at \(\sqrt{s_{NN}} = 5.02\,\text {TeV} \) for comparison, the evolution of event characteristics from RHIC (\(\sqrt{s_{NN}} = 0.2\,\text {TeV} \), [3, 4, 6]) to LHC energies [33] suggests that yield ratios should remain similar, while \(\langle p_{\mathrm {T}} \rangle \) values will increase by about 5 % when going from \(\sqrt{s_{NN}} =\) 2.76\(\,\text {TeV}\) to 5.02\(\,\text {TeV}\).

For low track multiplicity (\(N_\mathrm{tracks} \lesssim 40\)), pPb collisions behave very similarly to pp collisions, while at higher multiplicities (\(N_\mathrm{tracks} \gtrsim 50\)) the \(\langle p_{\mathrm {T}} \rangle \) is lower for pPb than in pp. The first observation can be explained since low-multiplicity events are peripheral pPb collisions in which only a few proton–nucleon collisions are present. Events with more particles are indicative of collisions in which the projectile proton strikes the thick disk of the lead nucleus. Interestingly, the pPb curves (Fig. 9, top panel) can be reasonably approximated by taking the pp values and multiplying their \(N_\mathrm{tracks}\) coordinate by a factor of 1.8, for all particle types. In other words, a pPb collision with a given \(N_\mathrm{tracks}\) is similar to a pp collision with \(0.55 \times N_\mathrm{tracks}\) for produced charged particles in the \(|\eta | < 2.4\) range. Both the highest-multiplicity pp and pPb interactions yield higher \(\langle p_{\mathrm {T}} \rangle \) than seen in central PbPb collisions. While in the PbPb case even the most central collisions possibly contain a mix of soft (lower-\(\langle p_{\mathrm {T}} \rangle \)) and hard (higher-\(\langle p_{\mathrm {T}} \rangle \)) nucleon-nucleon interactions, for pp or pPb collisions the most violent interaction or sequence of interactions are selected.

The transverse momentum spectra could also be successfully fitted (\(\chi ^2/\hbox {ndf}\) in the range 0.7–1.8) with a functional form proportional to \(p_{\mathrm {T}} \exp (-m_{\mathrm {T}}/T')\), where \(T'\) is called the inverse slope parameter, motivated by the success of Boltzmann-type distributions in nucleus–nucleus collisions [34]. In the case of pions, the fitted range was restricted to \(m_{\mathrm {T}} > 0.4{\,\text {GeV/}c} \) in order to exclude the region where resonance decays would significantly contribute to the measured spectra. The inverse slope parameter as a function of hadron mass is shown in Fig. 10, for a selection of event classes, both for pPb data and for MC event generators (Ampt, Epos Lhc, and Hijing). While the data display a linear dependence on mass with a slope that increases with particle multiplicity, the models predict a flat or slowly rising behavior versus mass and only limited changes with track multiplicity. This is to be compared with pp results [8], where both data and several tunes of the pythia 6 [35] and pythia 8 event generators show features very similar to those in pPb data. A similar trend is also observed in nucleus–nucleus collisions [3, 6], which is attributed to the effect of radial flow velocity boost [1].

Fig. 10
figure 10

Inverse slope parameters \(T'\) from fits of pion, kaon, and proton spectra (both charges) with a form proportional to \(p_{\mathrm {T}} \exp (-m_{\mathrm {T}}/T')\). Results for a selection of multiplicity classes, with different \(\langle N_\mathrm{tracks} \rangle \) as indicated, are plotted for pPb data (top) and for MC event generators Ampt, Epos Lhc, and Hijing (bottom). The curves are drawn to guide the eye

Average rapidity densities \(\langle \text {d}N/\text {d}y\rangle \) and average transverse momenta \(\langle p_{\mathrm {T}} \rangle \) of charge-averaged pions, kaons, and protons as a function of center-of-mass energy are shown in Fig. 11 for pp and pPb collisions, both corrected to the DS selection. To allow comparison at the pPb energy, a parabolic (linear) interpolation of the pp collision values at \(\sqrt{s}= 0.9\), 2.76, and 7\(\,\text {TeV}\) is shown for \(\text {d}N/\text {d}y\) \((\langle p_{\mathrm {T}} \rangle )\). The rapidity densities are generally about three times greater than in pp interactions at the same energy, while the average transverse momentum increases by about 20, 10, and 30 % for pions, kaons, and protons, respectively. The factor of three difference in the yields for pPb as compared to pp can be compared with the estimated number of projectile collisions \(N_\mathrm{coll}/2 = 3.5 \pm 0.3\) or with the number of nucleons participating in the collision \(N_\mathrm{part}/2 = 4.0 \pm 0.3\), based on the ratio of preliminary pPb and pp cross-section measurements, that have proven to be good scaling variables in proton–nucleus collisions at lower energies [36].

Fig. 11
figure 11

Average rapidity densities \(\langle \text {d}N/\text {d}y\rangle \) (top) and average transverse momenta \(\langle p_{\mathrm {T}} \rangle \) (bottom) as a function of center-of-mass energy for pp [8] and pPb collisions, for charge-averaged pions, kaons, and protons. Error bars indicate the uncorrelated combined uncertainties, while boxes show the uncorrelated systematic uncertainties. The curves show parabolic (for \(\langle \text {d}N/\text {d}y\rangle \)) or linear (for \(\langle p_{\mathrm {T}} \rangle \)) interpolation on a log-log scale. The pp and pPb data are for laboratory rapidity \(|y |<1\), which is the same as the center-of-mass rapidity only for the pp data

6 Conclusions

Measurements of identified charged hadron spectra produced in pPb collisions at \(\sqrt{s_{NN}} =5.02\,\text {TeV} \) have been presented, normalized to events with simultaneous hadronic activity at pseudorapidities \(-5 < \eta < -3\) and \(3 < \eta < 5\). Charged pions, kaons, and protons were identified from the energy deposited in the silicon tracker and other track information. In the present analysis, the yield and spectra of identified hadrons for laboratory rapidity \(|y | \!<\! 1\) have been studied as a function of the event charged particle multiplicity in the range \(|\eta |\!<\!2.4\). The \(p_{\mathrm {T}} \) spectra are well described by fits with the Tsallis–Pareto parametrization. The ratios of the yields of oppositely charged particles are close to one, as expected at mid-rapidity for collisions of this energy. The average \(p_{\mathrm {T}} \) is found to increase with particle mass and the event multiplicity. These results are valid under the assumption that the particle yield distributions follow the Tsallis–Pareto function in the unmeasured \(p_{\mathrm {T}} \) regions.

The results can be used to further constrain models of hadron production and contribute to the understanding of basic non-perturbative dynamics in hadron collisions. The Epos Lhc event generator reproduces several features of the measured distributions, a significant improvement from the previous version, attributed to a new viscous hydrodynamic treatment of the produced particles. Other studied generators (Ampt, Hijing) predict steeper \(p_{\mathrm {T}}\) distributions and much smaller \(\langle p_{\mathrm {T}} \rangle \) than found in data, as well as substantial deviations in the \(\mathrm {p}/{\pi }\) ratios.

Combined with similar results from pp collisions, the track multiplicity dependence of the average transverse momentum and particle ratios indicate that particle production at LHC energies is strongly correlated with event particle multiplicity in both pp and pPb interactions. For low track multiplicity, pPb collisions appear similar to pp collisions. At high multiplicities, the average \(p_{\mathrm {T}} \) of particles from pPb collisions with a charged particle multiplicity of \(N_\mathrm{tracks}\) (in \(|\eta |<2.4\)) is similar to that for pp collisions with \(0.55 \times N_\mathrm{tracks}\). Both the highest-multiplicity pp and pPb interactions yield higher \(\langle p_{\mathrm {T}} \rangle \) than seen in central PbPb collisions.