Tutorials

Initialisation and finalisation of the physics

Prior to any usage physics tables must be initialised. This can be done in two ways:

1. PUMAS physics can be created from a Materials Description File (MDF), using the pumas_physics_create function, e.g. as:

   struct pumas_physics * physics;
   pumas_physics_create(&physics, PUMAS_PARTICLE_MUON,
       "examples/data/materials.xml", NULL, NULL);
   

   The 3^nd argument to pumas_physics_create specifies the path to the MDF.

   Note

   The 4^th and 5^th arguments to pumas_physics_create are optionnal. The 4^th argument allows one to indicate a folder where stopping power tables are looked-for and/or generated to. Relative paths w.r.t. the MDF location can be given by prepending the path with an @. Ordinary paths are assumed otherwise. This argument can also be left NULL in which case the base directory of the MDF is used.

   The 5^th argument can be left empty for most applications (NULL) resulting in PUMAS using its default physics settings. See below for more advanced usage.

2. PUMAS physics can also be loaded from a binary dump of a a previously generated physics instance. Loading is done using the pumas_physics_load function, e.g. as:

   struct pumas_physics * physics;
   pumas_physics_load(&physics, "materials/dump");
   

When the physics is initialised from a MDF, PUMAS pre-computes and tabulates properties for all the materials defined in the MDF file, should there be latter used or not. During this operation intermediary stopping power tables are generated in the PDG format. These tables can be browsed using any text editor. Fully tabulating the physics can last several seconds depending on the number of materials defined and on your machine performances.

Once the physics has been created, it can be saved in binary format with the pumas_physics_dump function. Initialising the physics from such a binary dump is much faster. However, the dump is a low level representation of the physics depending on the machine and OS used. Thus, it is safer to use MDFs as an exchange format rather than binary dumps.

When the physics ressources are no more needed they can be released using the pumas_physics_destroy function. From there on it is possible to re-initialise (load) the physics again.

To sum up, the following code is an elementary program using the PUMAS library. It computes the low level material tables used by PUMAS and dumps them to a binary file.

#include "pumas.h"

int main()
{
    /* Initialise PUMAS physics from a MDF */
    struct pumas_physics * physics;
    pumas_physics_create(&physics, PUMAS_PARTICLE_MUON,
        "examples/data/materials.xml", NULL, NULL);

    /* Dump the materials data to a binary file */
    FILE * fid = fopen("examples/data/materials.pumas", "wb+");
    pumas_physics_dump(physics, fid);
    fclose(fid);

    /* Release PUMAS memory */
    pumas_physics_destroy(&physics);

    return 0;
}

Error handling

By default the PUMAS library is configured to print out to stderr whenever an error occurs and to resume back to the OS. While this is practical for short programs it is anoying when integrating PUMAS in a larger project. Therefore, the user can override the default error handling by providing its own handler callback. This is done with the pumas_error_handler_set function. A brief human readable description of the error is forwarded to the handler. The following illustrates the corresponding mechanism with a basic example.

#include "pumas.h"
#include <stdlib.h>

/* Error handler for PUMAS with a hard exit. Note that this is actually
 * the behaviour of the default PUMAS error handler
 */
static void error_handler(
    enum pumas_return rc, pumas_function_t * caller, const char * message)
{
        /* Dump the error summary */
        fputs("pumas: library error. See details below\n", stderr);
        fprintf(stderr, "error: %s\n", message);

        /* Exit to the OS */
        exit(EXIT_FAILURE);
}

int main()
{
    /* Set the error handler callback or disable it at all by providing
     *  a `NULL` pointer instead
     */
    pumas_error_handler_set(&error_handler);

    /* Do something with PUMAS */
    ...
}

Simulation context

Monte-Carlo simulations with PUMAS are performed within isolated contexts. This allows the transport engine to handle multiple simulation streams simultaneously. A simulation stream is created using the pumas_context_create function. The allocated memory is released using the pumas_context_destroy function. The pumas_context_create function returns a handle to a semi-opaque pumas_context structure. Note that it must not be re-allocated since it refers to internal (hidden) data required by the transport engine. If extra (user) data must be attached to the simulation context, this can be done at the context creation by reserving additional memory for them.

A simulation stream can be configured directly by modifying the pumas_context data. In particular the following flags allow one to customise the transport:

By default each simulation context ships with its own pseudo random stream based on a Mersenne Twister algorithm seeded from the OS, e.g. from /dev/urandom on UNIX. A specific random seed can be provided with the pumas_context_random_seed_set function. The default pseudo random engine can also be overriden by providing an alternative pumas_random_cb to the simulation context generating a pseudo random number in [0,1] uniformly. Then, it is the user's responsibility to ensure that the random stream is thread safe if multiple simulation contexts are used simultaneously.

Below is an example of context creation with the built-in pseudo random engine. The context is configured with a mixed Monte Carlo for energy loss but with scattering disabled. This is simular e.g. to MUM. For the sake of clarity errors are handled by the default PUMAS error handler.

#include "pumas.h"
#include <stdlib.h>

int main()
{
    /* Initialise the physics from a MDF */
    struct pumas_physics * physics;
    pumas_physics_create(&physics, PUMAS_PARTICLE_MUON,
        "examples/data/materials.xml", NULL, NULL);

    /* Create a new simulation context */
    struct pumas_context * context = NULL;
    pumas_context_create(&context, physics, 0);

    /* Configure the context */
    context->mode.energy_loss = PUMAS_MODE_MIXED;
    context->mode.scattering = PUMAS_MODE_DISABLED;

    /* Release the allocated memory */
    pumas_context_destroy(&context);
    pumas_physics_destroy(&physics);

    exit(EXIT_SUCCESS);
}

Monte-Carlo state and transport

An elementary Monte-Carlo (particle) state in PUMAS is summarised by a pumas_state structure. It contains the minimal set of data required by PUMAS for performing a Monte-Carlo transport. Extra user data can be added by embedding this structure in a larger one.

Prior to any usage the particle state must be initialised. One has to provide an electric charge number, a kinetic energy, a Monte-Carlo weight and a starting position and momentum direction. For example, below is an example of static initialisation:

struct pumas_state state = {
    .charge = -1.,
    .energy = 1E+01,
    .weight = 1.,
    .position = { 0., 0., 10. }
    .direction = { 0., 0., -1. } };

In addition, during the simulation, the Monte-Carlo state carries information on the traveled distance, on the corresponding grammage, and on the spent proper time.

A Monte-Carlo state is transported (propagated) using the pumas_contex_transport function. This requires a simulation stream, i.e. providing a properly configured pumas_context. Once everything is configured, transporting the Monte-Carlo particle is as simple as:

enum pumas_event event;
struct pumas_medium * medium[2];
pumas_context_transport(context, &state, &event, medium);

The two last arguments are optionnal and can be set to NULL. The pumas_context_transport call returns if:

• a user specified event occurs, e.g. a user supplied limit has been reached,

• the particle escapes the simulation area, i.e. a NULL pumas_medium is reached (see geometry below),

• the particle decays (if enabled),

• an error occurs, e.g. a non unit direction was provided.

Limits can be set on the particle traveled distance or grammage, on its spent proper time or on its kinetic energy. This is done per simulation stream by setting the event flag of the pumas_context to one or more of PUMAS_EVENT_LIMIT_DISTANCE, PUMAS_EVENT_LIMIT_GRAMMAGE, PUMAS_EVENT_LIMIT_ENERGY or PUMAS_EVENT_LIMIT_TIME. Combinations of limits are specified with the binary or operator (|). In addition a strictly positive value must be provided for the corresponding limit field(s) of the simulation context. For example the code below sets a limit of 1 km on the particle traveled distance:

/* Enable limits on the traveled distance */
context->event |= PUMAS_EVENT_LIMIT_DISTANCE;

/* Set a limit of 1 km on the traveled distance */
context->limit.distance = 1E+03;

Below is an example of forward transport where the simulation mode is adapted depending on the particle kinetic energy.

#include <float.h>
#include "pumas.h"

int main()
{
    /* Initialise the physics, create a simulation stream and a
     * Monte-Carlo state
     */
    ...

    /* Enable kinetic energy limits */
    context->event |= PUMAS_EVENT_LIMIT_ENERGY;

    /* Transport the Monte-Carlo state */
    const double energy_min = 1E-03;
    while (state.energy > energy_min) {
        if (state.energy < 1E+02 + FLT_EPSILON) {
            /* Below 100 GeV do a detailed simulation
             * à la Geant4 including scattering
             */
            context->mode.energy_loss = PUMAS_MODE_STRAGGLED;
            context->mode.scattering = PUMAS_MODE_MIXED;
            context->limit.energy = energy_min;
        } else {
            /* Do a fast simulation à la MUM */
            context->mode.energy_loss = PUMAS_MODE_MIXED;
            context->mode.scattering = PUMAS_MODE_DISABLED;
            context->limit.energy = 1E+02;
        }
        pumas_context_transport(context, &state, NULL, NULL);
    }

    /* Finalise PUMAS */
    ...
}

Specifying a Geometry

PUMAS is a pure transport engine. It does not provide geometric primitives, e.g. boxes, orbs, tubes, etc ... like Geant4. Instead it relies on a simple callback mechanism for describing the propagation media. This mechanism is decribed in the following.

In PUMAS, a geometry is defined by a user supplied pumas_medium_cb callback. This callback must map a Monte Carlo location to a pumas_medium. A NULL medium indicates the outer of the simulation geometry. In addition, the pumas_medium_cb callback must also indicate a maximum geometric stepping distance. Typically, this is the distance to the next medium assuming a straight line propagation. Finally the callback must return a pumas_step enum indicating if the proposed step needs to be cross-checked by PUMAS or if it should be used as is. Managing steps that end on a geometry boundary can be tricky numerically. Therefore it is recommended to return PUMAS_STEP_CHECK if you are unsure of what to do since it is more robust. The raw mode is usefull if your geometry engine already performs those checks in order to avoid double work.

A pumas_medium is described by a material index and a set of local properties. The material must be the same over the whole medium. The material index can be mapped to an explicit name, as defined in the MDF, with the pumas_physics_material_index function. CPU wise, this is best done during the initialisation and configuration stage since the mapping implies strings comparisons.

The medium local properties are its bulk density and any magnetic field. Those are set by a pumas_locals_cb callback. If not constant, they must vary continuously over the medium. The pumas_locals_cb callback can be NULL in which case the material default density is used as defined in the MDF. In addition, the magnetic field is set to zero in this case. If a pumas_locals_cb is provided then it must return a length, in m, consistent with the size of the propagation medium inhomogeneities, e. g. $|\frac{\rho}{\nabla \rho}|$ for a density gradient. Returning a value of zero or less indicates a uniform medium but e.g. with a non standard density or with a magnetic field.

Below is an example of a uniform medium of infinite extension in PUMAS. Although this simple geometry might seem rather complicated to implement with PUMAS, this callback mechanism is however flexible enough in order to deal with more complicated geometries, e.g. a non uniform atmosphere as in the examples. Furthermore, using this callback mechanism existing geometry engines can be integrated with PUMAS.

/* A uniform medium without magnetic field. Note that we could also set the
 * locals callback to `NULL` instead if using the material's default density.
 */
double locals(struct pumas_medium * medium,
    struct pumas_state * state, struct pumas_locals * locals)
{
        /* Set the medium density in kg/m^3. Set this to zero or less in
         * order to use the material's default density.
         */
        locals->density = 2.65E+03;

        /* Propose a maximum stepping distance. Returning zero or less indicates
         * a uniform medium
         */
        return 0.;
}

/* A simple medium of infinite extension */
enum pumas_step medium(struct pumas_context * context,
    struct pumas_state * state, struct pumas_medium ** where, double * step)
{
    static struct pumas_medium medium_ = { .material = 0, .locals = &locals };

    if (where != NULL) /* Beware that `where` is `NULL` if not requested */
        *where = &medium_;

    /* Propose a maximum stepping distance. Providing zero or less indicates
     * no boundaries.
     */
    if (step != NULL) /* Beware that `step` is `NULL` if not requested */
        *step = 0;

    /* Notify PUMAS that the returned step length should be used raw without
     * further cross-checks. Note that for more complex geometries you might
     * want to use `PUMAS_STEP_CHECK` instead which is more robust.
     */
    return PUMAS_STEP_RAW;
}

Practical Backward Monte-Carlo and weights

In Backward Monte-Carlo (BMC) mode PUMAS allows one to compute a point estimate of the muon (tau) flux for a given final state, see e.g. lemma 1 of the Backward Monte Carlo paper. In order to sample a distribution of final states, one relies on corollary 4. I.e. the final states are generated over a bias distribution as for classical Importance Sampling. In practice, the procedure can be described as below:

1. At generation the particle weight is initialised to $\omega = \frac{1}{ \text{PDF}_\text{gen}}$, where $\text{PDF}_\text{gen}$ is the generation PDF. Let say that one generates final states over a surface of interest and with some solid angle of aperture but at a fixed energy. Then $\text{PDF}_\text{gen}$ could be in units $\text{m}^{-2}\text{sr}^{-1}$. Following the weight would start with unit $\text{m}^2\text{sr}$.

2. The backward propagation modifies the weight by a unitless Jacobian factor, $\text{J}_{i,f}$, due to the change in the integration variable for the flux from initial state to final state. The particle’s weight is updated by PUMAS step by step.

3. When reaching the primary flux surface, the weight must be multiplied by the corresponding flux, $\Phi_0$. The units used for this flux must be consistent with the one used for $\text{PDF}_\text{gen}$. So in the present case it could be $\text{GeV}^{-1}\text{m}^{-2}\text{s}^{-1}\text{sr}^{-1}$. Following, the final weight would have unit $\text{GeV}^{-1}\text{s}^{-1}$. But, if the final state energy is randomised as well, the weight's unit would be $\text{s}^{-1}$, i.e. a rate of events.

To sum up the particle weight is:

This result is similar to a classical Importance Sampling procedure except for the extra Jacobian factor, $\text{J}_{i,f}$, due to the backward transport.

For illustration, below is a simplified example of backward Monte-Carlo integration over the kinetic energy. A 1 / E bias distribution is used, i.e. a log uniform sampling.

#include <float.h>
#include <math.h>
#include "pumas.h"

int main()
{
    /* Initialise PUMAS, create a simulation stream and a
     * Monte-Carlo state
     */
    ...

    /* Randomise the final state kinetic energy */
    const double energy_min = 1E-03;
    const double energy_max = 1E+06;
    const double r = log(energy_max / energy_min);
    state.energy = energy_min * exp(r * context->random(context));
    state.weight = state.energy * r; /* 1 / PDF_{gen}(k_f) */

    /* Backward transport the Monte-Carlo state. Note that at exit the
     * Monte-Carlo weight is updated by the Jacobian BMC factor, J_{i,f}.
     */
    pumas_context_transport(context, &state, NULL, NULL);

    /* Sample the primary flux */
    state.weight *= flux(state.energy); /* Phi(k_i) */

    /* Finalise PUMAS */
    ...
}

Tuning the physics

By default PUMAS is configured in order to deliver accurate yet fast results for muography applications. But for specific use cases or in order to estimate systematics it can be useful to modify PUMAS default settings. To do so there are three relevant parameters.

• First, one can provide alternative differential cross-sections (DCS) models for radiative energy loss processes, i.e. Bremsstrahlung, $e^+e^-$ pair production and photonuclear interactions. This is done when creating physics tables, e.g. with the pumas_physics_create function, by providing a struct pumas_physics_settings as 5^th argument. By default PUMAS uses Sandrock, Soedingrekso and Rhode's parametrizations for the Bremsstrahlung and $e^+e^-$ pair production processes. For photonuclear interactions the DRSS cross-section (Dutta et al.) is used. Alternative models available in PUMAS are described in the API documentation.

• When creating the physics tables with a struct pumas_physics_settings as 5^th argument one can also specify a custom cutoff value, $x_\text{cut}$, between continuous energy loss or discrete processes. By default PUMAS uses $x_\text{cut} = 5\%$ which is a good compromise between speed and accuracy for the transport of a continuous flux of muons (see e.g. Sokalski et al. or Koehne et al.).

  Note

  Modifying $x_\text{cut}$ only affects the PUMAS_MODE_MIXED and PUMAS_MODE_STRAGGLED energy loss modes. If CSDA is used then $x_\text{cut}$ is forced to 100%, i.e. all losses are continuous.

  Warning

  In backward mode, with mixed or straggled energy loss, cutoff values lower than 1% ($x_\text{cut} < 1\%$) are not currently supported.

• Finally, the accuracy of the Monte Carlo stepping can be modified with the accuracy parameter of the simulation context. PUMAS uses a mixed Monte Carlo algorithm (Fernandez-Varea et al.) for rendering the multiple scattering, the magnetic bending etc. Lowering the accuracy parameter results in shorter (more accurate) steps but at the cost of extra CPU time.