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 heavyion 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, Helium3, and Uranium, as well as the treatment
of deformed nuclei and GlauberGribov 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”)
I Introduction
In heavyion 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 eventbyevent 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 eventbyevent 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, Helium3, and Uranium, as well as the treatment of deformed nuclei and GlauberGribov 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 singleparticle 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 internucleon 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 lowenergy electron scattering experiments (6). The nuclear charge density is usually parameterized by a Fermi distribution with three parameters
(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), Helium3 (He3, He), and Sulfur (S) nuclei. For Sulfur, a three parameter Gaussian form is used
(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:

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

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
(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.
Name  R [fm]  a [fm]  w 

p 
0.234  
d 
0.010  0.5882  0.0000 
H3 

He3  
O  2.608  0.5130  0.5100 
Si  3.340  0.5800  0.2330 
S  2.540  2.1910  0.1600 
Ca  3.766  0.5860  0.1610 
Ni  4.309  0.5170  0.1308 
Cu  4.200  0.5960  0.0000 
CuHN  4.280  0.5000  0.0000 
Xe 
5.360  0.5900  0.0000 
Xes 
5.420  0.5700  0.0000 
W  6.580  0.4800  0.0000 
Au  6.380  0.5350  0.0000 
AuHN  6.420  0.4400  0.0000 
Pb 
6.620  0.5460  0.0000 
PbHN  6.650  0.4600  0.0000 
Pb*  6.624  0.5490  0.0000 
Name  R [fm]  a [fm]  

Al  3.340  0.5800  0.448  0.239 
Si2  3.340  0.5800  0.478  0.250 
Cu2  4.200  0.5960  0.162  0.006 
Xe2 
5.360  0.5900  0.161  0.003 
Xe2a 
5.360  0.5900  0.18  0 
Au2  6.380  0.5350  0.131  0.031 
U2  6.670  0.4400  0.280  0.093 
U 
6.670  0.4400  0.280  0.093 
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 nucleonnucleon 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
(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 eventbyevent 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 nucleonnucleon 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 GlauberGribov 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
(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 nucleusnucleus 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 eventbyevent quantities. It takes as parameters, the number of events to be generated, the collision system, the nucleonnucleon cross section, its width to simulate GlauberGribov 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 eventbyevent quantities:

Npart: Number of participating nucleons.

Npart0: Number of singlywounded 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 eventbyevent 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 nucleonnucleon 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 eventbyevent 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 nucleonnucleon 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.
Iv Summary
Glauber Model calculations are commonly used by heavyion 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, Helium3, and Uranium, as well as the treatment of deformed nuclei and GlauberGribov 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
 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.
 As explained in the text, also two further parameterizations (“dh” and “dhh”) for the Deuterium (H) are implemented. Please refer to the code.
 The configurations for H and He were obtained from Green’s function MC calculations as in LABEL:Nagle:2013lja.
 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.
 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.
 These values are usually also used for Pb, since the BesselFourier 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.
 Same parameters as for Xe with and from LABEL:xedef.
 Same parameters as for Xe with from interpolation between measured deformation parameters for the evenA Xe isotopes (19), and set to zero.
 Implementation as done in LABEL:Heinz:2004ir.
References
 B. Alver, M. Baker, C. Loizides and P. Steinberg (2008), arXiv:0805.4411.
 M.L. Miller, K. Reygers, S.J. Sanders, and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57 (2007) 205, arXiv:nuclex/0701025.
 B.Alver et al. [PHOBOS Collaboration], Phys. Rev. C 77 (2008) 014906, arXiv:0711.3724.
 M. Rybczynski, G. Stefanek, W. Broniowski and P. Bozek, Comput. Phys. Commun. 185, 1759 (2014), arXiv:1310.5475.
 C. Loizides, J. Kamin and D. d’Enterria, arXiv:1710.07098.
 H. De Vries, C.W. De Jager, and C. De Vries, Atom. Data Nucl. Data Tabl. 36 495 (1987).
 C. W. De Jager, H. De Vries and C. De Vries, Atom. Data Nucl. Data Tabl. 14, 479 (1974).
 L. Hulthen and M. Sugawara, Handbuch der Physik 39 1 (1957).
 S.S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 91 072303 (2003), arXiv:nuclex/0306021.
 S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. C 74, 024904 (2006), arXiv:nuclex/0603010.
 J. L. Nagle et al., Phys. Rev. Lett. 113 (2014) 112301, arXiv:1312.4565 [nuclth].
 T. Hirano and Y. Nara, Phys. Rev. C 79, 064904 (2009), arXiv:0904.4080.
 Q. Y. Shou et al. Phys. Lett. B 749 (2015) 215, arXiv:1409.8375 [nuclth].
 R. Hofstadter et al., Rev. Mod. Phys. 28 (1956) 214.
 P. Möller, J.R. Nix, W.D. Myers and W.J. Swiatecki, Atomic Data and Nuclear Data Tables 59 185 (1995).
 R. Corke and T. Sjostrand, JHEP 1105, 009 (2011), arXiv:1101.5953.
 K. Tsukada et al., Phys. Rev. Lett. 118 (2017) 262501 [arXiv:1703.04278].
 W. Fischer et al. Z. Physik 270 (1974) 270 113
 ALICE Collaboration, public note, ALICEPUBLIC2018003 (2018).
 U. W. Heinz and A. Kuhlman, Phys. Rev. Lett. 94, 132301 (2005), arXiv:nuclth/0411054
 See https://twiki.cern.ch/twiki/bin/view/Main/LHCGlauberBaseline.
 See http://root.cern.ch for installation files and documentation.
 See the TGlauberMC page on HepForge (http://www.hepforge.org/downloads/tglaubermc) for the most recent TGlauberMC release (latest version v2.6).
 M. Alvioli and M. Strikman, Phys. Lett. B 722, 347 (2013) arXiv:1301.0728.
 G. Aad et al. [ATLAS Collaboration], Eur. Phys. J. C 76 (2016) 199, arXiv:1508.00848.