Improved version of the PHOBOS Glauber Monte Carlo

# Improved version of the PHOBOS Glauber Monte Carlo

## Abstract

“Glauber” models are used to calculate geometric quantities in the initial state of heavy ion collisions, such as impact parameter, number of participating nucleons and initial eccentricity. Experimental heavy-ion collaborations, in particular at RHIC and LHC, use Glauber Model calculations for various geometric observables for determination of the collision centrality. In this document, we describe the assumptions inherent to the approach, and provide an updated implementation (v2) of the Monte Carlo based Glauber Model calculation, which originally was used by the PHOBOS collaboration. The main improvements w.r.t. the earlier version (v1) (1) are the inclusion of Tritium, Helium-3, and Uranium, as well as the treatment of deformed nuclei and Glauber-Gribov fluctuations of the proton in p+A collisions. A users’ guide (updated to reflect changes in v2) is provided for running various calculations.

Revisions && changes of the arXiv document and code:
v1
, 11 Aug 2014: initial document, code v2.0
v2, 20 Dec 2014: updated Table 1 with “p” and “PbHN”, code v2.1 (fixes bug for Hulthen)
v3, 01 Sep 2016: published version, fixed Eq. 3, code v2.3 (added , and various proton profiles
v3, 01 Sep 2016: from Pythia, changes in calculation of )
v4, 07 Jun 2017: updated Table 1 with “Xe” and corrected typo in proton radius, code v2.4
v4, 07 Jun 2017: (fixed drawing function)
v5, 16 Oct 2017: updated Table 1 with “P*”
v6, 13 Nov 2017: updated Table 1 with “Xes”
v7, 01 Mar 2018: split previous Table 1 into two tables: Table 1 for spherical and Table 2 for
v4, 01 Mar 2018: deformed nuclei parameters, added “Xe2” and “Xe2a” to include deformation of
v7, 01 Mar 2018: Xenon, code v2.5 (included core/corona calculation, see “Npart0”)
v8, 30 Apr 2018: fixed typo in , fixed in “Si2”, added “Al”, describes new function
v8, 25 Apr 2018: runAndCalcDens, latest code v2.6 (includes fixes and “runAndCalcDens”)

Glauber Monte Carlo model, participant, eccentricity, centrality

## I Introduction

In heavy-ion collisions, initial geometric quantities such as impact parameter and shape of the collision region cannot be directly determined experimentally. However, it is possible to relate the number of observed particles and number of spectator neutrons to the centrality of the collision. Using the percentile centrality of a collision, the initial geometric configuration can be estimated with models assuming the configuration nuclei from nucleons and the collision process of two nuclei.

These models fall in two main classes. (For a review, see LABEL:glaubreview). In the so called “optical” Glauber calculations, a smooth matter density is assumed, typically described by a Fermi distribution in the radial direction and uniform over solid angle. In the Monte Carlo (MC) based models, individual nucleons are stochastically distributed event-by-event and collision properties are calculated by averaging over multiple events. As discussed in LABEL:glaubreview and LABEL:eccentricity, the two type of models lead to mostly similar results for simple quantities such as the number of participating nucleons () and impact parameters (), but give different results in quantities where event-by-event fluctuations are significant, such as participant frame eccentricity ().

In this paper, we discuss several improvements of the Monte Carlo Glauber calculation originally implemented by the PHOBOS collaboration (1). The main improvements w.r.t. the earlier version (v1) are the inclusion of Tritium, Helium-3, and Uranium, as well as the treatment of deformed nuclei and Glauber-Gribov fluctuations of the proton in p+A collisions. Another widely used implementation of a MC Glauber calculation can be found in LABEL:Rybczynski:2013yba. In section II, the MC approach is outlined and the assumptions that go into the calculation are introduced. In section III, we discuss the implementation and the tutorial functions provided. Section IV summarizes the paper. The latest version 2 code described in this document is v2.6. For further improvements, mainly focusing on nuclear collisions involving Pb at high collision energies, see LABEL:Loizides:2017ack.

## Ii The model

The MC Glauber model calculation is performed in two steps. First, the nucleon positions in each nucleus are stochastically determined. Then, the two nuclei are “collided”, assuming the nucleons travel in a straight line along the beam axis (eikonal approximation) such that nucleons are tagged as wounded (participating) or spectator.

### ii.1 Makeup of nuclei

The position of each nucleon in the nucleus is determined according to a probability density function. In a quantum mechanical picture, the probability density function can be thought of as the single-particle probability density and the position as the result of a position measurement. In the context of the MC Glauber model, one commonly requires a minimum inter-nucleon separation () between the centers of the nucleons for the determination of the nucleon positions inside a given nucleus.

The probability distribution for spherical nuclei is taken to be uniform in azimuthal and polar angles. The radial probability function is modeled from nuclear charge densities extracted in low-energy electron scattering experiments (6). The nuclear charge density is usually parameterized by a Fermi distribution with three parameters

 ρ(r)=ρ01+w(r/R)21+exp(r−Ra), (1)

where is the nucleon density, is the nuclear radius, is the skin depth and corresponds to deviations from a spherical shape. The overall normalization () is not relevant for this calculation. Values of the other parameters used for different nuclei are listed in Table 1.

Exceptions from Eq. 1 are the Deuteron (H), Tritium (H3, H), Helium-3 (He3, He), and Sulfur (S) nuclei. For Sulfur, a three parameter Gaussian form is used

 ρ(r)=ρ01+w(r/R)21+exp(r2−R2a2). (2)

The values of , and for Sulfur are also given in Table 1. For Deuteron, three options are supported with the option being the one used in PHOBOS analyses:

1. The three parameter Fermi form with the values given in Table 1.

2. The Hulthén form

 ρ(r)=ρ0(e−ar−e−brr)2, (3)

where  fm and  fm (8); (9); (10).

3. The proton is taken from the Hulthén form given above with the neutron placed opposite to it.

For H and He, the configurations were computed (and stored in a database) from Green’s function MC calculations using the AV18 + UIX model interactions, which correctly sample the position of the three nucleons, including correlations, as in LABEL:Nagle:2013lja. Finally, note that we also provide the rescaled values for the radius and skin depth of Cu (“CuHN”), Au (“AuHN”) and Pb (“PbHN”) to take into account the finite nucleon profile as derived in (12) and (13). For deformed nuclei we use

 ρ(x,y,z)=ρ011+exp(r−R(1+β2Y20+β4Y40))a, (4)

where , with the deformation parameters and taken from Ref. (15). The values used for different nuclei (Si, Cu, Au and U) are listed in Table 2.

### ii.2 Collision Process

The impact parameter of the collision is chosen randomly from a distribution up to some large maximum with fm. The centers of the nuclei are calculated and shifted to and . The reaction plane, defined by the impact parameter and the beam direction, is given by the - and -axes, while the transverse plane is given by the - and -axes. It is assumed that the nucleons move along a straight trajectory along the beam axis. The longitudinal coordinate does not play a role in the calculation.

The inelastic nucleon-nucleon cross section (), which is only a function of the collision energy, is extracted from p+p collisions. A value of  mb is used at top RHIC energy of  GeV, while at the LHC  mb for  TeV and  mb for  TeV are used (21). When estimating the systematic uncertainties on calculated quantities like one varies the  (along with the other parameters) by typically 3 mb and 5 mb at RHIC and LHC, respectively.

The “ball diameter” defined as

 D=√σNN/π (5)

parameterizes the interaction strength of two nucleons: Two nucleons from different nuclei are assumed to collide if their relative transverse distance is less than the ball diameter. If no such nucleon–nucleon collision is registered for any pair of nucleons, then no nucleus–nucleus collision occurred. Counters for determination of the total (geometric) cross section are updated accordingly.

## Iii Users’ guide

The MC Glauber code works within the ROOT framework (ROOT 4.00/08 or higher (22)). The code is contained in the macro runglauber_vX.Y.C (23). The version described here is v2.6. Three classes, TGlauNucleon, TGlauNucleus and TGlauberMC and three functions (macros) runAndSaveNtuple(), runAndSaveNucleons(), and runAndSmearNtuple() are defined in the provided macro. Note that for generating events with H or He you will need the text files called “h3_plaintext.dat” or “he3_plaintext.dat” in the current path. While the functionality is essentially complete for known applications of the Glauber approach, users are encouraged to write their own functions to access results of the Glauber simulation or to modify the code.

The main classes are the following:

• TGlauNucleon is used to store information about a single nucleon. The stored quantities are the position of the nucleon, the number of binary collisions that the nucleon has had and which nucleus the nucleon is in, “A” or “B”. For every simulated event, the user can obtain an array containing all nucleons (via TGlauberMC::GetNucleons()).

• TGlauNucleus is used to generate and store information about a single nucleus. The user is not expected to interact with this class.

• TGlauberMC is the main steering class used to generate events and calculate event-by-event quantities such as the number of participating nucleons.

The steering class TGlauberMC has one constructor

TGlauberMC::TGlauberMC(Text_t* NA,
Text_t* NB,
Double_t xsect,
Double_t xsectsigma)


where NA and NB are the names of the colliding nuclei and xsect and xsectsigma are the nucleon-nucleon cross section mean and width given in mb. The defined nuclei names are given in Tables 1 and 2. For Deuteron, the names “d”, “dhh” and “dh” correspond to the three options described in section II.1 respectively. Units are generally given in fm for distances, while in mb for cross sections. If a finite value () for the width of the cross section is given, the calculation will use the Glauber-Gribov model for fluctuations of the proton. This option currently is only implemented for p+A collisions. In this case, following (24) as described in (25) will be distributed according to with

 P=ρσσ+σ0exp−(σ−σ0σ0Ω)2 (6)

where denotes the mean and the width. The normalization and the rescaling parameter are computed from the provided input (mean and ) making use of with for . Example distributions for p+Pb collisions at the LHC are given in Fig. 2.

### iii.1 Running the Code

To generate Au+Au collisions at GeV ( mb) one has to construct a TGlauberMC object by issuing at the “root” prompt the commands:

root [0] gSystem->Load("libMathMore")
root [1] .L runglauber_X.Y.C+
root [2] TGlauberMC glauber("Au","Au",42);


where the first ROOT command loads the needed library libMathMore, the second compiles, links and loads the compiled macro (where you have to replace X.Y with the current version number of the code, for example v2.6) including the Glauber code as explained in chapter 2 of the ROOT users’ guide, while the third commands sets up the Glauber simulation for Au+Au collisions at  GeV.

Events can be generated interactively using the two functions

• TGlauberMC::NextEvent(Double_t bgen), which is used to run an event at a specified impact parameter, or over a range of impact parameters (if bgen=-1, the default value) as described in section II.2).

• TGlauberMC::Run(Int_t nevents) which is used to run a large event sample by invoking NextEvent many times.

Other important public member functions are:

• TGlauberMC::SetMinDistance(Double_t d), which is used to set minimum nucleon separation within a nucleus, (default is  fm)

• TGlauberMC::SetBmin(Double_t bmin) and TGlauberMC::SetBmax(Double_t bmax), which can be used to set the range of impact parameter values generated in Run().

• TGlauberMC::GetTotXSect() which returns the total nucleus-nucleus cross section, calculated when the function Run() is called.

• TGlauberMC::Draw() which draws the current event in the current pad.

### iii.2 Provided functions

Three functions (macros) are provided to demonstrate how to run the model:

#### runAndSaveNtuple

The macro runAndSaveNtuple() generates a number of Monte Carlo events and saves some event-by-event quantities. It takes as parameters, the number of events to be generated, the collision system, the nucleon-nucleon cross section, its width to simulate Glauber-Gribov fluctuations (for p+A), the minimum separation distance and the output file name. It creates and stores an ntuple in the output file with the following event-by-event quantities:

• Npart: Number of participating nucleons.

• Npart0: Number of singly-wounded participating nucleons.

• Ncoll: Number of binary collisions.

• B: Generated impact parameter.

• BNN: Mean nucleon–nucleon impact parameter

• MeanX: Mean of for wounded nucleons, .

• MeanY: Mean of for wounded nucleons, .

• MeanX2: Mean of for wounded nucleons, .

• MeanY2: Mean of for wounded nucleons, .

• MeanXY: Mean of for wounded nucleons, .

• VarX: Variance of for wounded nucleons, .

• VarY: Variance of for wounded nucleons, .

• VarXY: Covariance of and for wounded nucleons, .

• MeanXSystem: Mean of for all nucleons.

• MeanYSystem: Mean of for all nucleons.

• MeanXA: Mean of for nucleons in nucleus A.

• MeanYA: Mean of for nucleons in nucleus A.

• MeanXB: Mean of for nucleons in nucleus B.

• MeanYB: Mean of for nucleons in nucleus B.

• NpartA: from nucleus A.

• NpartB: from nucleus B.

• PhiA: Azimuthal angle for nucleus A.

• ThetaA: Polar angle for nucleus A.

• PhiB: Azimuthal angle for nucleus B.

• ThetaB: Polar angle for nucleus B.

• PsiN: Event plane angle of Nth harmonic.

• EccN: Participant eccentricity for Nth harmonic.

It is important to note that for each of these event-by-event quantities a “getter” function is implemented providing the users the option to write their own event loop (using TGlauberMC::NextEvent()).

#### runAndSaveNucleons

The macro runAndSaveNucleons() generates a number of Monte Carlo events and saves an array of TGlauNucleon objects for each event. It is also possible to use this function to print out the values stored in the nucleons by setting the verbosity parameter. The function takes as parameters the number of events to be generated, the collision system, the nucleon-nucleon cross section, the minimum separation distance, the verbosity flag and the output file name.

#### runAndSmearNtuple

The macro runAndSmearNtuple() generates a number of Monte Carlo events and saves a reduced set of event-by-event quantities, where the position of the nucleon is smeared similar to what is done in Ref. (11). It takes the same arguments as the runAndSaveNtuple() macro.

#### runAndCalcDens

The macro runAndCalcDens() generates a number of Monte Carlo events and saves for every event a normalized density profile in and calculated from smeared participant and binary collision positions with a relative weight. As suitable input for hydrodynamical calculations, a scale factor needs to be applied to adjust to the multiplicity as desribed in Ref. (11). The function takes as parameters the number of events to be generated, the weight, the collision system, the nucleon-nucleon cross section, the minimum separation distance, and the output file name.

### iii.3 Sample Results

To provide example applications of this code, 100k events were generated for Cu+Cu, Au+Au and U+U at RHIC energies ( mb), and Pb+Pb at LHC beam energy ( mb) using the runAndSaveNtuple() function. The resulting ntuples were used to plot the distributions of and , shown in Fig. 3.

Using the event-by-event quantities, one can construct combinations of moments like (3):

• Reaction-plane eccentricity

 ϵRP=VarY−VarXVarY+% VarX (7)
• Participant-plane eccentricity

 ϵpart=√(VarY−VarX)2+4% VarXY2VarY+VarX (8)

which are shown in Fig. 4 for different systems.

## Iv Summary

Glauber Model calculations are commonly used by heavy-ion physics experiments to study the initial state configurations of nuclear matter. This document describes the updated implementation (v2) of the Monte Carlo based Glauber Model calculation, which originally was used by the PHOBOS collaboration. The main improvement w.r.t. the earlier version (1) are the inclusion of Tritium, Helium-3, and Uranium, as well as the treatment of deformed nuclei and Glauber-Gribov fluctuations of the proton in p+A collisions. The code, accessible online, can be used within user code or in a standalone mode allowing analysis of various distributions. The authors welcome comments on the code and suggestions on how to make it more useful to both experimentalists and theorists.

### Footnotes

1. An exponential is used with  fm based on the rms charge radius of the proton (14). Also single (“pg”) and double Gaussian (“pdg”) profiles are implemented (16). See code.
2. As explained in the text, also two further parameterizations (“dh” and “dhh”) for the Deuterium (H) are implemented. Please refer to the code.
3. The configurations for H and He were obtained from Green’s function MC calculations as in LABEL:Nagle:2013lja.
4. Parameters for Xe and fm from LABEL:Tsukada:2017llu were used. The radius is scaled down by and was reduced by fm to symmetrize the uncertainty and to approximate the smaller neutron skin of Xe.
5. Parameters for Sb (Antimony) and fm from LABEL:atomicdata were used. Natural Sb has two isotopes: (%) and (%). The radius of Sb is scaled up by using the average of . The resulting parameters are consistent with those obtained from Xe.
6. These values are usually also used for Pb, since the Bessel-Fourier coefficients for the two nuclei are similar (6). Indeed, an earlier publication (7) gives and  fm for Pb (“Pb*”), consistent within uncertainties with the values for Pb.
7. Same parameters as for Xe with and from LABEL:xedef.
8. Same parameters as for Xe with from interpolation between measured deformation parameters for the even-A Xe isotopes (19), and set to zero.
9. Implementation as done in LABEL:Heinz:2004ir.

### References

1. B. Alver, M. Baker, C. Loizides and P. Steinberg (2008), arXiv:0805.4411.
2. M.L. Miller, K. Reygers, S.J. Sanders, and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57 (2007) 205, arXiv:nucl-ex/0701025.
3. B.Alver et al. [PHOBOS Collaboration], Phys. Rev. C 77 (2008) 014906, arXiv:0711.3724.
4. M. Rybczynski, G. Stefanek, W. Broniowski and P. Bozek, Comput. Phys. Commun. 185, 1759 (2014), arXiv:1310.5475.
5. C. Loizides, J. Kamin and D. d’Enterria, arXiv:1710.07098.
6. H. De Vries, C.W. De Jager, and C. De Vries, Atom. Data Nucl. Data Tabl. 36 495 (1987).
7. C. W. De Jager, H. De Vries and C. De Vries, Atom. Data Nucl. Data Tabl. 14, 479 (1974).
8. L. Hulthen and M. Sugawara, Handbuch der Physik 39 1 (1957).
9. S.S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 91 072303 (2003), arXiv:nucl-ex/0306021.
10. S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. C 74, 024904 (2006), arXiv:nucl-ex/0603010.
11. J. L. Nagle et al., Phys. Rev. Lett. 113 (2014) 112301, arXiv:1312.4565 [nucl-th].
12. T. Hirano and Y. Nara, Phys. Rev. C 79, 064904 (2009), arXiv:0904.4080.
13. Q. Y. Shou et al. Phys. Lett. B 749 (2015) 215, arXiv:1409.8375 [nucl-th].
14. R. Hofstadter et al., Rev. Mod. Phys. 28 (1956) 214.
15. P. Möller, J.R. Nix, W.D. Myers and W.J. Swiatecki, Atomic Data and Nuclear Data Tables 59 185 (1995).
16. R. Corke and T. Sjostrand, JHEP 1105, 009 (2011), arXiv:1101.5953.
17. K. Tsukada et al., Phys. Rev. Lett. 118 (2017) 262501 [arXiv:1703.04278].
18. W. Fischer et al. Z. Physik 270 (1974) 270 113
19. ALICE Collaboration, public note, ALICE-PUBLIC-2018-003 (2018).
20. U. W. Heinz and A. Kuhlman, Phys. Rev. Lett. 94, 132301 (2005), arXiv:nucl-th/0411054
21. See http://root.cern.ch for installation files and documentation.
23. M. Alvioli and M. Strikman, Phys. Lett. B 722, 347 (2013) arXiv:1301.0728.
24. G. Aad et al. [ATLAS Collaboration], Eur. Phys. J. C 76 (2016) 199, arXiv:1508.00848.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters