Tests & internals¶

Unit tests¶

https://travis-ci.org/Duke-QCD/trento.svg?branch=master

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

Note

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.

Collider¶

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.

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

Output¶

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

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");
double radius = lead_nucleus->radius();
lead_nucleus->sample_nucleons(0);
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 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(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.\]

http://inspirehep.net/record/1261055

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.

NucleonData¶

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

NucleonData()¶: 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.

NucleonCommon¶

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

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