PUMAS v1.2

This release improves the sampling efficiency of backward collisions, for radiative processes. It also introduces a breaking change to the API.

Full details on PUMAS v1.2 physics are published in CPC (arXiv:2206.01457).

Physics update

• The power-law exponent for sampling backward decays has been updated from 1.4 (previously) to 2.0. The new value fits radiative DCSs for fractional energy losses approximately in [1, 10]%. As a result, the Monte Carlo efficiency of backward collisions is improved.

• For very long distances and small energy cutoffs, the current sampling method for backward collisions is numerically problematic. Therefore, for now PUMAS forbids using energy cutoffs values lower than 1% in backward mixed or straggled modes. The latter value has been checked to be numerically robust at least up to 10 km of standard rock. Note that in forward mode, arbitrary cutoffs values can still be used.

API update

• The pumas_version function now has an additional patch argument. I.e., full version number uses the traditonal MAJOR.MINOR.PATCH semantic. However, GitHub releases are still limited to minor version changes.

Bug patches

• Patch 1 corrects a memory alignment issue with some compilers / versions combinations when using agressive optimisation settings, typically -O3 with gcc. Note that this bug appeared in v1.1 already.

• Patch 3 should slightly improve the accuracy for small steps, at high energy and in CSDA mode.

v1.1

This is a major release improving the physics, the numeric accuracy and the quality of life. It introduces breaking changes to the API.

Full details on PUMAS v1.1 physics are published in [CPC][cpc] ([arXiv:2206.01457][arxiv]).

Warning : due to a typo in the simulation of hard elastic collisions, this version has been patched since it initial release (b16e7fe). Please check your git hash in order to ensure that you are using the latest patch, i.e. b52648f . The typo resulted in too large scattering, the higher the energy of the muon. Note that it only impacts muons with kinetic energy higher than 100 GeV.

Physics update

• The density effect in electronic collisions is now computed directly from binding energies following [Fano (1963)][Fano63]. In the case of standard rock the electronic structure of CaCO3 is assumed, see e.g. [Groom et al. (2001)][Groom01]. Previously the density effect was parametrized according to [Sternheimer et al. (1984)][Sternheimer84]. As a result, at GeV energies the stopping power in standard rock given by PUMAS is now 0.4% higher than the PDG one i.e. [Groom et al. (2001)][Groom01]. Note that both models are consistent with available experimental data on the density effect.

• Concerning elastic collision, the description of atomic charge densities has been improved following [Salvat et al. (1987)][Salvat87]. For the nuclear charge a more accurate model is used as well with radii from [De Vries et al. (1987)][DeVries87].

• Concerning radiative process, several models have been implemented in PUMAS following the work done for [PROPOSAL][PROPOSAL]. By default, the Sandrock, Soedingrekso and Rhode (SSR) parametrizations are used for Bremsstrahlung and pair production, see e.g. [Soedingrekso et al. (2019)][Soedingrekso19] . Previously, the Kelner, Kokoulin and Petrukhin (KKP) parameterizations were used. Those are still available, but they have been updated with improved radiation logarithms from [Kelner et al. (1999)][Kelner99]. This leads to a 2% decrease of the cross-section in standard rock.

  A description of available parametrizations can be found in the API documentation (or in include/pumas.h).

• The angular distribution in radiative collisions is now described directly from the Doubly Differential Cross-Sections (DDCs). For bremsstrahlung and pair production the DDCS of [Tsai (1974)][Tsai74] is used, rescaled to the selected parametrizations. In the case of photonuclear interactions, when available the DDCS in Q² is used. Otherwise, the DRSS DDCS is used rescaled to the selected parametrization.

  The angular deflection in hard collisions is randomised from the DDCS using rejection sampling. In addition, the contribution of all soft radiative collisions has been properly included in the multiple scattering. Previously only soft photonuclear interactions were considered with an approximate model.

Numeric improvements

• Physics tables are now interpolated using cubic Hermite polynomials following the method of [Fritsch and Butland (1984)][Fritsch84]. This method preserves the local monotonicity of the tabulations. Previously, a linear interpolation was used instead. This could lead to unpleasant results due to the discontinuity of the derivative at nodes. In addition, the new method reduces the interpolation error below 0.1%. Before, errors of ~1% could occur. However, the new interpolation method is longer than before resulting in a ~20% slow down at runtime.

• The integration accuracy of integral quantities has been increased by various means, e.g. by exploiting the new cubic interpolation or by increasing the order of Gaussian quadratures. Overall, the integration errors have been reduced below 0.1%. As a result, the initialisation of physics tables is ~10 times longer than before. Yet, the runtime speed is not impacted by this specific accuracy increase.

• Several inconsistencies related to rounding errors or overflows have been patched. This concerns extreme cases. E.g. previously, in rare cases where a random number was generated very close to 1, due to a numeric rounding error a particle could cross infinite depths of rocks without any hard collision. This was leading to wrong values of the transmitted flux for very large depths, ~10 km of standard rock.

• The sampling in forward mode is now 100% un-biased. A rejection sampling method is used with pre-computed envelopes for speed-up. Thus, PUMAS Monte-Carlo is now analog in forward mode provided that the projectile decay is randomised instead of being weighted.

API changes

pumas_physics

• The initialisation of the physics and the generation of stopping power tables has been simplified. The pumas_physics_create function now only requires a MDF. Stopping power tables are automatically generated from the MDF when creating the physics. By default, if a table already exists then it is used as so when creating the physics, i.e. it is not re-computed. Thus, by default the previous behaviour is recovered in this case.

  Optional settings can be given when creating the physics as an additional pumas_physics_settings structure argument to pumas_physics_create. For example, the energy grid for tables can be specified. See the API documentation for a complete list of options.

  In particular, the cutoff between soft and hard collisions is now configurable. The elastic_ratio parameter controlling the cutoff for elastic scattering is also configurable. Its default value has been tuned to 5% resulting in a good compromise between speed and accuracy.

• For muons, the default energy grid for stopping power tables has been extended to 1 EeV. Previously, the PDG maximum value was used, i.e. 1 PeV. This was not enough for muography applications when considering large rock depths of several km.

• Specifying a filename in the MDF for stopping power tables is now optionnal. If not provided a default filename is used according to the material name.

• Some physics properties have been renamed in order to be more canonical. E.g. The stopping power is now given by pumas_physics_property_stopping_power instead of what was previously called energy_loss.

• The pumas_physics_table_value now accepts negative values for row indicating an index relative to the last entry, e.g. as in Python. Thus, a row index of -1 corresponds to the last entry.

pumas_context

• An accuracy parameter was added to the pumas_context structure. This parameter directly controls the accuracy of the Monte Carlo stepping. The lower the more accurate but the longer the simulation. The default value of 1% has not been changed.

  Warning: the return value for pumas_locals_cb callback has changed. Now, the user should return the characteristic length of inhomogeneities, L, no more the stepping distance. The stepping distance is set internally by PUMAS given L and the value of the accuracy parameter.

• A default pseudo random number generator (PRNG) was added to pumas_context. It can be overriden by setting the random field of the context, as previously.

  A Mersenne Twister algorithm is used. Note that each context manages its own random stream. The random seed can be set (get) with pumas_context_random_seed_set (pumas_context_random_seed_get) function. If no random seed is specified it is initialised from the OS, e.g. using /dev/urandom on Unix. The state of the PRNG can be saved (loaded) with the pumas_context_random_dump (pumas_context_random_load) function.

• The Monte Carlo modes have been renamed in order to be consistent with the canonical terminology, i.e. following [Berger (1963)][Berger63]. E.g., PUMAS_MODE_MIXED and PUMAS_MODE_STRAGGLED should now be used for what was previously the HYBRID and DETAILED modes.

• The pumas_context_transport function now checks that the state direction has unit norm. If not an error is returned. PUMAS stepping algorithm requires the direction to be of norm one in order to yield correct results. Yet, this was not checked previously.

pumas_state

• The kinetic field of a Monte Carlo state has been renamed to energy. The usage is unchanged. This field still refers to the particle kinetic energy in GeV.

pumas_medium_cb

• The return flags of the pumas_medium_cb have been renamed to PUMAS_STEP_CHECK and PUMAS_STEP_RAW in order to be more explicit.

pumas_handler_cb

• The PUMAS library is now initialised with a default error handler. Previously no error handler was set. The default error handler prints the library error message to stderr and then it resumes back to the OS.

pumas_dcs_t

• Various parametrizations have been implemented for radiative processes as discussed above. Please, check the pumas_dcs API for a detailed description.

• A pumas_dcs_range function was added in order to get the range of valid energy transfer values for a radiative process.

• Radiative DCSs now check their inputs. If an invalid value is provided then the DCS returns 0. Thus, NaN should no more occur.

• A specific API has been added for elastic and electronic collisions, under pumas_elastic and pumas_electronic. These processes use built-in models in PUMAS. Related properties, e.g. the DCS or stopping power, can be accessed with the new API functions.

pumas_constant

• Additional physical constants used by PUMAS have been added as pumas_constant enumeration values.

Installation

• Building and installing PUMAS with CMake should now properly work on Windows and on OSX.

• The CMake options specific to PUMAS have been simplified and better documented. In particular, the Test build type has been removed.

Documentation...

PUMAS v1.0

This is a major release with several improvements resulting in breaking changes to the API and MDF files.

Add an explicit pumas_physics object

Previously physics data were loaded at the library initialisation and shared by all simulation contexts. Now the physics must be explicitly instanciated using the pumas_physics_create function. The syntax is similar to the previous pumas_initialise function which has been retired. The resulting pumas_physics object must be provided when creating a simulation context. Multiple physics can be used simultaneously, e.g. both μ and τ physics.

Explicit mode flags and limits

The simulation context now has explicit mode and limit attributes. Configuration options are similar but with slightly different names. E.g. PUMAS_SCHEME_DETAILED is now part of the generic enum pumas_mode flags renamed as PUMAS_MODE_DETAILED. The energy loss mode is set as context.mode.energy_loss = PUMAS_MODE_DETAILED. External limits are now collected under context.limit. The purpose of these changes is to make configuration options more explicit.

User provided DCS models

The Differential cross-sections (DCS) for radiate processes can now be overridden by the user. This can be done on a per pumas_physics basis using the pumas_physics_dcs_get and pumas_physics_dcs_set functions. Note that it is the user responsibility to use consistent energy loss tables if changing the DCS.

Accessor for physical constants

The pumas_constants function allows to get the values of some hard-coded physics constants used by PUMAS, e.g. the μ rest mass.

Provide default densities for base materials

The MDF file now requires to specify the default density for base materials. Consequently one no more needs to specify a component density for composites. This change is consistent with the fact that energy loss tables actually depend on the material density due to the density effect in ionization loss. While PUMAS allows to specify a different bulk density for the material at runtime this is actually an approximation, e.g. of a porous material filled with air. Note that it is the user's responsibility to use this feature consitively.

Accessors for elements and materials properties

The pumas_physics_element_properties and pumas_physics_material_properties allow to access the properties of atomic elements and materials used by PUMAS, e.g. the atomic mass number, A, or the material default density.

Optional locals callback

The pumas_locals_cb is now optional (can be NULL). In this case the material's default density is used. In addition, setting a null or negative density in the locals callback result in the material's default density being used.

Bug patches

• Patched a bug when looking for tabulated elements (c5ed32d).
• Patched a memory bug when using NO_LOSS scheme, now MODE_VIRTUAL (18bad78).
• Patched a sign error for non uniform end steps (8dde4a9).

Minor features

• Changed pumas_tag to pumas_version.

PUMAS v0.14

This release improves the geometry callback mechanism, e.g. allowing more efficient integration of ray tracers. A bug has been corrected for non uniform magnetic fields as well.

Stepping algorithm update

The signature of the pumas_medium_cb has changed. The step length is now an optional parameter provided by address. Old codes can be migrated as following:

double medium(struct pumas_context * context, struct pumas_state * state,
    struct pumas_medium ** medium_p)
{
    ...
    return step
}

must be changed to

void medium(struct pumas_context * context, struct pumas_state * state,
    struct pumas_medium ** medium_p, double * step_p)
{
    ...
    if (step_p != NULL) *step_p = step;
}

Note however that the new stepping algorithms allows better optimisations to be implemented, e.g. when using a ray tracer. The geometry callback gets called with 4 different cases, as following.

1. medium_p and step_p are both non NULL.
   This is an initialisation call. The stepping algorithm requires the start medium and an initial step length.

2. medium_p is NULL but step_p is not.
   This is a tracing call. The stepping algorithm requires a geometric step length.

3. medium_p is not NULL but step_p is NULL.
   This is a location call. It always occurs along the direction of the last preceding tracing call and not farther than the provided step length, within a few μm. Those call are used by PUMAS in order to do an approximate stepping, e.g. using a binary search. If an exact geometry is use, one can cache the result of the last tracing cal and re-use it at this step. Note that due to Physics interactions the used stepping distance might be lower than the one proposed at the last tracing call.

4. medium_p and step_p are NULL.
   This is a compatibility case for non uniform media. It ensures that the medium callback is always called before querying the local properties of the medium. If you don't need this feature, it is safe to add the following at the begining of your medium callback:

   if ((medium_p == NULL) && (step_p == NULL)) return;

Bug patches

• Non uniform magnetic fields should be properly handled now.

Minor features

• Versioning was added to the CMake build. The PUMAS library should now be generated with its version number, e.g. libpumas.so.0.14 on Linux.

PUMAS v0.13

This new release adds unit tests and correct some bugs, mostly in the API. Note that the signature of the pumas_recorder_create function has been changed as well, in order to match the same convention than for other objects.

List of changes

• Patch several bugs in the API mostly related to some corner cases not being triggered properly. Standard usage should not be affected though.
• Patch some error messages not being triggered properly, and polish the error reports.
• Update the doc to match the current behaviour.
• Patch a bug for direction vectors loosing unitarity for strong magnetic fields.
• Add unit tests with coverage reports using gcov & lcov.
• Expand the CMake build with examples and tests. Note that those are not built by default but must be explicitly enabled, e.g. as cmake -DPUMAS_BUILD_EXAMPLES.
• Changed the signature of the pumas_recorder_create function. Now the recorder is the 1st argument, as for other
  API functions manipulating their base object.

PUMAS v0.12-p1

This is a patch for v0.12. It corrects some bugs related to the transport:

• PUMAS_EVENT_MEDIUM should now interrupt pumas_transport as expected.
• Particles should now be transported again in uniform & infinite media.
• Straight steps have been optimised when navigating through a geometry.

In addition, the signature of the pumas_context_create function has been changed to match the rest of the API, i.e. the context goes as 1st argument.

PUMAS v0.12

This new release adds several improvements for controlling and monitoring the transport. Extra verbosity has been added to error reports as well. Note that some changes are not backward compatible with v0.11. In particular the API signature of pumas_initialise, pumas_transport, pumas_recorder_create and pumas_error_handler_cb have changed. In addition the pumas_error_print and pumas_error_string functions have been removed.

List of changes

1. Optimize the stepping and patch some bugs.

2. Add more verbosity to error reports. When using an error handler a brief description of the error is provided, including some contextual information, e.g. file and line number.

3. Energy loss path is now absolute by default. In order to specify a relative path, w.r.t. the MDF location, as previously, one must now prefix the energy loss path with a @. E.g. dedx_path = "@dedx".

4. The pumas_transport function can now return an event flag indicating its stop condition. In addition it can also return the start and end media of the transport. Furthermore, stop conditions can be configured by modifying the new event field of the simulation context.

5. The pumas_recorder API has been cleaned. In addition it also accepts a user callback now for custom recording. If no callback is provided the previous in memory recorder is used. Extra PUMAS_EVENT_START and PUMAS_EVENT_STOP events are provided to the recorder callback in order to notify the start or stop of a pumas_transport.

6. The primary muon flux in the examples has been changed for a more realistic model, i.e. Guan et al..

PUMAS v0.11

• Add some examples and expand the documentation.

• Initialise the materials locals to zero when entering a new medium.

  • No need to overwrite the magnetic field anymore, when null.

• Patch a bug with the maximum distance in backward mode note being properly handled.

  • Introduced by a previous commit in v0.10.

PUMAS v0.10

• Patch several bugs:

  • Protect against a numeric error in the stepping, for tiny density changes.
  • Patch a critical bug where the kinetic limit was not properly updated at a change of medium.
  • Patch a bug in the deallocation of dedx file names.
  • Patch the bias weight for using x = q / K fraction instead of nu = q / E

• Improve the DEL algorithms at low energy for ionisation losses:

  • Add an analytical form for the integral of the ionisation DCS, and a direct (unbiased) sampling at low energies for forward MC.
  • Tune the backward sampling at low energies.

PUMAS v0.9

• Optimise the tracking for the case where exact steps to the next boundary are provided.
• Correct an underflow in the Pair Creation DCS, at UHE, e.g. for τ.
• Tune the detailed stepping at UHE.