Charged-particle pseudorapidity density at mid-rapidity in p–Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pmb {\sqrt{s_{\scriptscriptstyle {\mathrm{NN}}}}}$$\end{document} = 8.16 TeV

The pseudorapidity density of charged particles, dNch/dη\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {d}N_{\mathrm{ch}}/\mathrm {d}\eta $$\end{document}, in p–Pb collisions has been measured at a centre-of-mass energy per nucleon–nucleon pair of sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\scriptscriptstyle {\mathrm{NN}}}}$$\end{document} = 8.16 TeV at mid-pseudorapidity for non-single-diffractive events. The results cover 3.6 units of pseudorapidity, |η|<1.8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\eta |<1.8$$\end{document}. The dNch/dη\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {d}N_{\mathrm{ch}}/\mathrm {d}\eta $$\end{document} value is 19.1±0.7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$19.1\pm 0.7$$\end{document} at |η|<0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\eta |<0.5$$\end{document}. This quantity divided by ⟨Npart⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle N_{\mathrm{part}} \rangle $$\end{document} / 2 is 4.73±0.20\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4.73\pm 0.20$$\end{document}, where ⟨Npart⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle N_{\mathrm{part}} \rangle $$\end{document}is the average number of participating nucleons, is 9.5% higher than the corresponding value for p–Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\scriptscriptstyle {\mathrm{NN}}}}$$\end{document} = 5.02 TeV. Measurements are compared with models based on different mechanisms for particle production. All models agree within uncertainties with data in the Pb-going side, while HIJING overestimates, showing a symmetric behaviour, and EPOS underestimates the p-going side of the dNch/dη\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {d}N_{\mathrm{ch}}/\mathrm {d}\eta $$\end{document} distribution. Saturation-based models reproduce the distributions well for η>-1.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta >-1.3$$\end{document}. The dNch/dη\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {d}N_{\mathrm{ch}}/\mathrm {d}\eta $$\end{document} is also measured for different centrality estimators, based both on the charged-particle multiplicity and on the energy deposited in the Zero-Degree Calorimeters. A study of the implications of the large multiplicity fluctuations due to the small number of participants for systems like p–Pb in the centrality calculation for multiplicity-based estimators is discussed, demonstrating the advantages of determining the centrality with energy deposited near beam rapidity.


Introduction
Particle production in proton-nucleus (pA) collisions is influenced by nuclear effects in the initial state. In particular, p-Pb collisions are a valuable tool to study initial-state effects, which are present as a consequence of the nucleons being bound into nuclei. Additionally, the particle multiplicity is an important tool to study the various theoretical models of gluon saturation, which contain different treatments of the upper limit in the growth of the parton density. Therefore, pseudorapidity density measurements can provide constraints to the modelling of the initial state at small Bjorkene-mail: alice-publications@cern.ch x. Moreover, evidence for collective phenomena have been observed in p-Pb collisions, with the magnitude of the effects increasing with event multiplicity [1][2][3][4][5][6][7][8][9]. Proton-nucleus collisions serve as a tool to study also final-state effects that are sensitive to the formation of a Quark-Gluon Plasma in heavyion collisions, under active scrutiny by the community [10]. For these reasons, it is important to understand the collision geometry and the global properties of the system produced in p-Pb collisions.
This paper presents a measurement of the primary chargedparticle density in p-Pb collisions, dN ch /dη lab , at a nucleonnucleon centre-of-mass energy of √ s NN = 8. 16 TeV for pseudorapidities |η lab | < 1.8 in the laboratory system. A primary charged particle is defined as a charged particle with a mean proper lifetime τ larger than 1 cm/c, which is either produced directly in the interaction, or from decays of particles with τ smaller than 1 cm/c, excluding particles produced in interactions with the beam pipe, material of the subdetectors, cables and support structures [11]. The dominant processes in p-Pb collisions are the non-diffractive ones. Diffractive events can be single-, double-or central-diffractive and results are presented for non-single-diffractive (NSD) events. Data are compared to other experimental measurements available in pp, p-Pb, d-Au and AA collisions. Results are compared also with simulations (performed with HIJING 2.1 [12,13], EPOS 3 [14][15][16] and EPOS LHC [17]) and calculations incorporating the saturation of the gluon density in the colliding hadrons (MC-rcBK [18,19] and KLN [20,21]).
The rest of this article is organised in the following way: Sect. 2 describes the experimental conditions and the detectors used to measure the centrality of the event and the pseudorapidity density of charged particles. In Sect. 3, the centrality determination methodologies are described, both the ones using the multiplicity distributions of charged particles and the alternative one that relies on the energy collected in the neutron Zero-Degree Calorimeters (ZDCs). Section 4 explains, in detail, the analysis procedure to measure the dN ch /dη. The systematic uncertainties are described in Sect. 5, and the results along with comparisons to models are presented in Sect. 6. A brief summary and conclusions are given in Sect. 7.

Experimental setup
The p-Pb data were provided by the Large Hadron Collider (LHC) in December 2016. There were two configurations that were exploited: in one, denoted by p-Pb below, the proton beam circulated towards the negative z direction in the ALICE laboratory system, while 208 Pb ions circulated in the opposite direction; in the second configuration, denoted by Pb-p, the direction of both beams was reversed. The total luminosity was 0.06 nb −1 , corresponding to around 120 million minimum-bias (MB) events in the p-Pb and Pb-p configurations. The beams in both rings have the same magnetic rigidity. The nucleon-nucleon centre-of-mass energy was √ s NN = 8. 16 TeV, with both p and Pb beams at 6.5 TeV per proton charge. Due to the asymmetric collision system, there is a shift in the centre-of-mass rapidity of y = 0.465 in the direction of the proton beam.
Full details of the ALICE detector are given elsewhere [22,23]. The main element used for the analysis was the Silicon Pixel Detector (SPD): the two innermost cylindrical layers of the ALICE Inner Tracking System [22], made of hybrid silicon pixel chips. The SPD is located inside a solenoidal magnet that provides a magnetic field of 0.5 T. The first layer covers |η lab | < 2.0 for collisions at the nominal Interaction Point (IP), while the second covers |η lab | < 1.4. The layers have full azimuthal coverage and radii of 3.9 cm and 7.6 cm, respectively. In total, the SPD has 9.8 × 10 6 silicon pixels, each of size 50 × 425 mm 2 .
The MB trigger signal is given by a hit in both the V0 hodoscopes [24]. The V0 detector is composed of two arrays of 32 scintillators positioned at 3.3 m (V0A) and -0.90 m (V0C) from the nominal IP along the beam axis. Each array has a ring structure segmented into 4 radial and 8 azimuthal sectors. The detector has full azimuthal coverage in the pseudorapidity ranges 2.8 < η lab < 5.1 and −3.7 < η lab < −1.7. The signal amplitudes and particle arrival times are recorded for each of the 64 scintillators. The V0 is well suited for triggering thanks to its good timing resolution (below 1 ns) and its large angular acceptance. The timing is used to discriminate the beambeam collisions from background events, like beam-gas and beam-halo events, produced outside the interaction region. The neutron ZDCs [25] are likewise utilised for background rejection. The neutron calorimeters, ZNs, are quartz-fibre spaghetti calorimeters placed at zero degrees with respect to the LHC beam axis, positioned at 112.5 m (ZNA) and −112.5 m (ZNC) from the nominal IP. ZNs detect neutral particles emitted at pseudorapidities |η lab | > 8.7 and have an energy resolution of around 18% for neutron energies of 2.56 TeV.
ALICE is equipped also with the proton calorimeters, ZPs, which are not used in the analysis.
A subsample of 6.8 million events is analysed for p-Pb collisions, with an average number of interactions per bunch crossing, μ of 0.004. A subsample of 2.7 million events is analysed for Pb-p collisions, with μ = 0.007. The comparison of p-Pb and Pb-p results is used to assess the systematic uncertainties. The hardware MB trigger is configured to have high efficiency for hadronic events, requiring a signal in both V0A and V0C. Beam-gas and beam-halo interactions are suppressed in the analysis by requiring offline the arrival time of particles in the V0 and ZN detectors to be compatible with collisions from the nominal IP. The contamination from background is estimated to be negligible through control triggers on non-colliding bunches.
The event sample after trigger and timing selection consisted of NSD, single-diffractive (SD), and electromagnetic (EM) interactions. The MB trigger efficiency for NSD events is estimated to be 99.2% using the DPMJet Monte Carlo event generator [26], and 99.5% using HIJING 1.36 [27]. HIJING 1.36 combines perturbative-QCD processes with soft interactions, and includes a strong impact parameter dependence of parton shadowing. DPMJet is based on the Gribov-Glauber approach and treats soft and hard scattering processes in a unified way. It includes incoherent SD collisions of the projectile proton with target nucleons; these interactions are concentrated mainly on the surface of the nucleus. The generated particles are transported through the experimental setup using the GEANT3 [28] software package. SD collisions are removed in DPMJet by requiring that at least one of the binary nucleon-nucleon interactions is NSD. The SD and EM contaminations are estimated from Monte Carlo simulation studies to be around 0.03% and below 0.3%, respectively.
Among the selected events in data, 99% had a primary interaction vertex. In DPMJet this fraction was 99.6% (99.8% for HIJING 1.36), with a trigger and selection efficiency for events without a primary vertex of 28% (23.1%). Taking into account the difference of the fraction of events without a vertex in the data and the simulation, the overall selection efficiency for NSD events in the analysis is estimated to be 97.0% (96.2%) according to DPMJet (HIJING 1.36).

Centrality determination
The Glauber model [29,30] is used to calculate the number of participating nucleons (participants), N part , and the corresponding number of nucleon-nucleon collisions, N coll , which depend on the collision impact parameter, b. Indeed, the number of produced particles changes with the variation of the amount of matter overlapping in the collision region; N part and N coll describe quantitatively this variation. In pA collisions, N coll = N part −1. Using the Glauber model, it is possible to calculate the probability distributions of the relevant parameters, N part and N coll , which for pA collisions are loosely correlated to b. Centrality classes are defined as percentile intervals of the visible cross section, which determines the event sample after the selections described in Sect. 2. The number of participating nucleons and nucleon-nucleon collisions are calculated, accordingly, for the visible cross section.
The centrality is determined for three different estimators, two of which are based on observables well separated in pseudorapidity to limit the effect of short-range correlations in the collision region. The method founded on multiplicity-based estimators is derived by fitting the measured charged-particle multiplicity distributions with an N coll distribution obtained from the Glauber model convoluted with a Negative Binomial Distribution (NBD) to model the multiplicity produced in a single collision. Multiplicity fluctuations play an important role in pA collisions. The range of multiplicities used to define a centrality class in the case of pA collisions is of the same order of magnitude as the multiplicity fluctuations width [31]. Therefore, a biased sample of nucleonnucleon collisions is selected using multiplicity. Samples of high-multiplicity events select not only a class with larger than average N part , but also one which is widely spread in N coll and that leads to deviations from the scaling of hard processes with Multiple Parton Interactions (MPI). These high-multiplicity nucleon-nucleon collisions have a higher particle mean transverse momentum p T , and are collisions where MPI are more likely [4]. The opposite happens for low-multiplicity events.
The centrality determined from the hybrid method, described in Sect. 3.2 using the energy deposited in the ZDCs, on the contrary, minimises biases on the binary scaling of hard processes. Indeed, the ZDCs detect, at large η separation from the central region, the nucleons produced in the interaction through the nuclear de-excitation process or knocked out by participants (called slow nucleons). A heuristic approach based on extrapolation from low-energy data is discussed in a previous publication [31].

Centrality from charged-particle distributions
In the method based on multiplicity estimators [31], the events are classified into centrality classes using either the number of clusters in the outer layer of the SPD (CL1 estimator) with acceptance η lab < 1.4, or the amplitude measured by the V0 in the Pb-remnant side, A-side, for p-Pb (V0A estimator) or in the C-side for Pb-p (V0C estimator) collisions. The amplitudes are fitted with a Monte Carlo implementation of the Glauber model assuming that the number of sources is given by the N part /2 convoluted with an NBD, which is the assumed particle production per source, parametrised with μ and k, where μ is the mean multiplicity per source and k controls the contribution at high multiplicity. The nuclear density for Pb is modelled by a Woods-Saxon distribution for a spherical nucleus with a radius of 6.62 ±0.06 fm and a skin thickness of 0.55 ± 0.01 fm [32]. The hard-sphere exclusion distance between nucleons is 0.40 ± 0.40 fm. For √ s NN = 8. 16 TeV collisions, an inelastic nucleon-nucleon cross section of 72.5 ± 0.5 mb is used, obtained by interpolation of cross section experimental values [32].
The measured V0A distribution with the NBD-Glauber fit is shown in Fig. 1. A similar fit has been performed for the CL1 estimator. The failure of the chosen fit function for amplitudes smaller than about 10 is due to trigger inefficiencies in peripheral collisions. The average number of participants, collisions and nuclear overlap function, T pPb , are calculated from the NBD-Glauber simulation for every defined centrality class. The values for the different estimators are given in Table 1. The systematic uncertainties are obtained by repeating the fit, varying the Glauber parameters (radius, skin thickness and hard-sphere exclusion) within their uncertainties. The number of participants for all selected events is on average N part = 8.09 ± 0.17. The increase in the average N part , when calculated for NSD collisions only, is of around 2% and within systematic uncertainties. The geometrical properties determined with the NBD-Glauber model are robust and approximately independent of the centrality estimator used, within the model assumptions of this approach.

Centrality from Zero degree Calorimeter and the hybrid method
The ZNs detect the slow neutrons produced in the interaction. The multiplicity of slow nucleons is monotonically related to N coll , and can, therefore, be used to determine the cen-  trality of the collision [31]. The ZPs are not used, since the uncertainty on N coll would be much larger. The experimental distribution of the neutron energy spectrum measured in the Pb-going side, E ZNA , is shown in Fig. 2 and it is used for the hybrid method, which aims to provide an unbiased centrality estimator. It is based on two assumptions, the first is that the event selection based on the energy deposited in the ZDCs is free from the multiplicity fluctuation biases in the particle production at mid-rapidity. The second assumption is that the wounded nucleon model holds [33] and that some observables, defined below, scale linearly with N coll and N part allowing one to establish a relationship to the collision geometry. Two sets of N coll are calculated: N mult coll and N Pb−side coll for each centrality bin i estimated using ZN. The first set is computed assuming that the charged-particle multiplicity at mid-rapidity is proportional to the N part : where N part MB is the average number of participating nucleons in MB collisions reported in Table 1, and, consequently: N coll The second set is calculated using the Pb-side multiplicity: N coll where S is the raw signal of the innermost ring of V0A for p-Pb (4.5 < η lab < 5.1) and V0C for Pb-p collisions (−3.7 < η lab < −3.2). A comparison of the N coll values obtained for the various estimators is reported in Table 2 for p-Pb collisions. The two different sets are consistent among each other and with the values calculated for Pb-p. The systematic uncertainties come from the uncertainty on the N coll for 0-100% in Table 1 summed with the maximum difference between the N mult coll and N Pb−side coll .

Analysis procedure
The technique for the dN ch /dη lab measurement is the same as the one employed at √ s NN = 5.02 TeV [31,34]. The pseudorapidity acceptance in the laboratory system depends on the position of the primary interaction vertex along the beamline, z vtx . The position of the primary vertex is obtained by correlating hits in the two silicon-pixel layers (SPD vertex). The selection of a reconstructed vertex within |z vtx | < 15 cm allows a range of |η lab | < 1.8 to be covered. In order to maximise the pseudorapidity coverage, instead of tracks we use tracklets (short track segments) formed using two hits in the SPD, one in the first and one in the second layer. In order to select combinations corresponding to charged particles, the angular difference in the azimuthal direction, ϕ, and in the polar direction, θ , of the inner and outer layer hit with respect to the reconstructed primary vertex is determined for each pair of hits. Afterwards, the sum of the squares of the weighted differences in azimuth and polar angles δ 2 = ( ϕ/σ ϕ ) 2 + ( θ/σ θ ) 2 is required to be less than 1.5, where σ ϕ = 60 mrad and σ θ = 25 sin 2 θ mrad, where the sin 2 factor takes the dependence of the pointing resolution on θ into account. With such a requirement, tracklets corresponding to charged particles with p T > 50 MeV/c are effectively selected. Particles with lower p T are mostly absorbed by the detector material or lost due to the bending in the magnetic field. A cross check utilising pp collisions [35] has shown full compatibility of analyses using tracklets and tracks, where the tracks have been reconstructed in the Time Projection Chamber matched with clusters in the Inner Tracking System. The raw multiplicity measured by tracklets needs to be corrected for (i) the acceptance and efficiency of a primary track to be reconstructed as a tracklet, (ii) the contribution from combinatorial tracklets, i.e. those whose two hits do not originate from the same primary particle, (iii) the difference between the fraction of events without a vertex in the data and in the simulation and (iv) the secondary-particle contamination. The first three corrections are computed using simulated data from the HIJING 1.36 or DPMJet event generators. The centrality definition in the simulated data is adjusted such that the particle density is similar to that in real data for the same centrality classes. The correction factors (i) and (ii), determined as a function of z and η lab , are on average around 1.5 for the acceptance and reconstruction efficiency, and around 0.02 for the combinatorial background removal in MB and centrality-dependent measurements at mid-rapidity, independently of the estimator selected and the centrality class. At |η lab | = 1.8 the combinatorial background contribution reaches a maximum value of 0.07. We further correct the measurement by the difference in the fraction of events without a vertex observed in data and simulation. The correction for MB dN ch /dη lab amounts to 2.2% (3.4%) when using DPMJet (HIJING 1.36). Since the centrality classes are defined as percentiles of the visible cross section, the centrality-dependent measurements are not corrected for the trigger inefficiencies. Differences in strange-particle content observed at lower beam energies [6,36] have been used for a data-driven correction applied to the generator output, giving rise to a correction factor of −0.6%, independent of centrality.

Systematic uncertainties
Several sources of systematic uncertainties were investigated. The uncertainty coming from the selection of the tracklet quality value δ 2 is negligible at mid-rapidity and amounts to 0.5% at |η lab | = 1.8. The other uncertainties associated to the MB dN ch /dη lab are independent of the pseudorapidity. The uncertainty resulting from the subtraction of the contamination from weak decays of strange hadrons is estimated to be about 1.3%. It is estimated by varying the amount of strange particles except kaons by ±50%. The uncertainty in detector acceptance and reconstruction efficiency is estimated to be 2.2% by carrying out the analysis for different slices of the z vtx position distribution and with subsamples in azimuth. The measurement for Pb-p collisions gives rise to an additional contribution of 1.8%, when reflected in η lab , for the most peripheral centrality bins (80-100%), and 1.1% for 60-80% at |η lab | = 1.8, and is added to the systematic uncertainty for acceptance. For the other centrality bins and the MB result the difference among p-Pb and Pb-p is negligible and already accounted for in the acceptance and reconstruction efficiency uncertainty. The uncertainty related to the trigger and event selection efficiency for NSD collisions is estimated to be 0.8% by taking into account the differences in the efficiency obtained with HIJING 1.36 and DPMJet. An additional 1.2% uncertainty comes from the difference in the scaling factors due to the events without vertex using the two event generators, as discussed in Sect. 4. A Monte Carlo test was also carried out with DPMJet to check the difference in the results obtained from NSD generated events and from selected events, resulting in a difference of 0.2% for the MB result, absorbed in the trigger efficiency uncertainty, and of 1.7% (0.2%) for 80-100% (60-80%) centrality bins. The contribution due to the subtraction of the background is studied using an alternative method where fake hits are injected into real events and it gives rise to a 0.3% uncertainty. The uncertainty from the material budget is 0.1%, while the uncer- tainty due to the particle composition amounts to 0.3%. The contributions from the extrapolation down to zero p T and from the pileup are found to be negligible. The final systematic uncertainties assigned to the measurements are the quadratic sums of the individual contributions. An overview of the systematic uncertainties is presented in Table 3. For MB dN ch /dη lab , they amount to 3.0%. For centrality-dependent measurements the total uncertainty for central events is 2.6%. For the most peripheral events it is 3.1% at mid-rapidity and 3.6% for |η lab | = 1.8. The difference in uncertainty between the MB and the centralitydependent measurement is mostly due to the contributions from the selection efficiency for NSD, which are not included in the centrality-dependent measurement, and to the difference among p-Pb and Pb-p collisions, which is more relevant for the most peripheral events at |η lab | = 1.8.

Results
The pseudorapidity density as a function of η lab is presented in Fig. 3 for |η lab | < 1.8. An asymmetry between the proton and the lead hemispheres is observed, and the number of charged particles is higher in the Pb-going side (positive η lab ). The ALICE measurement is compared with the pseudorapidity density measured by CMS [37] showing very good agreement within systematic uncertainties, although CMS results exclude prompt leptons. The result is also compared with several models with different descriptions of particle production, all shifted by η lab = 0.465 to take into account the shift to the laboratory system. In the improved HIJING 2.1 [12,13] version the Cronin effect is included, as well as a strong nuclear shadowing effect (sg = 0.28) in order to explain the global properties of the final hadron system in p-Pb collisions [34].  16 TeV in ALICE, with total systematic uncertainties shown as bands, compared with CMS results [37] and theoretical predictions shifted to the laboratory system [12, 14,17,18,20]. The bottom panel shows the ratio to ALICE data The model describes well both the normalisation and the shape of the distribution for the Pb-going side, while it overestimates the p-going side, showing a symmetric behaviour, as for the p-Pb collisions at 5.02 TeV. The dN ch /dη lab versus η lab is compared with two different versions of EPOS. EPOS LHC [17] is a tune of EPOS 1.99 based on LHC data. It is designed to describe all bulk properties of hadronic interactions and based on Gribov-Regge theory for partons. It incorporates collective effects with a separation of the initial state into a core and a corona. EPOS LHC reproduces the Pb-going side, although it underestimates the p-going side of the distribution, showing a stronger asymmetry than data. EPOS 1.99 contains collective flow parametrised at freezeout, while EPOS 3 [14][15][16] includes a full viscous hydrodynamical simulation. It starts from flux tube initial conditions, which are generated in the Gribov-Regge multiple scattering framework. It reproduces the most forward part of the distribution in the Pb-going side, but underestimates both the normalisation, the mid-rapidity part and the p-going side of the dN ch /dη lab distribution. Finally, the distribution is compared with two saturation-based models: MC-rcBK [18,19] and KLN [20,21], which contain a mechanism to limit the number of partons and particles produced. The MC-rcBK results are obtained using the McLerran-Venugopalan model (γ = 1) [59] for the Albacete-Armesto-Milhano-Quiroga-Salgado initial conditions [60]. Saturation-based models are the ones which perform better, underlining the necessity of a mechanism to limit the number of partons produced. Indeed, both MC-rcBK and KLN reproduce the distribution well, within the uncertainties of data, and start to deviate in the region η lab < −1.3. The MC-rcBK model better predicts the p-Pb collisions at 8. 16 TeV than the distribution at 5.02 TeV. The shadowing mechanism used by HIJING is not sufficient to limit the partons produced in the p-going side. Both EPOS and HIJING contain final-state effects, and the performance is worse than for models based on initial-state effects only, like MC-rcBK and KLN. This means that for the dN ch /dη observable final-state effects do not play a role, for the models considered. Nevertheless, all models lie within about 10% when compared with data, and reproduce within systematic uncertainties the Pb-going side.
The charged-particle pseudorapidity density in the laboratory system for |η lab | < 0.5 is dN ch /dη lab = 20.08 ± 0.01 (stat.) ± 0.61 (syst.). In the following, the statistical uncertainty is considered to be negligible. The data are integrated in the range −0.965 < η lab < 0.035 and corrected for the effect of the rapidity shift to retrieve the dN ch /dη in the centre-of-mass system. The correction for the pseudorapidity shift is estimated from HIJING 1.36 [27] to be −3.7% ± 1.9%. The resulting pseudorapidity density in the centre of mass is dN ch /dη = 19.1 ± 0.7.
The charged-particle production is scaled by N part /2, calculated with a Glauber model as explained in Sect. 3, in order to compare the bulk particle production in different collision systems. The number of participants for MB events is 8.09 ± 0.17. The value normalised to the number of participants divided by 2 gives dN ch /dη ×(2/N part ) = 4.73 ± 0.20. In Fig. 4, this quantity is compared with lower energy p-Pb measurements by ALICE [34] as well as by CMS [37] and d-Au measurements at RHIC [38], showing that the values overlap with dN ch /dη measurements

Fig. 4 Values of 2
Npart dN ch /dη for pA [34,37,38], pp and pp [35,[39][40][41][42][43][44][45][46][47] 16 TeV confirms the trend established by lower energy data since the exponent β is not significantly different when the new point is excluded from the fit. The values for p-Pb and d-Au collisions fall on the inelastic pp curve, indicating that the strong rise in AA might not be solely related to the multiple collisions undergone by the participants since the proton in pA collisions also encounters multiple nucleons. As the contribution of diffractive processes to the selected p-Pb sample is negligible, it is expected that the NSD and inelastic selection belong to the same curve for p-Pb, and that this slope corresponds to the one obtained from the inelastic pp curve.
The pseudorapidity density as a function of η lab is presented in Fig. 5 for |η lab | < 1.8 for different centrality intervals, from most central 0-5% to most peripheral 80-100% events. The results for the CL1 estimator have a strong bias due to the complete overlap with the tracking region. V0A has a small multiplicity fluctuation bias due to the enhanced con- It is worth noting that for all the estimators used to select centrality the asymmetry is evident for most central events, while the results for 60-80% and 80-100% classes, where the N part are around 4.5 and 3, respectively, are symmetric.
The left panel of Fig. 6 shows 2 N part dN ch /dη lab as a function of N part for various centrality estimators. For CL1 and V0A the N part from the Glauber model are used and the resulting 2 N part dN ch /dη lab has a steep increase for most central events (higher N part ) due to the strong multiplicity bias discussed in Sect. 3. The rise is steeper for CL1, where the overlap of the centrality selection region with the tracking region is maximal. For the ZNA estimator, two sets of N part are used corresponding to the two different hybrid method selections. For both N mult part and N Pb−side part the trend is similar and extrapolates to the pp point at √ s = 8 TeV. The overall N part dependence of 2 N part dN ch /dη lab for the ZNA estimator is flat and the N part range is more limited when the selection is made in a well separated pseudorapidity region, rather than for multiplicity-based estimators (CL1 and V0A).
A Glauber Monte Carlo calculation based on single quark scattering is also performed [61,62], as it was done for AA collisions [48,49]. Quark constituents are located around the nucleon centre, where the proton density is modelled by a function of the proton radius. To account for effective partonic degrees of freedom, N c = 5 quark constituents have been selected, since this number of constituents was tested for AA collisions and resulted in a constant charged-particle production rate per constituent quark. The effective inelastic cross section for constituent-quark collisions is set to 11.0 mb for 5 constituent quarks to match the 72.5 mb nucleon cross section for p-Pb interactions at 8. 16 TeV [30]. The effective cross sections are constrained using nuclear reaction cross sections [62]. The right panel of Fig. 6 shows the μ N q−part dN ch /dη lab scaled by the average number of participating quarks, μ, in pp collisions, which is 4.44 out of 10 participating quarks for N c = 5, as a function of N part (open points). For the multiplicity-based estimators, CL1 and V0A, there is an increase for the most central and decrease for the most peripheral events with a trend that resembles the one for N part scaling (full points) but with decreased slope. This fact suggests that nuclear-geometrical effects are represented in terms of constituent participant quarks, but not as well as observed for AA collisions [48,49,63], meaning that the multiplicity-fluctuation bias might influence also the quark participants scaling. The μ N q−part dN ch /dη lab has been measured also for 3 constituent quarks, with an inelastic cross section of 22.5 mb and μ = 3.54, showing a distribution in between the N part and N q−part points.

Summary and conclusions
Summarising, the charged-particle pseudorapidity density in |η lab | < 1.8 in NSD p-Pb collisions at √ s NN = 8.16 TeV is presented. A value of dN ch /dη = 19.1 ± 0.7 is measured at mid-rapidity, corresponding to 4.73 ± 0.20 charged particles per unit of pseudorapidity per participant pair, N part /2, calculated with the Glauber model. The new measurement is 9.5% higher than the value at √ s NN = 5.02 TeV. The dependence of dN ch /dη on the centre-of-mass energy is fitted with a power-law function, which gives a much weaker sdependence than for AA collisions. The MB dN ch /dη lab distribution as a function of η lab is compared with CMS results, showing good agreement within uncertainties, and to different models: HIJING 2.1, EPOS (versions LHC and 3) and two saturation-based models, MC-rcBK and KLN. All models can reproduce the data within about 10%, which is a sound achievement given the complexity in describing soft-QCD processes. The best performance comes from saturationbased models, and final-state effects seem not to improve the description of dN ch /dη. Nevertheless, the results provide further constraints for models describing high-energy hadron collisions. The pseudorapidity density for various centrality estimators has been shown and the asymmetry, typical of asymmetric collision systems like p-Pb, is evident for most central events, while results for 60-80% and 80-100% centrality classes are symmetric. The methods to select centrality in p-Pb collisions based on multiplicity measurements have been presented and they induce a multiplicity-fluctuation bias. Results with a selection based on multiplicity estimators at mid-rapidity or within a few units of pseudorapidity and N part from the Glauber model are lower for peripheral values of 2 N part dN ch /dη and higher for most central collisions than the pp value. On the contrary, with centrality selected by the energy deposited in the ZDC, and assuming that the multiplicity in the Pb-going direction is proportional to N Pb−side part , the overall behaviour of 2 N part dN ch /dη as a function of N part is flat, and agrees with the pp measurement at 8 TeV.