Usage¶

frzout is a particlization model (Cooper-Frye sampler) for relativistic heavy-ion collisions. Its primary purpose is to sample particles in a multistage collision model, after hydrodynamics and before a hadronic afterburner. It can also be used to compute hadron resonance gas thermodynamic quantities (e.g. to calculate an equation of state).

Physics references:

The original paper by Cooper & Frye
“Particlization in hybrid models” by Huovinen & Petersen
The author’s dissertation, section 3.4

Sampling particles¶

In typical heavy-ion collision models, the hydro code outputs a spacetime hypersurface as a collection of finite volume elements. The surface is usually defined by a switching temperature.

The basic procedure for sampling particles:

Create a Surface object from the data output by the hydro code.
Create an HRG object, representing the hadron resonance gas from which particles will be sampled.
Call the sample function.

1. Create a Surface¶

Suppose a 3D ideal hydro code outputs surface data as a text table with the following columns:

Column index:  0  1  2  3  4        5        6        7        8   9   10
Quantity:      t  x  y  z  sigma_t  sigma_x  sigma_y  sigma_z  vx  vy  vz

The first four columns (0–3) are the components of the spacetime position vector, the next four (4–7) are the covariant normal vector σ_μ, and the last three (8–10) are the velocity.

Each row of the table represents a volume element.

We will now load the data into arrays and create a Surface object:

import numpy as np
import frzout

# read surface data from text file
surface_data = np.loadtxt('surface.dat')

# separate into sub-arrays
# (there are of course other ways to do this)
x, sigma, v = np.hsplit(surface_data, [4, 8])

# create Surface object
surface = frzout.Surface(x, sigma, v)

Note

The components of the vectors must be in standard order (t, x, y, z). If the surface data columns are in a different order, they must be rearranged in the arrays x, sigma, v.

The units of the spacetime position must be fm.

The units of the normal vector must be fm³. It must be covariant, i.e. σ_μ, not σ^μ. Most surface finding algorithms output the covariant vector. If for some reason it is contravariant, simply negate the spatial components.

The velocity must be dimensionless (relative to the speed of light). It is the three-vector v, not the four-velocity u^μ.

Boost invariance¶

Suppose a 2D (boost-invariant) ideal hydro code outputs a similar format, but without the z components:

Column index:  0  1  2  3        4        5        6   7
Quantity:      t  x  y  sigma_t  sigma_x  sigma_y  vx  vy

Analogously, we can load the data as:

x, sigma, v = np.hsplit(np.loadtxt('surface.dat'), [3, 6])
surface = frzout.Surface(x, sigma, v, ymax=1)

The Surface class infers whether the surface is 2D (boost-invariant) or 3D from the shapes of the input arrays. After creating the object, attribute surface.boost_invariant is a boolean indicating whether the surface is boost-invariant or not (2D or 3D).

Parameter ymax sets the maximum momentum rapidity: The dN/dy of sampled particles will be flat from −ymax to +ymax. The default value of ymax is 0.5. It has no effect for 3D surfaces.

Note

The units of the normal vector are formally fm² for boost-invariant surfaces. This is what most surface finding algorithms output. The full differential volume element is

\[d^3\sigma_\mu = \tau \, \Delta y \, d^2\sigma_\mu\]

The Surface class internally scales 2D normal vectors by 2 × ymax × τ.

Viscous pressures¶

If the hydro code is viscous, we must also load the viscous pressures and input them to the Surface class.

Suppose we are using the author’s version of the OSU viscous hydro code, which outputs surface data in the format:

Column index:  0  1  2  3        4        5        6   7
Quantity:      t  x  y  sigma_t  sigma_x  sigma_y  vx  vy

Column index:  8      9      10     11     12     13     14     15
Quantity:      pi^tt  pi^tx  pi^ty  pi^xx  pi^xy  pi^yy  pi^zz  Pi

The first 8 columns are the same as above, columns 8–14 are the components of the shear pressure tensor π^μν, and column 15 is the bulk pressure Π.

We can load the data as follows:

# data file is binary, not text
surface_data = np.fromfile('surface.dat', dtype='float64').reshape(-1, 16)

# extract usual sub-arrays
x, sigma, v, _ = np.hsplit(surface_data, [3, 6, 8])

# create mapping of pi components
pi = dict(
   xx=surface_data[:, 11],
   xy=surface_data[:, 12],
   yy=surface_data[:, 13]
)

# extract bulk pressure
Pi = surface_data[:, 15]

# create Surface object
surface = frzout.Surface(x, sigma, v, pi=pi, Pi=Pi)

As shown in the example, the shear tensor pi must be provided as a mapping (dict-like object).

For 2D surfaces, only components (xx, xy, yy) are required.
For 3D surfaces, components (xz, yz) are additionally required.

The remaining components are computed internally by using that the shear tensor is traceless and orthogonal to the velocity.

When pi and/or Pi are provided, viscous corrections are applied to particle sampling. See Viscous corrections for more information.

Note

The units of the viscous pressures must be GeV/fm³.

Warning

Surface data files can be very large, especially for central events and/or 3D hydro. The Surface class stores a copy of the data internally, effectively doubling memory usage. Loading data and creating a Surface may exceed system memory. Currently, the only way around this is to break the surface into smaller chunks and load them one at a time.

One way to ensure that the extra memory is freed after creating a Surface is to write a function that loads the data and returns the object:

def load_surface():
   x, sigma, v = ...
   return frzout.Surface(x, sigma, v, ...)

The local variables will go out of scope when the function returns and their memory will be released.

Text files exacerbate the problem, since typical methods of reading text files in Python (such as np.loadtxt) create a temporary list which uses even more memory than the final array. Further, reading and writing text files is slow. I recommend using a binary format if possible. The example above shows how to read the binary file output by osu-hydro.

2. Create an HRG¶

The HRG class represents a hadron resonance gas. The only required parameter is the temperature of the gas, which should be set to the switching temperature of the hydro surface. For example, to create an HRG with temperature 150 MeV:

hrg = frzout.HRG(.150)  # temperature units are GeV

Optional parameters control the composition of the gas and the behavior of resonances. See the reference section (below) for details, and Resonance mass distributions for information on the physical treatment of resonances.

Examples of creating HRG with more options:

# for use with UrQMD
hrg = frzout.HRG(.150, species='urqmd', res_width=True)

# a pion gas
hrg = frzout.HRG(.150, species=[111, 211])

Tip

It is possible, and in fact advisable, to reuse an HRG object for multiple Surface (for example if running several events in a row, all with the same particlization temperature). There are some preliminary calculations in constructing an HRG that are not repeated when the object is reused.

3. Call the sample function¶

After creating a Surface and HRG, sample particles by:

parts = frzout.sample(surface, hrg)

The function returns a numpy structured array with fields 'ID', 'x', 'p', containing the ID number, spacetime position four-vector (t, x, y, z), and energy-momentum four-vector (E, p_x, p_y, p_z), respectively, for each sampled particle. These can be extracted with a dict-like interface:

ID = parts['ID']  # 1D array of integers, shape (nparts,)
x = parts['x']  # 2D array of floats, shape (nparts, 4)
p = parts['p']  # 2D array of floats, shape (nparts, 4)

The particle array can be converted to a record array so that the fields are accessible by attribute rather than dict key:

parts = frzout.sample(surface, hrg).view(np.recarray)
x = parts.x
# etc...

We can loop over particles:

for p in parts:
   x = p['x']  # a single particle's position vector
   # ...

# or unpack every particle's attributes
for ID, x, p in parts:
   # ...

Examples of oversampling a surface:

# use samples one at a time
for _ in range(nsamples):
   parts = frzout.sample(surface, hrg)
   # do something with particles

# alternatively, to keep all samples at once
# (this could use a lot of memory!)
samples = [frzout.sample(surface, hrg) for _ in range(nsamples)]

The number of sampled particles is Poisson distributed. The average number may be computed as surface.volume * hrg.density(). Note that if the surface has any negative contributions (most physical surfaces do), the actual average will be somewhat larger than the computed average. The number of particles in a particular sample is given by parts.size.

Computing hadron gas quantities¶

The HRG class can also be used to compute thermodynamic quantities of a hadron resonance gas, e.g. to calculate the equation of state. For example, to compute the energy density of a hadron gas at 150 MeV:

frzout.HRG(.15).energy_density()

The calculations use the actual composition of the gas and take resonance width into account:

# pressure of 150 MeV gas of UrQMD species
frzout.HRG(.15, species='urqmd').pressure()

# same, but neglecting resonance width (result will be slightly different)
frzout.HRG(.15, species='urqmd', res_width=False).pressure()

See the reference section for all available quantities.

By calculating energy density, pressure, etc over a range of temperatures, we can construct an HRG equation of state. The script eos.py in the osu-hydro repository does this and connects the result to the HotQCD lattice equation of state at high temperature.

Note, the HRG class emits a warning about inaccurate momentum sampling for high temperatures. This may be safely ignored when only computing thermodynamic quantities.

Reference¶

Main classes¶

class frzout.Surface(x, sigma, v, pi=None, Pi=None, ymax=None)¶

A particlization hypersurface.

A surface consists of a collection of volume elements. This class is initialized with arrays containing data for each volume element.

This class handles both 2D (boost-invariant) and 3D surfaces.

Required parameters

These must be array-like objects of shape (n_volume_elements, n_components), where the number of volume elements must be the same for all arrays.

x – spacetime position
- components: (t, x, y, z)
- units: fm
sigma – normal vector (covariant, i.e. σ_μ, not σ^μ)
- components: (σ_t, σ_x, σ_y, σ_z)
- units: fm³
v – velocity
- components: (v_x, v_y, v_z)
- units: dimensionless (relative to speed of light)

The components listed above are for 3D surfaces. For 2D (boost-invariant) surfaces, there is no longitudinal (z) component and time (t) is replaced by proper time (τ). For example, the components of x for a 2D surface are (τ, x, y). In addition, the units of sigma are fm² for a 2D surface (normal vectors are internally multiplied by τ). The dimensionality of a surface is inferred from the shapes of these arrays.

Optional parameters

pi – shear viscous pressure tensor

A mapping (dict-like object) whose keys are two-character strings 'ij' representing components of π^ij, e.g. 'xx' for π^xx, and values are array-like objects of shape (n_volume_elements,) containing the appropriate components in units of GeV/fm³.

2D surfaces require components 'xx', 'xy', 'yy'. 3D surfaces additionally require 'xz', 'yz'. The remaining components are computed internally by tracelessness and orthogonality to the velocity.
Pi – bulk viscous pressure

Array-like object of shape (n_volume_elements,) containing the bulk pressure in units of GeV/fm³.
ymax – maximum momentum rapidity

For 2D surfaces, momentum rapidity is sampled uniformly from −ymax to +ymax. The default is 0.5. Note that the number of particles per unit rapidity is a physical property, so the average number of produced particles is proportional to ymax.

Has no effect for 3D surfaces.

Class attributes

boost_invariant¶: Whether the surface is boost-invariant or not (2D or 3D).

volume¶: The total volume of all elements [fm³].

class frzout.HRG(T, species='all', res_width=True, decay_f500=False)¶

A hadron resonance gas.

T sets the temperature of the gas in GeV.

species sets the composition of the gas, i.e. the list of particle species. It may be one of:

'all' (default) – all known particle species
'urqmd' – species in UrQMD
'id' – identified particles (pions, kaons, and protons)
a sequence of particle IDs (do not include antiparticle IDs — they will be added automatically)

If res_width is True (default), the mass width of resonances is taken into account: Mass distributions are integrated when computing thermodynamic quantities, and the masses of resonances are randomly chosen during sampling. Otherwise, mass distributions are ignored, and every sampled resonance is assigned its pole mass. Enabling this option increases computation time.

If decay_f500 is True, f₀(500) resonances are decayed into pion pairs. The default is False. This is a workaround for hadronic afterburners that are not aware of the f₀(500), such as UrQMD. This option is automatically enabled if species is set to 'urqmd'. If for some reason this is not desirable, set species to frzout.species.urqmd, which does not enable decay_f500.

Thermodynamic quantities

The following methods compute thermodynamic quantities of the HRG, given its temperature, composition, and whether or not resonance width is enabled.

density()¶: Particle density [fm^-3].

energy_density()¶: Energy density [GeV/fm³].

pressure()¶: Pressure [GeV/fm³].

entropy_density()¶: Entropy density [fm^-3].

mean_momentum()¶: Average magnitude of momentum [GeV].

cs2()¶: Speed of sound squared [relative to speed of light].

eta_over_tau()¶: Shear viscosity / relaxation time [GeV/fm³].

zeta_over_tau()¶: Bulk viscosity / relaxation time [GeV/fm³].

Bulk viscous corrections

The following methods relate to the parametric bulk correction method (see Viscous corrections).

bulk_scale_factors(Pi)¶: Returns a tuple (nscale, pscale) of the density and momentum scale factors at the given bulk pressure Pi [GeV/fm³].

Pi_lim()¶: Returns a tuple (Pi_min, Pi_max) of the limiting values of bulk pressure [GeV/fm³].

Sampling function¶

frzout.sample(surface, hrg)¶

Sample particles using the given Surface and HRG.

Returns particle data as a numpy structured array with fields

'ID' – particle ID number
'x' – position four-vector (t, x, y, z)
'p' – momentum four-vector (E, p_x, p_y, p_z)

Species information¶

Species information is read from the table http://pdg.lbl.gov/2017/mcdata/mass_width_2017.mcd, which is updated annually by the Particle Data Group (PDG). The file is included with frzout.

Particle IDs follow the PDG numbering scheme (http://pdg.lbl.gov/2017/reviews/rpp2017-rev-monte-carlo-numbering.pdf).

frzout.species_dict¶

Dictionary of species data.

The keys are particle IDs; values are subdicts containing the following keys:

'name' – name of the species
'mass' – mass in GeV
'width' – width in GeV
'degen' – spin degeneracy
'boson' – whether the species is a boson or fermion
'charge' – electric charge
'has_anti' – whether or not the species has a corresponding antiparticle

For example, frzout.species_dict[211]['mass'] is the mass of a charged pion (ID = 211).

Use list(frzout.species_dict) to obtain the list of all particle IDs.

frzout.species.identified¶: List of IDs of the standard identified particles (pions, kaons, protons).

frzout.species.urqmd¶: List of particle IDs in UrQMD.