Centrality determination in heavy-ion collisions with the LHCb detector

The centrality of heavy-ion collisions is directly related to the created medium in these interactions. A procedure to determine the centrality of collisions with the LHCb detector is implemented for lead-lead collisions at √s NN = 5 TeV and lead-neon fixed-target collisions at √s NN = 69 GeV. The energy deposits in the electromagnetic calorimeter are used to determine and define the centrality classes. The correspondence between the number of participants and the centrality for the lead-lead collisions is in good agreement with the correspondence found in other experiments, and the centrality measurements for the lead-neon collisions presented here are performed for the first time in fixed-target collisions at the LHC.


Introduction
In the context of heavy-ion collisions, centrality is a quantity of relevance since it is directly related to the medium formed by the colliding nuclei, and measures the overlap region between the two nuclei in a collision. The centrality of a collision is characterised by the impact parameter (b) between the two nuclei, i.e. the distance between their centres in the plane transverse to the beam axis. The impact parameter defines the overlap region of the nuclei and thus determines also the size and shape of the resulting medium. A schematic view of a heavy-ion collision is shown in Fig. 1.
The geometry of the collision is related to the number of nucleons that participate in it and the number of nucleon-nucleon collisions. These quantities are not directly accessible and hence need to be derived from the data recorded during the collisions by making use of other quantities that scale approximately with the number of participating nucleons, such as the outgoing particle multiplicity. For this purpose, a Glauber model is often used [1].
This paper presents the centrality determination with the LHCb detector for lead-lead (PbPb) collisions at a centre-of-mass energy √ s NN = 5 TeV, which is in agreement with results obtained by the ALICE [2,3], ATLAS [4] and CMS [5] collaborations, and the first centrality determination in fixed-target mode at the LHC, for lead-neon (PbNe) collisions at √ s NN = 69 GeV. After introducing the Glauber model (section 3) and the assumptions it relies on, the datasets (section 4), the centrality determination procedure with its results (section 5) and the study of systematic uncertainties (section 6), for PbPb and for PbNe collisions are described.

The LHCb detector
The LHCb detector [6,7] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector (VELO) surrounding the beam interaction region with two pile-up (PU) stations upstream from the interaction point, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex (PV), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad (SPD) and pre-shower detectors, an electromagnetic (ECAL) and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. A gas-injection system for beam-gas interactions (SMOG) [8] is installed in the LHCb detector, which gives the unique possibility of injecting a low pressure noble gas and collecting fixed-target collisions (proton-nucleus or nucleus-nucleus).

Glauber model
The centrality of a nucleus-nucleus collision is related to the overlap region between the nuclei where the nucleons are colliding. In practice, the particles produced by the collisions are not originating purely from hadronic interactions between the nuclei, but also from electromagnetic processes. Therefore, a model is needed to isolate the hadronic part and subsequently define the centrality classes. The most common approach in heavy-ion physics to model the collisions of two nuclei is to consider the transverse shapes of the nuclei. This shape, i.e. the nuclear density, is described by a two-parameter Fermi distribution (2pF), also known as Woods-Saxon distribution [1], for each nuclear species considered, defined as where r stands for the radial distance from the centre of the nucleus. The constant ρ 0 corresponds to the density at the centre of the nucleus, and R to the nuclear radius, which is approximately the radial extension of the bulk of the nucleus. The diffusivity a describes how abruptly the density falls at the edge of the nucleus. The last parameter, w, is used to describe nuclei whose maximum density is reached at a radius r > 0. The values of these parameters are taken from other experiments, typically involving lepton-nucleus collisions and other types of nuclear spectroscopy [9,10]. The 2pF distribution can be seen in Fig. 2 with parameters R = 6 fm, a = 0.5 fm and w = 0 for illustration. The Glauber model is generally approached in two ways, the optical Glauber model and the Monte Carlo (MC) Glauber model. The two colliding nuclei are labelled A and B. In the optical model, it is considered that nucleons from projectile A see the target B as a continuous distribution, which is described by an analytical function, and vice versa. This is also called the optical limit approximation. Subsequently, the overlap area, the number of participating nucleons (N part ) and the number of binary nucleon-nucleon collisions (N coll ) can be obtained analytically.
On the other hand, in the MC Glauber model, the calculation is performed through a MC method where nucleons from each nucleus A and B are generated as hard spheres 1 and are placed around the respective centres of the nuclei following the 2pF distributions. Then a random impact parameter b is sampled from the distribution dσ/db = 2πb. Finally the nuclei are made to collide, with the following assumptions: • nucleus-nucleus collisions are considered to be a superposition of several independent nucleon-nucleon collisions; • nucleons are treated as hard spheres moving in straight lines all along the process, even if they have undergone a collision; • nucleons have a geometrical transverse cross-section (σ inel NN ) and two nucleons collide if the transverse distance between their centres is d < σ inel NN /π.
The average values of the number of participating nucleons N part , of binary collisions N coll and other quantities, are obtained with a Monte Carlo simulation. The distributions of N part , N coll and other quantities of interest can then be obtained for any centrality interval.
The Glauber model has two relevant external inputs, the nucleon-nucleon inelastic crosssection σ inel NN and the spatial distribution given by the 2pF distribution with its parameters R, a and w. The cross-section is obtained from a phenomenological parametrisation tuned on data, using measurements from a broad range of energies from ∼ 20 GeV to ∼ 60 TeV, given by σ inel NN (s) = A + B ln 2 (s), with A = 25.0 ± 0.9 and B = 0.146 ± 0.004 [12]. The The distribution of the VELO clusters exhibits a peak structure with a sharp fall at 45 000 clusters. This is related to the total number of readout channels in the VELO, leading to saturation for high occupancy events.
optical approach describes fairly well the collision process but does not completely capture the physics of the total cross-section and leads to distortions in the estimation of N part and N coll compared to the estimation made with the MC approach [1]. Therefore, in the following, the MC Glauber model is used.

Data
For the centrality determination a minimum bias (MB) data sample is needed, that is, data that have the minimum possible number of selections applied, to not bias the sample.

PbPb collisions
For PbPb collisions, the sample used for this analysis corresponds to the data recorded in a special run of the 2018 PbPb data-taking period, at a centre-of-mass energy per nucleon of √ s NN = 5 TeV. Throughout the data-taking period, PbPb collisions were recorded in parallel with fixed-target PbNe collisions, between Pb beams and atoms of Ne injected in the VELO with the SMOG system. However, for this particular MB data sample, no gas was injected. The resulting distribution of the number of VELO clusters (nVeloClusters), which are clustered energy deposits in the VELO stations, and of the energy deposited in the ECAL for this MB sample can be seen in Fig. 3.

PbNe collisions
The PbNe data sample corresponds to a MB sample from the whole data-taking period when gas was injected. The PbNe collisions occur at a centre-of-mass energy per nucleon of √ s NN = 69 GeV where the Ne atoms act as a fixed target. In order to avoid contamination from PbPb collisions that were recorded simultaneously, only events for those bunch crossings where a filled bunch from the incoming Pb beam crosses an empty bunch of the outgoing Pb beam are selected. There is some residual contamination in the data sample and extra selections are applied to increase the fraction of PbNe collisions.
The different topology of the PbPb and PbNe events allows to disentangle these two types of events by setting an upper limit on the number of clusters in the PU stations, which are located upstream from the nominal interaction point. Since the PbNe collisions are all boosted downstream, i.e. towards the detector, naturally a low number of clusters in the PU stations is expected. On the other hand, since PbPb collisions are symmetric, a larger number of clusters in the PU stations is observed for these collisions.
Another source of contamination are the beam-gas collisions that take place far upstream. Since the injected gas can travel up to 20 m in either direction from the nominal interaction point, the incoming Pb beam can undergo interactions with the gas before arriving into the VELO tank. These events can produce forward particles, hitting the PU stations and depositing energy in the detector.
To ensure a high enough purity of PbNe collisions in the sample, events with clusters in the PU stations are rejected. The effect of this requirement can be seen in Fig. 4. On the left, three different populations can be seen, the high-slope population which corresponds to the very upstream events (indicated by the red line), the middle-slope population which corresponds to ghost PbPb collisions (indicated by the green line), 2 and finally the continuum which corresponds to the PbNe collisions of interest which present no clusters in the PU stations (enclosed by the black lines).
Only the central-region PbNe events are used in what follows, i.e. events whose primary vertex is located in the range z PV ∈ [−200, 200] mm. The number of VELO clusters and distribution of ECAL energy are shown in Fig. 5.

Centrality determination
In theory, any observable that scales monotonically with impact parameter could be used for classification according to centrality. In practice, the reach in centrality possible with the tracking detectors of LHCb is limited by the performance of these detectors at high track multiplicities. In the case of PbPb collisions this means that the VELO information cannot be used for this purpose, since its sensors saturate under these conditions as seen in the rightmost part of the plot on the left of Fig. 3. In contrast, the ECAL has the advantage of not saturating even in the most central collisions, as can be seen in the right plot of Fig. 3. In the case of PbNe collisions, the VELO does not saturate but still cannot be used for the centrality determination, since the relevant events take place all along the length of the VELO. This means that the measured VELO multiplicity depends on the position of the collisions along the beam axis. For this reason the energy deposited in the ECAL is used for the multiplicity determination of both PbPb and PbNe collisions.
Centrality classes are defined as quantiles of the inelastic PbPb or PbNe cross-section. The data contain contributions to the deposited energy in the ECAL from both hadronic and electromagnetic origins (the latter originating from peripheral collisions where the electromagnetic interactions dominate and ultra-peripheral collisions (UPC)). Thus, the energy spectrum cannot be used directly to define the desired quantiles for centrality. To estimate the hadronic component, a MC Glauber model [12] is used to simulate the colliding nuclei, and from the resulting quantities such as b, N coll or N part , the expected observable can be constructed, which is in this case the energy deposited in the ECAL. The parameters of the model are then tuned to fit the ECAL energy distribution from the data. Finally, the centrality quantiles are defined from the simulated distribution that corresponds only to the hadronic part of the interaction. Geometric quantities from the Glauber MC can then be mapped to the data for each centrality class.

Methodology
In this section the simulation of the events is described first, then the generation of the simulated ECAL energy distribution and the steps to fit it to the data are explained. Once the fit has been performed, the simulated distribution is split into centrality classes based  on the fraction of the total hadronic distribution integral, and the geometric variables of each class are mapped to the measured events falling in the same class.

Events simulation
The first step is to simulate the collisions using the TGlauberMC software from Ref. [12]. 3 The parameters for the 2pF density function are listed in Table 1. In all cases, the parameter w is assumed to be equal to 0.
One million PbPb collisions are simulated using the corresponding nucleon-nucleon cross-section σ inel NN = 67.6 ± 0.6 mb for a centre-of-mass energy of √ s NN = 5 TeV, and one million PbNe collisions were simulated using the corresponding nucleon-nucleon cross-section σ inel NN = 35.4 ± 0.9 mb for a centre-of-mass energy of √ s NN = 69 GeV.
The N part and N coll values of every simulated collision is then computed. With these numbers, the number of ancestors N anc is defined as which effectively scales with the number of sources of particle production, with a relative weight for the participating nucleons and the number of collisions. This is motivated by the fact that the particle multiplicity is expected to scale with N part when soft processes dominate and to scale with N coll when hard processes dominate [16][17][18][19][20]. Below a centreof-mass energy of 100 GeV soft processes are expected to dominate. The parameter f determines the fraction of soft processes that contribute to the particle production and has to be determined with a fit. In Fig. 6, the distributions of N part , N coll and N anc in simulation are displayed with, as an example, f = 0.751 for the PbPb case. To get the distribution of particles originating from the collision, N anc is convoluted with a negative binomial distribution (NBD) which has been extensively used to model particle production and has been shown to be a reasonable approach at diverse energy and rapidity regimes [21][22][23][24][25]. The NBD is given in its discrete form by with p = µ k + 1 −1 , where µ and k are parameters related to the mean and spread of the NBD respectively, and n is the number of particles that are produced and deposited energy in the ECAL. Since there are N anc particle sources for each nucleus-nucleus collision, each producing particles following an NBD, the NBD is sampled N anc times to get the particle multiplicity distribution. Figure 7 illustrates the NBD function and the result after sampling it N anc times for each event to obtain the distribution of the number of outgoing particles (N out ) which deposit energy in the ECAL. The mean energy per particle in the ECAL is assumed to come mainly from π 0 decays. The π 0 spectrum is approximated with the charged pion spectrum seen in data in MB pp collisions at 5 TeV. The value for the mean energy deposited per particle in the ECAL is found to be E PbPb = 10.4 GeV. The simulated ECAL energy distribution can be seen in Fig. 8.
The same procedure is repeated for the PbNe case. However, since there are no pp data at √ s NN = 69 GeV, pNe collisions from the 2015 data-taking period at this centre-of-mass energy are used. The value for the mean energy found is E PbNe = 10.4 GeV. This value is consistent with the one found for PbPb, but this is compensated by a much lower number of particles produced in the collisions.

Fit model
The parameter k is linked to the width of the NBD distribution. The ALICE collaboration uses a value k = 1. 6 [26] in their analysis. A comparison is made between PbPb simulated distributions varying k between 1.0 and 2.0 while everything else is kept constant. It is found that there is no significant dependence on this parameter. The resulting distributions  can be seen in Fig. 9. For this reason the parameter k was fixed to k = 1.5 to be in the middle of the explored range. Finally the model has only two free parameters, f from Eq. 2 and µ from the NBD.

Fit to the PbPb collision data
To fit the simulated distribution to the data a χ 2 function on both distributions is minimised. It is defined as where E i and O i are the expected and observed values for the i th bin, i.e. the simulated and measured values for a given energy bin, assuming the values are counts with Poissonian errors. As previously mentioned, in order to avoid a possible contamination at low energy of electromagnetic origin, the fitting range is chosen to be from 2 to 52 TeV. The MC Glauber energy distribution is normalised to the data in the energy range of 5 to 15 TeV to avoid the extremes of the distributions. The fit procedure is delicate since the variables f and µ are highly correlated as they both modulate the horizontal reach of the distribution: µ is related to the mean value of the NBD, and thus the high µ then the higher the number of particles produced per collision and consequently the more energy is deposited in the ECAL. On the other hand, f controls how alike the final distribution is to the distributions of N part or N coll , which have a different reach on the x-axis. For the same reason f controls the shape of the right shoulder, which helps disentangle f from µ.
In the first step µ is fixed to a test value, µ = 3.85, and the similarity between the right shoulder of the data and the simulated distributions is evaluated for 1000 different values of f ranging from 0 to 1. For this procedure the simulated distribution is horizontally scaled to match the data by a factor H s defined as where E Data 300 is the energy of the last bin containing more than 300 events in the data (bin centre), and E Glauber 300 is the same in the simulated distribution. This scaling procedure is shown in Fig. 10, where the data have been normalised to 1, meaning that 300 events now correspond to roughly 1.6 × 10 −4 . For this step, the χ 2 is computed from 35 to 52 TeV to only consider the right shoulder. The resulting values for the χ 2 as a function of f can be seen in Fig. 11.
The minimum of the χ 2 is found to be at f = 0.83. This result is taken as a reference to reduce the range in f for the subsequent grid search. The allowed range for f is set to be [0.60, 0.93], from which the range for µ is chosen to be [3.7, 9.0]. A grid of 100 × 100 is defined and the χ 2 is computed from 2 to 52 TeV at every point of the grid. The result of this grid search can be seen in Fig. 12.
The best fit is not associated to the point of the grid with the lowest value, since this is prone to be affected by the fluctuations from the random NBD sampling. Instead the χ 2 map from Fig. 12 is considered and the minimum parametrised by f and by µ separately. In order to get the parametrisation as a function of f , the minimum of the χ 2 as a function of µ in bins of f is found. Like this, the value of the minimum in each slice along µ (at fixed f ) is assigned to the corresponding value of f . The same is done for the µ dependence. Consequently, for all values of f and µ the f -parametrised minimum and µ-parametrised minimum are constructed. In Fig. 13 the parametrisations as a function of f and µ are shown and an example slice is displayed to illustrate the process. From this procedure, the best fit is found at (f, µ) = (0.866, 6.778) with a χ 2 /ndf = 3.025, which is the one shown in Fig. 12. The number of points on the grid is not limiting the precision on the result.
Finally another grid search is performed with the same amount of points but on a narrower range, namely f ∈ [0.79, 0.92] and µ ∈ [5.7, 7.9]. The resulting χ 2 map is shown in Fig. 14. Two results are shown as a best fit. The star is the nominal value, whereas the cross is obtained by using an alternative method where the slice is fitted to obtain the minimum of this fit. The best fits found are (f, µ) = (0.869, 6.814) with χ 2 /ndf = 2.82 and (f, µ) = (0.869, 6.853) with χ 2 /ndf = 2.83. These best fits correspond to the star and the cross on Fig. 14 respectively. Since the goodness of fit is virtually the same for the two best fits, the one which is kept is (f, µ) = (0.869, 6.814).
The final result of the fit can be seen in Fig. 15. On the right plot of Fig. 15 a zoom for all values of f and µ, that is, the f -parametrised minimum (left) and the µ-parametrised minimum (right) for the PbPb case using the coarse grid. These are fitted by a 5 th and 6 th degree polynomial respectively whose minima are at f = 0.866 and µ = 6.778.  of the low-energy part of the distribution is displayed, where the discrepancy between the MC Glauber and the data, due to the presence of events of electromagnetic origin, becomes clear. This will be addressed in more detail in Sect. 5.2. This region below 0.5 TeV is well outside the fitting range, which starts at 2 TeV. The Glauber model, with its parameters obtained from a fit to the data, can thus be used to define the centrality classes in PbPb collision data.

Fit to the PbNe collision data
The same χ 2 function from Eq. 4 is used to evaluate the goodness of fit of the Glauber MC to the data for PbNe collisions. The fitting range is chosen to be from 0.5 to 3.9 TeV in order to avoid possible contamination from electromagnetic origin present at low energy. The Glauber MC energy distribution is normalised to the data in the energy range of 0.5 to 2 TeV to not consider the tails of the distributions, even in the case where f takes on the extreme values. Since in this scenario the shape of the energy distribution at high energy is not as characteristic as it is for PbPb, the approach of trying to fit f first cannot be applied, and the allowed range remains f ∈ [0.0, 1.0], since there is still sensitivity to the contributions of N part or N coll . The range in µ is chosen accordingly to be µ ∈ [1.0, 3.4]. A grid of 200 × 200 is defined in the previously mentioned ranges and the χ 2 is computed at every point of the grid. The result of this grid search can be seen in Fig. 16.
To find the best fits shown in Fig. 16, the same approach described for the PbPb case is used. The values of f and µ-parametrised minima are obtained by taking the minimum value for every slice (first method, star in Fig. 16), and by fitting each slice and getting the minimum of the fit (second method, cross in Fig. 16). However, in the following step only the µ-parametrised minima are fit to get the optimal µ whereas for the f -parametrised minima, the f value where the χ 2 is minimum is picked without fit. This is because the f -parametrised distributions show a steep decrease as f approaches one, and fitting a function around that region is difficult. The best fit is found at (f, µ) = (0.980, 3.156)  with a χ 2 /ndf = 1.039, and at (f, µ) = (0.995, 3.174) with a χ 2 /ndf = 1.031. These best fits are shown in Fig. 16 by a star and a cross, respectively.
Finally another grid search is performed where the same amount of points is used but on a narrower range, namely f ∈ [0.8, 1.0] and µ ∈ [2.9, 3.4]. From this grid, the best fits are found at (f, µ) = (0.996, 3.157) with a χ 2 /ndf = 1.026, and at (f, µ) = (0.992, 3.173) with a χ 2 /ndf = 1.031, computed using the same methods described above. The resulting χ 2 map and the best fits are shown in Fig. 17 by a star and a cross, respectively. Since the goodness of fit is virtually the same for all best fits in the coarse and fine grid, the one kept is the one which results in the smallest χ 2 value, that is, (f, µ) = (0.996, 3.157).
The final result of the fit can be seen in Fig. 18. On the right plot, a zoom of the low-energy part of the distribution is shown, where the discrepancy between the Glauber MC and the data, due to the presence of events of electromagnetic origin, becomes clear. This region below 0.1 TeV is well below the fitting range, which starts at 0.5 TeV.
The fact that f is close to one (hence N anc ∼ N part ) is due to the fact that, below 100 GeV, the particle production is dominated by soft processes, which scale geometrically, as mentioned in Sect. 5.1.

Centrality classes
After obtaining the simulated distribution of energy deposited in the ECAL, which corresponds only to the hadronic contribution, the distribution can be divided into centrality classes. To determine the ECAL energy boundary values for each class, the simulated distribution is integrated from a value of deposited energy to infinity, until a starting value is found giving a percentage of the total integral. Defining I T as the total integral of the energy distribution, the ECAL energy requirement for any percentage p of centrality, would be the value of E p such that Similarly, as an example, the centrality class (10 − 20)% would correspond to the events depositing an energy E such that E 20 < E < E 10 .

Results
The centrality classification of the MB dataset of PbPb collisions in percentile intervals of 10%, with the cuts in energy obtained as in the previously described procedure, is shown in Fig. 19 as well as the b, N part and N coll distributions for each class obtained from the Glauber MC model. For each class a mean number is estimated for each of the quantities of interest, together with their corresponding standard deviations. The same distributions for the PbNe case can be seen in Fig. 20.
In this way, one can define as many classes as desired and of arbitrary width in percentiles. The values of the geometric quantities for each class, as well as the corresponding energy requirements, can be seen in Table 2 for PbPb and in Table 3 for PbNe. Ten classes are shown for each case.   N coll and b) of PbPb collisions for centrality classes defined from a Glauber MC model fitted to the data. The classes correspond to sharp cuts in the energy deposited in the ECAL. Here σ stands for the standard deviation of the corresponding distributions.   N coll and b) of PbNe collisions for centrality classes defined from a MC Glauber model fitted to the data. The classes correspond to sharp cuts in the energy deposited in the ECAL. Here σ stands for the standard deviation of the corresponding distributions. The PbNe results, when compared to the PbPb case, exhibit a larger uncertainty in the values of the geometrical quantities. This limited precision is an effect of the system size which is much smaller in the case of PbNe, and not a result of the use of the ECAL to estimate centrality. If the centrality classes were defined purely from the Glauber model by applying sharp cuts in the b distribution, the same large overlaps in the N coll and N part distributions would be found, without ever including the ECAL in the procedure.

Centrality % E [ GeV ]
One important caveat is that at low energy the dominating events are of electromagnetic nature or from UPC. Because of this, it is important to exclude in the analyses the energy region where there is a sizeable contamination from these events. If no further selection has been applied to reject UPC events, its contamination in PbPb events will be below 5% at energies higher than 585 GeV, that is at centralities lower than 84% (more central than 84%), and for PbNe at energies higher than 98.9 GeV, that is at centralities lower than 89%.
To determine this threshold, the data are compared to the fitted Glauber MC as in the right plot of Fig. 15. The point from which the two distributions match is found by computing a centred mean of the data over MC ratio around each bin. When this ratio is below a chosen tolerance of 1.05 (meaning 5% contamination of UPC events) for three consecutive bins, the centre of the bin of lower energy is chosen as the energy threshold. 4 If UPC events were identified and rejected, then this limit, of 84% for PbPb or 89% for PbNe, could be increased to include more peripheral events.
The results in the PbPb case are in very good agreement with the results obtained by the ALICE [2,3], ATLAS [4] and CMS [5] collaborations at the same centre-of-mass energy. The PbNe results correspond to the first centrality measurements in fixed-target collisions at the LHC.

Systematic uncertainties
In the given centrality classes with fixed deposited energy boundaries, the uncertainties on the geometric quantities, like the mean values of N part , N coll and b, are assessed. The previously found energy selections are kept, but the systematic uncertainties are quantified on the geometric properties. In what follows, the systematic uncertainties are reported in tables of ten classes of ten percentiles each.

Bin-width dependence
To find the boundary values of the ECAL energy bins, an integration procedure is performed on the histogram of the simulated energy deposition per event. The larger the bins, the less precise the percentile of events in the energy bins can be determined.
A binning scheme of 6000 bins is used, which has an average miss percentage of 0.04% for the PbPb case and 0.02% for PbNe. Since each percentile corresponds to a 1% interval, an average miss of 0.04% means that on average 4% of the events of one percentile "migrate" to the class immediately below in energy. The effect on the geometric uncertainties is estimated and the results are stored for every percentile and N part , N coll and b.

Hadronic cross-section uncertainty
One of the main ingredients for the Glauber MC model is the nucleon-nucleon crosssection. For PbPb at a centre-of-mass energy of √ s NN = 5 TeV, it corresponds to σ inel NN = 67.6 ± 0.6 mb, where the uncertainty comes from the data driven parametrisation described in Ref. [12]. For PbNe at a centre-of-mass energy of √ s NN = 69 GeV, the nucleon-nucleon cross-section corresponds to σ inel NN = 35.4 ± 0.9 mb. To quantify the effect this uncertainty has in the geometric properties, the simulated energy distribution is generated from a Glauber MC simulation made with σ inel NN + 1σ and with σ inel NN − 1σ for each collision system. Then the centrality classes are defined with these new distributions and the effect on the mean values for N part , N coll and b is taken as the associated systematic uncertainty.

Fit uncertainty
One of the most important steps in the process of determining centrality is the choice of the parameters f and µ. From Sect. 5.1, for PbPb there are three best fits that are found, and four for PbNe. These best fits are used to compute the systematic uncertainties due to the choice of a given set of (f, µ) for PbPb and then for PbNe.
In order to compute the uncertainty the Glauber MC energy distribution is generated with the different sets of values. For both PbPb and PbNe separately, the centrality classes are defined for each set of best fits and finally the resulting mean values of the geometric quantities (N part , N coll and b) are compared between the best fits that were not kept and the one that was kept as the definitive best fit.
The precision at which the mean energy deposited per particle is determined in the ECAL does not affect the final result, since any variation of this value would be compensated by the value of µ found. Thus the uncertainty on this energy is absorbed into the uncertainty due the choice of a given set of parameters (f, µ).

NBD uncertainty
The NBD sampling introduces statistical fluctuations that affect directly the observed χ 2 value when comparing the Glauber MC simulation and the data. This effect is noticeable even when looking at the same point of the (f, µ) parameter space. To estimate how these fluctuations affect the final computed geometric quantities, ten simulated energy distributions are simulated with the same best fit parameters and the geometric quantities are computed for all of them.
For each percentile, the standard deviation for N part , N coll and b is computed and used as uncertainty. As before, this is done separately for PbPb and for PbNe.

Total systematic uncertainties
These uncertainties are added together in quadrature to obtain the total uncertainty for each centrality class. The result can be seen in Table 4 for ten classes for the PbPb case,   Table 5 for the PbNe case. The uncertainties on the geometric quantities in both cases are dominated by the systematic uncertainties, as expected. In the PbPb case, the dominant one is the uncertainty due to the binning effect, while in the PbNe case, the dominant one is the uncertainty due to the binning effect in more peripheral collisions (centrality higher than 50%) and the uncertainty due to the hadronic cross-section uncertainty for more central events (centrality lower than 50%).

Conclusions
A procedure to determine the centrality in PbPb collisions at √ s NN = 5 TeV and in PbNe collisions at √ s NN = 69 GeV with the LHCb detector is implemented. The distributions of measured energy deposits in the ECAL are fitted to obtain the parameters of the Glauber model for the simulation. After the fit is performed, the simulated distribution is divided in percentiles, which are delimited by sharp energy boundaries obtained by integrating the distribution. These energy selections allows to classify the data into the same percentiles and subsequently the geometric quantities from the Glauber MC model can be mapped to the real data. The obtained centrality classification is limited to the 84% (89%) most central PbPb (PbNe) events, avoiding the region with large contamination from ultra-peripheral collisions. The correspondence between the results obtained for the PbPb collisions is in good agreement with the results from the ALICE, ATLAS and CMS experiments, and the centrality measurements for the PbNe collisions presented here are the first performed in fixed-target collisions at the LHC.