Measurement of charged pion, kaon, and proton production in proton-proton collisions at ﬃﬃ s p = 13 TeV

Transverse momentum spectra of charged pions, kaons, and protons are measured in proton-proton collisions at ﬃﬃﬃ s p ¼ 13 TeV with the CMS detector at the LHC. The particles, identified via their energy loss in the silicon tracker, are measured in the transverse momentum range of p T ≈ 0 . 1 – 1 . 7 GeV =c and rapidities j y j < 1 . The p T spectra and integrated yields are compared to previous results at smaller ﬃﬃﬃ s p and to predictions of Monte Carlo event generators. The average p T increases with particle mass and charged particle multiplicity of the event. Comparisons with previous CMS results at ﬃﬃﬃ s p ¼ 0 . 9 , 2.76, and 7 TeV show that the average p T and the ratios of hadron yields feature very similar dependences on the particle multiplicity in the event, independently of the center-of-mass energy of the pp collision.

Transverse momentum spectra of charged pions, kaons, and protons are measured in proton-proton collisions at ffiffi ffi s p ¼ 13 TeV with the CMS detector at the LHC. The particles, identified via their energy loss in the silicon tracker, are measured in the transverse momentum range of p T ≈ 0.1-1.7 GeV=c and rapidities jyj < 1. The p T spectra and integrated yields are compared to previous results at smaller ffiffi ffi s p and to predictions of Monte Carlo event generators. The average p T increases with particle mass and charged particle multiplicity of the event. Comparisons with previous CMS results at ffiffi ffi s p ¼ 0.9, 2.76, and 7 TeV show that the average p T and the ratios of hadron yields feature very similar dependences on the particle multiplicity in the event, independently of the center-of-mass energy of the pp collision. DOI: 10.1103/PhysRevD.96.112003

I. INTRODUCTION
The study of hadron production has a long history in high-energy particle, nuclear, and cosmic ray physics. The absolute yields and the transverse momentum (p T ) spectra of identified hadrons in high-energy hadron-hadron collisions are among the most basic physical observables. They can be used to improve the modeling of various key ingredients of Monte Carlo (MC) hadronic event generators, such as multiparton interactions, parton hadronization, and final-state effects (such as parton correlations in color, p T , spin, baryon and strangeness number, and collective flow) [1]. The dependence of the hadron spectra and yields on the impact parameter of the proton-proton (pp) collision provides additional valuable information to tune the corresponding MC parameters. Indeed, parton hadronization and final-state effects are mostly constrained from elementary e þ e − collisions, whose final states are largely dominated by simple qq final states, whereas lowp T hadrons at the LHC issue from the fragmentation of multiple gluon "minijets" [1]. Such large differences have a particularly important impact on baryons and strange hadrons, whose production in pp collisions is not well reproduced by the existing models [2,3], and also affect the modeling of hadronic interactions of ultrahigh-energy cosmic rays with Earth's atmosphere [4]. Spectra of identified particles in pp collisions also constitute an important reference for high-energy heavy ion studies, where various final-state effects are known to modify the spectral shape and yields of different hadron species [5][6][7][8][9].
The present analysis uses pp collisions collected by the CMS experiment at the CERN LHC at ffiffi ffi s p ¼ 13 TeV and focuses on the measurement of the p T spectra of charged hadrons, identified primarily via their energy depositions in the silicon detectors. The analysis adopts the same methods as used in previous CMS measurements of pion, kaon, and proton production in pp and pPb collisions at ffiffi ffi s p of 0.9, 2.76, and 7 TeV [2,10], as well as those performed by the ALICE Collaboration at 2.76 and 7 TeV [3,11].

II. THE CMS DETECTOR AND EVENT GENERATORS
A detailed description of the CMS detector can be found in Ref. [12]. The CMS experiment uses a right-handed coordinate system, with the origin at the nominal interaction point (IP) and the z axis along the counterclockwisebeam direction. The pseudorapidity η 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 η ¼ − ln½tanðθ=2Þ, where θ is the polar angle with respect to the z axis and y ¼ 1 2 ln½ðE þ p z Þ=ðE − p z Þ. The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter. Within the 3.8 T field volume are the silicon pixel and strip tracker, the crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter. The tracker measures charged particles within the range jηj < 2.4. It has 1440 silicon pixel and 15 148 silicon strip detector modules with thicknesses of either 300 or 500 μm, assembled in 13 detection layers in the central region. Beam pick-up timing for the experiment (BPTX) devices were used to trigger the detector readout. They are located around the beam pipe at a distance of 175 m from the IP on either side, and are designed to provide precise information on the bunch structure and timing of the incoming beams of the LHC.
In this paper, distributions of identified hadrons produced in inelastic pp collisions are compared to predictions from MC event generators based on two different theoretical frameworks: perturbative QCD (PYTHIA6.426 [13] and PYTHIA 8.208 [14]) and Reggeon field theory (EPOS v3400 [15]). On the one hand, the basic ingredients of PYTHIA 6 and PYTHIA8 are (multiple) leading-order perturbative QCD 2 → 2 matrix elements, complemented with initial-and final-state parton radiation (ISR and FSR), folded with parton distribution functions in the proton, and the Lund string model for parton hadronization. Two different "tunes" of the parameters governing the nonperturbative and semihard dynamics (ISR and FSR showering, multiple parton interactions, beam-remnants, final-state color-reconnection, and hadronization) are used: the PYTHIA6 Z2* [13,16] and PYTHIA8 CUETP8M1 [16] tunings, based on fits to recent minimum bias and underlying event measurements at the LHC. On the other hand, EPOS starts off from the basic quantum field-theory principles of unitarity and analyticity of scattering amplitudes as implemented in Gribov's Reggeon field theory [17], extended to include (multiple) parton scatterings via "cut (hard) Pomerons." The latter objects correspond to color flux tubes that are finally hadronized also via the Lund string model. The version of EPOS used here is run with the LHC tune [18] which includes collective final-state string interactions resulting in an extra radial flow of the final hadrons produced in more central pp collisions.

III. EVENT SELECTION AND RECONSTRUCTION
The data used for the measurements presented in this paper were taken during a special low luminosity run where the average number of pp interactions in each bunch crossing was 1.0. A total of 7 × 10 6 collisions were recorded, corresponding to an integrated luminosity of approximately 0.1 nb −1 .
The event selection consisted of the following requirements: (i) at trigger level, the coincidence of signals from both BPTX devices, indicating the presence of both proton bunches crossing the interaction point; (ii) offline, to have at least one reconstructed interaction vertex; (iii) beam-halo and beam-induced background events, which usually produce an anomalously large number of pixel hits, were identified [19] and rejected. The event selection efficiency as well as the tracking and vertexing acceptance and efficiency are evaluated using simulated event samples produced with the PYTHIA8 (tune CUETP8M1) MC event generator, followed by the CMS detector response simulation based on GEANT4 [20]. Simulated events are reconstructed and analyzed in the same way as collision data events. The final results are given for an event selection corresponding to inelastic pp collisions, which will be presented in Sec. VI. According to the three MC event generators considered, the fraction of inelastic pp collisions not resulting in a reconstructed pp interaction amounts to about 14% AE 3%, where the uncertainty is based on the variance of the predictions coming from the event generators. These events are mostly diffractive ones with negligible central activity.
The reconstruction of charged particles in CMS is limited by the acceptance of the tracker (jηj < 2.4) and by the decreasing tracking efficiency at low momentum caused by multiple scattering and energy loss. The identification of particle species using specific ionization (Sec. IV) is restricted to p < 0.15 GeV=c for electrons, p < 1.20 GeV=c for pions, p < 1.05 GeV=c for kaons, and p < 1.70 GeV=c for protons [2,10]. Pions are measured up to a higher momentum than kaons because of their larger relative abundance. In order to have a common kinematic region where pions, kaons, and protons can all be identified, the range jyj < 1 is chosen for this measurement.
The extrapolation of particle spectra into unmeasured ðy; p T Þ regions is model dependent, particularly at low p T . A precise measurement therefore requires reliable track reconstruction down to the lowest possible p T values. Special tracking algorithms [21], already used in previous studies [2,10,19,22], made it possible to extend the present analysis to p T ≈ 0.1 GeV=c with high reconstruction efficiency and low background. Compared to the standard tracking algorithm used in CMS, these algorithms feature special track seeding and cleaning, hit cluster shape filtering, modified trajectory propagation, and track quality requirements. The charged-pion mass is assumed when fitting particle momenta.
The acceptance of the tracker (C a ) is defined as the fraction of primary charged particles leaving at least two hits in the pixel detector. Based on MC studies, it is flat in the region jηj < 2 and p T > 0.4 GeV=c, and at values of 96%-98% as can be seen in Fig. 1. The loss of acceptance at p T < 0.4 GeV=c is caused by energy loss and multiple scattering, which are both functions of particle mass. The reconstruction efficiency (C e ), which is defined as the fraction of accepted charged particles that result in a successfully reconstructed trajectory, is usually in the range 80%-90%. It decreases at low p T , also in a mass-dependent way. The misreconstructed-track rate (C f ), defined as the fraction of reconstructed primary charged tracks without a corresponding genuine primary charged particle, is very small, reaching 1% for p T < 0.2 GeV=c. The probability of reconstructing multiple tracks (C m ) from a single charged particle is about 0.1%, mostly from particles spiralling in the strong magnetic field of the CMS solenoid. The efficiencies and background rates (misreconstruction, multiple reconstruction) are found not to depend on the charged-particle multiplicity of the event in the range of multiplicities of interest for this analysis. They largely factorize in η and p T , but for the final corrections (Sec. V) an ðη; p T Þ matrix is used.
The region where pp collisions occur (beam spot) is measured from the distribution of reconstructed interaction vertices. Since the bunches are very narrow in the plane transverse to the beam direction (with a width of about 50 μm for this special run), the x-y 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 that have p T > 0.1 GeV=c and originate from the vicinity of the beam spot, i.e. their transverse impact parameters d T (with respect to the center of the beam spot) satisfy the condition d T < 3σ T . Here σ T is the quadratic sum of the uncertainty in the value of d T and the root mean square of the beam spot distribution in the transverse plane. In order to reach higher efficiency in special-topology low-multiplicity events, an agglomerative vertex reconstruction algorithm [23] is used, with the z coordinates of the tracks (and their uncertainties) at the point of closest approach to the beam axis as input. The distance distributions of reconstructed vertex pairs in data indicates that the fraction of merged vertices (with tracks from two or more true vertices) and split vertices (two or more reconstructed vertices with tracks from a single true vertex) is about 1%. For single-vertex events, there is no minimum requirement on the number of tracks associated with the vertex (those assigned to it during vertex finding), and even one-track vertices, which are defined as the point of closest approach of the track to the beam line, are allowed. The fraction of events with more than one (three) reconstructed primary vertices is about 26% (1.8%). Only events with three or fewer reconstructed primary vertices were considered and only tracks associated with a primary vertex are used in the analysis.
The vertex resolution in the z direction is a strong function of the number of reconstructed tracks and is always less than 0.1 cm. The distribution of the z coordinates of the reconstructed primary vertices is Gaussian with a width of σ ¼ 4.2 cm. Simulated events are reweighted in order to have the same vertex z coordinate distribution as in collision data.
The contribution to the hadron spectra from particles of nonprimary origin arising from the decay of particles with proper lifetime τ > 10 −12 s was subtracted. The main sources of these secondary particles are weakly decaying particles, mostly K 0 S , Λ=Λ, and Σ þ =Σ − . According to the simulations, this correction (C s ) is approximately 1% for pions and rises to 15% for protons with p T ≈ 0.2 GeV=c. Because none of these particles decay weakly into kaons, the correction for kaons is less than 0.1%. Charged particles from interactions of primary particles or their decay products with detector material are suppressed by the impact parameter cuts described above. For p < 0.15 GeV=c, electrons can be clearly identified based on their energy loss (Fig. 2, left) and their contamination of the hadron yields is below 0.2%. Although muons cannot be distinguished from pions, according to MC predictions their fraction is below 0.05%. Since both contaminations are negligible with respect to the final uncertainties, no corrections are applied.

IV. ESTIMATION OF ENERGY LOSS RATE AND YIELD EXTRACTION
For this paper an analytical parametrization [25] is used to model the energy loss of charged particles in the silicon detectors. It provides the probability density PðΔjε; lÞ of finding an energy deposit Δ, if the most probable energy loss rate ε at a reference path length l 0 ¼ 450 μm and the path length l are known. The choice of 450 μm is motivated by being the approximate average path length traversed in the silicon detectors. The value of ε depends on the momentum and mass m of the charged particle. The parametrization is used in conjunction with a maximum likelihood fit for the estimate of ε. All details of the methods described below are given in Ref. [2].
Using the cluster shape filtering mentioned in Sec. III, only hit clusters compatible with the particle trajectory are used. For clusters in the pixel detector, the energy deposits are calculated based on the individual pixel deposits. In the case of clusters in the strip detector, the energy deposits are corrected for truncation performed by the readout electronics and for losses due to deposits below threshold because of capacitive coupling and cross-talk between neighboring strips. The readout threshold, the strength of coupling, and the standard deviation of the Gaussian noise for strips are determined from data. The response of all readout chips is calibrated with multiplicative gain correction factors.
After the readout chip calibration, the measured energy deposit spectra for each silicon subdetector are compared to the expectations of the energy loss model as a function of p=m and l using particles satisfying tight identification criteria. These comparisons allow the computation of hitlevel corrections to the energy loss model that is used to estimate the particle energy loss rate ε and its associated distribution.
The best value of ε for each track is calculated from the measured energy deposits by minimizing the negative loglikelihood function of the combined energy deposit for all hits (index i) associated with the particle trajectory, where the probability density functions include the hit-level corrections mentioned above. Hits with incompatible energy deposits (contributing more than 12 units to the combined χ 2 ) are excluded.
For the determination of ε, removal of at most one hit per track is allowed; this affected about 1.5% of the tracks.
Low-momentum particles can be identified unambiguously and can therefore be counted (Fig. 2). Conversely, at high momenta (above about 0.5 GeV=c for pions and kaons and above 1.2 GeV=c for protons) the ε bands overlap. Therefore the particle yields need to be determined by means of a series of template fits in ε, in bins of η and p T (Fig. 2, right panel). Fit templates with the expected ε distributions for all particle species (electrons, pions, kaons, and protons) are obtained from reconstructed tracks in data. All track parameters and hit-related quantities are kept but, in order to populate the distributions, the energy deposits are regenerated by sampling from the hit-level corrected analytical parametrization assuming a given particle type. Possible residual discrepancies between the observed and expected ε distributions, present in some regions of the parameter space (mostly at low p T ), are taken into account by means of the track-level corrections consisting, as for the hit-level corrections, of a linear transformation of the parametrization using scale factors and shifts. For a less biased determination of these track-level residual corrections, enriched samples of each particle type are employed for determining starting values of the parameters to be fitted. For electrons and positrons, photon conversions in the beam-pipe and in the innermost pixel layer are used. For high-purity pion and enriched proton samples, weakly decaying hadrons are selected (K 0 S , Λ=Λ). The following criteria and methods described in Ref.
[2] are also exploited to better constrain the parameters of the fits: fitting the ε distributions in slices of number of hits (n hits ) and track fit χ 2 =ndf (where ndf is number of degrees of freedom) simultaneously; setting constraints on the n hits distribution for specific particle species; imposing the expected continuity of track-level residual corrections in adjacent ðη; p T Þ bins; and using the expected convergence of track-level residual corrections as the ε values of two particle species approach each other at large momentum.
Distributions of ε as a function of total momentum p for positive particles are plotted in the left panel of Fig. 2 and compared to the predictions of the energy loss parametrization [25] for electrons, pions, kaons, and protons. The results of the (iterative) ε fits are the yields for each particle species and charge in bins of ðη; p T Þ or ðy; p T Þ, both inclusive and divided into classes of reconstructed primary charged-track multiplicity. Although pion and kaon yields could not be determined for p > 1.30 GeV=c, their sum is measured. This information is an important constraint when fitting the p T spectra.

V. YIELD EXTRACTION AND SYSTEMATIC UNCERTAINTIES
The measured yields in each ðη; p T Þ bin, ΔN measured , are first corrected for the misreconstructed-track rate C f and the fraction of secondary particles C s : The bin widths are Δη ¼ 0.1 and Δp T ¼ 0.05 GeV=c. The distributions are then unfolded to take into account bin migrations due to the finite η and p T resolutions. The η distribution of the tracks is almost flat and the η resolution is significantly smaller than the bin width. At the same time the p T distribution is steep in the low-momentum region and separate p T -dependent corrections in each η slice are necessary. For that, an unfolding procedure with a linear regularization method (Tikhonov regularization [26]) is used, based on response matrices obtained from PYTHIA 8 MC samples separately for each particle species. This procedure guarantees that the uncertainties associated with the assumption of the pion mass in the track fitting step are taken into account. The bin purities of the matrices are above 80%-90%. The chosen regularization term reflects that the original distribution changes only slowly, but that particular choice has negligible influence on the results.
Further corrections for acceptance, efficiency, and multiple track reconstruction probability are applied: where N ev is the corrected number of inelastic pp collisions in the data sample. Bins that meet at least one of the following criteria are not used in order to ensure robustness of the fits described below and to minimize the impact on the systematic uncertainties: acceptance less than 50%; efficiency less than 50%; multiple-track rate greater than 10%; multiplicity below 80 tracks. Finally, the η-differential yields d 2 N=dηdp T are transformed into d 2 N=dydp T yields by multiplying with the Jacobian of the η to y transformation (E=p), and the ðη; p T Þ bins are mapped onto a ðy; p T Þ grid. The differential yields exhibit a slight (5%-10%) dependence on y in the narrow region considered (jyj < 1), an effect that decreases with the event multiplicity. The yields as a function of p T are obtained averaged over the rapidity window.
The p T distributions are fit using a Tsallis-Pareto-type function, which empirically describes both the low-p T exponential and the high-p T power-law behaviors while employing only a few parameters. Based on the good reproduction of previous measurements of unidentified and identified particle spectra [2,10,19,27], the following form of the distribution [28,29] is used: where and m T ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The free parameters are the integrated yield dN=dy, the exponent n, and the parameter T. According to some models of particle production based on nonextensive thermodynamics [29], the parameter T is connected with the average particle energy, while n characterizes the "nonextensivity" of the process, i.e. the departure of the spectra from a Boltzmann distribution (n ¼ ∞). Equation (3) is useful for extrapolating the spectra down to zero and up to high p T , and thereby extracting hp T i and hdN=dyi. Its validity for different multiplicity bins is cross-checked by fitting MC spectra in the p T ranges where there are data points, and verifying that TABLE I. Summary of the systematic uncertainties affecting the p T spectra. Values in parentheses indicate uncertainties in the hp T i measurement. Representative, particle-specific uncertainties (π, K, p) are given for p T ¼ 0.6 GeV=c in the third group of systematic uncertainties.

Source
Uncertainty the fitted values of hp T i and hdN=dyi are consistent with the generated values. Nevertheless, for a more robust estimation of both hp T i and hdN=dyi, the unfolded binby-bin yield values and their uncertainties are used in the measured range while the fitted functions are employed for the extrapolation into the unmeasured regions. As discussed earlier, pions and kaons cannot be unambiguously distinguished at high momenta. For this reason the pion-only, the kaon-only, and the joint pion and kaon d 2 N=dydp T distributions are fitted for jyj < 1 and p < 1.20 GeV=c, jyj < 1 and p < 1.05 GeV=c, and jηj < 1 and 1.05 < p < 1.7 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 T =m T at η ≈ 0, Eq. (3) becomes Moreover, below p T values of 0.1-0.3 GeV the detector acceptance and the tracking efficiency significantly decrease. The Tsallis-Pareto function is used to extrapolate the measured yields both into this latter region and to the region at high momenta such that the integrated yield (dN=dy) and the average transverse momentum (hp T i) can be reported for the full p T range. This choice allows measurements performed by different experiments in various collision systems and center-of-mass energies to be compared.
The fractions of particles outside the measured p T range are 15%-30% for pions, 40%-50% for kaons, and 20%-35% for protons, depending on the track multiplicity of the event.
The systematic uncertainties are very similar to those in Ref.
[2] and are summarized in Table I. They are obtained from the comparison of different MC event generators, differences between data and simulation, or based on previous studies (hit inefficiency, misalignment). The uncertainties in the corrections C a , C e , C f and C m , which are related to the event selection, and the effects of pileup, are fully or mostly correlated and are treated as normalization uncertainties: altogether they propagate to a 3.0% uncertainty in the yields and a 1.0% uncertainty in the average p T . In order to study the influence of the high-p T extrapolation on the hdN=dyi and hp T i averages, the reciprocal of the exponent (1=n) of the fitted Tsallis-Pareto function was increased and decreased by AE0.05 only in the region above the highest measured p T ; in this same region both the function and its first derivative were required to fit continuously the data points. The choice of the magnitude for the variation is 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 (though they are mostly indistinguishable from the nominal fit curves). The high-p T extrapolation introduces systematic uncertainties of 1%-3% for hdN=dyi, and 4-8% for hp T i. The systematic uncertainty related to the low p T extrapolation is small compared to the contributions from other sources and therefore is not included in the combined systematic uncertainty of the measurement.  (5) are superimposed. For the π þ K fit, only the region corresponding to the range jηj < 1 and 1.05 < p < 1.7 GeV=c is plotted. Boxes show the uncorrelated systematic uncertainties, while error bars indicate the uncorrelated statistical uncertainties (barely visible). The fully correlated normalization uncertainty (not shown) is 3.0%. Dotted lines (mostly indistinguishable from the nominal fit curves) illustrate the effect of varying the inverse exponent (1=n) of the Tsallis-Pareto function by AE0.05 beyond the highest-p T measured point.
The tracker acceptance and the track reconstruction efficiency generally have small uncertainties (1% and 3%, respectively), but at very low p T they reach 6%. For the multiple-track and misreconstructed-track rate corrections, the uncertainty is assumed to be 50% of the correction, while for the correction for secondary particles it is estimated to be 25% based on the differences between predictions of MC event generators and data. These bin-bybin, largely uncorrelated uncertainties are caused by the imperfect modeling of the detector: regions with incorrectly modeled tracking efficiency, alignment uncertainties, and channel-by-channel varying hit efficiency. All these effects are taken as uncorrelated.
The statistical uncertainties in the extracted yields are given by the fit uncertainties. Variations of the track-level correction parameters, incompatible with statistical fluctuations, are observed. They are used to estimate the systematic uncertainties in the fitted scale factors and shifts and are at the level of 10 −2 and 2 × 10 −3 , respectively. The systematic uncertainties in the yields in each bin are thus obtained by refitting the histograms with the parameters changed by these amounts. For the present measurement, systematic uncertainties dominate over the statistical ones.
The systematic uncertainties originating from the unfolding procedure are also studied. Since the p T response matrices are close to diagonal, the unfolding of the p T distributions does not introduce substantial uncertainties. The correlations between neighboring p T bins are neglected, and therefore statistical uncertainties are regarded as uncorrelated. The systematic uncertainty of the fitted yields is in the range 1%-10%, depending primarily on total momentum.

VI. RESULTS
The results discussed in the following are averaged over the rapidity range jyj < 1. In all cases, error bars in the figures indicate the uncorrelated statistical uncertainties, while boxes show the uncorrelated systematic uncertainties. The fully correlated normalization uncertainty is not shown. For the p T spectra, the average transverse momentum hp T i, and the ratios of particle yields, the data are compared to the predictions of PYTHIA 8, EPOS, and PYTHIA 6.

A. 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 χ 2 =ndf values in the range 0.4-1.2 (Table II). Figure 4 presents the same data compared to the PYTHIA 8, EPOS, and PYTHIA 6 predictions. While pions are described well by all three generators, kaons are best modelled by PYTHIA 8 and EPOS. For protons and very low p T pions only PYTHIA 8 gives a good description of the data.
Ratios of particle yields as a function of the transverse momentum are plotted in Fig. 5. Only PYTHIA 8 is able to predict both the K=π and p=π ratios as a function of p T . The ratios of the yields for oppositely charged particles are close to one (Fig. 5, right), as expected at this center-ofmass energy in the central rapidity region.

B. Multiplicity-dependent measurements
The study of the p T spectra as a function of the event track multiplicity is motivated partly by the intriguing hadron correlations measured in pp and pPb collisions at high track multiplicities [30][31][32][33], suggesting possible collective effects in "central" collisions at the LHC. We have also observed that in pp collisions at LHC energies [2,10], the characteristics of particle production (hp T i, ratios of yields) are strongly correlated with the particle multiplicity in the event, which is in itself closely related to the number of underlying parton-parton interactions, independently of the concrete center-of-mass energy of the pp collision.
The event track multiplicity, N rec , is defined as the number of tracks with jηj < 2.4 reconstructed using the same algorithm as for the identified charged hadrons [21]. The event multiplicity is divided into 18 classes as defined in Table III. To facilitate comparisons with models, the event charged-particle multiplicity over jηj < 2.4 (N tracks ) is determined for each multiplicity class by correcting N rec for the track reconstruction efficiency, which is estimated with the PYTHIA 8 simulation in ðη; p T Þ bins. The corrected yields are then integrated over  (3) and (5)], associated goodness-of-fit values, and extracted hdN=dyi and hp T i averages, for charged pion, kaon, and proton spectra measured in the range jyj < 1 in inelastic pp collisions at 13 TeV. Combined statistical and systematic uncertainties are given. p T , down to zero yield at p T ¼ 0 (with a linear extrapolation below p T ¼ 0.1 GeV=c). Finally, the integrals for each η slice are summed up. The average corrected charged-particle multiplicity hN tracks i is shown in Table III for each event multiplicity class. The value of hN tracks i is used to identify the multiplicity class in Figs. 6-9.
Transverse-momentum distributions of pions, kaons, and protons, measured over jyj < 1 and normalized such that the fit integral is unity, are shown in Fig. 6 for various multiplicity classes. The distributions of negatively and positively charged particles are summed. The Tsallis-Pareto parametrization is fitted to the distributions with χ 2 =ndf values in the range 0.3-2.3 for pions, 0.2-2.6 for kaons, and 0.1-0.8 for protons. It is observed that for kaons and protons, the parameter T increases with multiplicity, while  III. Relationship between the number of reconstructed tracks (N rec ) and the average number of corrected tracks (hN tracks i) in the region jηj < 2.4 in the 18 multiplicity classes considered.  7  16  28  40  51  63  74  85  97  108  119  130  141  151  162  172  183  for pions T slightly increases and the exponent n slightly decreases with multiplicity. The ratios of particle yields are displayed as functions of track multiplicity in Fig. 7. The K=π and p=π ratios are relatively flat as a function of hN tracks i, and none of the models is able to accurately reproduce the track multiplicity dependence. The ratios of yields of oppositely charged particles are independent of hN tracks i as shown in the lower panel of Fig. 7. The average transverse momentum hp T i is shown as a function of multiplicity in Fig. 8. Although PYTHIA 8 gives a good description of the (multiplicity integrated) inelastic p T spectra (Fig. 4), none of the MC event generators reproduces well the multiplicity dependence of hp T i for all particle species. In particular, all generators overestimate the measured values for kaons. Pions are well described by PYTHIA 6 and EPOS, while protons are best described by PYTHIA 8.
In the lower multiplicity events, with fewer than 50 tracks, we observe a reasonable agreement between the data and the MC generator predictions for the different particle yields. However in higher multiplicity events, the measured kaon (proton) yield is smaller (higher) than predicted by the models. This indicates that the MC parameters that control the strangeness and baryon production as a function of parton multiplicity, need additional fine-tuning.

C. Comparisons with lower energy pp data
The comparison of these results with lower-energy pp data taken at various center-of-mass energies (0.9, 2.76, and 7 TeV) [2] is presented in Fig. 9, where the trackmultiplicity dependence of hp T i (left) and the particle yield ratios (K=π and p=π, right) are shown. In the previous publication [2], the final results are corrected to a particlelevel selection that requires at least one particle (with proper lifetime τ > 10 −18 s) with E > 3 GeV in the range   (pions, kaons, protons) in the range jyj < 1, as functions of the corrected track multiplicity for jηj < 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 PYTHIA 8, EPOS, and PYTHIA 6.
−5 < η < −3 and at least one in the range 3 < η < 5. This selection is referred to as the "double-sided" (DS) selection. Average rapidity densities hdN=dyi and average transverse momenta hp T i of charge-averaged pions, kaons, and protons as a function of center-of-mass energy are shown in Fig. 10 corrected to the DS selection (pp DS'). Based on the predictions of the three MC event generators studied, the inelastic hdN=dyi result is corrected upwards by 28%, with an additional systematic uncertainty of about 2%. No such correction is applied in the case of hp T i, since the inelastic value is close to the DS' one, with a difference of about 1%. The average p T increases with particle mass and event multiplicity at all ffiffi ffi s p , as predicted by all considered event generators. We note that both hp T i and ratios of hadron yields show very similar dependences on the particle multiplicity in the event, independently of the center-ofmass energy of the pp collisions. The ffiffi ffi s p evolution of the average hadron p T provides useful information on the socalled "saturation scale" (Q sat ) of the gluons in the proton [34]. Minijet-based models such as PYTHIA have an energydependent infrared p T cutoff of the perturbative multiparton cross sections that mimics the power-law evolution of Q sat characteristic of gluon saturation models [35]. In addition, the latter saturation models consistently connect Q sat to the impact parameter of the hadronic collision, thereby providing a natural dependence of hp T i on the final particle multiplicity in the event.

VII. SUMMARY
Transverse momentum spectra have been measured for different charged hadron species produced in inelastic pp collisions at ffiffi ffi s p ¼ 13 TeV. Charged pions, kaons, and protons are identified from the energy deposited in the silicon tracker and the reconstructed particle trajectory. The yields of such hadrons at rapidities jyj < 1 are studied as a function of the event charged particle multiplicity measured in the pseudorapidity range jηj < 2.4. The transverse momentum (p T ) spectra are well described by fits using the Tsallis-Pareto parametrization. The ratios of the yields of oppositely charged particles are close to unity, as expected in the central rapidity region for collisions at this center-of-mass energy. The average p T is found to increase with particle mass and event multiplicity, and shows features a slow (logarithmiclike or power-law) dependence on ffiffi ffi s p . As observed in lower-energy data, the hp T i and the ratios of particle yields are strongly correlated with event particle multiplicity. The PYTHIA 8 CUETP8M1 event generator reproduces most features of the measured distributions, which represents a success of the preceding tuning of this model, and EPOS LHC also gives a satisfactory description of several aspects of the data. Although soft QCD effects are intertwined with other effects, the present results could be used to further constrain models of hadron production and to contribute to a better understanding of multiparton interactions, parton hadronization, and final-state effects in high-energy hadron collisions.

ACKNOWLEDGMENTS
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Austria); FNRS and FWO (