GLISSANDO 3: GLauber Initial-State Simulation AND mOre..., ver. 3

TL;DR

GLISSANDO 3 is an advanced Monte-Carlo generator for modeling initial states in ultra-relativistic heavy-ion collisions, incorporating new features like the wounded quark model and deformation effects, with user-friendly interface and ROOT integration.

Contribution

This version introduces the wounded quark model and deformation effects, enhancing the simulation of initial collision states with simplified input and ROOT compatibility.

Findings

Successfully reproduces particle multiplicities at RHIC and LHC

Allows study of nucleon substructure fluctuations in p-A collisions

Supports simulations of light nuclei and deformed heavy nuclei

Abstract

We present ver. 3 of GLISSANDO, a versatile Monte-Carlo generator for Glauber-like models of the initial stages of ultra-relativistic heavy-ion collisions. The present version incorporates the wounded quark model, which is phenomenologically successful in reproducing multiplicities of particle production at RHIC and the LHC. Within this model, one can study the nucleon substructure fluctuation effects, recently explored in p-A collisions. In addition, the code includes the possibility of investigating collisions of light nuclei, such as $^{3}$He and $^{3}$H, or the $\alpha$-clustered $^{7,9}$Be, $^{12}$C, and $^{16}$O, where the deformation of the intrinsic wave function influences the transverse shape of the initial state. The current version, being down-compatible, retains the functionality of the previous releases, such as incorporation of various variants of Glauber-like models, a smooth $NN$ inelasticity profile in the impact parameter obtained from a parametrization of experimental data, fluctuating strength of the entropy deposition, or realistic nuclear distributions of heavy nuclei with deformation. The code can provide output in the format containing the event-by-event source location, which may be further used in modeling the intermediate evolution phase, e.g., with hydrodynamics or transport models. The interface is simplified, such that in the control input file the user may supply only the very basic information, such as the collision energy, the mass numbers of the colliding nuclei, and the model type. GLISSANDO 3 is integrated with the CERN ROOT platform. The package includes numerous illustrative and useful ROOT scripts to compute and plot various results.

Figures

Tables

Table 1. Table 1: Parameters of the distribution of light clustered nuclei used in GLISSANDO 3 .

Nucleus	$l$ [fm]	$r_{\alpha }$ [fm]	$r_{{}_{}{}^3\mathrm{He}}$ [fm]	$r_n$ [fm]
7Be	3.2	1.2	1.4	-
9Be	3.6	1.1	-	1.9
12C	2.8	1.1	-	-
16O	3.2	1.1	-	-

Table 2. Table 3: Parameters in the input file.

name	default	description
ECM	-1	center of mass energy of the colliding system [GeV]
NUMA	208	mass number of nucleus $A$
NUMB	208	mass number of nucleus $B$
EVENTS	10000	number of generated events
BMIN	0.	minimum impact parameter [fm]
BMAX	25.	maximum impact parameter [fm]
W0	2	minimum allowed number of wounded nucleons
W1	100000	maximum allowed number of wounded nucleons
RDS0	0	minimum allowed RDS
RDS1	100000	maximum allowed RDS
MODEL	0	0 - constant superimposed weight, 1 - Poisson, 2 - Gamma, 3 - Negative Binomial
Uw	1.	Poisson, Gamma, or Negative Binomial mean value for wounded objects
Ubin	1.	Poisson, Gamma, or Negative Binomial mean value for binary collisions
Vw	2.	Gamma or Negative Binomial variance for wounded objects
Vbin	2.	Gamma or Negative Binomial variance for binary collisions
NCS	3	number of partons in the nucleon
NBIN	40	number of bins for histograms in $\rho$, $x$, or $y$
RWSA	-1	Woods-Saxon radius for the distribution of centers, nucleus $A$ [fm]
AWSA	0.459	Woods-Saxon width, nucleus $A$ [fm]
BETA2A	-1	deformation parameter ${\beta }_2$, nucleus $A$
BETA4A	-1	deformation parameter ${\beta }_4$, nucleus $A$
ROTA_THETA	-1	rotation parameter (angle $\theta$), -1 - random rotation, nucleus $A$
ROTA_PHI	-1	rotation parameter (angle $\varphi$), -1 - random rotation, nucleus $A$
RWSB	-1	Woods-Saxon radius for the distribution of centers, nucleus $B$ [fm]
AWSB	0.459	Woods-Saxon width, nucleus $B$ [fm]
BETA2B	-1	deformation parameter ${\beta }_2$, nucleus $B$
BETA4B	-1	deformation parameter ${\beta }_4$, nucleus $B$
ROTB_THETA	-1	rotation parameter (angle $\theta$), -1 - random rotation, nucleus $B$
ROTB_PHI	-1	rotation parameter (angle $\varphi$), -1 - random rotation, nucleus $B$
SNN	-1	NN (parton-parton) “wounding” cross section [mb]
SBIN	-1	NN (parton-parton) binary cross section [mb]
ALPHA	0.12	mixed model parameter, 0 - wounded, 1 - binary, 0.12 - LHC@2.76 TeV/nucleon
CD	0.9	closest allowed distance between centers of nucleons [fm]
ISEED	0	seed for the random number generator, if 0 a random seed is generated
BTOT		range parameter for histograms [fm], default= max(RWSA,RWSB)+AWSA+AWSB
SHIFT	1	1 - shift the coordinates of the fireball to the c.m. in the fixed-axes case, 0 - do not shift
DW	0.	dispersion of the location of the source for wounded objectss [fm]
DBIN	0.	dispersion of the location of the source for binary collisions [fm]
WFA	0	the w parameter of the Fermi distribution, nucleus $A$
WFB	0	the w parameter of the Fermi distribution, nucleus $B$
FBIN	72	number of bins for histograms in the azimuthal angle
DOBIN	0	1 - compute the binary collisions also for the case ALPHA=0, 0 - do not
GA	1.0	central value of the Gaussian wounding profile
PP	-1	power of the transverse radius in the Fourier moments
RO	0	rank of the rotation axes (0 - rotation rank = rank of the Fourier moment)
RCHA	5.66	harmonic oscillator shell model density mean squared charge radii of nucleus $A$ (12C-nucleus)
RCHB	5.66	harmonic oscillator shell model density mean squared charge radii of nucleus $B$ (12C-nucleus)
RCHP	0.7714	harmonic oscillator shell model density mean squared charge radii of proton
OMEGA	-1	relative variance of cross-section fluctuations for the Gamma wounding profile
GAMA	-1	central value of the Gamma wounding profile
SCALEA	1.	scale parameter for the size of the nucleus (cluster version) [fm]
SIGMAA	1.	standard deviation of x, y, z coordinates of nucleons in the alpha cluster [fm]
SIGMABISA	1.	standard deviation of x, y, z coordinates of nucleons in the 3He cluster or nucleon no. 9 [fm]
SCALEB	1.	scale parameter for the size of the nucleus (cluster version) [fm]
SIGMAB	1.	standard deviation of x, y, z coordinates of nucleons in the alpha cluster [fm]
SIGMABISB	1.	standard deviation of x, y, z coordinates of nucleons in the 3He cluster or nucleon no. 9 [fm]
QSCALE	-1	scale parameter in the parton distribution function [fm]
DS	-1	source smearing parameter

Equations

$$p_{\mathrm{in}}(b,s)={\sigma }_{\mathrm{tot}}(b,s)-{\sigma }_{\mathrm{el}}(b,s)=4p_{\mathrm{CM}}\mathrm{Im}h(b,s)-4p_{\mathrm{CM}}^2|h(b,s){|}^2,$$

$$2p_{\mathrm{CM}}h(b,s)=2{\int }_0^{\infty }qdq{J}_0(bq)f(s,-q^2),$$

$$p_{\mathrm{in}}(s,b)=G\Gamma \left(\frac{1}{\omega (s)},\frac{\pi G(s)b^2}{\omega (s){\sigma }_{\mathrm{in}}(s)}\right)/\Gamma \left(\frac{1}{\omega (s)}\right),$$

$${\sigma }_{\mathrm{in}}(s)$$

$$G(s)$$

$$\omega (s)$$

$$\frac{d{N}_{\mathrm{ch}}}{d\eta }\sim Q_{\mathrm{W}}$$

$$\rho (r)=\mathrm{const}{r}^2\exp \left[-\sqrt{\frac{2}{3}\left(1-\frac{1}{p}\right)}\frac{r}{r_0(s;p)}\right].$$

$$p_{\mathrm{in}}^{qq}(s,b)=\exp \left[-\frac{\pi b^2}{{\sigma }_{\mathrm{in}}^{qq}(s;p)}\right],$$

$${\sigma }_{\mathrm{in}}^{qq}(s;p)=[A_1+A_2(\sqrt{s}/\mathrm{GeV}{)}^{A_3}+A_4/(\sqrt{s}/\mathrm{GeV})]\mathrm{mb},$$

$$A_1=-1.66+8.73/p+0.11p,A_2=-1.58+9.03/p+0.08p,$$

$$A_3=0.22-0.07/p+0.01p,A_4=-1.02+23.10/p-0.05p,$$

$$r_0(s;p)=[B_1+B_2(\sqrt{s}/\mathrm{GeV}{)}^{B_3}]\mathrm{fm},$$

$$B_1=0.37-0.95/p^2+0.11,B_2=-0.21-0.39/p^2+0.22/p,$$

$$B_3=-0.28-5.21/p^2+1.36p.$$

$$f_i(\vec{r})=\mathrm{const}\exp \left(-\frac{3}{2}\frac{(\vec{r}-\overrightarrow{c_i}{)}^2}{r_c^2}\right),$$

$$f_n(\vec{r})=\mathrm{const}{r}^2\exp \left(-\frac{3}{2}\frac{r^2}{r_n^2}\right),$$

$$\Delta (x,y;x_i,y_i)=\frac{1}{\pi {\sigma }^2}\exp \left[-\frac{(x-x_i{)}^2+(y-y_i{)}^2}{{\sigma }^2}\right],$$

$${ϵ}_n{e}^{i{\Psi }_n}\equiv -\frac{\int dxdys(x,y){\rho }^n{e}^{in\varphi }}{\int dxdys(x,y){\rho }^n},$$

$${ϵ}_n{e}^{i{\Psi }_n}=-\frac{{\sum }_i{\rho }_i^n{e}^{in{\varphi }_i}}{{\sum }_i{\rho }_i^n},$$

$${ϵ}_n{e}^{i{\Psi }_n}=-\frac{{\sum }_i{N}_i^{(n)}}{{\sum }_i{D}_i^{(n)}},$$

$$N_i^{(n)}=\int dxdy\Delta (x,y;x_i,y_i){\rho }^n{e}^{in\varphi },D_i^{(n)}=\int dxdy\Delta (x,y;x_i,y_i){\rho }^n.$$

$$N_i^{(n)}=\frac{1}{\pi \sigma }\int \rho d\rho \int d\varphi \exp \left(-\frac{{\rho }^2+{\rho }_i^2-2\rho {\rho }_i\cos (\varphi -{\varphi }_i)}{{\sigma }^2}\right){\rho }^n{e}^{in(\varphi -{\varphi }_i)}{e}^{in{\varphi }_i}.$$

$$N_i^{(n)}=\frac{2}{{\sigma }^2}\int \rho d\rho {\rho }^{n+1}{e}^{-\frac{{\rho }^2+{\rho }_i^2}{{\sigma }^2}}{I}_n\left(\frac{2\rho {\rho }_i}{{\sigma }^2}\right)e^{in{\varphi }_i},$$

$$N_i^{(n)}={\rho }_i^n{e}^{in{\varphi }_i},$$

$$D_i^{(n)}=\frac{1}{\pi \sigma }\int \rho d\rho \int d\varphi \exp \left(-\frac{{\rho }^2+{\rho }_i^2-2\rho {\rho }_i\cos (\varphi -{\varphi }_i)}{{\sigma }^2}\right){\rho }^n$$

$$=\frac{2}{{\sigma }^2}\int \rho d\rho {\rho }^{n+1}{e}^{-\frac{{\rho }^2+{\rho }_i^2}{{\sigma }^2}}{I}_0\left(\frac{2\rho {\rho }_i}{{\sigma }^2}\right)={\sigma }^n\Gamma \left(\frac{n}{2}+1\right)L_{\frac{n}{2}}\left(-\frac{{\rho }_i^2}{{\sigma }^2}\right),$$

$$D_i^{(2)}={\rho }_i^2+{\sigma }^2,$$

$$D_i^{(4)}={\rho }^4+4{\rho }_i^2{\sigma }^2+2{\sigma }^4,$$

$$D_i^{(6)}={\rho }_i^6+9{\rho }_i^4{\sigma }^2+18{\rho }_i^2{\sigma }^4+6{\sigma }^6,$$

$$\dots$$

$$D_i^{(3)}=\frac{1}{4\sigma }\sqrt{\pi }{e}^{-\frac{{\rho }_i^2}{2{\sigma }^2}}[2{\rho }_i^2\left({\rho }_i^2+2{\sigma }^2\right)I_1\left(\frac{{\rho }_i^2}{2{\sigma }^2}\right)$$

$$+\left(2{\rho }_i^4+6{\rho }^2{\sigma }^2+3{\sigma }^4\right)I_0\left(\frac{{\rho }_i^2}{2{\sigma }^2}\right)],$$

$$D_i^{(5)}=\frac{1}{8\sigma }\sqrt{\pi }{e}^{-\frac{{\rho }_i^2}{2{\sigma }^2}}[{\rho }_i^2\left(4{\rho }_i^4+24{\rho }_i^2{\sigma }^2+23{\sigma }^4\right)I_1\left(\frac{{\rho }_i^2}{2{\sigma }^2}\right)$$

$$+\left(4{\rho }_i^6+28{\rho }_i^4{\sigma }^2+45{\rho }_i^2{\sigma }^4+15{\sigma }^6\right)I_0\left(\frac{{\rho }_i^2}{2{\sigma }^2}\right)],$$

$$\dots$$

$$\mathtt{RDS}=\frac{1-\alpha }{2}\mathtt{NwAB}+\alpha \mathtt{Nbin},$$

Full text

PROGRAM SUMMARY

Program Title: GLISSANDO 3 ver. 3.42

Licensing provisions: CC By 4.0

Programming language: C++ with the ROOT libraries

*Nature of problem: The code implements in a versatile way the Glauber modeling of the initial stages of ultra-relativistic nuclear collisions, including the wounded nucleon and wounded quark models, with possible admixture of binary collisions. A state-of-the art inelastic nucleon-nucleon collision profile is implemented. A statistical distribution of the strength of the sources can be overlaid. The $\alpha$ clustered structure of light nuclei is built in.

*Solution method: Monte-Carlo simulation of nuclear collisions, analyzed off-line with numerous ROOT scripts. The software allows for a straightforward event-by-event analysis of eccentricity coefficients and their correlations, size fluctuations, or multiplicity distributions.

*Additional comments including Restrictions and Unusual features: The input may consist only of the model type, the mass numbers of the nuclei and the collision energy. The output can also be used as initial conditions for further hydrodynamic studies.

1 Introduction

We present an extended version of the package GLISSANDO – GLauber Initial-State Simulation AND mOre, previously described in [1] (ver. 1) and [2] (ver. 2), where the user is directed for a more detailed description of the underlying physics and the features which remain unchanged in the present version.

With the on-going efforts to understand the nature of the ultra-relativistic nuclear collisions, more ideas and models are being explored. Consequently, the dedicated open-access software should follow these efforts and supply useful analysis tools. Concerning the earliest stage of such collisions, modeling based on the so-called Glauber approach [3, 4, 5] has been commonly used, as it works well phenomenologically.

Another dedicated publicly available code implementing the Glauber approach is the PHOBOS Glauber Monte Carlo [6, 7, 8], provided as a ROOT package. The Glauber initial condition is also an option in the parameterizations of TRENTo code [9]. Many other ultra-relativistic nuclear collision codes, e.g., HIJING [10], AMPT [11], URQMD [12], or EPOS [13] use internally the Glauber approach to model the early stage of the collision.

A very successful variant of the Glauber model is based on wounded quarks [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26] acting as the particle production sources, rather than the commonly used wounded nucleons [27, 28] amended with binary collisions [29, 30]). Implementation of the wounded quark (or, in general, wounded parton) model is the first major extension in GLISSANDO 3.

Second, collisions of small nuclei with non-uniform distribution of nucleons may offer further stringent tests of the geometry-flow transmutation mechanism, including the event-by-event fluctuations [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 13, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65]. Moreover, the heavy-light reactions may provide insight into the ground-state clusterization structure in light nuclei, as first proposed in [66, 67] and further explored in [68, 69, 70]. We thus offer a possibility of simulating reactions with 3He (studied experimentally in [71]) or 3H nuclei, which may be thought of as small triangles, as well as collisions involving light $\alpha$-clustered nuclei, such as 7,9Be (dumbbell), 12C (triangle), and 16O (tetrahedron).

The nuclear configurations for 3He and 3H from [72, 73], as supplied with the PHOBOS Glauber Monte Carlo input files [7], are used in GLISSANDO 3 for simulations with $A=3$ nuclei. Similarly, the user has an option to use the configurations for 12C and 16O from [72, 73] as provided in [7].

The user interface to the code has been greatly simplified. In GLISSANDO 3, it suffices to provide the type of the model used (wounded nucleons or wounded partons), the energy of the collision, and the mass numbers of the colliding nuclei. The appropriate $NN$ inelastic cross section and its profile in the impact parameter are then automatically generated from interpolation formulas based on the experimental $pp$ scattering data. Specifically, we use here the COMPETE model parametrization implemented in the Particle Data Group review [74], which offers statistically the best description on the $pp$ and $p\bar{p}$ scattering data over a very broad range of the collision energies.

As in the previous releases, the output of GLISSANDO 3 is generated in two ways. The key features of the formed initial states are stored in a ROOT tree containing ready to use results in a compact form. In addition, a full event-by-event output can be generated for use in event-by-event studies of further evolution stages, such as hydrodynamics [35, 75, 76, 77, 59, 78, 79, 80, 81] or transport.

The new features implemented in GLISSANDO 3 include:

• •

State-of-the art nucleon-nucleon inelastic cross sections and their impact-parameter profiles.

• •

Implementation of the wounded quark (in general, wounded parton) model.

• •

Possibility of colliding 3He and 3H, with the distributions from external files.

• •

Inclusion of $\alpha$-clustered structure of 7,9Be, 12C, and 16O nuclei.

• •

Possibility of studying the effect of the proton fluctuations in the framework of the wounded parton model.

• •

Simplified user interface.

The main retained features of the previous releases incorporate:

• •

Parametrization of distributions of all popular nuclei, in particular those used in the energy - system-size scan of the NA61/SHINE [82].

• •

Inclusion of the nuclear deformation [83, 84, 85].

• •

Possibility of using externally-generated (correlated) nuclear distributions, e.g., [86, 87] for 208Pb, 197Au or 40Ca.

• •

Possibility of studying proton-nucleus and deuteron-nucleus collisions.

• •

Possibility of overlaying weights over the distribution of sources (Poisson, Gamma, negative binomial).

• •

Study of the core-corona effect [88, 89, 90, 91].

• •

Output of the event-by-event data with location of the sources to a text file, to be used for initialization of hydrodynamics of transport codes.

• •

A reference manual generated with doxygen [92] may be useful for those who wish to tailor the code.

2 New features in GLISSANDO 3

2.1 Inelastic cross sections and inelasticity profiles

With new measurements of the total, elastic, and differential elastic cross section, in particular at the LHC [93, 94, 95] and in cosmic rays [96], new parameterizations of the $pp$ and $p\bar{p}$ scattering amplitudes became available. A standard description is provided by the COMPAS group within the COMPETE model, presented in Sec. 52 of the 2016 edition of the Review of Particle Physics [74], which parametrizes the $pp$ and $p\bar{p}$ scattering data in a very accurate way. In GLISSANDO 3 we use this parameterization of the data as the most complete one.

The inelastic profile function, determining the probability of an inelastic $NN$ collision at impact parameter $b$ and collision energy $\sqrt{s}$ is defined as

[TABLE]

where $h(b,s)$ is the Fourier-Bessel transform of the $pp$ elastic scattering amplitude $f(s,-q^{2})$,

[TABLE]

with $p_{\rm CM}=\sqrt{s/4-m_{N}^{2}}$ denoting the CM momentum of the nucleon. Using the COMPETE model for $f(s,-q^{2})$ with parametrization from [74], we thus obtain $\sigma_{\rm in}(b,s)$. In the code, rather than using these long formulas involving a numerical integration in Eq. (2), we apply a simple parametrization of the inelastic profile of the form [97],

[TABLE]

where $\Gamma(a,z)$ denotes the incomplete Euler $\Gamma$ function, $\sigma_{\rm in}(s)$ is the inelastic cross section, and $G(s)$ and $\omega(s)$ are suitably adjusted parameters such that the COMPETE results are accurately reproduced. The values obtained from our fit can be efficiently represented with interpolating functions (which we use for $\sqrt{s}\geq 5$ GeV):

[TABLE]

Note that the applied expressions are for shear numerical convenience and do not reflect physical mechanisms of high-energy $pp$ scattering.

The quality of the interpolating formulas is illustrated in Fig. 1, where we compare the inelastic profile functions from the COMPETE model to Eq. (3) with parameters (4) for three selected collision energies. In the left panel of Fig. 2 we show the results of the COMPETE model for the total inelastic $pp$ cross section $\sigma_{\rm in}(s)$ (solid line) and our simple fit of Eq. (4). We also show (dotted line) the parametrization $A+B\,{\rm ln}^{2}(s/{\rm GeV}^{2})$, with $A=25$ mb and $B=0.146$ mb, used in [8]. We note that whereas at higher collision energies ($\sqrt{s}>100$ GeV) this parametrization agrees with the COMPETE model (and thus with the data) accurately, for $\sqrt{s}<100$ GeV the COMPETE parametrization, or our Eq. (4), is clearly advantageous. The interpolated dependence of parameters $G$ and $\omega$ on the collision energy, as given by formulas (4), is displayed in the right panel Fig. 2.

We recall that the use of a realistic $NN$ inelasticity profile (as opposed to sometimes used hard-sphere profile in older software, e.g., [10]) is phenomenologically important. First, it leads to proper elastic differential cross section in $pp$ collisions, as shown above. Second, it yields sizable effects in standard observables in ultra-relativistic nuclear collisions. In particular, the effect on ellipticity of the fireball was discussed in detail in [98], where the use of a realistic (smooth) wounding profile led to a significant (10-20%) reduction for peripheral A-A collisions, compared to the case with a hard-sphere profile.

2.2 Wounded partons

In the wounded quark model [14, 15, 17, 16] the valence quarks play the role of elementary scatterers. The approach is implemented in GLISSANDO 3, as it has proved to be phenomenologically successful in describing the multiplicity of produced hadrons in a variety of reactions, with the linear scaling

[TABLE]

well satisfied. Here $N_{\rm ch}$ denotes the observed charged hadrons, and $Q_{\rm W}$ the number of wounded quarks in a given reaction and centrality class. It was first noticed that the RHIC data follow the scaling [18], further explored at RHIC by the PHENIX Collaboration [20, 21], as well as at the SPS [19]. See also [7, 20, 21, 23, 22, 24, 25, 26] for further developments.

Our implementation of the model is described in detail in [23]. The nucleon consists of $p$ partons ($p=3$ for the wounded quark model) distributed according to a radial density

[TABLE]

The scale parameter $r_{0}(s;p)$ controls the size of the nucleon built of partons. The factor $1-1/p$ accounts for the center-of-mass corrections; a shift to the center of mass of the nucleon is carried out after randomly generating the positions of the centers of the partons. As $r_{0}$ is adjusted phenomenologically (see the following), the presence of the $1-\tfrac{1}{p}$ factor in formula (6) is conventional.

In the wounded parton approach one needs to choose the parton-parton inelasticity profile. We take it in the Gaussian form

[TABLE]

where $\sigma_{\rm in}^{qq}(s)$ is the parton-parton inelastic cross section.

Thus our wounded parton model implementation brings in two parameters: $r_{0}(s;p)$ and $\sigma_{\rm in}^{qq}(s;p)$. Their values are chosen in such a way that the resulting $NN$ collision profile reproduces the COMPETE parametrization discussed in Sec. 2.1. In other words, we impose the desired feature that the $NN$ inelasticity profile is the same if one uses nucleons or partons as elementary scatterers. With our chosen parameterizations of Eq. (6) and (7) this cannot be accomplished in an exact manner, but optimization leads to a close agreement (within a few percent). This can be seen in Fig. 3, where we compare the $NN$ elasticity profiles (multiplied with $b$) obtained in the nucleon and parton models for $\sqrt{s}=7$ TeV. The quality of agreement is similar at other collision energies or values of $p$.

The optimum values of the wounded parton model parameters are interpolated as follows for the parton-parton cross section and the parton distribution parameter:

[TABLE]

[TABLE]

The above parameterizations work for $3\leq p\leq 10$. The dependence of $\sigma_{\rm in}^{qq}(s;p)$ and $r_{0}(s;p)$ on $\sqrt{s}$ is visualized in Fig. 4. We note that in our model the parton-parton cross section decreases with $p$, whereas the dependence of the size parameter $r_{0}$ exhibits a rather weak growth with $p$. We also note that an approximate scaling $\sigma_{\rm in}^{qq}(s;p)\sim 1/p^{2}$ holds, which is accurate to a few percent for $p>3$.

As in previous releases of our code, centers of nucleons in a nucleus are distributed according to a Woods-Saxon density with additional $NN$ repulsion [1], or are taken from external calculations, e.g., [86] of [7]. Nuclear shape deformation [99, 83, 84, 100] is implemented for deformed nuclei, such as 63Cu, 197Au, or 238U, in the same manner as in GLISSANDO 2.

2.3 Reactions with light clustered nuclei

GLISSANDO 3 offers a possibility to collide nuclei whose structure exhibits clustering. Specifically, we implement 7Be as $\alpha+^{3}$He, 9Be as $2\alpha+$neutron, 12C as $3\alpha$, and 16O as $4\alpha$ (see Fig. 5). The details of our fixing of geometric parameters of the distributions are presented in [69]. The parameters are adjusted in such way that the one-body nuclear distributions [72, 73] are matched. We first arrange the positions of clusters as in Fig. 5, separating the centers by the distance $l$. The distribution of the centers of nucleons in each cluster is randomly generated according to a Gaussian

[TABLE]

where $\vec{r}$ denotes the coordinate of the nucleon and $\vec{c_{i}}$ is the position of the center of the cluster $i$. The width of the cluster is controlled with $r_{c}$, the rms radius of the cluster. Then the positions of the nucleons are generated sequentially, switching between clusters 1, 2,…, 1, 2,…, until all the needed nucleons are placed.

For the case of 9Be we place the extra neutron on top of the two $\alpha$ clusters. The distribution is

[TABLE]

and it exhibits a hole in the middle.

The values of parameters used in the code for the cluster model are collected in Table 1.

One may do many physical studies with light clustered nuclei, in particular looking for clusterization signatures in harmonic flow patterns [66, 67, 69, 101]. In Fig. 6 we present an example of a quantity, the symmetric cumulant of elliptic and triangular flows, which displays different behavior on the model with clusters compared to the case without clusters (uniform).

2.4 Smearing of sources

The Monte Carlo Glauber model provides, event-by event, the positions of point-like sources (wounded objects) in the transverse plane. Physically, the sources must have transverse size, which is accomplished by smoothing them, typically with a Gaussian form

[TABLE]

where $(x,y)$ are the transverse coordinates and $(x_{i},y_{i})$ gives the location of the $i$th source, and $\sigma$ is the smearing width. Smearing is necessary, i.a., in forming the initial condition for subsequent hydrodynamic evolution. Moreover, smearing increases the size of the fireball and reduces the eccentricities, which are its basic characteristics. The effect is particularly relevant for small systems. The effects of Gaussian smearing (12) are simple to implement, as the corresponding expressions are analytic.

The eccentricity of rank $n$ ($n\geq 2$) in a given event, $\epsilon_{n}$, and the corresponding angle of the principal axis (the event-plane angle), $\Psi_{n}$, are defined in the standard way as

[TABLE]

where $s(x,y)$ is the transverse entropy distribution, and we have introduced the polar coordinates via $x=\rho\cos\phi$, $y=\rho\sin\phi$. When the sources are point-like, then $s(x,y)=\sum_{i}\delta(x-x_{i})\delta(u-y_{i})$ and

[TABLE]

where $i$ runs over all the sources in the event, and $x_{i}=\rho_{i}\cos\phi_{i}$, $y_{i}=\rho_{i}\sin\phi_{i}$. With Gaussian smearing (12), we have $s(x,y)=\sum_{i}\Delta(x,y;x_{i},y_{i})$ and, correspondingly,

[TABLE]

In polar coordinates we may write explicitly

[TABLE]

Integration over $\phi$ yields

[TABLE]

and integration over $\rho$

[TABLE]

which means that the numerator with smearing in Eq. (2.4) is exactly the same as for point-like sources in Eq. (14).

Repeating the above steps for the denominator yields

[TABLE]

where $\Gamma(z)$ is the Euler Gamma function and $L_{a}(z)$ is the Laguerre polynomial. The forms of $D_{i}^{(n)}$ for the first few even values of $n$ are

[TABLE]

For odd values of $n$

[TABLE]

where $I_{0,1}(z)$ are the modified Bessel functions of the first kind. It is apparent from the above formulas that smearing increases the denominators, therefore $\epsilon_{n}$ is quenched, in agreement with intuition. At the same time we observe that the principal axes angles $\Psi_{n}$ are not altered by smearing.

The reduction of the contribution of a single source placed at $\rho_{i}$ for first few eccentricities, quantified with $\rho_{i}^{2}/D_{i}^{(n)}$, is shown in Fig. 7. We note that, as expected, the effect increases with $n$ and with the value of the smearing parameter $\sigma$.

Taking for definiteness the ellipticity $(n=2)$ and collecting the contributions from all sources, we find the smearing reduction factor $\langle r^{2}\rangle/(\langle r^{2}\rangle+\sigma^{2})$, which becomes increasingly important with decreasing mean squared radius, $\langle r^{2}\rangle$, of the fireball.

Smearing is incorporated in GLISSANDO 3 in the evaluation of eccentricities and the mean squared radius of the fireball. The amount of smearing is controlled with the parameter DS$=\sigma/2$, whereby it has the interpretation of the Gaussian width in one dimension.

2.5 Rapidity modeling

The modeling in rapidity, present to some extent in ver. 2 of the code, has been removed. This is because, unlike the modeling in the transverse plane where the Glauber approach became one of the standards, the longitudinal modeling has much more freedom and is under development. If the user wishes to model the rapidity dependence of the fireball, she/he may use the transverse distributions (generated event-by-event and written to an output file) in an external application.

3 Installation and running

The user must first install the CERN ROOT package [102]. After downloading and unpacking GLISSANDO 3, the command

make

should be run, which creates the executable binary glissando3.

To optionally recreate the doxygen manual (provided as /doc/latex/refman.pdf with the distribution), the doxygen [92] package should be installed and the following commands executed:

make cleandoc make doc

With the code installed, the user may look at an instructive presentation of the main features of the present version by executing the shell script

./demo_ver_3.sh

To see the created pdf files evince should be installed prior to the run. To use another pdf viewer, the line VIEW=evince should be edited in the *.sh files. The features of ver. 2, retained in the present release, may be explored by executing

./demo_ver_2.sh

For collisions of two nuclei, the standard running command has the syntax

./glissando3 [input_file] [output_file]

where the nuclei can be any, including the proton and the deuteron. When the input or output arguments are absent, the defaults are

input.dat - default input glissando.root - default output

The input parameters with their defaults are collected in Appendix A.

3.1 Makefile

The Makefile contains commands for compilation and linking. It can be used for profiling the code by modifying the preprocessor options. Most importantly, the user should decide if she/he uses the wounded nucleon or the wounded parton version of the model. In the nucleon case the standard preprocessor options should read

PREPROCESS := -D_partons_=0 -D_nnwp_=2 -D_bindep_=1 -D_files_=0 -D_profile_=0 -D_weight_=0 -D_evout_=0 -D_clusters_=0 -D_uncluster_=0 -D_pardis_=0 -D_rdsconv_=1

whereas for the wounded parton model

PREPROCESS := -D_partons_=1 -D_nnwp_=1 -D_bindep_=1 -D_files_=0 -D_profile_=0 -D_weight_=0 -D_evout_=0 -D_clusters_=0 -D_uncluster_=0 -D_pardis_=0 -D_rdsconv_=1

The meaning of all preprocessor parameters is as follows:

nnwp=0 - hard-sphere NN wounding profile, =1 - Gaussian, =2 - gamma bindep=1 - binary and wounding depend on each other, =0 - do not files=1 - read the nuclear distributions from external files, =0 - generate randomly profile=1 - generate the nucleon profile and NN correlation data, =0 - do not weight=1 - generate data for the NN collision profiles and for the weight (RDS) ditributions, 0 - do not evout=0 - do not generate the text event data, =1 - short, =2 -long clusters=1 - generate clusters in light nuclei, =0 do not uncluster=1 - smoothe out the alpha clusters in nucleus A, =2 - in nucleus B, =3 in both nuclei, =0 - do not smoothe partons=1 - wounded parton model, =0 - wounded nucleon model pardis=0 - exponential parton distribution in nucleus, =1 - Gaussian rdsconv=2 - convention for RDS: 1 - Nw/2, 2 - Nw

The user may modify the Makefile file to suit his needs, or, alternatively, run make for example as follows:

make ’PREPROCESS = -D_partons_=1 -D_nnwp_=1 -D_evout_=1 ’

to produce the binary code for the wounded parton model with the Gaussian wounding profile and with generation of some results to an external ASCII file.

Another functionality is the storage of the current version of the package,

make package

as well as the cleaning options:

make clean make cleandoc make cleanoutput

3.2 Input

The input file is a standard ASCII file. Every line contains the name of the parameter separated with space from the assigned value. The sequence of the lines with parameters is flexible. When a parameter is missing in the input file, or a line containing it is commented out with #, the default value is used. See Appendix A for details concerning the input parameters.

3.3 Output

A typical output from a run of GLISSANDO 3 to the console is shown in Table 2.

The results of the simulation are stored in the ROOT output file or, optionally, also in text files when the preprocessor parameter evout=1 or 2 is used.

To see the physical results, the user should enter the ROOT environment

root

and execute one of the supplied scripts *.C

.x macro/demo_3/<script name>.C(<optional parameters>)

or directly examine the output ROOT file, e.g. with a

TBrowser a

command (executed within ROOT environment). An alternative method of executing the ROOT scripts is provided in the example shell demo_ver_3.sh. Proceeding this way does not require a greater familiarity with ROOT.

The preprocessor option evout controls the possible event-by-event output to an additional external ASCII file with the name of the ROOT output file appended with .points. If evout=2, then the file stores the transverse positions event-by-event. Its content consists of the blocks

RDS NwAB specA atan2(ycmA,xcmA) specB atan2(ycmB,xcmB)

followed with NwAB lines of the format

x y z c w

The second part of the block is

Nbin

followed with Nbin lines of the format

x y z c w

Above we use the notation RDS for the relative deposited strength in the event (defined below), NwAB for the number of wounded objects (nucleons or partons), specA for the number of spectators from nucleus A and atan2(ycmA,xcmA) for the azimuthal angle of their center of mass relative to the $x$ axis, and specB and atan2(ycmB,xcmB) for analogous quantities for nucleus B, and Nbin for the number of binary collisions. Next, x and y are the transverse coordinates of the sources in fm, z is the longitudinal coordinate, |c| indicates how many times a wounded nucleon collided, with positive (negative) c corresponding to nucleus A (B), while c=0 indicates the binary collisions. The last entry is the weight w of the given source (by source we customarily mean a wounded nucleon or a binary collision). The number of blocks is equals to the number of events.

If evout=1, then a simplified output is generated, with the format

NwA NwB Nbin RDS

where NwA and NwB are the numbers of wounded objects in nuclei A and B, respectively. The number of output lines equals to the number of events.

The relative deposited strength (term introduced in [1]) is equal to

[TABLE]

in accordance to the Glauber model with the binary collision admixture [29, 30]. Correspondingly, the weight of a given source is $(1-\alpha)/{2}$ for the wounded objects and $\alpha$ for the binary collisions.

3.4 Reading external nuclear distributions

The user may provide external files with nuclear distributions. We have used two sets of such files. The first one comes from Alvioli et al. [86] for the case of 16O, 40Ca, and 208Pb. They can be obtained from http://sites.psu.edu/color/. The user must create from the downloaded files a singe file, for instance running

cat ca40-1.dat ca40-2.dat ca40-3.dat [more files] > ca40.dat

The resulting file o16.dat should be placed in the relative subdirectory nucl.

The files have the format

x y z k

where x, y, and z denote the Cartesian coordinates of the centers on nucleons in fm, while $k=0$ for neutrons and $k=1$ for protons. Groups of A lines, where A is the mass number of the nucleus, correspond to a single nuclear configuration.

The second case concerns 3He or 3H distributions, which are taken from the Green’s function Monte Carlo simulations [103] as provided in [7]. For the user’s convenience the files he3_plaintext.dat and h3_plaintext.dat are also provided at the GLISSANDO web site. The format of these files is

x1 y1 z1 x2 y2 z2 x3 y3 z3 foo foo foo foo

where (xi, yi, zi) are the Cartesian coordinates of the $i$th nucleon in a given configuration.

Similarly, the configurations of 12C or 16O from [103] as provided in [7] may be used. They are provided at the GLISSANDO web site as carbon_plaintext.dat and oxygen_plaintext.dat.

To use the nuclear configurations from external files, the code must be compiled with

make ’PREPROCESS = -D_files_=1’

and executed as

./glissando3 [input_file] [output_file] [nucleus_A_file] [nucleus_B_file]

When the syntax

./glissando3 [input_file] [output_file] [nucleus_A_file]

is used, then the distribution of the nucleons in nucleus A is read from an external file, whereas for nucleus B the positions of nucleons are generated randomly. This syntax must also be used for the collisions of nucleus A with the proton or the deuteron.

3.5 Fixing centrality cuts

We explain how to impose centrality cuts on GLISSANDO 3 simulations and how to properly choose the range for the impact parameter $b$ (a too large range slows down the simulations). The supplied shell script centrality.sh implements the following procedure. First, a minimum-bias simulation must be run, with no (or broad-range) values for the W0, W1, RDS0, and RDS1 parameters, as well as BMIN=0 and BMAX set to a larger value than the sum of the radii of the two colliding nuclei. We typically use BMAX as about twice this sum to account for the tails in the Woods-Saxon distributions. Note, however, that large values for BMAX lead to a larger Monte Carlo rejection rate,slowing down the simulation. Next, the macro/demo_3/cent.C script should be executed in root. The values of W0, W1, or RDS0, RDS1 corresponding to the desired centrality classes can be read off from the generated file output/centrality.dat. Next, the input file must be modified with the proper values for W0, W1 supplied (if centrality is determined by the multiplicity of the wounded objects), or RDS0, RDS1 (if centrality is given by the relative deposited strength RDS), and the code must be rerun (with low statistics). Running /macro/demo_3/b_dist.C and looking at the plot b_dist.pdf allows us to determine the optimum values for BMIN=0 and BMAX. We thus should update these parameters in the input file and rerun simulations with full statistics.

The above steps can be simply traced by executing centrality.sh.

4 Structure of the code

The structure of the code is provided in the supplied reference manual created with doxygen.

5 Summary

We have presented GLISSANDO 3, with hopes it will continue to be a useful and versatile tool for the heavy-ion community. The authors most welcome suggestions and questions from the users.

We would like to thank dr Maksym Deliyergiyev for his help to make our scripts compatible with ROOT6.

This research was carried out in laboratories created under the project ”Development of research base of specialized laboratories of public universities in Swietokrzyskie region” no POIG 02.2.00-26-023/08 dated May 19, 2009.

Appendix A Input and output

The input parameters are collected in Table LABEL:tab:input. The hash sign # at the beginning of the line in the input file comments it out, in which case a default value of the parameter as set in the code (file functions.h) is used.