Charged-particle multiplicities in proton-proton collisions at $\sqrt{s}$ = 0.9 to 8 TeV

A detailed study of pseudorapidity densities and multiplicity distributions of primary charged particles produced in proton-proton collisions, at $\sqrt{s} =$ 0.9, 2.36, 2.76, 7 and 8 TeV, in the pseudorapidity range $|\eta|<2$, was carried out using the ALICE detector. Measurements were obtained for three event classes: inelastic, non-single diffractive and events with at least one charged particle in the pseudorapidity interval $|\eta|<1$. The use of an improved track-counting algorithm combined with ALICE's measurements of diffractive processes allows a higher precision compared to our previous publications. A KNO scaling study was performed in the pseudorapidity intervals $|\eta|<$ 0.5, 1.0 and 1.5. The data are compared to other experimental results and to models as implemented in Monte Carlo event generators PHOJET and recent tunes of PYTHIA6, PYTHIA8 and EPOS.


Introduction
The multiplicity of emitted charged particles is one of the most basic characteristics of high-energy hadron collisions and has been the subject of longstanding experimental and theoretical studies, which have shaped the understanding of the strong interaction. Following on from earlier ALICE studies of global properties of proton-proton (pp) collisions [1][2][3][4][5][6][7][8], this publication presents a comprehensive set of measurements of the pseudorapidity density ( dN ch / dη) of primary 1 charged particles and of their multiplicity distributions over the energy range covered by the LHC, from 0.9 to 8 TeV. The pseudorapidity density of primary charged particles was studied over the pseudorapidity range |η| < 2, and their multiplicity distributions in three intervals: |η| < 0.5, 1.0 and 1.5. Results are given for three conventional event classes: (a) inelastic (INEL) events, (b) non-single diffractive (NSD) events and (c) events with at least one charged particle in |η| < 1 (INEL>0). At LHC energies, particle production is still dominated by soft processes but receives significant contributions from hard scattering, thus multiplicity and other global event properties measurements allow to explore both components. As these properties are used as input in Glauber inspired models [9][10][11][12], such studies are also contributing to a better modelling of Pb-Pb collisions. Already at 8 TeV, high multiplicity proton-proton collisions provide energy densities comparable, for instance, to energy densities in Au-Au central collisions at RHIC, allowing a comparison of nuclear matter properties in strongly interacting systems with similar energy densities but with volumes orders of magnitude smaller.
It is worth noting that, already at √ s = 2.36 TeV, hadron collision models tuned to pre-LHC data failed to reproduce basic characteristics of proton-proton collisions at the LHC, such as pseudorapidity density of charged particles, multiplicity distributions, particle composition, strangeness content, transverse momentum distributions and sphericity (see for instance [2], [3], [4] and [13]). Therefore, a more precise measurement of charged-particle multiplicity distributions and a study of their energy dependence contribute to a better understanding of particle production mechanisms and serve to improve models. In turn, a better simulation of collision properties improves the determination of the detector response and background estimates of underlying event properties relevant to the study of high-p T phenomena.
In the Regge theory [14], one of the most successful models for describing soft hadronic interactions, the asymptotic behaviour of cross-sections for elastic scattering and multiple production of hadrons is determined by the properties of the Pomeron, the t-channel right-most pole, in the elastic scattering amplitude. In QCD, the Pomeron, which has vacuum quantum numbers, is usually related to gluonic exchanges in the t-channel. The experimentally observed increase of the total cross-section with increasing collision energy made it necessary to consider a Pomeron as a Regge trajectory with t = 0 intercept: α P (0) = 1 + ∆ > 1 [14]. The energy dependence of the particle (pseudo-)rapidity density provides information about the Pomeron trajectory intercept parameter, ∆. If interactions between Pomerons are neglected, the inclusive particle production cross-section, σ Incl. , is determined only by the contribution of the single (cut-)Pomeron exchange diagram. In this approximation, dσ Incl. / dy (∼ dσ Incl. / dη) at midrapidity is proportional to s ∆ [15]. Thus, the energy dependence of the inclusive cross-section gives more reliable information about the value of ∆ than the energy dependence of the total interaction crosssection, for which contributions from multi-Pomeron exchanges strongly modify the energy dependence of the single Pomeron exchange diagram. In the same approximation, the energy dependence of the particle (pseudo-)rapidity density in the central rapidity region is given by dN/ dy ∝ s ∆ σ Int. , where σ Int. is the interaction cross-section (see for instance [16][17][18][19]). Up to LHC energies, σ Int. is well represented by a power law of s. However, for reasons of unitarity [20], it is expected that this power law should be broken at sufficiently high energy, although well above LHC energies. Therefore, the energy dependence of the particle (pseudo-) rapidity density in the central region at LHC, dN/ dy ≈ dN/ dη, should follow the same power law trend. In this publication, this relationship is explored further for three event classes and using 5 ALICE data points.
It was more than 40 years ago that A. M. Polyakov [21] and then Z. Koba, H.B. Nielsen and P. Olesen [22] proposed that the probability distribution of producing n particles in a collision, P (n), when expressed as a function of the average multiplicity, n , should reach an asymptotic shape at sufficiently high energy where Ψ is a function supposed to describe the energy-invariant shape of the multiplicity distribution. Such scaling behaviour is a property of particle multiplicity distributions known today as Koba-Nielsen-Olesen (KNO) scaling.
One well identified mechanism for KNO scaling violation is the increasing probability of multi-parton scattering with increasing √ s. Moreover, since the topologies and multiplicities of diffractive and nondiffractive (ND) events are different, their KNO behavior may be different. Even if KNO scaling were to be valid for each, it might not be valid for their sum. Nevertheless, KNO scaling is expected to be violated for both diffractive and non-diffractive processes [23,24] at sufficiently high collision energies and the LHC provides the best opportunity to study the extent of these scaling violations.
Indeed, deviation from KNO scaling was already observed long ago at ISR energies (proton-proton collisions at √ s from 30.4 to 62.2 GeV), in the full phase space, for inelastic events [25]. On the other hand, for NSD collisions, scaling was still found to be present [25], suggesting that diffractive processes might also play a role in KNO scaling violations. In e + e − collisions, at √ s from 5 to 34 GeV, KNO scaling was found to hold within ±20% [26]. In proton-antiproton collisions at the CERN collider ( √ s = 200, 546 and 900 GeV), KNO scaling was found to be violated for NSD collisions in full phase space [27], [28], [29]. Nevertheless, for NSD collisions, in limited central pseudorapidity intervals, KNO scaling was still found to hold up to 900 GeV, and at √ s = 546 GeV, KNO scaling was found to hold in the pseudorapidity interval |η| < 3.5 [30,31]. In NSD proton-proton collisions at the LHC, at √ s = 2.36 and 7 TeV and in |η| < 0.5, ALICE [2] and CMS [32] observed no significant deviation from KNO scaling.
This publication presents a study of KNO scaling, at √ s from 0.9 to 8 TeV, in three pseudorapidity intervals (|η| < 0.5, 1.0 and 1.5) and for a higher multiplicity reach compared to previous ALICE publications, quantified with KNO variables (moments) [22] as well as with the parameters of Negative Binomial Distributions (NBD) used to fit measured multiplicity distributions.
With respect to previous ALICE publications, the analysis reported here makes use of improved tracking and track-counting algorithms; better knowledge and improved simulation of diffraction processes; an expanded pseudorapidity range for dN ch / dη studies and better statistical precision at √ s = 0.9 and 7 TeV, extending by a factor of 2 the previously published multiplicity distribution reach. Results at √ s = 2.76 and 8 TeV are presented for the first time in this publication.
Previous measurements of both dN ch / dη and multiplicity distributions from CMS [33,34] and UA5 [27] allow a direct comparison to our data. Others by ATLAS [35] and LHCb [36] use different definitions (η and p T ranges) making direct comparison impossible.
This publication is organized as follows: Section 2 describes the ALICE sub-detectors relevant to this study; Section 3 provides the details of the experimental conditions and of the collection of data; Section 4 explains the event selection; Section 5 describes the track selection criteria and the three track counting algorithms; Sections 6 and 7 report the analyses for the measurement of the pseudorapidity density and of multiplicity distributions, respectively; Section 8 discusses systematic uncertainties; Section 9 presents the multiplicity measurements, NBD fits of the multiplicity distributions, KNO scaling and q-moment studies. Finally, in Section 10, the results are summarized and conclusions are given.

ALICE subdetectors
The ALICE detector is fully described in [37]. Only the main properties of subdetectors used in this analysis are summarized here. Charged-particle tracking and momentum measurement are based on data recorded with the Inner Tracking System (ITS) combined with the Time Projection Chamber (TPC) [38], all located in the central barrel of the ALICE detector and operated inside a large solenoid magnet providing a uniform 0.5 T magnetic field parallel to the beam line.
The V0 detector [39] consists of two scintillator hodoscopes, each one placed at either side of the interaction region, at z = 3.3 m (V0A) and at z = −0.9 m (V0C) (z is the coordinate along the beam line, with its origin at the centre of the ALICE barrel detectors), covering the pseudorapidity ranges 2.8 < η < 5.1 and −3.7 < η < −1.7, respectively. The time resolution of each hodoscope is better than 0.5 ns.
The ITS is composed of high resolution silicon tracking detectors, arranged in six cylindrical layers at radial distances to the beam line from 3.9 to 43 cm. Three different technologies are employed. For the two innermost layers, silicon pixels (SPD [40]) are used, covering pseudorapidity ranges |η| < 2 and |η| < 1.4, respectively. The SPD is followed by two Silicon Drift Detector layers (SDD, [41]). The Silicon Strip Detector (SSD, [42]) constitutes the two outmost layers consisting of double-sided silicon microstrip sensors. The intrinsic spatial resolution (σ rϕ × σ z ) of the ITS subdetectors is: 12×100 µm 2 for SPD, 35×25 µm 2 for SDD, and 20×830 µm 2 for SSD, where ϕ is the azimuthal angle and r the distance to the beam line. The ITS sensors were aligned using survey measurements, cosmic muons and collision data [43]. The estimated alignment accuracy is 8 µm for SPD and 15 µm for SSD in the most precise coordinate (rϕ). For the SDD, the intrinsic space point resolution is σ z = 30 µm in the z direction and σ rϕ = 40 to 60 µm, depending on the sensor, along rϕ (drift). Because of some anomalous drift field distributions, in the reconstruction, a systematic uncertainty up to 50 µm in z and 500 µm in rϕ was added to account for differences between data and simulation. The ITS resolution in the determination of the transverse impact parameter measured with respect to the primary vertex is typically 70 µm for tracks with p T = 1 GeV/c, including the contribution from the primary vertex position resolution.
The SPD and the V0 scintillator hodoscopes provided triggers for collecting data.
The TPC [38] is a large cylindrical drift detector with a central high voltage membrane at z = 0, maintained at +100 kV and two readout planes at the end-caps. The material budget between the interaction point and the active volume of the TPC corresponds to 11% of a radiation length, when averaged over |η| < 0.8.
The TPC and the ITS were aligned relative to each other within a few hundred micrometers using cosmicray and proton collision data [43].
The momentum measurement is not explicitly used in this study, however, the simulation of the detector response is sensitive to the particle momentum spectrum. Since event generators used in Monte Carlo simulations do not reproduce the observed momentum distributions, the difference between data and Monte Carlo simulation is taken into account when evaluating systematic errors. For momenta lower than 2 GeV/c, representing the bulk of the data, the p T resolution for tracks measured in the TPC and in the ITS, is about 0.80% at p T = 1 GeV/c, it increases to 0.85% at p T = 2 GeV/c and to 3% at p T = 0.1 GeV/c. Charged-particle multiplicities were measured using information from the TPC in |η| < 0.9 and from the ITS in |η| < 1.3. At larger pseudorapidities, the SPD alone was used to expand the range of dN ch / dη measurements to |η| < 2.0.

Proton beam characteristics
Data were selected during LHC collision periods at a luminosity low enough to allow the minimum bias trigger rate not to exceed 1 kHz. At √ s = 0.9 TeV, the number of protons per colliding bunch varied from 9×10 9 to 3.4×10 11 , while the number of colliding bunches was either 1 or 8. At √ s = 2.76 TeV, the number of protons per colliding bunch varied from 5×10 12 to 7×10 12 , while the number of colliding bunches was either 48 or 64. At √ s = 7 TeV, the number of protons per colliding bunch varied from 8.6×10 9 to 1.4×10 12 , resulting in a luminosity between 10 27 and 10 30 cm −2 s −1 . There were up to 36 bunches per beam colliding at the ALICE interaction point. When needed, the luminosity was kept below 10 30 cm −2 s −1 by a transverse displacement of the beams with respect to one another. At √ s = 8 TeV, there were 3 proton bunches colliding at the ALICE interaction point each containing about 1.6×10 11 protons.
Data used for this study were collected at low beam currents, so that beam-induced backgrounds (beamgas or beam-halo events) were low and could be removed offline using V0 and SPD detector information, as discussed in Section 4.1.

Triggers
The ALICE trigger system is described in [44]. Data were collected with a minimum bias trigger, MB OR , requiring a hit in the SPD or in either one of the V0 hodoscopes; i.e. essentially at least one charged particle anywhere in the 8 units of pseudorapidity covered by these detectors. Triggers were required to be in time coincidence with a bunch crossing the ALICE interaction point. Control triggers, taken for various combinations of beam and empty-beam buckets, were used to measure beam-induced and accidental backgrounds.

Characteristics of data samples used in this study
General characteristics of the data samples used are given in Table 1.
The data at √ s = 0.9 TeV were collected in May 2010, with one polarity of the ALICE solenoid magnet (solenoid magnet field pointing in the positive z direction).
The first LHC data above Tevatron energy were collected in 2009, at √ s = 2.36 TeV, in a run with unstable LHC beams, during which only the SPD was turned on. Therefore, in this case, the charged-particle multiplicity was measured using exclusively the SPD information. In this publication, the previously published results at √ s = 2.36 TeV [2] are used for comparison.
Proton-proton data were collected at √ s = 2.76 TeV, an energy that matches the nucleon-nucleon centreof-mass energy in the first Pb-Pb collisions provided by the LHC, in 2011.

Data at
√ s = 7 TeV were collected in 2010. About 20% of the data were taken with a magnet polarity opposite (solenoid field pointing in the negative z direction) to that of √ s = 0.9 TeV data. A sample of 12.3×10 6 events, collected without magnetic field, was used to check some of the systematic biases in track reconstruction.

At
√ s = 8 TeV only a subset of runs was collected with the MB OR as a minimum bias trigger in 2012, 10 were selected for this analysis. At 0.9 and 7 TeV, data samples are substantially larger than those available in previous ALICE publications on charged-particle multiplicities [1][2][3]. For the charged-particle multiplicity analysis, the event sample at √ s = 0.9 and 7 TeV increased by a factor of 50 and 2000, respectively, giving significant extension of the multiplicity reach and better statistical precision. The precision of dN ch / dη is not substantially limited by event sample size. However, the large number of runs available made it possible to The main sources of event background are beam gas and beam halo collisions. Such events were removed by requiring that the timing signals from the V0 hodoscopes, if present, be compatible with the arrival time of particles produced in collision events. In addition, because of the different topology of beam background events, the ratio between the number of SPD clusters and the number of SPD tracklets 2 is much higher in beam background events, therefore a cut on this ratio was applied. The remaining fraction of beam background events in the data, estimated by analysing special triggers taken with non-colliding bunches or empty beam buckets, does not exceed 10 −4 for all centre-of-mass energies. The track beam background is mostly significant in the last η bins (|η| ≈ 2) where it reaches 4×10 −3 in the worst case.

Event pileup
The other type of potential event background comes from multiple collision overlap. For the data used in this publication, the proton bunch spacing was 50 ns or longer, the luminosity did not exceed 10 30 cm −2 s −1 , and the probability to have collisions from different bunch crossings in the 300 ns integration time of the SPD was negligible. However, multiple collisions in the same bunch crossing, also referred to as event pileup or overlap, have to be considered in case their vertices are not distinguishable. In order to avoid or minimize corrections for event pileup, runs with a low number of interactions per bunch crossing, µ ≤ 0.061, were selected resulting in an average µ, µ ≤ 0.04, for all data samples (Table 1). This corresponds to at most 2% probability of more than one interaction per event.
The identification of pileup events relies on multiple vertex reconstruction in the SPD, with algorithms using three basic parameters: (a) The distance of closest approach (DCA) to the main vertex for a SPD tracklet to be included in the search for an additional interaction: DCA > 1 mm; (b) The distance between an additional vertex and the main vertex, ∆z > 8 mm; (c) The number of SPD tracklets (N trk ) used to determine an additional vertex (number of contributors to the vertex): N trk ≥ 3.
With this choice of parameters, and with the relatively broad z vertex distribution at the LHC (FWHM ≥ 12 cm), typically only 10% to 15% of multiple collisions are missed, and the fraction of fake multiple collisions due to SPD vertex splitting from a single interaction is low (typically a few times 10 −5 ).
The pileup detection efficiency was studied both by overlapping two Monte Carlo proton-proton collisions and by measuring pileup in the data. The pileup fraction, estimated from identified pileup events in the data, is found to be consistent with what is expected from the µ values derived from trigger information (Table 1).
In multiplicity measurements, pileup affects the data mainly when two vertices are not distinguishable. When they are distinguishable, the multiplicity is taken from the vertex with the highest number of tracks. The small bias induced by choosing systematically the highest multiplicity vertex is negligible in our low pileup data samples.
Comparing dN ch / dη measurements, for different runs, no correlation is found between dN ch / dη values at η = 0 and µ values. Comparing data with and without identified pileup rejection, the change in dN ch / dη values is smaller than 0.5%, which is smaller than systematic uncertainties. Note that the requirements for track association to the main vertex reject a further fraction of the tracks coming from the 10% to 15% of unidentified pileup collisions. The conclusion is that event pileup corrections to dN ch / dη are negligible in these low pileup data samples.
For multiplicity distributions, even though data were selected with a low pileup probability, it is important to verify that the pileup does not distort the distributions, as the relative pileup fraction increases with multiplicity. The fraction of pileup events, which the ALICE pileup detection algorithm identifies after the event selection, is about 10 −2 , with no significant differences between the four centre-of-mass energies. Moreover, tight DCA cuts allow tracks originating from the main vertex to be distinguished from those coming from a pileup vertex even when the vertices are closer than 0.8 cm in z. This was confirmed by simulating events, where two Monte Carlo pp collisions were superimposed, demonstrating that only 5% of the events passing the selection had extra tracks from the secondary vertex. In 90% of such cases, the distance along the beam line between the two vertices was ∆z < 0.5 cm. In the data samples with a pileup fraction of order µ/ 2 ≤ 0.02, the residual average fractions of events with pileup is at most 0.4%. Furthermore, the simulation shows that the pileup that does affect the multiplicity of an event is rather broadly distributed across events with different multiplicity, but becomes significant only outside the multiplicity range studied here. The multiplicity at which the pileup contribution reaches 10% of the measured multiplicity at √ s = 7 TeV is N ch = 105, 170 and 310, for |η| < 0.5, 1.0, and 1.5, respectively, which is beyond multiplicity ranges covered in this publication. Therefore, no pileup corrections were applied. Other background contributions from cosmic muons or electronics noise are also negligible.

Offline trigger requirement
Both for the INEL and INEL>0 normalizations, the online MB OR trigger was used. However, for the NSD analysis, a subset of the total sample was selected offline by requiring a coincidence (MB AND ) between the two V0 hodoscope arrays. This corresponds to the detection of at least one charged particle in both hemispheres, in the V0 hodoscope arrays separated by 4.5 units of pseudorapidity, a topology that tends to suppress single-diffraction (SD) events; therefore, model dependent corrections and associated systematic errors are minimized.

Vertex requirement
The position of the interaction vertex is obtained either by correlating hits in the two silicon-pixel layers (SPD vertex), or from the distribution of the impact parameters of reconstructed global tracks 3 (global track vertex) [37,38,43,45,46]. The next step in the event selection consists of requiring the existence of a reconstructed vertex.  Two SPD vertex algorithms were used: a three-dimensional vertexer (3D-vertexer) that reconstructs the x, y and z positions of the vertex, or a one-dimensional vertexer (1D-vertexer) that reconstructs the z position of the vertex. The vertex position resolution achieved depends on the track multiplicity. For the 3D-vertexer it is typically 0.3 mm both in the longitudinal (z) direction and in the plane perpendicular to the beam direction. The 1D-vertexer resolution in the z direction is on average 30 µm. If the 3D-vertexer algorithm does not find a vertex (typically 47% of the cases at √ s = 7 TeV), then the simpler 1D-vertexer is used to determine the z position of the vertex, and the x and y coordinates are taken from the average x and y vertex positions of the run. The 3D-vertexer efficiency is strongly multiplicity dependent. As the bulk of the events have a low multiplicity, this explains the relatively low average vertex finding efficiency. For the z coordinate, if no reliable vertex is found (typically 14% of the cases), either because the 1D-vertexer did not find a vertex or the 1D-vertex quality was not sufficient (the dispersion of the difference of azimuthal angles between the two hits, one in each SPD layer, of tracklets contributing to the vertex is required to be smaller than 0.02 rad), the event is rejected. For the global track vertex, the resolution is typically 0.1 mm in the longitudinal (z) direction and 0.05 mm in the direction transverse to the beam line.
Both SPD and global track vertices have to be present and consistent by requiring that the difference between the two z positions be smaller than 0.5 cm. If not, in 3 to 4% of the cases, the event is rejected. The cut was chosen to be compatible with DCA z cut applied to tracks to ensure that we combine tracklets and tracks from the same collision (see Section 5.2). This condition removes mainly non-Gaussian tails in the columns of the detector response matrix 4 at low multiplicity, coming from the fact that SPD and track vertices, when separated, tend to have different multiplicities associated to them. In the data, this requirement also removes 80% of pileup events with well-separated vertices.
The η acceptance is correlated with the vertex z position (z vtx ) (Fig. 1). For multiplicity distribution measurements, in order for tracks to remain within the acceptance of the SPD in the η versus z vtx plane, the following requirements were imposed on the vertex position along the z axis: |z vtx | < 10, 5.5 and 1.5 cm for |η| < 0.5, 1 and 1.5, respectively. In the measurement of dN ch / dη, the requirement on the vertex was relaxed to |z vtx | < 30 cm, in order to allow extending the η range to |η| < 2.

Event selection efficiency
As described in [47], PYTHIA6 [48][49][50] and PHOJET [51,52] event generators used by ALICE were adjusted to reproduce the measured diffraction cross-sections and the shapes of the diffracted mass (M X ) distributions extracted from the Kaidalov-Poghosyan model [53]. These modified versions of event generators are referred to as "tuned for diffraction". Typically, σ SD / σ INEL ≈ 0.20, where σ INEL is the inelastic cross-section, σ SD is the SD cross-section for M X < 200 GeV/c 2 , and σ DD / σ INEL ≈ 0.11, where σ DD is the double diffraction cross-section for ∆η > 3 (∆η is the size of the particle gap in the pseudorapidity distribution). These fractions have insignificant energy dependence between 0.9 and 7 TeV [47], and the values at 7 TeV were used for 8 TeV data. Table 1 shows the number of events selected at each centre-of-mass energy prior to the z vtx requirement. Selection efficiencies using criteria defined above in this section, were estimated for INEL, NSD and SD events (classified at generator level by event generator flags) as a function of the number of generated charged particles (shown on Fig. 2 for the case |η| < 1 and the various centre-of-mass energies considered). The particular selection is designated by the offline trigger used to construct it, MB OR or MB AND . Note that for dN ch / dη measurement selection efficiencies are defined in a separate way (see Section 6.1). At √ s ≥ 7 TeV the INEL event selection efficiency based on the MB OR trigger reaches 100% for a charged-particle multiplicity above 8.
For SD events, the efficiency of the MB AND selection reduces significantly when going to higher energies ( Fig. 2), because the Lorentz boost of the diffracted system increases with increasing centre-of-mass energies. This implies that in the normalization to the NSD event class, corrections for the remaining SD contribution become smaller when going to higher energies. The MB AND trigger selects 84%, 86%, 87% and 87% of the MB OR triggers, and 13%, 4%, 1% and 1% of the SD events satisfy the MB AND selection, at √ s = 0.9, 2.76, 7 and 8 TeV, respectively.
5 Track selection and multiplicity algorithms 5.1 Track quality requirements The following criteria were used to select reconstructed tracks associated to the main event vertex: -for tracks reconstructed from both ITS and TPC information (global tracks), the selection requires at least 70 pad hit clusters in the TPC, a good track quality ( χ 2 dof < 4), a distance of closest approach (DCA) along the z direction (DCA z ) < 0.5 cm, and a p T -dependent transverse DCA (DCA T ) requirement, which corresponds to a 7 sigma selection. DCA T conditions are relaxed by a factor 1.5 for tracks lacking SPD hits.
-for tracks reconstructed with ITS information only (ITS-only tracks) the number of ITS hit clusters associated to the track must be larger than 3, among the 6 layers of the ITS, and χ 2 dof < 2.5. The DCA z and DCA T requirements are the same as for global tracks.
-for SPD tracklets, the association to the vertex is ensured through a χ 2 requirement. Using the SPD vertex as the origin, differences in azimuthal (∆ϕ = ϕ 2 −ϕ 1 , bending plane) and polar (∆θ = θ 2 −θ 1 , nonbending direction) angles are calculated between hits in the inner (layer 1) and in the outer (layer 2) SPD layers. Hit combinations, called tracklets, are selected with the following condition where σ ϕ = 0.08 rad, σ θ = 0.025 rad and the sin 2 factor takes into account the θ dependence of ∆θ. reproduced by the simulation. The cut imposed on the difference in azimuthal angles rejects charged particles with a transverse momentum below 30 MeV/c; however, the effective transverse-momentum cut-off is determined mostly by particle absorption in the material and is approximately 50 MeV/c, in |η| < 1. If more than one hit in an SPD layer matches a hit in the other layer, only the hit combination with the smallest χ 2 value is used.
Some of the SPD elements had to be turned off, resulting in lower efficiency in some regions of the η versus azimuthal angle plane. In order to reach the best possible precision in the measurement of dN ch / dη, fiducial cuts were applied to both tracks and tracklets, excluding azimuthal regions where the tracking efficiency corrections are relatively large. These fiducial cuts vary with data taking periods, following the evolution of the SPD acceptance. At √ s = 0.9, 2.76, 7, and 8 TeV, the fractions of the acceptance removed were 64%, 68%, 65%, and 35%, respectively. Some of the SPD elements could be recovered before collecting 8 TeV data, explaining the improvement.
For multiplicity distribution studies, fiducial cuts were not applied because they increase statistical uncertainty, hence limiting the high multiplicity reach.

Track counting algorithms
In previous ALICE publications [1][2][3], the charged-particle multiplicity was measured in |η| < 1.3 using only SPD tracklets built from SPD pixel hits. In order to extend the pseudorapidity range to |η| < 2, an improved tracklet algorithm, initially used in [54], was introduced to take into account the θ dependence of the uncertainty in the χ 2 (Eq. (2)). With this improvement, the efficiency for detecting SPD tracklets became uniform as a function of pseudorapidity and z position of the vertex, which allowed vertices further away from the nominal interaction point along the beam direction to be used, thereby extending significantly the pseudorapidity range.
To be less sensitive to the SPD acceptance, track counting algorithms were developed, that make use of tracking information from other ALICE detectors, the SDD, the SSD and the TPC. Each track is counted as primary if it fulfills the transverse DCA requirements listed in Section 5.1 and it is not associated to a secondary vertex identified by a dedicated algorithm [45] tuned to tag γ-conversions, K 0 and Λ decays.
Three multiplicity estimators were developed by ALICE using three different samples of tracks: -SPD tracklets, with |η| < 2 (referred to as Tracklet algorithm) 5 . The Tracklet algorithm stores, for each tracklet, references to ITS or global track candidates using at least one of its pixel clusters.
-ITS-only tracks, with |η| < 1.3, obtained using all hit clusters in this detector, plus tracklets (|η| < 2) built out of SPD pixel clusters not matched to any ITS track (referred to as ITS+ algorithm).
-TPC tracks, with |η| < 0.9, matched to hits in the ITS, plus ITS-only tracks (up to |η| < 1.3) built out of silicon hit clusters not matched to any TPC track, plus tracklets (|η| < 2) built out of SPD pixel clusters not matched to any ITS or TPC track (referred to as ITSTPC+ algorithm).   [48] combined with a simulation of the ALICE detector, at √ s = 7 TeV, for three pseudorapidity intervals (|η| < 0.5, 1.0, and 1.5 from left to right, respectively), and for the three track counting algorithms, Tracklet, ITS+ and ITSTPC+, from top to bottom, respectively. Horizontal axes show generated primary charged-particle multiplicities and vertical axes measured multiplicities.
In order to keep away from the edges of the detectors, where the acceptance is less precisely known, ITS and TPC tracks used in this study are limited to |η| < 1.3 and |η| < 0.9, respectively. 5 Potentially |η| 3 can be reached using event vertices displaced from the detector center at distances |z vtx | 30 cm (see Fig. 1), however the sample size of such events is too small.   Properties of the three track counting algorithms are compared in Fig. 3, showing that, going from Tracklet to ITS+ and to ITSTPC+ algorithms, the detector response matrix becomes narrower and has a topology closer to that of a diagonal matrix. When going from |η| < 0.5 to |η| < 1.5, the response matrix becomes broader and has a less diagonal topology, as geometrical acceptance effects become more important, and dominated by the SPD with significant inefficiency due to some missing modules. Note that by restricting the azimuthal angle to good regions of the SPD, the difference between algorithms in dN ch / dη measurements is of order ±1% in the central region (Fig. 4). However, the result with the Tracklet algorithm is not sensitive to this cut, and, as it is needed to measure multiplicities beyond |η| = 1.3, it is used alone for dN ch / dη measurement. For multiplicity distribution measurements all three algorithms are used without the ϕ region restrictions with a corresponding systematic uncertainty contribution.
In the pseudorapidity region |η| < 0.9, the TPC accounts for 90% of the tracks, the ITS complement 9% and the SPD complement 1%. These fractions vary with the η range. Outside |η| < 1.3, SPD tracklets are the only contribution. The small fluctuations between points in |η| > 1.3 come from the slightly different number of events used for averaging between algorithms, after efficiency corrections in each η bin.
6 Pseudorapidity density of primary charged particles: analysis Raw dN ch / dη distributions have to be corrected for detector and trigger acceptance and efficiency, and for contamination from daughters of strange particles. Note that this section only describes the particularities of dN ch / dη measurement, unless specifically stated otherwise. For charged particle multiplicity distribution measurement see Section 7.

Acceptance and efficiency corrections
Three types of corrections have to be applied to the raw data: (a) a track-to-particle correction to take into account the difference between measured tracks and "true" charged primary particles. This correction mainly depends on acceptance effects and on detector and reconstruction efficiencies; (b) corrections for the bias coming from the vertex reconstruction requirement, at both track and event levels (vertex reconstruction correction). This bias exists on both the number of tracks and the events used, since events without a reconstructed vertex are not selected, and tracks from those events therefore do not contribute; (c) corrections at both track and event levels, to take into account the bias due to the MB OR trigger required for INEL and INEL>0 event classes or the MB AND offline selection for the NSD event class.
In practice, the number of tracks is corrected as a function of η and z vtx and the number of events is corrected as a function of reconstructed track multiplicity and z vtx . The number of events without trigger or without reconstructed vertex is estimated from the simulation and included in the corrected number of events. Finally, the quantity dN ch / dη, averaged over all events, is obtained for each η bin. The range of z vtx contributing to the multiplicity varies with η ( Fig. 1). For instance, at η = 2, tracks originate mostly from vertices in the range: -30 cm < z vtx < -5 cm. Therefore, for each η bin, a z vtx acceptance correction is applied. See [55] for details of the procedure.

Strangeness correction
Since ALICE's definition of primary charged particles excludes particles originating from the weak decays of strange particles, data have to be corrected for cases when daughter particles of such decays pass the track selection. Current Monte Carlo event generators have a strangeness content which differs from data by a factor approaching 2. Therefore, the strangeness content in the Monte Carlo simulation was normalized to data using ALICE's K 0 and Λ measurements in |η| < 0.9 [4], that were extrapolated to |η| ≤ 2 using the shape from simulation. The ratios of strangeness contents between data and Monte Carlo generators are slightly centre-of-mass energy dependent. For √ s varying from 0.9 to 8 TeV they increase from 1.6 to 1.85 according to PYTHIA6, and from 1.4 to 1.6 according to PHOJET. The uncertainty on these ratios coming from the uncertainty in ALICE measurements of strange particle production [4], is estimated to be 5%. The strangeness contamination is slightly η dependent, and varies from 1.7% at η = 0 to 2.5% at η = 2 at √ s = 7 TeV. The strangeness correction is about 1%, has no significant η variation in |η| < 2 and no significant energy dependence between √ s = 0.9 and 8 TeV. This correction is explained in more detail in the Section 8.1.3, where the corresponding systematic uncertainty is discussed.

Event class normalization
The final correction applied to the data is the normalization to one of the three event classes defined in this study: NSD, INEL and INEL>0. In the normalization to NSD, corrections have to be made for the fraction of SD events remaining in the selection and for the fraction of double-diffraction (DD) events not included in the selection. In the normalization of results to the INEL event class, corrections have to be made for the fraction of single-and double-diffractive events not included in the selection. The INEL>0 class is of interest because it minimizes diffractive corrections. In addition, ALICE measurements of SD and DD cross-sections [47] reduced the systematic uncertainties coming from diffraction. Corrections for higher order diffractive processes associated with events with two or more pseudorapidity gaps (regions devoid of particles) are neglected in the normalization to INEL, NSD and INEL>0 classes, as their contribution to inelastic collisions is expected to be smaller than 1% [14,53]. Furthermore, such events tend to have a high trigger efficiency, which makes corresponding corrections even smaller.   Table 2: MB AND trigger efficiencies for NSD events and MB OR trigger efficiencies for inelastic events at four centre-of-mass energies, obtained from diffraction-tuned versions of PYTHIA6 Perugia0 [48] and PHOJET [51]. Uncertainties listed are total uncertainties. Statistical errors are negligible. The asymmetry of the MB OR errors is due to the asymmetric uncertainties in the diffraction efficiencies.
To normalize measurements to a given event class, trigger biases must be corrected for, both at event and track levels. For the INEL and INEL>0 classes, the correction is straightforward using the MB OR trigger efficiency (Table 2).
For the NSD event class, contamination of the event sample by SD events must be taken into account. The measured quantity may be re-written as: ∝ ε NSD MB AND σ NSD , where ε and σ are efficiencies and cross-sections, respectively, for SD or NSD events [47], one obtains: The coefficient in front of the single diffraction term in Eq. (4), varies from 0.04 at √ s = 0.9 TeV to 0.003 at √ s = 8 TeV. As the single diffraction term is not measured, but corresponds to a relatively small correction, this term was calculated using the simulation. The corresponding uncertainty was estimated by varying the single diffraction term conservatively between extreme cases, assuming either no SD, or assuming that all events are from SD. The last step consists of correcting for the MB AND trigger efficiency to obtain the desired quantity, The DD event content of the MB OR and MB AND data samples, is small, of the order of 5.5% and 4.5%, respectively. These fractions do not vary significantly between 0.9 and 8 TeV. The corrections for DD efficiency are included in the general efficiency correction. For the INEL and INEL>0 event classes, the MB OR trigger efficiency for DD events as a function of multiplicity is the same as for the other inelastic events. The MB AND selection, which is used for the NSD event sample, has an efficiency for DD events that is lower than that of the other inelastic events. However, we checked in the simulation that the average efficiency correction for the NSD event class gives the same result as separate efficiency corrections implemented for DD and ND events. The data samples used in these measurements are described in Table 1. The next step in the analysis consists of correcting the raw distributions for detector acceptance and efficiencies, using an unfolding method.
The unfolding procedure follows the same approach as in [2], i.e. the corrected distribution is constructed by finding the vector U, which minimizes a χ 2 given by where M represents the raw multiplicity distribution vector with uncertainty vector s, U the unfolded multiplicity distribution vector, and R the detector response matrix. Indices m and t run from 0 to the maximum number of multiplicity bins, in raw and corrected distributions respectively. The regularization term β × F(U) is used to decrease the sensitivity of the unfolding to statistical fluctuations. For F(U) a usual Tikhonov-type of function [56] was used, which has a smoothing effect on the unfolded distribution where N is the number of unfolded multiplicity bins, evaluated with the help of the response matrix, from the maximum raw multiplicity.
The weight β (Table 3) was chosen to minimize the mean squared error [56]. The solution is found to be stable over a broad range of β values (±50%), and the correct minimum was ensured in each case by scanning β over few orders of magnitude. The particular values of optimal weights depend on many features of the unfolding problem, such as distribution size, a pattern of fluctuations in the input raw data, properties of the response matrix and the regularization term. The most obvious dependence was eliminated by factorizing N in Eq. (6).  Table 3: Values of the weight parameter β used in the regularization term (Eq. (5)), for each centre-of-mass energy and for each pseudorapidity range.
For each generated multiplicity bin N gen = t, the response matrix column R mt consists of the distribution of the probability to measure multiplicity N ch = m. To extend the response matrix to the highest multiplicities encountered in this study, beyond the reach of the available simulation, probability distributions were parameterized and extrapolated towards high multiplicities (Fig. 5). In the low-N gen region (N gen < 10 to 20, depending on the η range) the response matrix was taken directly from the simulation. In the large N gen region (N gen ≥ 10 to 20), the column R mt is well described by a Gaussian distribution and mean values follow a linear trend (Fig. 5). Widths were parameterized using two different functions, a Padé function and a power law C 0 , C 1 , C 2 , C 3 and γ are constants to be fitted. These functions have different asymptotic behaviours (Fig. 5), however, using either function makes a difference only for multiplicities above 100 (in |η| < 1.5).
The switch to parameterization occurs at N gen = 10, 15 and 20, for |η| < 0.5, 1 and 1.5, respectively, for all energies. These values ensure that using the parameterized response matrix introduces no distortions in the low multiplicity region. The range of multiplicities in the final unfolded distribution was further restricted by requiring that the bias (an estimate of how far is the result from the true solution [56]) is less than 10% in each bin. As unfolding is performed for each correction scenario (see Section 8.2 on systematic uncertainties), in the end the multiplicity range is limited by the unfolding resulting in the shortest range. The quality of the unfolding was verified by comparing the raw distribution M to the product R ⊗ U.

Event class normalization
After the unfolding step, distributions have to be corrected for event selection efficiency (including trigger efficiency and vertex reconstruction efficiency).
For the INEL and INEL>0 event classes this is straightforward, given that where ε t is the selection efficiency for true multiplicity t. Thus the unfolded distribution can be normalized to the inelastic event class by dividing the contents of each multiplicity bin by the corresponding efficiency.
For the NSD event class the procedure is different as the unfolded distribution, U * , includes a contamination by SD events where upper indices denote the event class. The overall fraction of SD in inelastic collisions (α t ) was measured by the ALICE Collaboration [47] The desired unfolded distribution normalized to the NSD event class is obtained by combining Eqs. (10) and (11) 8 Study of systematic uncertainties 8.1 Common sources of systematic uncertainties

Material budget
The material budget in the ALICE central barrel was checked in the range |η| < 0.9, by comparing measured and simulated gamma conversion maps. The conclusion is that, in this pseudorapidity range, the material budget is known with a precision of 5% [46]. The corresponding systematic uncertainty was obtained by varying the material budget in the simulation, conservatively over the whole pseudorapidity range by ±10%, which induces a variation of dN ch / dη of ±0.2% for all event classes. For multiplicity distributions, the η range considered does not exceed ±1.5, making the effect of the higher material budget uncertainty outside |η| < 0.9 relatively small. The systematic uncertainty from material budget is negligible compared to other sources of uncertainty.

Magnetic field
The magnetic field map was measured and simulated with finite precision. To check the sensitivity of the detector response to the precision of the simulation of the ALICE solenoidal magnetic field, data samples collected at √ s = 7 TeV with opposite polarities were compared. The differences in dN ch / dη values are consistent with observed fluctuations between runs within the data-taking period at this energy. Therefore, the contribution from systematic uncertainties associated with the magnetic field are smaller than, or of the order of, the run-to-run fluctuations, and have been neglected.

Strangeness correction and particle composition
The main sources of uncertainty associated to the correction for strange particles originate from: (a) the difference (≤ 5%) in K 0 and Λ detection efficiency between data and simulation [4]; (b) the difference in p T distributions of strange particles in data compared to simulation, which implies a difference in the fractions of daughter particles meeting the vertex association condition; and (c) the uncertainty in the simulation of the strange particle content.
For dN ch / dη measurements, the systematic uncertainty from strangeness correction is found to have a small η variation, it is slightly larger at η = 0 compared to |η| = 2. These uncertainties have a small energy dependence, they increase slightly with increasing √ s, from 0.14% at √ s = 0.9 TeV to 0.16% at √ s = 8 TeV. The uncertainties at η = 0 are listed in Table 4.
For multiplicity distributions, the strangeness contamination was studied with the Monte-Carlo simulation by evaluating the survival probability of strange particle decay products for the track selection used in the analysis. The probability that a track from strange particle decay passes the track requirements is less than 0.1% on average, leading to a negligible contribution to the uncertainty on multiplicity distributions.
The particle composition affects the efficiency estimate, because different particle species have different efficiencies and effective p T cut off. The influence of the uncertainty in particle composition was estimated by varying, in the simulation, the relative fractions of charged kaons and protons with respect to charged pions by ±30%, which covers conservatively the uncertainties in the measured particle composition at the LHC [57], and found to range, for dN ch / dη, between 0.1% for the INEL events class and 0.2% for the NSD and INEL>0 event classes ( Table 4). The effect is negligible for multiplicity distributions.

Detector Simulation
The systematic uncertainty related to the limited precision with which the detector performance is simulated was evaluated by varying the threshold on parameters used to select the various types of tracks,  Table 4: Relative contributions, in percent, to systematic uncertainties in the measurement of dN ch / dη, at η = 0, for centre-of-mass energies considered in this study and for the three event classes defined in the text; the p T uncertainties combine the contributions from the difference between MC models and data with the uncertainty on the p T distribution below 50 MeV/c. Run-to-run fluctuation contributions indicated here for comparison are not included in total uncertainties. over a range obtained from the observed difference in the distributions of these parameters between simulation and data. For dN ch / dη measurements based on tracklets, the χ 2 cut was varied between 1.3 and 4. The spread of dN ch / dη values over the range of z-positions of the vertex covered by a given η bin was used as a measure of the bias introduced by the z-dependence of the tracklet reconstruction efficiency. The corresponding uncertainty ranges from 0.6% to 1.1% depending on the data sample (Table 4).
Since for multiplicity distribution measurements all three track-counting algorithms are used, track parameters were also varied for global tracks, and for ITS-only tracks (details may be found in [46]): -for global tracks, the minimum number of TPC clusters was varied between 60 and 80 and the χ 2 cut between 3 and 5.
-for ITS-only tracks, the χ 2 cut was varied between 2 and 3.
The uncertainty on the DCA distribution is not considered here, as it is included implicitly in the uncertainty on the strangeness content, which affects the DCA distribution.
We find that, in the case of multiplicity distributions, uncertainties in the detector simulation are negligi-ble compared to other sources of uncertainties.

Model dependence
The remaining SD fraction in the sample selected with the MB AND trigger in view of the normalization to the NSD event class is 3%, 1% and negligible at √ s = 0.9, 2.76 and ≥ 7 TeV, respectively. Uncertainties coming from diffraction contributions are included in the trigger efficiency uncertainties (Table 2) obtained in [47], for √ s ≤ 7 TeV. At √ s = 8 TeV, the efficiency values were taken to be the same as for √ s = 7 TeV. In addition, the model uncertainties at √ s = 0.9, 2.76, 7 and 8 TeV were obtained from the difference between PYTHIA6 Perugia0 and PHOJET. A test of the efficiency evaluation was obtained by comparing simulated MB AND to MB OR trigger efficiency ratios to the measured values. Excellent agreement was found at all energies [47].
For multiplicity distributions, the systematic uncertainty from the model dependence is included in both efficiency correction and p T dependence uncertainties, as different event generators and tunes are used to estimate independently efficiencies and response matrices (see Section 8.2).

p T dependence
None of the MC generators used in the detector simulation reproduces correctly the p T distribution of charged particles observed in the data [58,59]. This introduces an uncertainty in the determination of the detector response, as it is integrated over transverse momenta, and the probability of detecting a particle decreases with decreasing p T . This affects in particular, together with the uncertainty on the material budget and the magnetic field, extrapolations of measurements to p T = 0.
In order to study p T spectrum effects, two different tunes of the PYTHIA6 event generator were used, ATLAS CSC and Perugia0, which give an average p T versus charged-particle multiplicity respectively below and above the data (Fig. 6 (left)), for most of the multiplicity range (N ch > 2). The difference between measurements obtained with response matrices corresponding to each of these Monte Carlo generators, is used as the corresponding systematic error contribution. Figure 6 (right) shows that this procedure introduces an uncertainty which is weakly dependent on η, and amounts to ±0.3%, when averaged over |η| ≤ 2.
For undetected particles, below a threshold of about 50 MeV/c, a value chosen to coincide with a track detection efficiency of 50%, the corresponding systematic uncertainty is obtained by varying their fraction by a conservative amount (−50% and +100%). (The fraction of the p T spectrum below 50 MeV/c is about 1% of the total, for both PYTHIA6 and PHOJET.) The resulting systematic uncertainty on dN ch / dη is η dependent, and found to range from −0.5% and +1.0%, at η = 0, to -0.75% and +1.5%, at |η| = 2.
The systematic uncertainty induced on P (N ch ) by the difference in p T between data and simulation is slightly sensitive to the tune of PYTHIA6 considered: for instance, at N ch = 90, varying the p T spectrum below 50 MeV/c by -50% and +100% induces a change of -5% and +9%, respectively for the ATLAS CSC tune and of -4% and +8%, respectively, for the Perugia0 tune.

Uncertainty evaluation
The results of the unfolding procedure and unfolding uncertainty estimate were cross-checked, in standard ways: -changing the regularization term by either varying β or the function F(U); -using two alternative unfolding procedures: Bayesian [60] and singular value decomposition [61].  The changes in unfolded distributions due to variations of the different elements used in the unfolding procedure, each having their own systematic uncertainty, were studied by considering all possible combinations: counting algorithm (3 cases), event generator used for efficiency correction (2 cases: PYTHIA6 and PHOJET tuned for diffraction), event generator used for the response matrix (2 cases PYTHIA6 CSC and PYTHIA6 Perugia0), low p T spectrum extrapolation below 50 MeV/c (3 cases: varying the integral below 50 MeV/c by −50%, 0% and +100%), response matrix parameterization (2 cases: Power law, and Padé), which correspond to 72 separate measurements that are correlated as a result of the unfolding procedure. To take these correlations into account, for a given energy and a given η range, the systematic uncertainty from all the sources considered was estimated as the overall spread between the resulting distributions. The average between these distributions (center of the band covered by the 72 curves) was used as the measurement. The total uncertainty originating from the unfolding procedure (evaluated for each multiplicity bin as a linear sum of the unfolding bias [56] and the unfolded distribution covariance calculated from the statistical uncertainty of the raw distribution) was added linearly to the systematic uncertainty defined as half the size of the band. This takes into account the fact that the unfolded distributions come from the same raw multiplicity histogram, hence the raw data statistical fluctuations are propagated to each of the unfolded distributions in a similar way and affect the spread of distributions uniformly. The resulting systematic uncertainty is expected to be highly correlated with multiplicity.
The systematic uncertainties from material budget, tracklet and track selection, detector alignment (evaluated by changing the geometry within alignment uncertainties), particle composition, strangeness corrections are found to be negligible. Among the non-negligible contributions to multiplicity distribution systematic uncertainties, the contribution from selection efficiency (Section 8.2.1) can be evaluated separately as it is a multiplicative correction. All the other contributions, from the uncertainty on p T , the extrapolation down to p T = 0, the counting algorithms, and the model dependence including the contribution from diffraction, are all mixed together through the procedure described above. For NSD and INEL event classes, the diffraction contribution is mainly significant in the zero multiplicity bin, which is absent for the INEL>0 event class.   Table 5: Systematic uncertainties, in percent, due to efficiency corrections including trigger and reconstruction efficiency, event generator dependence and diffraction, for INEL, NSD and INEL>0 event classes, in pseudorapidity intervals |η| < 0.5, 1.0 and 1.5, at √ s = 0.9, 2.76, 7 and 8 TeV. Numbers are given, from left to right, for multiplicities 0, 1 and 2. For the INEL>0 event class, the efficiencies are the same as for the INEL class, except that there is no bin with multiplicity 0. Note that the uncertainty listed here for bin 0 is not the total uncertainty in that bin. Bin 0 values are recalculated separately from simulation, which adds a significant spread of values that is seen in the total uncertainty estimate.  Table 6: Total systematic uncertainties in charged-particle multiplicity distribution measurements, in percent, for INEL, NSD and INEL>0 event classes, in pseudorapidity intervals |η| < 0.5, 1.0 and 1.5, at √ s = 0.9, 2.76, 7 and 8 TeV. Values for INEL and NSD event classes are given for three characteristic multiplicities, from left to right: zero multiplicity -mean multiplicity -five times the mean multiplicity. For the INEL>0 class, the first value is for N ch = 1, as there is no entry for N ch = 0.

Bin-to-bin correlations in systematic uncertainty
Due to specific nature of the correction procedure, the final multiplicity distributions contain bin-to-bin correlations coming from various sources with different properties. These sources can be categorized by their effect: -statistical correlations, resulting from the propagation of the raw distribution statistical uncertainties through unfolding process, their characteristics largely depend on the response matrix structure; -fully correlated shift of the distribution as a result of the uncertainty in normalization; -long-range correlations in systematic uncertainties as a result of multiplicity scale change (that is determined by the position of the response matrix bulk relative to the diagonal, see Section 7.1) in the unfolded distribution.
We found that first category correlations (statistical) are negligible compared to the last category (scaling), while the second category correlations can be easily factorized for fitting the multiplicity distribution. The three main sources of correlated uncertainty, falling into the third category, are the change of counting algorithm, the change of event generator tune used to produce the response matrix (see Fig. 6) and the variation of p T distribution under the detection threshold. In order to evaluate the effect on the final distributions we construct 18 intermediate distributions corresponding to all possible combinations of counting algorithm, event generator tunes and variations of p T spectra. All other sources of systematics are included (except the uncertainty corresponding to bin N ch = 0 renormalization, as it is not applied) with the same procedure described for the final distribution previously in this section. These 18 distributions are then treated independently in Section 9.5 to evaluate the effect of correlated uncertainties on NBD fits.

Pseudorapidity density
The various contributions to systematic uncertainties in dN ch / dη are summarized in Table 4 for the three event classes and the four centre-of-mass energies studied in this publication. For the INEL>0 event class, a precision of 1.5% is achieved, as the sensitivity to diffraction is negligible and the η range is reduced in the definition of this event class.
In the η range covered in this study (|η| < 2), we find that systematic uncertainties show essentially no η variation (Fig. 7), and are therefore strongly correlated bin-to-bin.

Multiplicity distributions
The efficiency uncertainties (Section 8.2.1) are only relevant at low multiplicity. For multiplicities above 8 to 9, the efficiency reaches 100% and the corresponding systematic uncertainty becomes negligible. Therefore, efficiency uncertainties are only given for a few characteristic low multiplicities. The total systematic uncertainties vary with multiplicity, therefore they are given in Section 8.2.1 only for a few characteristic values of the multiplicity.

Consistency checks
In the measurement of dN ch / dη, statistical errors are negligible. Therefore, the study of run-to-run fluctuations (measured RMS of dN ch / dη results in different runs) provides a check that all run dependent corrections are properly handled.
For the INEL and NSD event classes, contributions from run-to-run fluctuations are significantly smaller than the total systematic error (Fig. 7). For the INEL>0 event class, for which the precision is highest, the relative importance of run-to-run fluctuations is larger than for the INEL and NSD event classes (Fig. 7), but reaches at most 5% of the total systematic uncertainty.
As data correction procedures are significantly different between dN ch / dη and multiplicity distribution measurements, a test was performed to verify the consistency of the two measurements. At the four centre-of-mass energies and for the three pseudorapidity intervals used in this study, integrals of the multiplicity distributions were found to be consistent with the direct measurements of dN ch / dη, within errors.  In their common η range, |η| < 0.9, the three track counting algorithms discussed in Section 5.2 achieve a similar precision of 1%, and were found to give consistent results. The main difference is that for ITS+ and ITSTPC+, there is a detector calibration contribution to systematics for the TPC and the SDD, not present for the SPD. To achieve the largest possible η range, in the measurement of dN ch / dη versus η, the Tracklet algorithm is used alone.
The measurement at √ s = 0.9 TeV is shown in Fig. 8 compared with previous results. At |η| > 0.9 the measurement for the INEL event class is slightly lower than in ALICE's previous publication [2]. The difference comes mainly from: (a) the tuning of the MC generators for diffraction, as larger pseudorapidities are more sensitive to SD; (b) the subtraction of particles coming from the decay of strange particles was improved, using ALICE's measurement of strangeness [4]; and (c) the improvement of the η dependence of the Tracklet algorithm.
The discrepancy with UA5 for the INEL event class at large η could perhaps be related to the fact that UA5 used a 1/ M X variation of the single-diffractive cross section (see [53]). Note also that UA5 data seem to be internally inconsistent (see discussion in [62]). The measurement at √ s = 2.76 TeV is shown  [2,3], UA5 [66] and CMS [33]. Note that to avoid overlap of data points on the figure, the INEL>0 data were displaced vertically, and for these data the scale is to be read off the right-hand side vertical axis. Systematic uncertainties on previous data are shown as error bars (except for UA5, with coloured bands), while they are shown as grey bands for the data from this publication. in Fig. 9 top, and found to be consistent with the ALICE measurement at √ s = 2.36 TeV [2], as expected because of the small change in center-of-mass energy. The new measurements of dN ch / dη, at √ s = 7 TeV (Fig. 9 middle), show agreement both with the previous ALICE results for INEL>0 [3], and with CMS NSD data [33]. The measurements of dN ch / dη at √ s = 8 TeV (Fig. 9 bottom) show the 3% increase with respect to the 7 TeV data, which corresponds to what is obtained in the extrapolation from lower energy data. Comparisons of the η distributions at the four centre-of-mass energies (Fig. 10), for the three events classes, show no significant change of shape and a smooth increase of the charged-particle density with increasing energy.
The data for the INEL event class at √ s = 0.9 and 7 TeV were compared to simulations with current event generators (Fig. 11). At √ s = 0.9 TeV, EPOS LHC [63] and PYTHIA8 4C [64,65] are consistent with the data. PHOJET overestimates the data, while PYTHIA6 Perugia0 and Perugia 2011 underestimate the data. At √ s = 7 TeV, EPOS LHC, PHOJET and PYTHIA6 Perugia 2011 are consistent with the data. PYTHIA8 4C overestimates the data, while PYTHIA6 Perugia0 underestimates the data. Note that PYTHIA6 Perugia 2011, PYTHIA8 4C and EPOS LHC were tuned using LHC data.

Energy dependence of dN ch / dη at η = 0
The traditional definition for dN ch / dη at η = 0 is an integral of the data over the pseudorapidity range |η| < 0.5 The results of the measurements of dN ch / dη at η = 0 are given in Table 7. The energy dependence of dN ch / dη at η = 0 is of interest not only because it provides information about the basic properties of pp collisions, but also because it is related to the average energy density achieved in the interaction of pro- . Systematic uncertainties are shown as error bars for the previous data and as grey bands for the data from this publication. The scale is to be read off the right-hand side axis for INEL>0.   tons, and constitutes a reference for the comparison with heavy ion collisions. At mid-rapidity, dN ch / dη can be parameterized as dN ch / dη ∼ s δ . Combining the ALICE data with other data at the LHC and at lower energies, we obtain δ = 0.102 ± 0.003, 0.114 ± 0.003 and 0.114 ± 0.0015 6 , for the INEL, NSD and INEL>0 event classes, respectively, to be compared to δ 0.15 for central Pb-Pb collisions [54]. This is clear evidence that the particle pseudorapidity density increases faster with energy in Pb-Pb collisions than in pp collisions. Fits are shown on Fig. 12 and Table 8 gives extrapolations to centre-of-mass energies of 13 and 14 TeV (LHC design energy). While this paper was being prepared, the first measurement at 13 TeV by CMS appeared [67], resulting in dN ch / dη| |η|<0.5 = 5.49 ± 0.01 (stat) ± 0.17 (syst) for inelastic events, which is consistent with our extrapolation of 5.30 ± 0.24. Over the LHC energy range, from 0.9 to 14 TeV, while the centre-of-mass energy increases by a factor 15.5, extrapolation of present data for dN ch / dη at η = 0 shows an increase by factors 1.75 ± 0.03, 1.87 ± 0.03 and 1.87 ± 0.01, respectively for the three event classes. The multiplicity increase is similar for NSD and INEL>0 classes but slightly lower for the INEL class.

Multiplicity distributions of primary charged particles: measurements
The results of ALICE measurements of multiplicity distributions of charged primary particles are displayed as probability distributions (P (N ch )) in Figures   5.37 ± 0.25 6.62 ± 0.20 6.98 ± 0.10 Table 8: Extrapolations of dN ch / dη, at η = 0, for the three event classes, to higher energies at the LHC ( √ s = 13 and 14 TeV), using the fits described in the text. first two event classes the measurements were obtained in three pseudorapidity intervals |η| < 0.5, 1 and 1.5, and for INEL>0 in |η| < 1. At √ s = 7 TeV, P (N ch ) varies over 6 to 7 orders of magnitude and the multiplicity range reaches up to 160 in |η| < 1.5 for both INEL and NSD event classes. In |η| < 0.5 and |η| < 1, the observed multiplicity reaches 10 times the mean multiplicity. It is expected that the average energy density in proton collisions at the LHC, at √ s = 14 TeV, is about 5 to 15 times smaller than energy densities reached in gold ions at RHIC [76]. Therefore, in proton-proton collisions of multiplicity exceeding 10 times the average multiplicity, energy densities should overlap with those of heavy ion collisions at RHIC, allowing to compare properties of systems with very different collision volumes (two to three orders of magnitude) but the same energy density. Future runs of the LHC should allow extending much further the range of multiplicities probed so far.
The high-multiplicity tail of the distribution increases as expected with increasing energy (Fig. 16). This behaviour is studied quantitatively in Section 9.6 on KNO scaling and q-moment analysis.
The measurements presented in this publication are consistent with previous ALICE data, for INEL [2] at √ s = 0.9 TeV and INEL>0 [3] at √ s = 7 TeV, in the multiplicity range where they overlap (Fig. 17). The wavy structure already observed by ALICE in [2,3], for multiplicities above N ch = 25, and previously by   UA5 [27], is still hardly significant, and it is not present in the raw data. This feature was also observed in a study of CMS data [77]. However, applying the same procedure to Monte Carlo data produces similar structures in the unfolded distribution, while no oscillation was present at particle level. We conclude that this structure is an unfolding procedure artifact. The period of the structure is related to response matrix width.

Comparison of multiplicity distributions with other experiments and models
CMS data are available for the NSD normalization only. At √ s = 0.9 and 7 TeV, in the three pseudorapidity intervals where they could be compared, the ALICE multiplicity measurements are generally in agreement with CMS data [34] (Fig. 17). However, at √ s = 0.9 TeV, above N ch ≈ 40 in |η| < 1 and above N ch ≈ 50 in |η| < 1.5, the ALICE multiplicity distributions tend to be higher than CMS data. The difference is not significant in each isolated bin, however the general trend seems to be different. As the bin-to-bin correlations in CMS's results are unknown, it is very difficult to obtain reliable unbiased quantitative comparison. A possible source of the discrepancy is the different treatment of single-diffractive events in the two analyses making the NSD event sample definitions not strictly compatible. More precisely, there was no diffraction tuning in the simulations used by CMS and, moreover, ALICE's criteria for events to be considered single-diffractive (see [47]) include a fixed cut on diffractive mass that differs from the value that CMS used. This explanation is supported by the fact that the difference is insignificant at √ s = 7 TeV, where single-diffractive contribution is negligible. The comparison with CMS data is further discussed in the Section 9.5, where NBD fits to the data are reported.
Measured multiplicity distributions are compared to models (Fig. 18), for two energies, the lowest one ( √ s = 0.9 TeV), and the energy at which there is the largest event sample ( √ s = 7 TeV), for the INEL event class, in the pseudorapidity range |η| < 1.   timate the high multiplicity part of the distribution. PYTHIA6 Perugia 2011, EPOS LHC and PYTHIA8 4C give a reasonable fit of the low multiplicity region but underestimate the data above N ch ≈ 60.

Parameterization of multiplicity distributions with NBDs
Single Negative Binomial Distributions (NBD) have been traditionally used to parameterize particle multiplicity distributions in hadron collisions P NBD (n, n , k) = Γ (n + k) Γ (k) Γ (n + 1) where n is the average multiplicity and the variance is given by The parameter k is related to the two-particle correlation function, in the pseudorapidity interval considered [78]. In the limit k → ∞, the NBD becomes a Poisson distribution.
In previous ALICE data, for the NSD event class, no strong deviation from a single NBD fit was observed at √ s = 0.9 and 2.36 TeV, for |η| ≤ 1, while a hint of a substructure appears at |η| < 1.3 [2]. At √ s = 7 TeV, for the INEL>0 event class, the single NBD fit slightly underestimated the data at low multiplicity (N ch < 5), and slightly overestimated the data at high multiplicities (N ch > 55) [3].
In the present data, for all event classes, already at √ s = 0.9 TeV, there is a hint that single NBD fits start diverging from the data at the higher multiplicity, for |η| < 0.5 and 1.0. More significant departure from   The appearance of substructures in multiplicity distributions attributed to the occurrence of several sources in the process of particle production [79][80][81][82], can be parameterized by fitting the data with two NBDs. Indeed, a much better fit to the data is obtained by using a weighted sum of two NBD functions P (n) = λ αP NBD (n, n 1 , k 1 ) + (1 − α) P NBD (n, n 2 , k 2 ) This type of function, however, is not meant to describe the value P(0) for INEL and NSD distributions, which occurs when the η acceptance is limited, therefore the bin n = 0 was excluded from the fit and an overall normalization factor (λ) was introduced, as a free parameter, to account for this. At √ s = 7 TeV, P. Gosh et al. in [77] perform a similar fit for CMS data, however without using an overall scale factor (λ) and with no account for the effect of bin-to-bin correlations.
As was described in Section 8.2.2, 18 multiplicity distributions, corresponding to main sources of correlated uncertainty, were fitted independently. For each of these distributions systematic uncertainties corresponding to the non-leading sources were found to have negligible correlations, as expected. Thus these residual systematic uncertainties were added to the diagonals of the statistical covariance matrices for these distributions, that were produced in the unfolding procedure. In Tables 9 to 11 we present the central values of the fit parameters, i.e. center of the maximum spread of each parameter between the 18 independent fits. Half of the spread is used to estimate the systematic uncertainty of the parameters (the statistical uncertainty being negligible). However these estimates should be treated with caution due to the significant correlations between parameters. The multiplicity distributions for each of the 18 cases differ by their slope and intersect at low-to-middle values of multiplicity (N ch ≈ 10 to 30 depending on η window and energy). There is no single point of intersection for a given set, but rather a small interval. Therefore extreme cases were selected as the curves with highest and lowest values at high multiplicities (and consequently lowest and highest values at low multiplicity). To indicate the change in shape of the parametrization allowed by the correlations present in systematic uncertainties, the parameters for these extreme cases are given in Tables 9 to 11 Table 9: Double-NBD fit parameters for charged-particle multiplicity distributions of INEL events. The shaded (middle) row for each entry contains central values of the parameters with their estimated uncertainty, and the surrounding rows represent the parameter values of the bounding fits: topmost in the first row and the bottommost in the last one.
The shape evolution is quantified by the parameter n 2 , which tends to increase with increasing η range and with increasing centre-of-mass energy. The observed relation n 2 ≈ 3 × n 1 , is consistent with the analysis of CMS data reported in [77], despite the fact that the scaling parameter λ was not used in their fit function.
In the bins |η| < 1 and |η| < 1.5, the CMS NSD data at √ s = 0.9 TeV (Fig. 17)   at high multiplicity as compared to the ALICE data (Section 9.4). We have made our own fit of CMS data excluding the bin N ch = 0 (Table 12) and using the overall scaling parameter (λ). As we cannot take into account correlations within CMS data, we compare these fits with fits of our final distributions, ignoring the correlations. We find that the relative weight of the second NBD component in the CMS data is smaller than in the ALICE data, while other parameters are compatible within their respective uncertainties. This confirms the different trends observed in distribution tails.
In conclusion, a double-NBD function provides a precise description of the entire set of multiplicity distributions measured in this experiment.

KNO studies
The KNO variable N ch / N ch provides another way to study the evolution of the shape of multiplicity distributions with varying centre-of-mass energies and varying pseudorapidity intervals. This study is carried out for the NSD event class only so that SD events, which may have a different behaviour, are not included in the data samples.
The quantities of interest are derived from the original set of 72 multiplicity distributions and the resulting spread is used to estimate the various uncertainties presented in this section.
9.6.1 Evolution of the shape of multiplicity distributions with √ s The KNO test in the range 0.9 to 8 TeV is limited by the multiplicity reach at 0.9 TeV. KNO-scaled distributions and their ratios were obtained for each of the available combinations of corrections with the same procedure used for multiplicity distribution measurements (averaging the 72 cases listed in Section 8.2 and using the spread as a measure of the systematic uncertainty). Bin to bin correlations were ignored when comparing KNO distributions and q-moments at various centre-of-mass energies. Consequently, the relative errors obtained on the ratios are somewhat overestimated. Ratios between the two highest energies and 0.9 TeV exceed the value 2 at N ch / N ch larger than 5.5, 5 and 4.5, for |η| < 0.5, |η| < 1 and |η| < 1.5, respectively (Fig. 19). This confirms that KNO scaling violation increases with increasing pseudorapidity intervals. The shape of the KNO scaling violation reflects the fact that the high-multiplicity tail of the distribution increases faster with increasing energy and with increasing pseudorapidity interval than the low (N ch ≤ 20) multiplicity part as already noted in Section 9.5.

Quantitative KNO study with normalized q-moments
Multiplicity distributions may be characterized by their normalized q-moments defined as where q is a positive integer studied here for values 2, 3, 4 and 5, for NSD events ( Fig. 20 and Table 13). For the three pseudorapidity intervals studied here (|η| < 0.5, 1 and 1.5), C 2 remains constant over the energy range, C 3 shows a small increase with increasing energy for the two largest η intervals, C 4 and C 5 show an increase with increasing energy, which becomes stronger for larger η intervals.   ∆η η η λ λ λ α α α n n n 1 1 1 k k k 1 1 1 n n n 2 2 2 k k k 2 2 2 χ χ χ 2 2 2 n n n dof ALICE 0.9 1 0.95 ± 0.02 0.56 ± 0.10 5.0 ± 0.8 2.6 ± 0. 8 Table 12: Comparison of double-NBD fit parameters for charged-particle multiplicity distributions of NSD events by ALICE and CMS, obtained without accounting for correlations using the function defined by Eq. (16). CMS data are taken from [34]. The last column gives the χ 2 value and the number of degrees of freedom (n dof ).

Discussion of results and conclusion
The ALICE Collaboration has carried out a detailed study of pseudorapidity densities and multiplicity distributions of primary charged particles produced in proton-proton collisions, at √ s = 0.9, 2.36, 2.76, 7 and 8 TeV, in the pseudorapidity range |η| < 2. A large increase of event sample size compared to previous ALICE publications, combined with improved measurement techniques, was used to study the evolution of charged-particle multiplicities over the whole centre-of-mass energy range covered so far by the LHC. The data at the highest energy appear as a smooth continuation of lower energy data, both in shape and in magnitude.
The pseudorapidity density of charged particles, dN ch / dη, was measured as a function of pseudorapidity   : KNO-scaled distribution N ch P (N ch ) versus the KNO variable N ch / N ch at √ s = 0.9, 2.76, 7 and 8 TeV, for three pseudorapidity intervals: |η| < 0.5 (top), 1.0 (middle) and 1.5 (bottom). In each case, ratios to the distribution at √ s = 0.9 TeV are shown, on the right-hand side parts of the figures. As N ch / N ch takes different values at different centre-of-mass energies, ratios were obtained by interpolating the KNO-scaled distributions, and uncertainties were taken from the nearest data point. Bands represent the total uncertainties. in the range |η| ≤ 2. The relative precision achieved at η = 0 for √ s = 7 TeV is 5.5%, 2.6% and 1.3% for INEL, NSD and INEL>0 event classes, respectively.
The power law parameterization of dN ch / dη at η = 0, s δ , provides a good description of the data from ISR to LHC energies: δ = 0.102 ± 0.003, 0.114 ± 0.003 and 0.114 ± 0.001, for the INEL, NSD and INEL>0 event classes, respectively, to be compared to δ 0.15 for Pb-Pb collisions [54]. The ALICE Collaboration has shown clearly that the particle pseudorapidity density increases faster with energy in Pb-Pb collisions than in pp collisions. The extrapolation of dN ch / dη at η = 0 to the nominal LHC energy ( √ s = 14 TeV) is obtained with a precision of 4.6%, 3.0% and 1.3% for INEL, NSD and INEL>0 event classes, respectively.
Multiplicity distributions of primary charged particles were measured in three pseudorapidity ranges: [34] Fig. 20: Centre-of-mass energy dependence of the q-moments (q = 2 to 4, left-hand scale, and q = 5, right-hand scale) of the multiplicity distributions for NSD events in three different pseudorapidity intervals (|η| < 1.5 top, |η| < 1.0 middle and |η| < 0.5 bottom). ALICE data (black) are compared to UA5 [27] (red) for |η| < 0.5 and |η| < 1, at √ s = 0.9 TeV, and with CMS [34] (blue) at √ s = 0.9 and 7 TeV for |η| < 0.5. The error bars represent the combined statistical and systematic uncertainties. The data at 0.9 and 7 TeV are slightly displaced horizontally for visibility. |η| < 0.5, |η| < 1.0 and |η| < 1.5. At √ s = 7 TeV, N ch reaches about 70, 120 and 150, in |η| < 0.5, 1.0 and 1.5, respectively, with the present statistics. These correspond to multiplicity densities 8 to 10 times the average multiplicity density. Based on the J. D. Bjorken formula [83], a characteristic collision energy density can be estimated, which increases by the same factor. For a qualitative estimate, assuming that the average energy density in pp collisions at √ s = 7 TeV is of the order of 1 GeV/fm 3 (see for example [76]), a density of 10 GeV/fm 3 should be reached with high multiplicity pp collisions, similar to the energy density of Au-Au central collisions at RHIC [84]. When LHC runs at its nominal centre-of-mass energy of 14 TeV, high multiplicity proton-proton collisions will provide further direct comparisons of nuclear matter properties for interacting systems with similar energy densities but very different volumes.
At LHC energies, above 2 TeV, multiplicity distributions can no longer be represented by a single NBD, but a double NBD gives a good representation of the data. The deviation from single NBD is already visible in the tail of the distributions at √ s = 0.9 TeV, but becomes increasingly large as the centre-of-mass energy increases.
A test of KNO scaling between √ s = 0.9 and 8 TeV confirms that KNO scaling violation increases with increasing √ s and, at a given centre-of-mass energy, with increasing width of pseudorapidity intervals.
Comparisons with models in the pseudorapidity range |η| ≤ 1 show that none of the event generators considered is able to describe all properties of charged particle production up to √ s = 8 TeV. PYTHIA6 Perugia 2011, PYTHIA8 4C and EPOS LHC describe the pseudorapidity densities fairly well at √ s = 7 TeV, as well as the multiplicity distributions, but not above N ch ∼ 60. The fact that PYTHIA6 Perugia 2011, PYTHIA8 4C and EPOS LHC are in reasonable agreement with the data presented in this publication can probably be attributed to the fact that these generators were adjusted using the first LHC data. Studies of charged-particle production are refining the understanding of global properties of protonproton collisions at the LHC. It was already demonstrated that, as with lower energy data, there is a strong correlation between multiplicity and mean transverse momentum p T [85]. However, there is also evidence [4] that the high multiplicity events have a topology more spherical than expected from current event generators PYTHIA6, PYTHIA8 and PHOJET, suggesting that single hard-parton collisions may not be the only contributors to high multiplicity events. The general picture that emerges from this study is that, from √ s = 0.9 to 8 TeV, multiplicity distributions and charged particle pseudorapidity densities follow a smooth evolution.

Acknowledgements
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector: State Committee of Science, World Federation of Scientists (WFS) and Swiss Fonds Kidagan, Armenia; Conselho Nacional de Desenvolvimento Científico e Tecnológico