Centrality Determination in Heavy-Ion Collisions Based on Monte-Carlo Sampling of Spectator Fragments

: The size and evolution of the matter created in a relativistic heavy-ion collision strongly depend on collision geometry, deﬁned by centrality. Experimentally the centrality of collisions can be characterized by the measured multiplicities of the produced particles at midrapidity or by the energy measured in the forward rapidity region, which is sensitive to the spectator fragments. This serves as a proxy for the true collision centrality, as deﬁned by the impact parameter in the models of collisions. In this work, the procedure for centrality determination based on Monte-Carlo sampling of spectator fragments has been proposed. The validity of the procedure has been checked using the fully reconstructed DCM-QGSM-SMM model events and published data from the NA61/SHINE experiment.


Introduction
A deconfined state of quarks and gluons, also known as Quark-Gluon Plasma (QGP) is produced in ultra-relativistic heavy-ion collisions at the Large Hadron Collider (LHC, CERN) and Relativistic Heavy Ion Collider (RHIC, BNL) [1,2]. Collisions at lower energies achievable at the SIS-18 (GSI) [3] or Nuclotron (JINR) [4] form matter with lower temperature and higher baryochemical potential, eventually reaching conditions at which the system is in the hadronic gas phase. The central goal of the existing Beam Energy Scan (BES) programs of the STAR experiment at RHIC ( √ s NN = 3-200 GeV) [5], NA61/SHINE experiment at SPS ( √ s NN = 5.1-17.3 GeV) [6], and BM@N experiment at Nuclotron ( √ s NN = 2.4-3.5 GeV) [4] is to understand the nature of the transition between QGP and hadronic matter. In the future, the MPD experiment at NICA ( √ s NN = 4-11 GeV) [7] and the CBM experiment at FAIR ( √ s NN = 2.7-4.9 GeV) [8] will further explore the phase diagram at a high baryon density region with high statistics data. The size and evolution of the matter created in a relativistic heavy-ion collision strongly depends on the collision geometry, defined by the collision centrality [9,10]. In theory, the collision centrality C b is defined as a fraction (expressed as a percentage) of a total inelastic nucleus-nucleus cross-section, σ AA inel : Here, b is the impact parameter and dσ/db is the differential cross-section of a nucleusnucleus collision. The impact parameter b, defined as the transverse distance between the centers of the two colliding nuclei, is a well-defined quantity and a key input to most theoretical calculations of heavy-ion collisions. However, one cannot directly measure the impact parameter in an experiment [11]. Experimentally, the heavy-ion collisions can be characterized by the measured multiplicity N ch of the produced charged particles around midrapidity or by the energy E sp measured in the forward rapidity region, which is sensitive to the spectator fragments. The centrality procedure is based on the correlation between measured N ch (E sp ) and b, which can be inferred by the comparison of experimental data with simulations of the collisions. High multiplicity events have a low average b and low multiplicity events have a large average b. The former are called central collisions and the latter are called peripheral collisions. The average charged-particle multiplicity N ch is assumed to decrease monotonically with increasing impact parameter b. Usually, the correlation between the impact parameter and the multiplicity N ch is determined using a Monte-Carlo version of the Glauber (MC-Glauber) model based on the eikonal approximation [9,12]. Geometrical properties of the collision, such as the impact parameter b, number of participating nucleons (N part ), number of spectator nucleons (N spec ) and number of binary nucleon-nucleon collisions (N coll ) can be calculated by simulating a large sample of minimum bias collision events. Charged hadron multiplicity can be built based on the MC-Glauber output and the simple phenomenological two-component model based on a negative binomial distribution (NBD) [13][14][15]. The produced particle multiplicity distribution can then be fitted to the experimentally measured one. Centrality classes are defined by cuts on N ch and corresponding parameter values (b, N part , N spec , . . . ) for each class are determined from the MC-Glauber model. The main disadvantages of this approach are limitations of the Glauber model, assumptions about the particle production mechanism, autocorrelations and large systematic uncertainties at low multiplicities [11,16]. Some of the problems can be avoided by using the centrality estimators, which are independent of particle production in the participant's zone: sum of energy, charge, or multiplicity of spectators. This approach might be beneficial for experimental studies that search for a critical point in the QCD phase diagram. Such measurements usually focus on baryon multiplicity fluctuations at midrapidity; hence, selecting events based on particles detected at forward rapidity may reduce the nontrivial autocorrelations [16,17]. The spectators are protons and neutrons as well as larger spectator fragments that continue to propagate in the same direction as the colliding ions before the interaction. The detection of all spectator nucleons and nuclear fragments in forward calorimeters would be the most straightforward method of estimating b from the total volume of spectator matter. Unfortunately, the number of event generators for nucleus-nucleus collisions that include a reliable model to describe the formation of nuclear fragments is very limited.
In this work, the new data-driven approach for centrality determination based on Monte-Carlo sampling of spectator fragments has been proposed. It is based on the output from the Monte-Carlo version of the Glauber model and hybrid model DCM-QGSM-SMM [18][19][20], which provides a reasonable treatment of spectator fragments. The validity of the procedure has been checked using the fully reconstructed DCM-QGSM-SMM model events and published data from the NA61/SHINE experiment (SPS, CERN) for Pb+Pb collisions at beam momentum p beam = 13 A GeV/c ( √ s NN = 5.12 GeV). The resulting values of the mean impact parameter b have been compared to the results of the MC-Glauber approach based on the multiplicity. The results would be useful for the upcoming Beam Energy Scan experiments: BM@N at Nuclotron, CBM at FAIR and MPD at NICA Collider. The paper is organized as follows. Section 2 briefly describes the used MC-Glauber and DCM-QGSM-SMM models, while Section 3 discusses the centrality determination from the multiplicity of charged particles based on the MC-Glauber approach. Section 4 discusses the proposed approach for centrality determination based on Monte-Carlo sampling of spectator fragments. The results of the check of the validity of the procedure are presented and discussed in Section 5. Finally, the summary and future outlook are given in Section 6.

Short Description of the Models
In this section, we briefly discuss the main features of the selected versions of the DCM-QGSM-SMM and MC-Glauber models, used in the present work.
DCM-QGSM-SMM is a hybrid heavy-ion event generator based on the Dubna Cascade Model (DCM) [18], the Quark-Gluon String Model (QGSM) and the Statistical Multifragmentation Model (SMM) [19,20]. The Dubna Cascade Model, DCM, is based on the Monte-Carlo solution of a set of the Boltzmann-Uehling-Uhlenbeck (BUU) relativistic kinetic equations with collision terms, including cascade-cascade interactions. For particle energies below 1 GeV, it considers only nucleons, pions and deltas. The model includes a proper description of pion and baryon dynamics for particle production and absorption processes. To make the DCM code applicable at higher energies (up to hundreds GeV/nucleon), it was merged with the Quark-Gluon String Model (QGSM). QGSM is used to describe elementary hadron collisions at energies higher than 5 GeV in the framework of independent quark-gluon strings. The QGSM considers the two lowest SU(3) multiplets in mesonic, baryonic and antibaryonic sectors, so interactions between almost 70 hadron species are treated on the same footing. The production of nuclear fragments is subdivided into three stages: (1) a dynamical stage leading to the formation of an equilibrated nuclear system, which is described by DCM; (2) disassembly of the system into individual primary fragments described by the Statistical Multifragmentation Model (SMM); (3) de-excitation of hot primary fragments according to evaporation/fission processes [20].
The purpose of Monte-Carlo implementations of the Glauber model is to compose two nuclei out of nucleons and simulate their collision process event by event. In the present work, we use the 3.2 version of the PHOBOS MC-Glauber model as described in refs. [9,12]. In this approach, the following steps are implemented: 1.
The nucleon position in each colliding nucleus is determined by the radial nuclear density function ρ(r), modeled by the modified Woods-Saxon distribution: where R is the radius of the nucleus (R = R Pb = 6.62 ± 0.06 fm for 208 Pb nucleus), ρ 0 is density in the center of the nucleus; a is the skin thickness of the nucleus, which defines how quickly the nuclear density falls off near the edge of the nucleus (a = 0.546 ± 0.010 fm) [13]. The additional parameter w is needed to describe nuclei whose maximum density is reached at radii r > 0 (w = 0 for Pb). The parameters are based on data from low-energy electron-nucleus scattering experiments. The radial coordinate of a nucleon is randomly selected from the distribution 4πr 2 ρ(r) with the condition, that no pair of nucleons inside the nucleus has a distance less than 0.4 fm.

2.
The impact parameter b of a given collision is selected randomly from the geometrical distribution dP/db The nucleus-nucleus collision is treated as a sequence of independent binary nucleonnucleon collisions, where the nucleons travel on straight-line trajectories and the inelastic nucleon-nucleon cross-section σ inel NN is assumed to be independent of the number of collisions and depends only on the collision energy. Two nucleons from different nuclei are assumed to collide if the relative transverse distance d between centers is less than the distance corresponding to the inelastic nucleon-nucleon crosssection: d < σ inel NN /π. For Pb+Pb collisions at beam momentum p beam = 13 A GeV/c ( √ s NN = 5.12 GeV), we use σ inel NN = 29.9 ± 3 mb [21]. Geometrical properties of the collision, such as the impact parameter b, number of participating nucleons (N part ), number of spectator nucleons (N spec ) and number of binary nucleon-nucleon collisions (N coll ), are calculated by simulating many nucleus-nucleus collisions. The relation between N part and the number of nucleons N spec remaining in the spectator fragments does not depend on the process of their formation: where A is the mass number of a colliding ion and the sum runs over all spectator fragments A i sp on both sides of the collision point.

Centrality Determination from the Multiplicity of Charged Particles
Below, we briefly discuss the well-established procedure of centrality (impact parameter) determination based on the charged hadron multiplicity and the MC-Glauber approach. In order to illustrate it, we used the DCM-QGSM-SMM model to simulate around 100 k minimum bias Pb+Pb collision events at p beam = 13A GeV/c. At the next step, the sample of events was made as an input for the full chain of realistic simulations of the NA61/SHINE detector subsystems based on the GEANT4 platform and realistic reconstruction algorithms [22]. The fully reconstructed events were used to generate the distribution of the multiplicity N ch of the produced charged particles with −0.9 < η < 3.1, see the left panel in Figure 1. In the MC-Glauber approach, the output from the 2M minimum-bias MC-Glauber Pb+Pb events was used to build the multiplicity of the produced particles N MC . The multiplicity N MC (N a , f , µ, k) is modeled as a sum of particles produced from a set of N a independent emitting sources ("ancestors") [13,14,23]. Each ancestor produces particles according to a negative binomial distribution (NBD) P µ,k with mean multiplicity per ancestor µ and width parameter k: The N a ( f ) parameterization is inspired by the phenomenological two-component model of Kharzeev and Nardi [15] for "soft" plus "hard" particle production processes. In this approach, the soft interactions contribute to the multiplicity dependence as N part , while hard interactions as N coll . The multiplicity distribution N MC is simulated for an ensemble of events and various values of the NBD parameters µ, k and the N a parameter f [13,14,23]. A minimization procedure is applied to find the optimal set of parameters that result in the smallest fitting criteria χ 2 . The MC-Glauber fit was performed for multiplicities above 40. A value of the multiplicity at which the fit starts to deviate from the multiplicity distribution defines the so-called "anchor point" below which the centrality determination is not reliable. The resulting fit parameters are: χ 2 /nd f = 1.26, µ = 0.87, k = 22, f = 1. The minimum of the χ 2 corresponds to a small contribution of hard processes ( f = 1). The resulting N MC multiplicity from the MC-Glauber fit is shown by blue solid triangles in the left panel of Figure 1. With the final set of parameters ( f , µ, k), the mean value of impact parameter b can be extracted for the centrality classes defined by the sharp cuts in the multiplicity distribution, see the solid vertical lines in the left panel of Figure 1. The right panel of Figure 1 shows the centrality dependence of the mean value of impact parameter b and its standard deviation extracted directly from the DCM-QGSM-SMM model (open symbols) and from the MC-Glauber approach (closed symbols). Overall good agreement is observed.
Centrality determination from the multiplicity of charged particles based on the MC-Glauber approach is a very well-established procedure, which was checked by different experiments and model calculations for heavy-ion (Au+Au and Pb+Pb) collisions from √ s NN = 2.2 GeV to 2.76 TeV [13,14,23]. Figure 2 shows a flowchart of different centrality determination approaches. The right-hand side shows the procedure based on the MC-Glauber model and NBD for sampling the multiplicity of produced particles, as described in this section.

Centrality Determination with Monte-Carlo Sampling of Spectator Fragments
The goal of the developed procedure is to build the path from geometry parameters from the output of the Monte-Glauber Model to the observable related to the spectator matter. In order to do this, several steps are needed to model the fragmentation process of nuclear remnants, to determine the parameters of spectator fragments and to model the detector's effects. The main steps of the developed procedure based on the Monte-Carlo sampling of spectator fragments are described below. The flowchart of this procedure is shown with the red boxes on the left-hand side of Figure 2. The main steps are the following:

0.
To provide a simple description of the fragmentation process, prepare a set of twodimensional maps (distributions) using a model with the formation of spectator fragments (e.g., hybrid DCM-QGSM-SMM model with statistical multifragmentation SMM [18][19][20] or dynamical multifragmentation implemented within the PHQMD model [24]): To provide a simple description of the detector effects, prepare a map (distribution) of spectator detector response (S f rag ) as a function of rapidity (y f rag : S f rag ) and energy (E f rag : S f rag ) of the spectator fragments (can be achieved, e.g., using a large sample of events from the full chain of realistic simulations of spectator detector based on the GEANT4 platform and realistic reconstruction algorithms). 1.
To model the initial geometry of the collisions, generate a large sample of MC-Glauber minimum bias events in order to obtain the information about the geometrical properties of collisions: impact parameter b and number of spectators N spec .

2.
To sample a total mass number of nuclei remnants A tot , use the map A tot : b. 3.
To model the process of spectator fragment formation based on the output of any generator with realistic fragmentation implication (step 0), generate an array of fragments (N A f rag [A f rag ]) using the following iterations: (a) Initialize all array elements with zeros: Use the map A f rag : N A f rag to generate a pair (A f rag , N A f rag ); (c) Verify if: • N A f rag [A f rag ] > 0, then redo the step (b); •

then go to the step (a); (d)
Write a pair (A f rag , N A f rag ), if checks above were passed;

then go to step (b). As a result of steps (a)-(e), a sample of
To generate values of energy and rapidity for each fragment in the event A 1 . . . A N f rag , use the maps A f rag : E f rag and A f rag : y f rag .

5.
To simulate a detector response S f rag (deposited energy of fragments for the hadronic calorimeter or charge for the scintillator detector, etc.), use the map from step 0. 6.
Calculate the total measured energy (or charge) of spectators, S tot .
As a result of steps (1)-(6) the centrality estimator based on the MC sampling of spectator fragments will be constructed. In order to optimize the description of the experimentally measured distribution, the maps A tot : b, A f rag : y f rag , A f rag : N A f rag , A f rag : N A f rag , y f rag : S f rag , E f rag : S f rag can be parameterized with a set of parameters, which can be later used to fit the experimental data. Then, iterating over all values of parameters, the optimal fit parameters can be found based on a minimum of χ 2 between the simulated and measured distributions of the total energy of spectator fragments or the total charge. Centrality classes are then determined using the distribution for the best-fit parameters. This provides a mapping between geometry parameters and centrality estimator.

Simplified Procedure with Gaussian Approximation for the Energy of Spectators
The procedure described in the previous section can be simplified based on the assumption that the distribution of the energy of spectators is Gaussian. This follows from the Monte-Carlo simulations based on the DCM-QGSM-SMM model, which realistically model the process of fragmentation of spectator fragments (see [18][19][20]). According to Figure 3, the energy distribution of a single spectator nucleon could be approximated with a Gaussian function. Distributions of heavier fragments are also close to Gaussian with mean values that are close to the product of the beam energy E beam and mass number of fragments A f rag . Hence, the energy of spectator fragments can be treated as a combination of energies of single spectator nucleons distributed according to the common Gaussian function: Here, P(E; µ, k) is the Gaussian function of energy E with mean µ and width k. Similar considerations can be applied for the detectors with Gaussian response to the energy of the fragments, e.g., hadronic calorimeters. According to experimental data obtained with hadronic calorimeters [25][26][27], measured detector signal S of the detected particle can be parameterized using the Gaussian distribution, which is proportional to the energy of the detected particle E and width σ S given by the following equation: Here, a = 0.56 GeV 1/2 -stochastic term, b = 0.02-constant term, c = 0.16 GeV 1/4 -leakage term.
The energy of spectator fragments (or single spectators) is distributed according to the convolution of two Gaussian distributions (which should also be a Gaussian distribution): the first one to sample the real energy of fragment E f rag (or single spectators E spec ) and the second one to sample the detector's response. Based on the considerations above, the following procedure, shown by the green box in Figure 2, is investigated. The energy of spectators S spec is sampled according to the Gaussian distribution G(µ, k) with free parameters µ and k. The measured detector signal, which is proportional to the energy of all detected spectators is calculated as S tot = f + ∑ N spec i=1 S i spec . The parameter f describes the offset at low energies due to the contribution of produced particles (participants) in the forward/backward rapidity region. The main limitation of this procedure is that the correlations due to the fragmentation of nuclei remnants are not taken into account. Parameters f , µ and k are fixed by searching for a minimum of χ 2 value between the simulated MC-Glauber energy distribution and the experimental one.

Results and Discussion
The validity of the procedure of centrality (impact parameter) determination based on the Monte-Carlo sampling of spectator fragments has been checked using the fully reconstructed DCM-QGSM-SMM model events and published data from the fixed target NA61/SHINE experiment (SPS, CERN) for Pb-Pb collisions at beam momentum p beam = 13 A GeV/c [22]. The NA61/SHINE fixed-target experiment uses a large acceptance hadron spectrometer located in the North Area H2 beamline of the CERN Super Proton Synchrotron accelerator [28]. A schematic view of the NA61/SHINE experiment is shown in the left part of Figure 4. The main detectors are five time projection chambers (TPC): two vertex TPCs mounted in the magnets, a gap TPC covering the beam area and two main TPCs (left and right) measuring tracks of charged particles. The Projectile Spectator Detector (PSD) is a segmented forward hadron calorimeter with the transverse and longitudinal segmentations assembled of sampling lead/scintillators modules [26]. The transverse layout of the PSD detector is shown in the right part of Figure 4. The central part of the PSD consists of 16 modules with 10 cm × 10 cm transverse size and the outer region has modules with 20 cm × 20 cm size. PSD is installed 18 m downstream of the target and covers the polar angle of the projectile nuclei fragments emission θ < 3.8 • . The acceptance of the PSD allows us to detect projectile spectators and produced particles emitted in forward rapidity region [26,27]. In order to test the developed procedure of centrality determination, a sample of 150 k minimum bias events was generated with the DCM-QGSM-SMM model and reconstructed with the GEANT4 transport package. The energy distribution of spectators was sampled from 2M minimum bias collisions simulated with the MC-Glauber model. Figure 5 shows the distribution of the total energy of spectator fragments E PSD in the PSD calorimeter (open boxes) for the fully reconstructed DCM-QGSM-SMM model events (left part) and NA61/SHINE experimental data (right part). The data from NA61/SHINE experiment are taken from [22]. The results of the approximation of PSD energy distribution using a Monte-Carlo sampling of spectator fragments with Gaussian approximation for the energy of spectators (see Section 4.1) are shown by closed triangles. With the final set of parameters the mean value of impact parameter b can be extracted for the centrality classes defined by the sharp cuts in the E PSD distribution, see the solid vertical lines in Figure 5. The events with a small value of E PSD have a small average b (central collisions) and events with a high value of E PSD have a large average b (peripheral collisions). In the case of fully reconstructed DCM-QGSM-SMM model events, the MC-Glauber fit for energy describes the E PSD distribution quite well, excluding very central collisions. This needs to be addressed by future studies. The situation is the opposite in the case of the real data. Here, the MC-Glauber fit for energy does not describe well the E PSD distribution for very peripheral collisions, see the right part of Figure 5. This difference is expected to be reduced in the future by accounting for effects due to spectator fragmentation in the full procedure described in Section 4.  Figure 6 shows the correlation between impact parameter b and the total energy in the PSD calorimeter E PSD for fully reconstructed DCM-QGSM-SMM model events and detector response obtained with the GEANT4 package (left part) and the results of the Monte-Carlo procedure based on Gaussian approximation for the spectator's energy (right part). Red lines show the regions used to map impact parameter b and energy E PSD in a certain centrality class. The shapes of the two distributions are very similar, except for very peripheral collisions, where it is broader for the detailed simulation with Geant4 and the DCM-QGSM-SMM model. This is due to a missing simulation of the fragmentation process in the simplified procedure. Figure 7 presents the resulting centrality dependence of the mean value of the impact parameter b for the MC-Glauber approach based on the multiplicity of the produced charged particles (left part) and MC-Glauber approach with Gaussian approximation for the energy of spectators (right part). The resulting values of b extracted from the Monte-Carlo approaches (closed symbols) are compared with the values used in the fully reconstructed DCM-QGSM-SMM model events and detector response obtained with the GEANT4 package (open symbols). The b values extracted from the MC-Glauber approach based on the multiplicity are in much better agreement with DCM-QGSM-SMM model results. For the case of impact parameter distributions obtained based on centrality classes determined by the energy distribution of spectators, there is a discrepancy, especially in the most central classes. This is again due to the bad quality of the fit in the region of the most central events and this needs to be addressed by future studies. At the same time, the width of the centrality classes determined using the PSD energy is slightly larger than in the case of track multiplicity. This is due to the fact that the width of the two-dimensional distribution of the impact parameter b vs energy E PSD is larger than the width of the distribution of the impact parameter b and multiplicity N ch . As can be seen in the comparison of open squares from the left and right plots of Figure 7, the centrality classes determined separately using the multiplicity of produced particles and spectators are slightly different. This is again because of the different shapes of two-dimensional distributions of impact parameters and corresponding centrality estimators. In order to understand the impact of this effect, an event-by-event comparison of centrality determined based on different centrality estimators should be performed during future investigations. As part of further work, it is planned to consider options for improving the fit of energy. For this, the applicability of Landau, Gamma and other distributions will be tested. Further, the developed methods can be applied to other heavy-ion experiments, such as MPD at NICA and BM@N at Nuclotron. These experiments are going to use the same type of forward hadronic calorimeters, similar to the PSD from the NA61/SHINE experiment.

Conclusions
In summary, a new procedure for centrality determination in heavy-ion collisions based on the Monte-Carlo sampling of spectator fragments is proposed. The connection between the averaged impact parameter and centrality classes is extracted using the energy of spectator fragments measured in the forward rapidity region by hadronic calorimeter. A simplified procedure based on Gaussian approximation for the energy of spectators has been implemented. The validity of the procedure has been checked using the distribution of energy of spectator fragments in the PSD forward calorimeter of the NA61/SHINE experiment for Pb+Pb collisions at beam momentum p beam = 13 A GeV/c ( √ s NN = 5.12 GeV) as well as for reconstructed DCM-QGSM-SMM model events. The comparison with the results of the standard centrality determination procedure based on produced particle multiplicity and the MC-Glauber approach is also provided. In the future, we plan to improve the procedure by taking into account the details of the spectator fragmentation and extending its application for the forward hadronic calorimeters of MPD at NICA and BM@N at Nuclotron heavy-ion experiments [4,7].