Measurement of prompt charged-particle production in pp collisions at a centre-of-mass energy of 13 TeV

The differential cross-section of prompt inclusive production of long-lived charged particles in proton-proton collisions is measured using a data sample recorded by the LHCb experiment at a centre-of-mass energy of ${\sqrt{s} = 13\,\mathrm{TeV}}$. The data sample, collected with an unbiased trigger, corresponds to an integrated luminosity of ${5.4\,\mathrm{nb}^{-1}}$. The differential cross-section is measured as a function of transverse momentum and pseudorapidity in the ranges ${p_\mathrm{T} \in [0.08, 10)\,\mathrm{GeV}\,c^{-1}}$ and ${\eta \in [2.0, 4.8)}$ and is determined separately for positively and negatively charged particles. The results are compared with predictions from various hadronic-interaction models.


Introduction
Hadron production in inelastic high-energy proton-proton (pp) collisions is dominated by soft processes described by quantum chromodynamics (QCD). These processes cannot be calculated from first principles in perturbative QCD due to the large coupling constant α s at small average momentum transfer. Instead, predictions are based on phenomenological models. The determination of the model parameters relies on input from experiments at particle accelerators. Monte Carlo event generators, in which these models are implemented, are used at the Large Hadron Collider (LHC) to simulate the final-state particles originating from the soft component of a collision. An introduction to soft-QCD theories is presented, for example, in Refs. [1][2][3][4][5].
In the field of cosmic-ray research, generators are employed to simulate interactions of ultra-relativistic nuclei with the atmosphere of the Earth, which induce extensive particle cascades, referred to as air showers. Although often used to predict interactions in a phase space that is not covered by the input from experiments, the generators are remarkably successful at describing many features of air showers. However, a long-standing excess is observed in the number of muons produced in high-energy air showers compared to simulations, termed the Muon Puzzle [6,7]. Precision measurements of the production of light hadrons at the TeV scale in the forward region are needed to guide further and constrain the models [8][9][10] and to extrapolate them safely to even higher collision energies that are of interest in astroparticle physics. The LHCb detector covers the forward pseudorapidity range 2 < η < 5, which is of particular interest for cosmic-ray research.
A suitable proxy for the production of light hadrons is the production of prompt long-lived charged particles. A particle is classified as long-lived if its lifetime is greater than 30 ps, and prompt if it is either produced directly in the primary interaction or does not have long-lived ancestors [11]. Long-lived charged particles are electrons, muons, pions, kaons and protons as well as Σ + , Σ − , Ξ − and Ω − baryons and their antiparticles.
In this paper, a measurement of the double-differential cross-section of inclusive production of prompt long-lived charged particles in pp collisions at a centre-of-mass energy of √ s = 13 TeV is presented. The data sample was recorded by the LHCb experiment with a zero-bias trigger in 2015 and corresponds to an integrated luminosity of 5.4 nb −1 . The double-differential cross-section is measured separately for positively and negatively charged particles as a function of transverse momentum and pseudorapidity in intervals that span the ranges p T ∈ [80, 10 000) MeV/c and η ∈ [2.0, 4.8). The minimum transverse momentum is a function of η due to the fiducial acceptance of the spectrometer. Both the charge-combined differential cross-section and the ratio of the differential cross-sections for the two charges are compared with predictions from four different hadronic-interaction models. The measurement goes a step further compared to the precursors in several ways. It is a double-differential measurement in the forward region performed separately for each charge, non-prompt production is removed and there is no trigger bias. High precision is achieved thanks to the high-performance tracker of the LHCb experiment and the detailed control measurements that can be performed with a general-purpose spectrometer.
The paper is structured as follows. In Sect. 2, the detector as well as the data and simulated samples used in this measurement are described. The analysis strategy is presented in Sect. 3. The efficiencies and the background contributions are detailed in Sects. 4 and 5, respectively. In Sect. 6, the results are discussed, and a summary is provided in Sect. 7.

Detector and data sample
The LHCb detector [40,41] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary pp collision vertex, the impact parameter, is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The online event selection for this measurement is performed by an unbiased trigger. Therefore, no trigger-related systematic uncertainty arises. At the hardware stage, events are accepted at a fixed rate. The software stage then restricts the data sample to collisions of leading bunches of the LHC bunch trains, which avoids background from previous events, while the bunch spacing of 50 ns in these low-intensity runs avoids contributions to the read-out from following events. The analysed data sample contains the events from two LHC fills, recorded with opposite magnetic-field configurations of the LHCb dipole magnet. The field configuration that bends positively (negatively) charged particles in the horizontal plane towards the centre of the LHC ring is referred to as upwards (downwards).
The fill recorded with the magnetic field pointing upwards comprises 226 × 10 6 events and corresponds to an integrated luminosity of 3.0 nb −1 , while the fill recorded with the magnetic field pointing downwards comprises 134 × 10 6 events and corresponds to an integrated luminosity of 2.4 nb −1 . The average numbers of collisions in a bunch crossing are 0.9 and 0.7 for these two fills, respectively. The data analysis is independent of the number of collisions per event, which means that events with multiple collisions are not treated differently. The results are obtained from the combined data sample, but as a cross-check, the analysis is also performed separately for each fill. To measure the background from interactions of the beams with residual gas in the beam pipe, beam-gas collisions are used, where only one of the two beams traverses the detector. Such collisions were also collected for each fill.
Simulation is required to model the effects of the imposed selection requirements and to study the background contributions. In the simulation, pp collisions are generated using Pythia 8.1 [17] with a specific LHCb configuration [42]. Decays of unstable particles are described by EvtGen [43], in which final-state radiation is generated using Photos [44]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [45] as described in Ref. [46]. For each of the two magnetic-field configurations, a trigger-unbiased sample containing 10 7 events was simulated.

Analysis strategy
The differential cross-section is computed as where n is the number of prompt long-lived charged particles produced in beam-beam collisions in a pseudorapidity interval with a width of ∆η and a transverse-momentum interval with a width of ∆p T that is obtained from a data set corresponding to an integrated luminosity of L. It is computed in finite intervals of η and p T . The count n cannot be obtained directly, it is computed from the observable number of recorded tracks after applying corrections, as further detailed below. There is no special treatment for events with multiple collisions per bunch crossing in this approach; every collision adds to the particle count and to the luminosity. Tracks that traverse the entire tracking system are selected as candidates for prompt long-lived charged particles. Among these tracks, some are fakes, which do not correspond to a real particle. Fake tracks mostly originate from random matches of unrelated track segments upstream and downstream of the magnet. Two kinds of fake tracks are distinguished, depending on whether the fake track forms in isolation or nearby a real track. The first kind is reduced by imposing a requirement on the fake-track probability, P fake , provided by a neural-network-based algorithm [47], but the remaining contribution is still non-negligible. The second kind is always accompanied by a nearby track and through that signature already detected. It is suppressed by the track-reconstruction software below the level of 0.1 % in all kinematic intervals and therefore negligible in this analysis. Non-prompt tracks passing the selection are another source of background. These tracks originate from interactions of particles with the detector material or from decays of long-lived particles. Interactions of the beams with residual gas are another source of background. These background contributions are further discussed in Sect. 5 as well as how much they contribute to the number of candidates.
Consequently, the number of candidate tracks, n cand , is related to the number of signal particles, n, according to where ε denotes the total efficiency, i.e. the product of the geometric acceptance of the detector, the track-reconstruction efficiency and the selection efficiency, and the sum includes the numbers of background tracks, n i , from source i. The values of ε and n i are taken from simulation, except for the background from beam-gas interactions, n gas , which is determined from a control sample of pure beam-gas events. Control measurements are performed to correct for imperfect modelling in the simulation. For this purpose, observables, P i , that are proportional to n i are chosen as proxies. The background contributions and the corresponding proxies are described in Sect. 5. The ratio of the background counts in data and simulation is assumed to be equal to the ratio of the proxies in data and simulation, which allows the background count in data to be estimated as n i = R i n i, sim . In summary, the number of signal particles in Eq. (1) can be expressed as where R ε is a correction to the total efficiency in simulation. Due to the slow variation of the differential cross-section in η, only six intervals in η are used with ∆η = 0.5. The width of the last interval, η ∈ [4.5, 4.8), is reduced to match the acceptance of the tracking system. Since the p T spectrum has a powerlaw shape, 50 logarithmic intervals are used in the range p T ∈ [10, 10 000) MeV/c with ∆log 10 (p T /(MeV/c)) = 0.06. The lower edge of this grid extends beyond the fiducial range of p T intervals that are used in the analysis, which starts at 80 MeV/c or higher, depending on the η interval. The lower limit is determined by the fact that a particle needs approximately 2 GeV/c to reach the tracking stations downstream of the magnet. In this analysis, an even tighter requirement p > 5 GeV/c is used, since a correction of the efficiency, which is described in Sect. 4, is based on the two control measurements presented in Refs. [33,48], which both applied this tighter requirement.
The track-reconstruction efficiency depends on the detector occupancy. To ensure that the simulation reproduces the occupancy observed in data, weights are assigned to the simulated events. The number of tracks that traverse the entire tracking system is used as a proxy for the occupancy. Simulated events are weighted by the ratio of the distributions of the number of these tracks per event in data and simulation. As a positive side effect, this also adjusts the normalisation of the simulated sample with respect to the data sample. Whenever the simulation is mentioned in the following, it is always the occupancy-weighted simulation.
The two simulated samples are each weighted to reproduce the occupancy of the LHC fill with the same magnetic-field configuration. As a cross-check, the chosen proxy for the occupancy is compared to an alternative proxy, the number of hits in the scintillating-pad detector, which is not affected by possible artefacts of the track reconstruction. A linear relation is observed between the number of tracks that traverse the entire tracking system and the number of hits in the scintillating-pad detector, confirming this choice of the proxy for the detector occupancy.
The detector resolution distorts the measured kinematic distributions, by inducing migration between different (η, p T ) intervals. This distortion is between 1 % and 2 % in the simulation, depending on the kinematic interval. The simulated efficiency, ε sim , implicitly takes these migration effects into account. It is calculated from simulation as the ratio of the number of candidate tracks that can be associated with signal particles and the respective number of signal particles, where the candidates are sorted into intervals using the reconstructed momenta, while the signal particles are sorted into intervals using the generated momenta. Deviations between the shapes of the kinematic distributions of charged particles in data and simulation can introduce a bias in this approach, but the bias in this case is estimated to be at the order of 0.01 % and therefore negligible.

Efficiencies
The total efficiency to observe prompt long-lived charged particles depends on the geometric acceptance of the detector, the track-reconstruction efficiency, and the particle loss due to decays and interactions with the detector material. The losses depend on the particle composition and the amount of traversed material. The composition of the prompt longlived charged particles has an impact on the total efficiency, because pions and protons have different interaction lengths, for example.
In this study, the composition is broken into four particle classes: charged pions, charged kaons, protons, and other long-lived charged particles (e − , µ − , Σ + , Σ − , Ξ − and Ω − ). The detector acceptance and track-reconstruction efficiency are taken from simulation, but corrections from control measurements are applied, which are detailed below. The simulated efficiencies, ε sim , for the four particle classes are shown in Fig. 1. The efficiency for the last class is generally lower, because Σ + , Σ − , Ξ − and Ω − baryons and their antiparticles have the shortest lifetimes. At low p T , the efficiency is higher, since the relative yield of these hadrons is small and leptons contribute more, which have high tracking efficiencies. At high p T , the efficiency increases again, because the lifetimes of the strange hadrons are Lorentz boosted. This particle class contributes little to the total efficiency, since the particles in the class are less abundantly produced compared to the particles of the other three classes. A dip in the efficiency can be seen in the η ∈ [3.0, 3.5) interval at p T ≈ 500 MeV/c, which becomes more pronounced in the following η intervals. This feature is due to an increased amount of material that particles with certain trajectories have to pass [63]. Efficiencies for particles and antiparticles are slightly different due to the different hadronic cross-sections. The simulated efficiency is validated LHCb simulation [4.5, 4.8) sim × q Figure 1: Efficiencies for different particle species in intervals of η and as a function of p T for the simulated sample generated with the magnetic field pointing upwards. The product of the efficiency and the particle charge in units of the elementary charge is shown to separate values for oppositely charged particles. The light-grey areas indicate the limit of the kinematic acceptance. against data. Small corrections for differences in the tracking efficiency and for differences in the hadron composition are applied in this analysis.
The track-reconstruction efficiency is corrected based on results of a separate control measurement [48], in which muon tracks from J/ψ → µ + µ − decays were studied to determine ratios of the track-reconstruction efficiencies in data and simulation in intervals of η and p. A linear transformation matrix is built to map the results of the control measurement to the η and p T intervals used in this analysis. Non-uniform track density is accounted for in this transformation. The obtained efficiency ratios are identical for both particle charges and are used as the first component of the correction to the efficiency in simulation. The uncertainty of this component is between 1 % and 5 % over the kinematic range covered by this analysis.
The simulated particle composition is then corrected by adjusting the relative yield of each particle class to match data. Since the composition has not yet been measured in pp collisions at √ s = 13 TeV, LHCb measurements of ratios of prompt hadron production are extrapolated from pp collisions at √ s = 0.9 TeV and 7 TeV [33]. The hadron ratios p/p, (K + + K − )/(π + + π − ) and (p + p)/(π + + π − ) depend on η and p T , and are separately extrapolated to 13 TeV in each η and p T interval with a linear function in ln √ s. The largest deviations of up to 40 % are observed between the extrapolated and simulated (p + p)/(π + + π − ) ratios.
The extrapolated ratios obtained from the first step cover only three p T intervals and do not cover the intervals η ∈ [2.0, 2.5) and η ∈ : Efficiency corrections, R ε , for positively and negatively charged particles in intervals of η and as a function of p T for the simulated sample generated with the magnetic field pointing upwards. The bands indicate the systematic uncertainty. The light-grey areas represent the limit of the kinematic acceptance. The gap between this limit and the correction is due to the tighter requirement of p > 5 GeV/c applied to the control measurements on which this correction is based.
in this simulation for each η and p T interval. Forming double ratios removes most of the η and p T dependence of the hadron ratios, which makes their extrapolation more robust. The double ratios for the intervals η ∈ [2.0, 2.5) and η ∈ [4.5, 4.8) are taken from the corresponding adjacent η intervals for further analysis.
In a final step, correction factors for the counts of long-lived particles in the four classes are introduced. These correction factors are parameterised as linear functions of ln p T separately for each η interval, particle class, and charge. The parameters of these correction functions are optimised in a least-squares fit so that the double ratios formed from the corrected and the original hadron ratios in this simulation approach the double ratios obtained in the previous step, under the constraint that the total number of particles per charge remains the same. Since this optimisation problem is underconstrained by the extrapolated ratios, Gaussian penalty terms are introduced as a regularisation to suppress deviations above 5 % of the corrected counts from their initial values. The composition-corrected efficiency in simulation is then computed by summing the products of the efficiencies for the four particle classes with their corrected fractions. Since this correction involves considerable modelling, the systematic uncertainty dominates over the statistical uncertainty. Half of the correction is assigned as a systematic uncertainty as a conservative estimate, which is at most 2.5 %.
The combined correction of the effects discussed so far for the simulated efficiency is shown in Fig. 2 and reduces the original efficiency in simulation by up to 5 % for negatively charged particles. This reduction arises primarily from a relative increase in the fourth particle class of the particle composition, which has a lower efficiency. The step-like structures arise from the comparably coarse (η, p) intervals over which the control measurement for the track-reconstruction efficiency is performed. The transfer to an (η, p T ) grid smoothens these steps, but they remain visible. Figure 2 also defines the usable p T range in each η interval. The p T intervals without an efficiency correction are not used to compute final results.
Material interactions contribute to the systematic uncertainty of the corrected efficiency as the simulated amount of material has an uncertainty of 10 %, as described in Ref. [48]. The hadronic losses are proportional to the amount of traversed material and therefore carry the same relative uncertainty. The interaction losses of charged pions and charged kaons are estimated to be 14 % and 11 % in Ref. [48], respectively. For protons, the loss is between 20 % and 30 % across the full kinematic range. This loss is estimated from the difference between the efficiency for protons in simulation and that for muons. Since protons are stable and muon decays inside LHCb are negligible, any additional loss of protons relative to muons is caused by material interactions. Hadronic loss for the species in the fourth class of long-lived particles is negligible, as they are either leptons or their loss is dominated by decays. Overall, the hadronic losses have uncertainties of up to 3 %, depending on the particle species. Yet, the final uncertainty on the efficiency due to material interactions amounts to only 1 %, because the uncertainties contribute proportional to the abundance of the corresponding particle class.
The total systematic uncertainty of the corrected efficiency is the sum of the contributions from the correction of the track-reconstruction efficiency, the particle-composition correction, material interactions and the statistical uncertainty of the simulated samples. The latter ranges between 0.06 % and 2 %, depending on the kinematic interval, while the total uncertainty ranges from 0.9 % to 5.1 %.

Background contributions
The background induced by beam-gas interactions is determined from the number of candidate tracks produced in a control sample of recorded pure beam-gas events. Both the configuration where the beam travels from the vertex detector towards the muon system and the opposite configuration are included. The contributions from these configurations are scaled to the corresponding number of recorded beam-beam events. The fraction of candidate tracks originating from beam-gas interactions is found to be 1 % for the LHC fill recorded with the magnetic field pointing downwards, independent of the kinematic interval. For the fill recorded with the magnetic field pointing upwards, the fraction is below the level of 0.1 %.
The origins of the candidate tracks in simulation are shown in Fig. 3. Prompt long-lived charged particles constitute more than 85 % of the candidate tracks around p T = 1 GeV/c. Towards lower or higher values of p T , the background contributions increase. At high p T , fake tracks are the largest background, while fake tracks and tracks from photon conversions dominate at low p T . These background contributions are quantified using simulation and adjusted with correction factors from proxies obtained from data and simulation. The background contributions that are sufficiently large to require a correction are: fake tracks, tracks originating from material interactions of charged pions and from LHCb simulation [4.5, 4.8) fraction of candidate tracks Figure 3: Origins of the candidate tracks in intervals of η and as a function of p T for the simulated sample generated with the magnetic field pointing upwards. The prompt category refers to the signal tracks, while the other categories correspond to non-prompt tracks originating from the listed parent particles and include both decays and material interactions. Fake tracks are indicated by the white areas above the stacked histograms. The light-grey areas represent the limit of the kinematic acceptance.
photon conversions, and tracks produced in decays of strange hadrons. In the following subsections, the proxies constructed for these sources of background are presented. The remaining sources, which are not adjusted, are combined and a conservative uncertainty of 50 % is assigned for their contribution. The uncertainty of the differential cross-section resulting from this assumption is negligible, as the contribution from the remaining sources of background is small.

Fake tracks
The largest background at the edges of the analysed p T range is due to isolated fake tracks, while the contribution of fake tracks accompanied by a nearby track is negligible, as discussed in Sect. 3. These are characterised by having a low quality of the track fit and missing hits in instrumented detector parts. This information is combined with further variables from the tracking system to estimate a fake-track probability, P fake . The imposed selection requirement, P fake < 0.3, rejects between 40 % and 80 % of the fake tracks and retains 99 % of the real tracks. The remaining background from fake tracks, shown in Fig. 3 as white areas above the stacked histograms, ranges between 1 % and 27 % in the effective p T range derived from Fig. 2.
For each kinematic interval, data and simulation are divided into ten equal-width intervals in the P fake variable. Tracks of both charges are combined in this analysis, since : Distributions of P fake in a kinematic interval (left) at low p T and (right) at high p T for the LHC fill recorded with the magnetic field pointing upwards and the simulated sample with the same magnetic-field configuration. The statistical uncertainty of the data is indicated by error bars, but they are not visible on these scales. The purity of the proxy is shown in the lower panels, where light-grey boxes represent the P fake interval selected to determine the proxy.
Vertical lines indicate the threshold P fake < 0.3 for candidate tracks. The normalisation of the simulated sample is described in Sect. 3.
the background from fake tracks is charge symmetric. Distributions in two different kinematic intervals are given as examples in Fig. 4. In simulation, a pure sample of fake tracks is present at high P fake values. Generally, good agreement is found in the shapes of the P fake distributions between data and simulation, but in simulation, fewer fake tracks are observed. This can be seen in Fig. 4, where the intervals with P fake > 0.5 contain more entries in data compared to simulation. Consequently, a correction is used to adjust the simulated contribution from fake tracks. The number of tracks in the first interval of the P fake distribution above 0.3, where the fake-track purity is greater than 80 %, is chosen as a proxy for fake tracks. This choice balances two sources of systematic uncertainty. A P fake interval with a high fake-track purity is required in order to select a proxy that is insensitive to possible differences in the number of signal tracks between data and simulation, which would favour the interval with the highest purity. However, selecting an interval close to the P fake region of interest, P fake < 0.3, minimises the impact of shape differences in the P fake distributions between data and simulation.
The uncertainty of the fake-track proxy is the quadratic sum of statistical and systematic contributions. The statistical uncertainty is propagated from the track counts in data and simulation. The systematic uncertainty is estimated by computing an alternative proxy. For this, the integral of the track counts in the selected P fake interval and all the intervals above is used. The alternative proxy is more affected by shape differences and is therefore an upper limit on the systematic error. This type of systematic variation is modelled with a uniform distribution of deviations of which the upper limit corresponds to the alternative proxy, while the centre is the value of the default proxy. The uncertainty of the proxy is then given by the standard deviation of this distribution.
The resulting proxy ratio is shown in Fig. 5. At low η, values of the proxy ratio of up to 2.8 are observed, but the p T intervals in which the proxy ratio deviates strongly from unity are also those where the background from fake tracks is very small, as it can be seen by comparing with Fig. 3. Where the background from fake tracks is large, the proxy ratio is between 1.2 and 1.4. In general, the ratio is smooth as a function of p T . Discontinuities in the range η ∈ [2.5, 3.5) are caused by changes of the P fake interval chosen to determine the proxy.

Material interactions
Non-prompt tracks produced in material interactions constitute the second-largest source of background. As only tracks that traverse the entire tracking system are used, these are interactions occurring inside the vertex detector. Electrons originating from photon conversions, which populate mainly the kinematic intervals at high η and low p T , contribute up to 20 % to the candidate tracks, while charged particles stemming from hadronic material interactions of charged pions contribute less than 5 %.
The number of tracks originating from interactions in the material is proportional to the product of the total particle flux, consisting primarily of charged pions, and the amount of traversed material. The number of tracks that originate from candidate vertices, in which three tracks converge and which pass certain selection criteria that are detailed below, is chosen as a proxy for this product. Vertices with three tracks are used because hadronic interactions frequently produce three or more charged particles, while decays into three or more charged particles are rare in the analysed data sample. This approach is similar to that presented in Ref. [64]. Instead of defining a separate proxy for photon conversions, their simulated yield is scaled with the proxy for charged-pion interactions in the material. This is motivated by the observation that most photons originate from decays of neutral pions, which are produced in a fixed ratio to charged pions due to isospin symmetry. Furthermore, the conversion probability of photons and the interaction probability of charged pions are both proportional to the amount of material.
In each event, all unique combinations of three candidate tracks are formed first. Their point of closest approach is computed with a vertex fit minimising the sum of the squared distances of the tracks. This point is the candidate vertex of an interaction. To reduce combinatorial background, where three tracks are randomly associated with a common vertex, a lower limit is imposed on the distance of the vertex from the z axis, which points from the radial centre of the vertex detector along the beam axis towards the muon system and has its origin close to the average primary vertex. This discards the region in which no material is present.
The purity of the obtained sample of vertices is further increased by applying requirements that are optimised using simulation. First, a simultaneous optimisation of a lower limit on the z coordinate of the vertex and an upper limit on the sum of the squared distances of the tracks to their associated vertex, d, is performed without subdividing the sample intervals in η and p T . The z requirement suppresses combinatorial background by selecting vertices downstream of the region in which most of the primary vertices are reconstructed. The upper limit on d rejects random combinations by requiring that the tracks are compatible with originating from a single point. Second, an optimisation of a lower limit on the mass of the three-track combination is performed separately for each kinematic interval. As the species of the produced particles are not known, a zero-mass hypothesis is assigned to each track. The mass requirement further suppresses the combinatorial background. For both optimisations, the figure of merit S/ √ S + B is used, where S denotes the number of tracks that represent the signal for this proxy, and B is the number of background tracks. Wider p T intervals are used in this study to reduce fluctuations in the proxy ratio to be determined, created by successively merging groups of three adjacent p T intervals. The requirements chosen are those that maximise the value of the figure of merit. This results in a lower limit on z of 290 mm and an upper limit on d of 0.1 mm 2 . The lower limits on the mass depend on the interval and vary around 200 MeV/c 2 .
After the application of all optimised requirements, the contribution from charged pions is approximately equal in size to that from all other particles combined, e.g. kaons or protons. For each kinematic interval and each particle charge, the number of tracks that are included in selected vertices is used as the proxy. In the first two η intervals and around p T = 500 MeV/c, the highest values of the purity of charged-pion interactions are obtained, which range from 40 % to 45 %. Towards higher values of η and lower or higher values of p T , the purity decreases.
The proxy ratio is independent of the purity for kinematic intervals with purities above 30 %. Therefore, the value of the proxy ratio, R mat , is only used if the purity in the kinematic interval under consideration is greater than the threshold value of 30 %. Moreover, an upper limit of 30 % is imposed on the statistical uncertainty of the purity to reject intervals that contain only a small number of tracks. As a substitute for the values of the proxy ratio in the intervals below the purity threshold, the median of the values of the proxy ratio in the intervals above the threshold is computed. The uncertainty of this median value is estimated assuming a uniform distribution covering the observed range of R mat . The systematic uncertainty of the proxy for material interactions is determined by loosening and tightening the optimised z and d requirements and computing alternative proxy ratios, R mat . The mass requirement is not varied, as a good agreement is found in the shapes of the distributions between data and simulation in all kinematic intervals. Deviations in some of the alternative proxies cannot be explained by statistical fluctuations and are therefore considered as a systematic uncertainty. The deviation of an alternative is considered significant if the condition is fulfilled, where σ 2 denotes the variance, which is a consequence of using the same data set to determine the default proxy ratio and the alternatives. Following Ref. [65], the systematic uncertainty is taken as the standard deviation of the variations that differ significantly and the initial ratio is replaced with the mean of the significant variations. This procedure is applied to the values of the proxy ratio in the kinematic intervals above the purity threshold and to the median value computed for the other intervals. The values

Strange-hadron decays
Up to 7 % of the candidate tracks are non-prompt tracks produced in decays of strange hadrons inside the vertex detector. Approximately 80 % of these tracks originate from K 0 S → π + π − , Λ → pπ − and Λ → pπ + decays. The number of K 0 S mesons and Λ and Λ baryons, of which both decay products create a candidate track, is chosen as a proxy for this background.
Candidate decays are obtained by forming pairs of oppositely charged candidate tracks. A small set of loose requirements is used to maintain high selection efficiency for strangehadron decays. The distance of closest approach of the two tracks is required to be less than 1 mm and the distance along the z axis of their point of closest approach from the average primary vertex is required to be greater than 150 mm. For K 0 S candidates, the mass is computed by assigning the pion-mass hypothesis to both tracks, while for Λ and Λ candidates, the proton-and pion-mass hypotheses are assigned to the tracks.  Figure 8: Ratio of the proxies for strange-hadron decays in data and simulation in intervals of η and as a function of p T for the LHC fill recorded with the magnetic field pointing upwards and the simulated sample with the same magnetic-field configuration. The V 0 model lines indicate the interpolated K 0 S , Λ and Λ yield ratios, with support points indicated by black dots. The lines labelled as products represent the proxy ratio of the decay products, with the bands representing the propagated systematic uncertainty. The light-grey areas indicate the limit of the kinematic acceptance. Fig. 7. Extended-maximum-likelihood fits [66] to the mass distributions of the K 0 S , Λ and Λ candidates in data and simulation are performed separately for each (η, p T ) interval of their kinematic distributions. As in the case of the proxy for material interactions, wider p T intervals are used in this study, obtained by merging groups of three adjacent p T intervals. In the fits, the K 0 S , Λ or Λ signal is modelled with a nonstandardised Student's t-distribution with free location and width parameters. The background, which is only combinatorial, is modelled with a third-degree Bernstein polynomial. The yields in each interval of the mass distributions are modelled with Poisson distributions for the unweighted histograms in data and with scaled Poisson distributions [67] for the occupancy-weighted histograms in simulation.

are shown as examples in
The ratios of the signal yields in data and simulation are only computed in those kinematic intervals where a signal peak is present with a significance of at least three standard deviations. The ratios are shown in Fig. 8. In all η intervals and for all three strange-hadron species, a pattern is observed in the yield ratios as a function of p T . This is interpreted as a general difference in the amount of produced strange quarks between data and simulation, since it is the same for a meson and a baryon. The pattern is modelled with a monotone cubic spline [68]. A least-squares fit of the model to the yield ratios is performed, in which a Gaussian penalty term is added in order to restrict the value of the spline at the upper limit of the p T range, where no data points are present.
The fitted model is then used to determine the proxy ratio of the strange hadrons, R parent , over the full kinematic range. However, the background for prompt long-lived charged particles is caused by the decay products of these hadrons. The proxy ratio of the products, R strange , in the (η, p T ) intervals of their kinematic distributions is computed as where: h is K 0 S , Λ, Λ; i and j iterate through the η and p T intervals of the kinematic distributions of the parent particles, respectively; k and iterate through the η and p T intervals of the kinematic distributions of the decay products, respectively; and n is the simulated yield of the decay products in the corresponding intervals. The ratio R strange is closer to unity than the ratio R parent , as the broad kinematic distributions of the decay products dilute deviations. The statistical uncertainty of the fitted model is negligible compared to the systematic uncertainty suggested by the deviations of the yield ratios from the model. A systematic uncertainty of ±15 % is assigned to the proxy ratio to cover these deviations.

Results
The differential cross-section of prompt inclusive production of long-lived charged particles is determined with Eqs. (1) and (4).
The statistical uncertainty of the data sample is subdominant in the full kinematic range, reaching 1.5 % at high p T . Non-Poissonian statistical fluctuations of the candidate counts in kinematic intervals are taken into account. Each event contributes entries to a number of kinematic intervals, introducing statistical correlations of up to 0.4 between intervals, and variances that are up to 50 % larger than the Poisson expectation. The covariance matrix of these statistical fluctuations is estimated by dividing the full sample into 100 equivalent subsets, computing the covariance of the subsets and extrapolating the result back to the full data set.
The uncertainty on the differential cross-section is computed from the uncertainties of the individual terms in Eqs. (1) and (4) using full error propagation of the respective covariance matrices. The contributions to the uncertainty of the number of prompt long-lived charged particles are shown in Fig. 9. The correction of the track-reconstruction efficiency is the largest contribution in most intervals. In the range η ∈ [2.0, 4.0) and at high p T , the uncertainty of the proxy for fake tracks dominates, while in the range η ∈ [4.0, 4.8) and at low p T , the uncertainty of the proxy for material interactions contributes most.
The final results are obtained by adding the fully corrected and background-subtracted numbers of prompt long-lived charged particles from both LHC fills and dividing by the combined integrated luminosity, after confirming that the two separate results from each fill are consistent with each other. As a further cross-check, data from six other fills, recorded at the same centre-of-mass energy and with the same trigger as the default data sample, are compared and found to be consistent with the main result. For the calculation of the covariance matrix of the combined result, statistical uncertainties are treated as uncorrelated, while all systematic uncertainties, including the uncertainty of the luminosity, are taken to be fully correlated between the fills.

Model or experiment Inelastic cross-section in mb
Pythia 8.3 [17] 78.05 EPOS-LHC [15] 78.98 QGSJet II-04 [16] 80.17 Sibyll 2.3d [12] 79 79.5 ± 1.8 The minimum and maximum values over all kinematic intervals for the uncertainties from each source that contribute to the final result are listed in Table 1. The total uncertainty of the differential cross-section is between 2.3 % and 15 %, which includes the fully correlated systematic uncertainty of 2.0 % from the integrated luminosity. The correlations are nonzero and positive between kinematic intervals and the two particle charges, since the systematic uncertainties dominate. The correlation matrix of the final differential cross-section is shown in Appendix B.
The measured differential cross-section of inclusive production of prompt long-lived charged particles is shown in Fig. 10. The values for positively and negatively charged particles are listed in Appendix A. In the figure, the measurement is compared with the predictions from the hadronic-interaction models Pythia 8.3 [17] with the default Monash tune [72], EPOS-LHC [15], QGSJet II-04 [16] and Sibyll 2.3d [12]. The latter three models are accessed through the CRMC package [73]. Also shown, but not comparable to the other models, is the occupancy-weighted LHCb tune of Pythia 8.1 that was used as the simulated sample in this analysis. The calculation of the differential cross-section, is based on: the inelastic cross-section, σ inel , that is implemented in each model; the number of generated inelastic events, N inel ; and the number of prompt long-lived charged particles, n, in each kinematic interval. The inelastic cross-sections of the models and recent measurements are listed in Table 2.
The deviations of the predictions from the measured values are between −26 % and +170 %. The models mostly overestimate the differential cross-section. These deviations are caused by differences in the kinematic distributions compared to the data and a potential difference in the overall multiplicity of produced particles. The latter would cause an overall change in the magnitude of the differential cross-section. Differences in the inelastic cross-section have a similar effect, but can be excluded as the driving factor since the implemented values of the inelastic cross-section in the models are compatible with the measured ones. The occupancy-weighted LHCb tune of Pythia 8.1 shows deviations between −17 % and +32 %. The comparably good agreement of this simulation is based on the occupancy weighting, which allows only deviations in the shape, while for the remaining models, deviations in the magnitude can occur. Out of the other models,  Figure 11: Ratios of the differential cross-sections of inclusive production of prompt longlived positively and negatively charged particles as a function of transverse momentum and pseudorapidity for the data and the models shown in Fig. 10.
EPOS-LHC shows the smallest differences overall. The shape of the p T distribution is well described; data and model differ mainly by an overall offset that weakly depends on η.
To address the question from the introduction, whether the pseudorapidity distribution is wide or narrow, the differential cross-section from LHCb needs to be integrated over the p T range and combined with the other measurements at lower and higher η values. This combination is complicated by the fact that other experiments use (minimum-)biased triggers, while this analysis is unbiased. Therefore, this combination is not attempted in this paper. The LHCb data do not cover the peak of the p T distribution in all η intervals. Since the peak contributes most to the integral, the integration of the LHCb data is model dependent for η < 3.5.
The ratio of the differential cross-sections for positively and negatively charged particles is shown in Fig. 11. The positively correlated components of the uncertainties cancel in the computation of the ratio. At η > 4.0 and low p T , the uncertainty of the ratio is large due to the subtraction of background from material interactions, which cannot be assumed to be charge symmetric and which is not well constrained by control measurements in that region. At high η and high p T , the production of positively charged particles increases as the initial state has a charge of +2, which transfers to the final state in the forward region. This effect is predicted by the models to a varying extent. The best description of the ratio is provided by Pythia 8.3, although the onset of the increase is shifted towards lower p T values. QGSJet II-04 predicts an onset at even lower p T values. EPOS-LHC and Sibyll 2.3d do not show a significant increase up to p T = 10 000 MeV/c, which is at tension with the measurement.

Summary
An unbiased measurement of the differential cross-section of prompt inclusive production of long-lived charged particles in pp collisions at √ s = 13 TeV is presented. The data sample was recorded by the LHCb experiment and corresponds to an integrated luminosity of 5.4 nb −1 . The differential cross-section is measured as a function of transverse momentum and pseudorapidity in the ranges p T ∈ [80, 10 000) MeV/c and η ∈ [2.0, 4.8) and is determined separately for positively and negatively charged particles. An uncertainty between 2.3 % and 15 %, depending on p T and η, is achieved. A comparison of the measured charge-combined differential cross-section with predictions from recent hadronic-interaction models shows that these models mostly overestimate the differential cross-section. The overall smallest deviations are observed for EPOS-LHC, while the ratio of the differential cross-sections for positively and negatively charged particles is best predicted by Pythia 8.3. The precision achieved in this measurement is essential for an improved simulation of the underlying event in collisions at the LHC and of interactions in the atmosphere of the Earth that cause air showers.

Acknowledgements
We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies:

Appendices A Differential cross-sections
The measured values of the differential cross-sections of prompt inclusive production of long-lived positively and negatively charged particles are listed in Table 3. The numerical values and the correlation matrix are provided in machine-readable form at [link]. Table 3: Differential cross-sections of prompt inclusive production of long-lived positively and negatively charged particles as a function of transverse momentum and pseudorapidity.

B Correlation matrix of differential cross-section
The correlation matrix of the measured differential cross-section is shown in Fig. 12. The correlations are positive since the systematic uncertainties dominate and most of the systematic uncertainties are positively correlated between kinematic intervals and the two particle charges. The positive correlations between neighbouring p T intervals are particularly large due to the correction of the tracking efficiency, which affects neighbouring intervals in the same way. The numerical values of the correlation matrix are provided in machine-readable form at [link].  Figure 12: Correlation matrix for the uncertainties of the differential cross-section of prompt inclusive production of long-lived charged particles. The four large quadrants correspond to the correlations between negatively and positively charged particles. The 36 cells within each quadrant correspond to the η intervals. In each cell, the correlations of the p T intervals are shown from low to high p T in logarithmic intervals.    [70] LHCb collaboration, R. Aaij et al., Measurement of the inelastic pp cross-section at a centre-of-mass energy of √ s =13 TeV, JHEP 06 (2018) 100, arXiv:1803.10974.