# 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),$

where

$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.

## 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.

### Collider¶

class trento::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.

Example:

Collider collider{var_map};
collider.run_events();


Public Functions

Collider(const VarMap &var_map)

Instantiate from the configuration.

void run_events()

Run events and output.

### Event¶

class trento::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.

Example:

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


Public Functions

Event(const VarMap &var_map)

Instantiate from the configuration.

void compute(const Nucleus &nucleusA, const Nucleus &nucleusB, NucleonProfile &profile)

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 entropy. May be interpreted as $$dS/dy$$ or $$dS/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.

### Output¶

class trento::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
• const Event& Event object

### Nucleus¶

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:

Example:

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


Return
a smart pointer std::unique_ptr<Nucleus>
Parameters
• 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)
Exceptions
• std::invalid_argument: for unknown species

class trento::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.

Parameters
• A: number of nucleons

void set_nucleon_position(iterator 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 trento::Proton

Trivial nucleus with a single nucleon.

Inherits from trento::Nucleus

class trento::Deuteron

Samples pairs of nucleons from the Hulthén wavefunction

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

http://inspirehep.net/record/1261055

Inherits from trento::Nucleus

class trento::WoodsSaxonNucleus

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.

Inherits from trento::MinDistNucleus

class trento::DeformedWoodsSaxonNucleus

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.

Inherits from trento::MinDistNucleus

class trento::ManualNucleus

Reads manual nuclear configurations from an HDF5 file.

Inherits from trento::Nucleus

### Nucleon¶

class trento::Nucleon

Represents a single nucleon. Stores its 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

Nucleon()

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

double x() const

The transverse x position.

double y() const

The transverse y position.

double z() const

The longitudinal z position.

bool is_participant() const

Whether or not this nucleon is a participant.

Private Functions

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

Set the position and reset participant status to false.

void set_participant()

Mark as a participant.

### Nucleon profile¶

class trento::NucleonProfile

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

Public Functions

NucleonProfile(const VarMap &var_map)

Instantiate from the configuration.

double radius() const

The radius at which the nucleon profile is truncated.

double max_impact() const

The maximum impact parameter for participation.

void fluctuate()

Randomly fluctuate the profile. Should be called prior to evaluating the thickness function for a new nucleon.

double thickness(double distance_sqr) const

Compute the thickness function at a (squared) distance from the profile center.

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

Randomly determine if a pair of nucleons participates.

### Fast exponential¶

template <typename T = double>
class trento::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)$$.

Example:

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).