Tests & internals

Unit tests


TRENTo has a thorough test suite. To run the tests, build the code in debug mode:

mkdir build-debug && cd build-debug
cmake .. -DCMAKE_BUILD_TYPE=Debug
make tests

Tests rely on the Catch C++ test framework, which CMake downloads automatically during the build process. If for some reason this fails, manually download the header catch.hpp and place it in the test directory.

If gcovr is installed a code coverage report will be generated.

If any sanitizers are available they will run. The sanitizers check for runtime errors such as memory leaks. They are included in recent versions of GCC and Clang.

Fitting the cross section

It is vital that the model reproduce the inelastic nucleon-nucleon cross section \(\sigma_{NN}\). This condition may be written

\[\sigma_{NN} = \int d^2b \, P_\text{coll}(b),\]


\[P_\text{coll}(b) = 1 - \exp[-\sigma_{gg}T_{AB}(b)]\]

is the probability of two nucleons colliding at impact parameter \(b\).

\(T_{AB}\) is the overlap integral of the two nucleons’ thickness functions, easily computed given Gaussian nucleons of width \(w\):

\[\begin{split}\begin{align*} T_A(x, y) &= \frac{1}{2\pi w^2} \exp\biggl( -\frac{x^2 + y^2}{2w^2} \biggr), \\[1em] T_{AB}(b) &= \int dx \, dy \, T_A(x - b/2, y) T_B(x + b/2, y) \\ &= \frac{1}{4\pi w^2} \exp\biggl( -\frac{b^2}{4w^2} \biggr). \end{align*}\end{split}\]

\(\sigma_{gg}\) is an effective parton-parton cross section now determined by the relation

\[\sigma_{NN} = \int_0^{b_\text{max}} 2\pi b \, db \, \biggl\{ 1 - \exp\biggl[ -\frac{\sigma_{gg}}{4\pi w^2} \exp\biggl( -\frac{b^2}{4w^2} \biggr) \biggr] \biggr\},\]

where \(b_\text{max}\) is the maximum impact parameter for a collision. Let \(b_\text{max} = Aw\), i.e. some number of nucleon widths (the actual code uses \(A = 6\)).

After appropriate change of variables, this relation may be written

\[\frac{\sigma_{NN}}{4\pi w^2} = \frac{A^2}{4} + \text{Ei}\biggl( -e^{-A^2/4} \frac{\sigma_{gg}}{4\pi w^2} \biggr) - \text{Ei}\biggl( -\frac{\sigma_{gg}}{4\pi w^2} \biggr)\]

where \(\text{Ei}\) is the exponential integral. This is still a transcendental equation but it can be quickly solved numerically for a given cross section and nucleon width.

The code determines \(\sigma_{gg}\) at runtime by solving this equation using a standard root finding algorithm. The relevant function is compute_cross_sec_param in nucleon.cxx. Note that it actually optimizes over the dimensionless variable \(\log({\sigma_{gg}/4\pi w^2})\).

The nucleon unit test verifies that the cross section is accurately reproduced.


The semi-analytic cross section procedure described above only works for a collision involving two Gaussian particles. If nucleon substructure is enabled, i.e. constituent number m > 1, the code uses a Monte Carlo algorithm to fit the cross section parameter \(\sigma_{gg}\) numerically. This procedure is somewhat slow, and it scales poorly with the number of constituents. For m > 10, it could easily take a minute or longer, so the code caches the value of \(\sigma_{gg}\) in the user’s XDG_DATA_HOME directory.

List of classes

This section is automatically generated from the source code by Doxygen and converted into Sphinx format by Breathe. In order to emphasize functionality—rather than implementation details—private class methods are not shown.


class Collider

Orchestrates event generation and output. Owns instances of several other TRENTO classes and knows how they work together. Responsible for sampling impact parameters. After instantiation, call run_events() to do everything.


Collider collider{var_map};

Public Functions

Collider(const VarMap &var_map)

Instantiate from the configuration.

void run_events()

Run events and output.


class Event

The primary computation class, responsible for constructing nuclear thickness functions and calculating event observables. Designed to be created once and used many times by repeatedly calling compute(). Stores its observables internally and provides inspector methods.


Event event{var_map};
for (int n = 0; n < nevents; ++n) {
  event.compute(nucleusA, nucleusB, nucleon_profile);

Public Functions

Event(const VarMap &var_map)

Instantiate from the configuration.

void compute(const Nucleus &nucleusA, const Nucleus &nucleusB, const NucleonCommon &nucleon_common)

Compute thickness functions and event observables for a pair of Nucleus objects and a NucleonProfile. The nuclei must have already sampled nucleon positions and participants before passing to this function.

const int &npart() const

Number of nucleon participants.

const double &multiplicity() const

Multiplicity—or more specifically, total integrated reduced thickness. May be interpreted as \(dS/d\eta\) or \(dE/d\eta\) at midrapidity.

const std::map<int, double> &eccentricity() const

Eccentricity harmonics \(\varepsilon_n\) for n = 2–5. Returns a map of \((n : \varepsilon_n)\) pairs, so e.g.:

double e2 = event.eccentricity().at(2);

const Grid &reduced_thickness_grid() const

The reduced thickness grid as a square two-dimensional array.


class Output

Simple interface for outputting event data. Determines which output formats to create based on the configuration and writes those formats when called.

Public Functions

Output(const VarMap &var_map)

Instantiate from the configuration.

template<typename ...Args>
void operator()(Args&&... args) const

Call the functor to output event data. Arguments are perfect-forwarded to each output function. The required arguments are

  • int event number

  • double impact parameter

  • int binary collisions

  • const Event& Event object


NucleusPtr trento::Nucleus::create(const std::string &species, double nucleon_dmin = 0)

The canonical way to create a Nucleus. Constructs the appropriate subclass for the given species. Sets the Woods-Saxon parameters for Au, Pb, U, etc; parameters copied from the PHOBOS Glauber model:


std::unique_ptr<Nucleus> lead_nucleus = Nucleus::create("Pb");
double radius = lead_nucleus->radius();
for (const auto& nucleon : *lead_nucleus)


a smart pointer std::unique_ptr<Nucleus>

  • species: standard symbol, e.g. “p” for proton or “Pb” for lead-208

  • nucleon_dmin: minimum nucleon-nucleon distance for Woods-Saxon nuclei (optional, default zero)

  • std::invalid_argument: for unknown species

class Nucleus

Interface class to all nucleus types. Stores an ensemble of nucleons and randomly samples their positions. Implements a standard iterator interface through begin() and end() functions. Iterating over a Nucleus means iterating over its Nucleon members.

Subclassed by trento::Deuteron, trento::ManualNucleus, trento::MinDistNucleus, trento::Proton

Public Functions

virtual ~Nucleus()

Default virtual destructor for abstract base class.

virtual double radius() const = 0

The “radius”, i.e. the maximum distance at which a nucleon could be placed.

void sample_nucleons(double offset)

Sample a new ensemble of nucleon positions with the given offset in the x-direction.

Protected Functions

Nucleus(std::size_t A)

Constructor only accessible by derived classes.

  • A: number of nucleons

void set_nucleon_position(NucleonData &nucleon, double x, double y, double z)

Set a Nucleon position. This function provides a consistent interface to derived classes and ensures the position is correctly offset. Nucleus is a friend of Nucleon and therefore able to set nucleon positions; the derived classes must use this function to set positions.

Nucleus types

class Proton : public trento::Nucleus

Trivial nucleus with a single nucleon.

class Deuteron : public trento::Nucleus

Samples pairs of nucleons from the Hulthén wavefunction

\[f(r) \propto \biggl( \frac{\exp(-ar) - \exp(-br)}{r} \biggr)^2.\]


class WoodsSaxonNucleus : public trento::MinDistNucleus

Samples nucleons from a spherically symmetric Woods-Saxon distribution

\[f(r) \propto \frac{1}{1 + \exp(\frac{r-R}{a})}.\]

For non-deformed heavy nuclei such as lead.

class DeformedWoodsSaxonNucleus : public trento::MinDistNucleus

Samples nucleons from a deformed spheroidal Woods-Saxon distribution

\[f(r, \theta) \propto \frac{1}{1 + \exp(\frac{r-R(1+\beta_2Y_{20}+\beta_4Y_{40})}{a})}.\]

For deformed heavy nuclei such as uranium.

class ManualNucleus : public trento::Nucleus

Reads manual nuclear configurations from an HDF5 file.


class NucleonData

Represents a single nucleon. Stores its transverse position and whether or not it’s a participant. These properties are globally readable, but can only be set through Nucleus and NucleonProfile.

Public Functions


Only a default constructor is necessary—the class is designed to be constructed once and repeatedly updated.

bool is_participant() const

Whether or not this nucleon is a participant.

bool constituents_exist() const

Whether or not its constituents have been sampled.

double x() const

The transverse x position.

double y() const

The transverse y position.

double z() const

The longitudinal z position.

Private Functions

void set_position(double x, double y, double z)

Set the transverse position and reset participant status to false.


class NucleonCommon

Encapsulates properties shared by all nucleons: transverse thickness profile, cross section, fluctuations. Responsible for sampling nucleon-nucleon participation with given \(\sigma_{NN}\).

Public Functions

NucleonCommon(const VarMap &var_map)

Instantiate from the configuration.

double max_impact() const

The maximum impact parameter for participation.

std::array<double, 4> boundary(const NucleonData &nucleon) const

Corners of the tile enclosing the nucleon thickness.

double thickness(const NucleonData &nucleon, double x, double y) const

Nucleon thickness as a function of transverse position.

bool participate(NucleonData &A, NucleonData &B) const

Randomly determine if a pair of nucleons participates.

Private Functions

void sample_constituent_positions(NucleonData &nucleon) const

Sample constituent positions inside the nucleon.

void set_participant(NucleonData &nucleon) const

Flag the nucleon as a participant.

Fast exponential

template<typename T = double>
class FastExp

Fast exponential approximation, to be used as a drop-in replacement for std::exp when it will be evaluated many times within a fixed range. Works by pre-tabulating exp() values and exploiting the leading-order Taylor expansion; for step size \(dx\) the error is \(\mathcal O(dx^2)\).


FastExp<double> fast_exp{0., 1., 11};  // tabulate at 0.0, 0.1, 0.2, ...
fast_exp(0.50);  // evaluate at table point -> exact result
fast_exp(0.55);  // midway between points -> worst-case error

Public Functions

FastExp(T xmin, T xmax, std::size_t nsteps)

Pre-tabulate exp() values from xmin to xmax in nsteps evenly-spaced intervals.

T operator()(T x) const

Evaluate the exponential at x (must be within range).