Celeritas#

The celeritas directory focuses on the physics and transport loop implementation for the Celeritas codebase, using components from the corecel and orange dependencies.

Fundamentals#

namespace units#

Units in Celeritas for macro-scale quantities.

The following units have numerical values of 1 in the default Celeritas system (Gaussian CGS):

cm for standard unit of length
s for standard unit of time
g for standard unit of mass
G for standard unit of field strength

Unless otherwise specified, this unit system is used for input and output numerical values. They are meant for macro-scale quantities coupling the different code components of Celeritas.

Celeritas	CLHEP	Notes
a0_bohr	Bohr_radius	Bohr radius
alpha_fine_structure	fine_structure_const
atomic_mass	amu	Not the same as 1/avogadro
eps_electric	epsilon0	Vacuum permittivity
h_planck	h_Planck
k_boltzmann	k_Boltzmann
mu_magnetic	mu0	Vacuum permeability
na_avogadro	Avogadro	[1/mol]
r_electron	classic_electr_radius	Classical electron radius
kcd_luminous	[none]	Lumens per Watt
lambdabar_electron	electron_Compton_length	Reduced Compton wavelength
stable_decay_constant	[none]	Decay for a stable particle

Problem definition#

class MaterialParams : public celeritas::ParamsDataInterface<MaterialParamsData>#

Data management for material, element, and nuclide properties.

Material metadata

inline MaterialId::size_type num_materials() const#: Get the label of a material.

Label const &id_to_label(MaterialId id) const#: Get the label of a material.

MaterialId find_material(std::string const &name) const#

Locate the material ID corresponding to a label.

If the label isn’t among the materials, a null ID will be returned.

SpanConstMaterialId find_materials(std::string const &name) const#

Get zero or more material IDs corresponding to a name.

This is useful for materials that are repeated with different uniquifying ‘extensions’.

Element metadata

inline ElementId::size_type num_elements() const#: Get the label of a element.

Label const &id_to_label(ElementId id) const#: Get the label of a element.

ElementId find_element(std::string const &name) const#

Locate the element ID corresponding to a label.

If the label isn’t among the elements, a null ID will be returned.

SpanConstElementId find_elements(std::string const &name) const#: Get zero or more element IDs corresponding to a name.

Isotope metadata

Number of distinct isotope definitions

inline IsotopeId::size_type num_isotopes() const#: Get the label of an isotope.

Label const &id_to_label(IsotopeId id) const#: Get the label of an isotope.

IsotopeId find_isotope(std::string const &name) const#

Locate the isotope ID corresponding to a label.

If the label isn’t among the isotopes, a null ID will be returned.

SpanConstIsotopeId find_isotopes(std::string const &name) const#: Get zero or more isotope IDs corresponding to a name.

Public Functions

explicit MaterialParams(Input const &inp)#: Construct from a vector of material definitions.

inline MaterialId::size_type size() const#: Number of material definitions.

inline MaterialView get(MaterialId id) const#: Get material properties for the given material.

inline ElementView get(ElementId id) const#: Get properties for the given element.

inline IsotopeView get(IsotopeId id) const#: Get properties for the given isotope.

inline IsotopeComponentId::size_type max_isotope_components() const#: Maximum number of isotopes in any one element.

inline ElementComponentId::size_type max_element_components() const#: Maximum number of elements in any one material.

inline bool is_missing_isotopes() const#: Whether isotopic data is present (may be false for EM-only physics)

inline virtual HostRef const &host_ref() const final#: Access material properties on the host.

inline virtual DeviceRef const &device_ref() const final#: Access material properties on the device.

Public Static Functions

static std::shared_ptr<MaterialParams> from_import(ImportData const &data)#: Construct with imported data.

struct ElementInput#

Define an element’s input data.

Public Members

AtomicNumber atomic_number#: Atomic number Z.

units::AmuMass atomic_mass#: Isotope-weighted average atomic mass.

std::vector<std::pair<IsotopeId, real_type>> isotopes_fractions#: Isotopic fractional abundance.

Label label#: Element name.

struct Input#: Input data to construct this class.

struct IsotopeInput#

Define an element’s isotope input data.

Public Members

AtomicNumber atomic_number#: Atomic number Z.

AtomicMassNumber atomic_mass_number#: Atomic number A.

units::MevMass nuclear_mass#: Nucleons’ mass + binding energy.

Label label#: Isotope name.

struct MaterialInput#

Define a material’s input data.

Public Members

real_type number_density#: Atomic number density [1/length^3].

real_type temperature#: Temperature [K].

MatterState matter_state#: Solid, liquid, gas.

std::vector<std::pair<ElementId, real_type>> elements_fractions#: Fraction of number density.

Label label#: Material name.

class ParticleParams : public celeritas::ParamsDataInterface<ParticleParamsData>#

Data management for Standard Model particle classifications.

This class represents “per-problem” shared data about standard model particles being used.

The ParticleParams is constructed on the host with a vector that combines metadata (used for debugging output and interfacing with physics setup) and data (used for on-device transport). Each entry in the construction is assigned a unique ParticleId used for runtime access.

The PDG Monte Carlo number is a unique “standard model” identifier for a particle. See “Monte Carlo Particle Numbering Scheme” in the “Review of

Particle Physics”:

https://pdg.lbl.gov/2020/reviews/rpp2020-rev-monte-carlo-numbering.pdf It should be used to identify particle types during construction time.

Public Types

using Input = std::vector<ParticleInput>#: Input data to construct this class.

Public Functions

explicit ParticleParams(Input const &defs)#: Construct with a vector of particle definitions.

inline ParticleId::size_type size() const#: Number of particle definitions.

inline std::string const &id_to_label(ParticleId id) const#: Get particle name.

inline PDGNumber id_to_pdg(ParticleId id) const#: Get PDG code for a particle ID.

inline ParticleId find(std::string const &name) const#: Find the ID from a name.

inline ParticleId find(PDGNumber pdg_code) const#: Find the ID from a PDG code.

ParticleView get(ParticleId id) const#: Get particle properties in host code.

inline virtual HostRef const &host_ref() const final#: Access material properties on the host.

inline virtual DeviceRef const &device_ref() const final#: Access material properties on the device.

Public Static Functions

static std::shared_ptr<ParticleParams> from_import(ImportData const &data)#: Construct with imported data.

struct ParticleInput#

Define a particle’s input data.

Public Members

std::string name#: Particle name.

PDGNumber pdg_code#: See “Review of Particle Physics”.

units::MevMass mass#: Rest mass [MeV / c^2].

units::ElementaryCharge charge#: Charge in units of [e].

real_type decay_constant#: Decay constant [1/s].

class PhysicsParams : public celeritas::ParamsDataInterface<PhysicsParamsData>#

Manage physics processes and models.

The physics params takes a vector of processes and sets up the processes and models. It constructs data and mappings of data:

particle type and process to tabulated values of cross sections etc,
particle type to applicable processes

During construction it constructs models and their corresponding list of ActionId values, as well as the tables of cross section data. Besides the individual interaction kernels, the physics parameters manage additional actions:

”pre-step”: calculate physics step limits
”along-step”: propagate, apply energy loss, multiple scatter
”range”: limit step by energy loss
”discrete-select”: sample a process for a discrete interaction, or reject due to integral cross sectionl
”integral-rejected”: do not apply a discrete interaction
”failure”: model failed to allocate secondaries

Public Functions

explicit PhysicsParams(Input)#: Construct with processes and helper classes.

inline ModelId::size_type num_models() const#: Number of models.

inline ProcessId::size_type num_processes() const#: Number of processes.

inline ParticleId::size_type num_particles() const#: Number of particle types.

inline ProcessId::size_type max_particle_processes() const#: Number of particle types.

inline SPConstModel const &model(ModelId) const#: Get a model.

inline SPConstProcess const &process(ProcessId) const#: Get a process.

inline ProcessId process_id(ModelId id) const#: Get the process ID of the given model.

inline ActionIdRange model_actions() const#: Get the action kernel IDs for all models.

SpanConstProcessId processes(ParticleId) const#: Get the list of process IDs that apply to a particle type.

inline virtual HostRef const &host_ref() const final#: Access physics properties on the host.

inline virtual DeviceRef const &device_ref() const final#: Access physics properties on the device.

struct Input#

Physics parameter construction arguments.

Public Members

SPConstRelaxation relaxation#: Optional atomic relaxation.

Transport interface#

template<MemSpace M> class Stepper : public celeritas::StepperInterface#

Manage a state vector and execute a single step on all of them.

Stepper<MemSpace::host> step(input);

// Transport primaries for the initial step
StepperResult alive_tracks = step(my_primaries);
while (alive_tracks)
{
    // Transport secondaries
    alive_tracks = step();
}

Public Functions

explicit Stepper(Input input)#: Construct with problem parameters and setup options.

~Stepper()#: Default destructor.

virtual StepperResult operator()() final#

Transport already-initialized states.

A single transport step is simply a loop over a toplogically sorted DAG of kernels.

virtual StepperResult operator()(SpanConstPrimary primaries) final#: Initialize new primaries and transport them for a single step.

virtual void reseed(EventId event_id) final#

Reseed the RNGs at the start of an event for “strong” reproducibility.

This reinitializes the RNG states using a single seed and unique subsequence for each thread. It ensures that given an event number, that event can be reproduced.

inline virtual ActionSequence const &actions() const final#: Get action sequence for timing diagnostics.

inline StateRef const &state_ref() const#: Access core data, primarily for debugging.

inline virtual CoreStateInterface const &state() const final#: Get the core state interface for diagnostic output.

Geant4 physics interfaces#

class GeantImporter#

Load problem data directly from Geant4.

This can be used to circumvent ROOT as a serialization tool, whether to simplify the toolchain or to integrate better with Acceleritas.

GeantImporter import(GeantSetup("blah.gdml"));
ImportData data = import();

or to import from an existing, initialized Geant4 state:

GeantImport import(world_volume);
ImportData data = import();

Public Functions

inline explicit GeantImporter(G4VPhysicalVolume const *world)#: Construct from an existing Geant4 geometry, assuming physics is loaded.

inline explicit GeantImporter(GeantSetup &&setup)#: Construct by capturing a GeantSetup object.

inline ImportData operator()(DataSelection const &selection)#: Load data from Geant4.

inline ImportData operator()()#: Fill all available data from Geant4.

Public Static Functions

static inline G4VPhysicalVolume const *get_world_volume()#

Get an externally loaded Geant4 top-level geometry element.

This is only defined if Geant4 has already been set up. It’s meant to be used in concert with GeantImporter or other Geant-importing classes.

class GeantSetup#

Construct a Geant 4 run manager and populate internal Geant4 physics.

This is usually passed directly into GeantImporter . It hides Geant4 implementation details (including header files) from the rest of the code. It is safe to include even when Geant4 is unavailable!

The setup is targeted specifically for physics that Celeritas supports.

Type aliases

using Options = GeantPhysicsOptions #: Prevent copying but allow moving

GeantSetup(GeantSetup const&) = delete#: Prevent copying but allow moving

GeantSetup &operator=(GeantSetup const&) = delete#: Prevent copying but allow moving

GeantSetup(GeantSetup&&) = default#: Prevent copying but allow moving

GeantSetup &operator=(GeantSetup&&) = default#: Prevent copying but allow moving

Public Functions

inline GeantSetup(std::string const &gdml_filename, Options options)#: Construct from a GDML file and physics options.

inline ~GeantSetup()#: Terminate the run manager on destruction.

inline G4VPhysicalVolume const *world() const#: Get the world detector volume.

inline explicit operator bool() const#: True if we own a run manager.

Geant4 physics options#

struct GeantPhysicsOptions#

Construction options for Geant physics.

These options attempt to default to our closest match to G4StandardEmPhysics.

Gamma physics

Enable discrete Coulomb

bool coulomb_scattering = {false}#: Enable Compton scattering.

bool compton_scattering = {true}#: Enable Compton scattering.

bool photoelectric = {true}#: Enable the photoelectric effect.

bool rayleigh_scattering = {true}#: Enable Rayleigh scattering.

bool gamma_conversion = {true}#: Enable electron pair production.

bool gamma_general = {false}#: Use G4GammaGeneral instead of individual gamma processes.

Electron and positron physics

Enable e- and e+ ionization

bool ionization = {true}#: Enable positron annihilation.

bool annihilation = {true}#: Enable positron annihilation.

BremsModelSelection brems = {BremsModelSelection::all}#: Enable bremsstrahlung and select a model.

MscModelSelection msc = {MscModelSelection::urban_extended}#: Enable multiple coulomb scattering and select a model.

RelaxationSelection relaxation = {RelaxationSelection::none}#: Enable atomic relaxation and select a model.

Physics options

Number of log-spaced bins per factor of 10 in energy

int em_bins_per_decade = {7}#: Enable universal energy fluctuations.

bool eloss_fluctuation = {true}#: Enable universal energy fluctuations.

bool lpm = {true}#: Apply relativistic corrections for select models.

bool integral_approach = {true}#: See PhysicsParamsOptions::disable_integral_xs.

Cutoff options

Lowest energy of any EM physics process

MevEnergy min_energy = {0.1 * 1e-3}#: Highest energy of any EM physics process.

MevEnergy max_energy = {100 * 1e6}#: Highest energy of any EM physics process.

double linear_loss_limit = {0.01}#: See PhysicsParamsOptions::linear_loss_limit.

MevEnergy lowest_electron_energy = {0.001}#: Tracking cutoff kinetic energy for e-/e+.

bool apply_cuts = {false}#: Kill secondaries below the production cut.

double default_cutoff = {0.1}#: Set the default production cut for all particle types [len].

Multiple scattering configuration

E-/e+ range factor for MSC models

double msc_range_factor = {0.04}#: Safety factor for MSC models.

double msc_safety_factor = {0.6}#: Safety factor for MSC models.

double msc_lambda_limit = {0.1 * units::centimeter}#: Lambda limit for MSC models [len].

Public Members

bool verbose = {false}#: Print detailed Geant4 output.

Geometry interfaces#

class VecgeomParams : public celeritas::GeoParamsInterface, public celeritas::ParamsDataInterface<VecgeomParamsData>#

Shared model parameters for a VecGeom geometry.

The model defines the shapes, volumes, etc.

Public Functions

explicit VecgeomParams(std::string const &gdml_filename)#: Construct from a GDML input.

explicit VecgeomParams(G4VPhysicalVolume const *world)#: Translate a geometry from Geant4.

~VecgeomParams()#: Clean up vecgeom on destruction.

inline virtual bool supports_safety() const final#: Whether safety distance calculations are accurate and precise.

inline virtual BBox const &bbox() const final#: Outer bounding box of geometry.

inline int max_depth() const#: Maximum nested geometry depth.

inline virtual VolumeId::size_type num_volumes() const final#: Number of volumes.

virtual Label const &id_to_label(VolumeId vol_id) const final#: Get the label for a placed volume ID.

virtual VolumeId find_volume(std::string const &name) const final#: Get the ID corresponding to a label.

virtual VolumeId find_volume(Label const &label) const final#

Locate the volume ID corresponding to a label.

If the label isn’t in the geometry, a null ID will be returned.

virtual VolumeId find_volume(G4LogicalVolume const *volume) const final#: Locate the volume ID corresponding to a Geant4 logical volume.

virtual SpanConstVolumeId find_volumes(std::string const &name) const final#

Get zero or more volume IDs corresponding to a name.

This is useful for volumes that are repeated in the geometry with different uniquifying ‘extensions’ from Geant4.

inline virtual HostRef const &host_ref() const final#: Access geometry data on host.

inline virtual DeviceRef const &device_ref() const final#: Access geometry data on device.

Public Static Functions

static bool use_surface_tracking()#: Whether surface tracking is being used.

static bool use_vgdml()#: Whether VecGeom GDML is used to load the geometry.

class GeantGeoParams : public celeritas::GeoParamsInterface, public celeritas::ParamsDataInterface<GeantGeoParamsData>#

Shared Geant4 geometry model wrapper.

This can be consructed directly by loading a GDML file, or in-memory using an existing physical volume. One “gotcha” is that due to persistent static variables in Geant4, the volume IDs will be offset if a geometry has been loaded and closed previously.

Public Functions

explicit GeantGeoParams(std::string const &gdml_filename)#

Construct from a GDML input.

This assumes that Celeritas is driving and will manage Geant4 exceptions etc.

explicit GeantGeoParams(G4VPhysicalVolume const *world)#: Use an existing loaded Geant4 geometry.

~GeantGeoParams()#: Clean up on destruction.

inline G4VPhysicalVolume const *world() const#: Access the world volume.

inline virtual bool supports_safety() const final#: Whether safety distance calculations are accurate and precise.

inline virtual BBox const &bbox() const final#: Outer bounding box of geometry.

inline virtual VolumeId::size_type num_volumes() const final#: Number of volumes.

virtual Label const &id_to_label(VolumeId vol_id) const final#: Get the label for a placed volume ID.

virtual VolumeId find_volume(std::string const &name) const final#: Get the ID corresponding to a label.

virtual VolumeId find_volume(Label const &label) const final#

Locate the volume ID corresponding to a label.

If the label isn’t in the geometry, a null ID will be returned.

virtual VolumeId find_volume(G4LogicalVolume const *volume) const final#: Locate the volume ID corresponding to a Geant4 logical volume.

virtual SpanConstVolumeId find_volumes(std::string const &name) const final#

Get zero or more volume IDs corresponding to a name.

This is useful for volumes that are repeated in the geometry with different uniquifying ‘extensions’ from Geant4.

G4LogicalVolume const *id_to_lv(VolumeId vol_id) const#

Get the Geant4 logical volume corresponding to a volume ID.

If the input volume ID is false, a null pointer will be returned.

inline virtual HostRef const &host_ref() const final#: Access geometry data on host.

inline virtual DeviceRef const &device_ref() const final#: No GPU support code.

On-device access#

class MaterialTrackView#

Read/write view to the material properties of a single particle track.

These functions should be used in each physics Process or Interactor or anything else that needs to access particle properties. Assume that all these functions are expensive: when using them as accessors, locally store the results rather than calling the function repeatedly. If any of the calculations prove to be hot spots we will experiment with cacheing some of the variables.

The element scratch space is “thread-private” data with a fixed size greater than or equal to the number of elemental components in the current material.

Public Functions

inline MaterialTrackView(MaterialParamsRef const &params, MaterialStateRef const &states, TrackSlotId tid)#: Construct from dynamic and static particle properties.

inline MaterialTrackView &operator=(Initializer_t const &other)#: Initialize the particle.

MaterialId material_id() const#: Current material identifier.

MaterialView make_material_view() const#: Get material properties for the current material.

inline Span<real_type> element_scratch()#: Access scratch space with at least one real per element component.

class ParticleTrackView#

Read/write view to the physical properties of a single particle track.

Public Functions

inline ParticleTrackView(ParticleParamsRef const &params, ParticleStateRef const &states, TrackSlotId id)#: Construct from dynamic and static particle properties.

inline ParticleTrackView &operator=(Initializer_t const &other)#: Initialize the particle.

inline void energy(Energy)#

Change the particle’s kinetic energy.

This should only be used when the particle is in a valid state. For HEP applications, the new energy should always be less than the starting energy.

inline void subtract_energy(Energy)#: Reduce the particle’s energy without undergoing a collision [MeV].

ParticleId particle_id() const#: Unique particle type identifier.

Energy energy() const#: Kinetic energy [MeV].

bool is_stopped() const#: Whether the track is stopped (zero kinetic energy).

ParticleView particle_view() const#: Get static particle properties for the current state.

units::MevMass mass() const#: Rest mass [MeV / c^2].

units::ElementaryCharge charge() const#: Elementary charge.

real_type decay_constant() const#: Decay constant.

bool is_antiparticle() const#: Whether it is an antiparticle.

bool is_stable() const#: Whether particle is stable.

inline real_type beta_sq() const#

Square of \( \beta \), which is the fraction of lightspeed [unitless].

Beta is calculated using the equality

\[ pc/E = \beta ; \quad \beta^2 = \frac{p^2 c^2}{E^2} \]

. Using

\[ E^2 = p^2 c^2 + m^2 c^4 \]

and

\[ E = K + mc^2 \]

the square of the fraction of lightspeed speed is

\[ \beta^2 = 1 - \left( \frac{mc^2}{K + mc^2} \right)^2 = 1 - \gamma^{-2} \]

where \( \gamma \) is the Lorentz factor (see below).

By choosing not to divide out the mass, this expression will work for massless particles.

Todo:: For nonrelativistic particles (low kinetic energy) this expression may be inaccurate due to rounding error. It is however guaranteed to be in the exact range [0, 1].

inline units::LightSpeed speed() const#

Speed [1/c].

Speed is calculated using beta so that the expression works for massless particles.

inline real_type lorentz_factor() const#

Lorentz factor [unitless].

The Lorentz factor can be viewed as a transformation from classical quantities to relativistic quantities. It’s defined as

\[ \gamma = \frac{1}{\sqrt{1 - v^2 / c^2}} \]

Its value is infinite for the massless photon, and greater than or equal to 1 otherwise.

Gamma can also be calculated from the total (rest + kinetic) energy

\[ E = \gamma mc^2 = K + mc^2 \]

which we ues here since K and m are the primary stored quantities of the particles:

\[ \gamma = 1 + \frac{K}{mc^2} \]

inline units::MevMomentum momentum() const#

Relativistic momentum [MeV / c].

This is calculated by taking the root of the square of the momentum.

inline units::MevMomentumSq momentum_sq() const#

Square of relativistic momentum [MeV^2 / c^2].

Total energy:

\[ E = K + mc^2 \]

Relation between energy and momentum:

\[ E^2 = p^2 c^2 + m^2 c^4 \]

therefore

\[ p^2 = \frac{E^2}{c^2} - m^2 c^2 \]

\[ p^2 = \frac{K^2}{c^2} + 2 * m * K \]

class PhysicsTrackView#

Physics data for a track.

The physics track view provides an interface for data and operations common to most processes and models.

Public Functions

inline PhysicsTrackView(PhysicsParamsRef const &params, PhysicsStateRef const &states, ParticleId particle, MaterialId material, TrackSlotId tid)#

Construct from shared and state data.

Particle and material IDs are derived from other class states.

inline PhysicsTrackView &operator=(Initializer_t const&)#: Initialize the track view.

inline void interaction_mfp(real_type)#

Set the distance to the next interaction, in mean free paths.

This value will be decremented at each step.

inline void reset_interaction_mfp()#

Set the distance to the next interaction, in mean free paths.

This value will be decremented at each step.

inline void dedx_range(real_type)#

Set the energy loss range for the current material and particle energy.

This value will be calculated once at the beginning of each step.

inline void msc_range(MscRange const&)#

Set the range properties for multiple scattering.

These values will be calculated at the first step in every tracking volume.

bool has_interaction_mfp() const#: Whether the remaining MFP has been calculated.

real_type interaction_mfp() const#: Remaining MFP to interaction.

real_type dedx_range() const#: Energy loss range.

MscRange const &msc_range() const#: Persistent range properties for multiple scattering within a same volume.

MaterialId material_id() const#: Current material identifier.

inline ParticleProcessId::size_type num_particle_processes() const#: Number of processes that apply to this track.

inline ProcessId process(ParticleProcessId) const#: Process ID for the given within-particle process index.

inline ValueGridId value_grid(ValueGridType table, ParticleProcessId) const#

Return value grid data for the given table type and process if available.

If the result is not null, it can be used to instantiate a grid Calculator.

If the result is null, it’s likely because the process doesn’t have the associated value (e.g. if the table type is “energy_loss” and the process is not a slowing-down process).

inline IntegralXsProcess const &integral_xs_process(ParticleProcessId ppid) const#

Get data for processes that use the integral approach.

Particles that have energy loss processes will have a different energy at the pre- and post-step points. This means the assumption that the cross section is constant along the step is no longer valid. Instead, Monte Carlo integration can be used to sample the interaction length for the discrete process with the correct probability from the exact distribution,

\[ p = 1 - \exp \left( -\int_{E_0}^{E_1} n \sigma(E) \dif s \right), \]

where \( E_0 \) is the pre-step energy, \( E_1 \) is the post-step energy, n is the atom density, and s is the interaction length.

At the start of the step, the maximum value of the cross section over the step \( \sigma_{max} \) is estimated and used as the macroscopic cross section for the process rather than \( \sigma_{E_0} \). After the step, the new value of the cross section \( \sigma(E_1) \) is calculated, and the discrete interaction for the process occurs with probability

\[ p = \frac{\sigma(E_1)}{\sigma_{\max}}. \]

See section 7.4 of the Geant4 Physics Reference (release 10.6) for details.

inline real_type calc_xs(ParticleProcessId ppid, MaterialView const &material, Energy energy) const#: Calculate macroscopic cross section for the process.

inline real_type calc_max_xs(IntegralXsProcess const &process, ParticleProcessId ppid, MaterialView const &material, Energy energy) const#

Estimate maximum macroscopic cross section for the process over the step.

If this is a particle with an energy loss process, this returns the estimate of the maximum cross section over the step. If the energy of the global maximum of the cross section (calculated at initialization) is in the interval \( [\xi E_0, E_0) \), where \( E_0 \) is the pre-step energy and \( \xi \) is min_eprime_over_e (defined by default as \( \xi = 1 - \alpha \), where \( \alpha \) is max_step_over_range), \( \sigma_{\max} \) is set to the global maximum. Otherwise, \( \sigma_{\max} = \max( \sigma(E_0), \sigma(\xi E_0) ) \). If the cross section is not monotonic in the interval \( [\xi E_0, E_0) \) and the interval does not contain the global maximum, the post-step cross section \( \sigma(E_1) \) may be larger than \( \sigma_{\max} \).

inline ModelFinder make_model_finder(ParticleProcessId) const#: Models that apply to the given process ID.

inline ValueTableId value_table(ParticleModelId) const#

Return value table data for the given particle/model/material.

A null result means either the model is material independent or the material only has one element, so no cross section CDF tables are stored.

inline TabulatedElementSelector make_element_selector(ValueTableId, Energy) const#: Construct an element selector to sample an element from tabulated xs data.

inline bool has_at_rest() const#: Whether the particle can have a discrete interaction at rest.

inline ModelId action_to_model(ActionId) const#: Convert an action to a model ID for diagnostics, false if not a model.

inline ActionId model_to_action(ModelId) const#: Convert a selected model ID into a simulation action ID.

inline ModelId model_id(ParticleModelId) const#: Get the model ID corresponding to the given ParticleModelId.

inline real_type range_to_step(real_type range) const#

Calculate scaled step range.

This is the updated step function given by Eq. 7.4 of Geant4 Physics Reference Manual, Release 10.6:

\[ s = \alpha r + \rho (1 - \alpha) (2 - \frac{\rho}{r}) \]

where alpha is max_step_over_range and rho is min_range .

Below min_range, no step scaling is applied, but the step can still be arbitrarily small.

PhysicsParamsScalars const &scalars() const#: Access scalar properties (options, IDs).

inline size_type num_particles() const#: Number of particle types.

template<class T> inline T make_calculator(ValueGridId) const#

Construct a grid calculator of the given type.

The calculator must take two arguments: a reference to XsGridRef, and a reference to the Values data structure.

inline ModelId hardwired_model(ParticleProcessId ppid, Energy energy) const#

Return the model ID that applies to the given process ID and energy if the process is hardwired to calculate macroscopic cross sections on the fly.

If the result is null, tables should be used for this process/energy.

inline ParticleProcessId eloss_ppid() const#: Particle-process ID of the process with the de/dx and range tables.

Propagation and magnetic field#

The propagation interface is built on top of the geometry to allow both curved and straight-line movement. Field propagation is based on a composition of:

Field: Maps a point in space and time to a field vector.
Equation of motion: Calculates the path derivative of position and momentum given their current state and the templated field.
Integrator: Numerically integrates a new position/momentum state given the start, path derivative, and step length.
Driver: Integrate path segments that satisfy certain error conditions, solving for the required segment length.
Propagator: Given a maximum physics step, advance the geometry state and momentum along the field lines, satisfying constraints (see field driver options) for the maximum geometry error.

Propagation#

template<class GTV> class LinearPropagator#

Propagate (move) a particle in a straight line.

Public Functions

inline LinearPropagator(GTV &&track)#: Construct from a geo track view.

inline result_type operator()()#: Move track to next volume boundary.

inline result_type operator()(real_type dist)#: Move track by a user-provided distance up to the next boundary.

Public Static Functions

static inline bool tracks_can_loop()#: Whether it’s possible to have tracks that are looping.

template<class DriverT, class GTV> class FieldPropagator#

Propagate a charged particle in a field.

For a given initial state (position, momentum), it propagates a charged particle along a curved trajectory up to an interaction length proposed by a chosen physics process for the step, possibly integrating sub-steps by an adaptive step control with a required accuracy of tracking in a field. It updates the final state (position, momentum, boundary) along with the step actually taken. If the final position is outside the current volume, it returns a geometry limited step and the state at the intersection between the curve trajectory and the first volume boundary using an iterative step control method within a tolerance error imposed on the closest distance between two positions by the field stepper and the linear projection to the volume boundary.

Note

This follows similar methods as in Geant4’s G4PropagatorInField class.

Public Functions

inline FieldPropagator(DriverT &&driver, ParticleTrackView const &particle, GTV &&geo)#: Construct with shared field parameters and the field driver.

inline result_type operator()()#: Propagate a charged particle until it hits a boundary.

inline result_type operator()(real_type dist)#

Propagate a charged particle in a field.

It utilises a field driver (based on an adaptive step control to limit the length traveled based on the magnetic field behavior and geometric tolerances) to track a charged particle along a curved trajectory for a given step length within a required accuracy or intersects with a new volume (geometry limited step).

The position of the internal OdeState state_ should be consistent with the geometry geo_’s position, but the geometry’s direction will be a series of “trial” directions that are the chords between the start and end points of a curved substep through the field. At the end of the propagation step, the geometry state’s direction is updated based on the actual value of the calculated momentum.

Caveats:

The physical (geometry track state) position may deviate from the exact curved propagation position up to a driver-based tolerance at every boundary crossing. The momentum will always be conserved, though.
In some unusual cases (e.g. a very small caller-requested step, or an unusual accumulation in the driver’s substeps) the distance returned may be slightly higher (again, up to a driver-based tolerance) than the physical distance travelled.

inline real_type delta_intersection() const#

Distance close enough to a boundary to mark as being on the boundary.

TODO: change delta intersection from property in FieldDriverOptions to another FieldPropagatorOptions

inline real_type bump_distance() const#: Distance to bump or to consider a “zero” movement.

inline real_type minimum_substep() const#: Distance to bump or to consider a “zero” movement.

Public Static Functions

static inline bool tracks_can_loop()#: Whether it’s possible to have tracks that are looping.

static inline short int max_substeps()#: Limit on substeps.

template<template<class EquationT> class StepperT, class FieldT, class GTV> decltype(auto) celeritas::make_mag_field_propagator(FieldT &&field, FieldDriverOptions const &options, ParticleTrackView const &particle, GTV &&geometry)#

Field data input and options#

struct FieldDriverOptions#

Configuration options for the field driver.

TODO: replace epsilon_rel_max with 1/epsilon_rel_max^2 TODO: replace safety with step_shrink_mul (or something to indicate that it’s a multiplicative factor for reducing the step, not anything with geometry) TODO: remove errcon

Public Functions

inline explicit operator bool() const#: Whether all data are assigned and valid.

Public Members

real_type minimum_step = 1.0e-5 * units::millimeter #: The minimum length of the field step.

real_type delta_chord = 0.25 * units::millimeter #: The maximum sagitta of each substep (“miss distance”)

real_type delta_intersection = 1.0e-4 * units::millimeter #: Accuracy of intersection of the boundary crossing.

real_type epsilon_step = 1.0e-5#: Discretization error tolerance for each field substep.

real_type epsilon_rel_max = 1.0e-3#: Targeted discretization error for “integrate step”.

real_type errcon = 1.0e-4#: UNUSED: Targeted discretization error for “one good step”.

real_type pgrow = -0.20#: Exponent to increase a step size.

real_type pshrink = -0.25#: Exponent to decrease a step size.

real_type safety = 0.9#: Scale factor for the predicted step size.

real_type max_stepping_increase = 5#: Largest allowable relative increase a step size.

real_type max_stepping_decrease = 0.1#: Smallest allowable relative decrease in step size.

short int max_nsteps = 100#: Maximum number of integrations (or trials)

Public Static Attributes

static constexpr real_type initial_step_tol = 1e-6#: Initial step tolerance.

static constexpr real_type dchord_tol = 1e-5 * units::millimeter #: Chord distance fudge factor.

static constexpr real_type min_chord_shrink = 0.5#: Lowest allowable scaling factor when searching for a chord.

Field data#

These classes correspond to JSON input files to the field setup.

struct UniformFieldParams#

Input data and options for a uniform field.

Public Members

Real3 field = {0, 0, 0}#: Field strength (native units)

struct RZMapFieldInput#

Input data for an magnetic R-Z vector field stored on an R-Z grid.

The magnetic field is discretized at nodes on an R-Z grid, and each point the field vector is approximated by a 2-D vector in R-Z. The input units of this field are in NATIVE UNITS (cm/gauss when CGS). An optional _units field in the input can specify whether the input is in tesla or native units, with allowable values of “T”, “tesla”, “gauss”, or “native”.

The field values are all indexed with R having stride 1: [Z][R]

Public Functions

inline explicit operator bool() const#: Whether all data are assigned and valid.

Public Members

double min_z = {}#: Lower z coordinate [len].

double min_r = {}#: Lower r coordinate [len].

double max_z = {}#: Last z coordinate [len].

double max_r = {}#: Last r coordinate [len].

std::vector<double> field_z#: Flattened Z field component [bfield].

std::vector<double> field_r#: Flattened R field component [bfield].

Random number distributions#

The 2011 ISO C++ standard defined a new functional paradigm for sampling from random number distributions. In this paradigm, random number engines generate a uniformly distributed stream of bits. Then, distributions use that entropy to sample a random number from a distribution. Distributions are function-like objects whose constructors take the parameters of the distribution: for example, a uniform distribution over the range \([a, b)\) takes the a and b parameters as constructor arguments. The templated call operator accepts a random engine as its sole argument.

Celeritas extends this paradigm to physics distributions. At a low level, it has random number distributions that result in single real values (such as uniform, exponential, gamma) and correlated three-vectors (such as sampling an isotropic direction).

class BernoulliDistribution#

Select one of two options with a given probability.

The constructor argument is the chance of returning true, with an optional second argument for normalizing with a fraction of false values.

BernoulliDistribution snake_eyes(1, 35);
if (snake_eyes(rng))
{
    // ...
}
BernoulliDistribution also_snake_eyes(1.0 / 36.0);

template<class T> class DeltaDistribution#: Distribution where the “sampled” value is just the given value.

template<class RealType = ::celeritas::real_type> class ExponentialDistribution#

Sample from an exponential distribution: -log(xi) / lambda.

This is simply an implementation of std::exponential_distribution that uses the Celeritas canonical generator and is independent of library implementation.

Note (for performance-critical sections of code) that if this class is constructed locally with the default value of lambda = 1.0, the inversion and multiplication will be optimized out (and the code will be exactly identical to -std::log(rng.ran()).

template<class RealType = ::celeritas::real_type> class GammaDistribution#

Sample from a gamma distribution.

The gamma distribution can be parameterized with a shape parameter \( \alpha \) and a scale parameter \( \beta \) and has the PDF:

\[ f(x; \alpha, \beta) = \frac{x^{\alpha - 1} e^{-x / \beta}}{\beta^\alpha \Gamma(\alpha)} \quad \mathrm{for}\ x > 0, \quad \alpha, \beta > 0 \]

The algorithm described in Marsaglia and Tsang [MT00] is used here to generate gamma-distributed random variables. The steps are:

Set \( d = \alpha - 1/3 \) and \( c = 1 / \sqrt{9d} \)
Generate random variates \( Z \sim N(0,1) \) and \( U \sim U(0,1) \)
Set \( v = (1 + cZ)^3 \)
If \( \log U < Z^2 / 2 + d(1 - v + \log v) \) return \( dv \). Otherwise, go to step 2. A squeeze function can be used to avoid the two logarithms in most cases by accepting early if \( U < 1 - 0.0331 Z^4 \).

Though this method is valid for \( \alpha \ge 1 \), it can easily be extended for \( \alpha < 1 \): if \( X \sim \Gamma(\alpha + 1) \) and \( U \sim U(0,1) \), then \( X U^{1/\alpha} \sim \Gamma(\alpha) \).

template<class RealType = ::celeritas::real_type> class IsotropicDistribution#: Sample points uniformly on the surface of a unit sphere.

template<class RealType = ::celeritas::real_type> class NormalDistribution#

Sample from a normal distribution.

This uses the Box-Muller transform to generate pairs of independent, normally distributed random numbers, returning them one at a time. Two random numbers uniformly distributed on [0, 1] are mapped to two independent, standard, normally distributed samples using the relations:

\[ x_1 = \sqrt{-2 \ln \xi_1} \cos(2 \pi \xi_2) x_2 = \sqrt{-2 \ln \xi_1} \sin(2 \pi \xi_2) \]

template<class RealType = ::celeritas::real_type> class PoissonDistribution#

Sample from a Poisson distribution.

The Poisson distribution describes the probability of \( k \) events occuring in a fixed interval given a mean rate of occurance \( \lambda \) and has the PMF:

\[ f(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}. \]

For small \( \lambda \), a direct method described in Knuth, Donald E., Seminumerical Algorithms, The Art of Computer Programming, Volume 2 can be used to generate samples from the Poisson distribution. Uniformly distributed random numbers are generated until the relation

\[ \prod_{k = 1}^n U_k \le e^{-\lambda} \]

is satisfied; then, the random variable \( X = n - 1 \) will have a Poisson distribution. On average this approach requires the generation of \( \lambda + 1 \) uniform random samples, so a different method should be used for large \( \lambda \).

Geant4 uses Knuth’s algorithm for \( \lambda \le 16 \) and a Gaussian approximation for \( \lambda > 16 \) (see G4Poisson), which is faster but less accurate than other methods. The same approach is used here.

template<class RealType = ::celeritas::real_type> class RadialDistribution#: Sample from a uniform radial distribution.

template<class RealType = ::celeritas::real_type> class ReciprocalDistribution#

Reciprocal or log-uniform distribution.

This distribution is defined on a positive range \( [a, b) \) and has the normalized PDF:

\[ f(x; a, b) = \frac{1}{x (\ln b - \ln a)} \quad \mathrm{for} a \le x < b \]

which integrated into a CDF and inverted gives a sample:

\[ x = a \left( \frac{b}{a} \right)^{\xi} = a \exp\!\left(\xi \log \frac{b}{a} \right) \]

template<class RealType = ::celeritas::real_type> class UniformBoxDistribution#: Sample a point uniformly in a box.

template<class RealType = ::celeritas::real_type> class UniformRealDistribution#

Sample from a uniform distribution.

This distribution is defined between two arbitrary real numbers a and b , and has a flat PDF between the two values. It is allowable for the two numbers to have reversed order.

Physics distributions#

At a higher level, Celeritas expresses many physics operations as distributions of updated track states based on original track states. For example, the Tsai-Urban distribution used for sampling exiting angles of bremsstrahlung and pair production has parameters of incident particle energy and mass, and it samples the exiting polar angle cosine.

class BhabhaEnergyDistribution#

Helper class for MollerBhabhaInteractor .

Sample the exiting energy for Bhabha scattering.

class EnergyLossGammaDistribution#

Sample energy loss from a gamma distribution.

This model is valid for heavy particles with small mean losses (less than 2 Bohr standard deviations). It is a special case of EnergyLossGaussianDistribution (see that class for more documentation).

Note that while this appears in G4UniversalFluctuation, the Geant4 documentation does not explain why the loss is sampled from a gamma distribution in this case.

class EnergyLossGaussianDistribution#

Sample energy loss from a Gaussian distribution.

In a thick absorber, the total energy transfer is a result of many small energy losses from a large number of collisions. The central limit theorem applies, and the energy loss fluctuations can be described by a Gaussian distribution. See section 7.3.1 of the Geant4 Physics Reference Manual and GEANT3 PHYS332 section 2.3.

The Gaussian approximation is valid for heavy particles and in the regime \( \kappa = \xi / T_\textrm{max} > 10 \). Fluctuations of the unrestricted energy loss follow a Gaussian distribution if \( \Delta E > \kappa T_{max} \), where \( T_{max} \) is the maximum energy transfer (PHYS332 section 2). For fluctuations of the restricted energy loss, the condition is modified to \( \Delta E > \kappa T_{c} \) and \( T_{max} \le 2 T_c \), where \( T_c \) is the delta ray cutoff energy (PRM Eq. 7.6-7.7).

class EnergyLossUrbanDistribution#

Sample from the Urban model of energy loss fluctuations in thin layers.

The Urban model is used to compute the energy loss fluctuation when \( \kappa \) is small. It assumes atoms have only two energy levels with binding energies \( E_1 \) and \( E_2 \) and that the particle-atom interaction will either be an excitation with energy loss \( E_1 \) or \( E_2 \) or an ionization with energy loss proportional to \( 1 / E^2 \). The number of collisions for each interaction type has a Poisson distribution with mean proportional to the interaction cross section. The energy loss in a step will be the sum of the energy loss contributions from excitation and ionization.

When the number of ionizations is larger than a given threshold, a fast sampling method can be used instead. The possible energy loss interval is divided into two parts: a lower range in which the number of collisions is large and the energy loss can be sampled from a Gaussian distribution, and an upper range in which the energy loss is sampled for each collision.

See section 7.3.2 of the Geant4 Physics Reference Manual and GEANT3 PHYS332 section 2.4 for details.

class MollerEnergyDistribution#

Helper class for MollerBhabhaInteractor .

Sample the exiting energy for Moller scattering.

class TsaiUrbanDistribution#

Polar angular distribution for pair-production and bremsstrahlung processes.

For pair-production, the polar angle of the electron (or positron) is defined with respect to the direction of the parent photon. The energy- angle distribution given by Tsai is quite complicated to sample and can be approximated by a density function suggested by Urban.

The angular distribution of the emitted photons is obtained from a simplified formula based on the Tsai cross-section, which is expected to become isotropic in the low energy limit.

Note

This performs the same sampling routine as in Geant4’s ModifiedTsai class, based on derivation from Tsai (Rev Mod Phys 49,421(1977) and documented in section 6.5.2 (pair-production), and 10.2.1 and 10.2.4 (bremsstrahlung) of the Geant4 Physics Reference (release 10.6).

Physics implementations#

Additional distributions are built on top of the helper distributions above. All discrete interactions (in Geant4 parlance, “post-step do-it”s) use distributions to sample an Interaction based on incident particle properties. The sampled result contains the updated particle direction and energy, as well as properties of any secondary particles produced.

class BetheHeitlerInteractor#

Relativistic model for electron-positron pair production.

The energies of the secondary electron and positron are sampled using the Bethe-Heitler cross sections with a Coulomb correction. The LPM effect is taken into account for incident gamma energies above 100 GeV.

Note

This performs the same sampling routine as in Geant4’s G4PairProductionRelModel, as documented in sections 6.5 (gamma conversion) and 10.2.2 (LPM effect) of the Geant4 Physics Reference Manual (release 10.7)

class EPlusGGInteractor#

Annihilate a positron to create two gammas.

This is a model for the discrete positron-electron annihilation process which simulates the in-flight annihilation of a positron with an atomic electron and produces into two photons. It is assumed that the atomic electron is initially free and at rest.

Note

This performs the same sampling routine as in Geant4’s G4eeToTwoGammaModel class, as documented in section 10.3 of the Geant4 Physics Reference (release 10.6). The maximum energy for the model is specified in Geant4.

class KleinNishinaInteractor#

Perform Compton scattering, neglecting atomic binding energy.

This is a model for the discrete Compton inelastic scattering process. Given an incident gamma, it adds a single secondary (electron) to the secondary stack and returns an interaction for the change to the incident gamma direction and energy. No cutoffs are performed for the incident energy or the exiting gamma energy. A secondary production cutoff is applied to the outgoing electron.

Note

This performs the same sampling routine as in Geant4’s G4KleinNishinaCompton, as documented in section 6.4.2 of the Geant4 Physics Reference (release 10.6).

class MollerBhabhaInteractor#

Perform Moller (e-e-) and Bhabha (e+e-) scattering.

This interaction, part of the ionization process, is when an incident electron or positron ejects an electron from the surrounding matter.

Note

This performs the same sampling routine as in Geant4’s G4MollerBhabhaModel class, as documented in section 10.1.4 of the Geant4 Physics Reference (release 10.6).

class LivermorePEInteractor#

Livermore model for the photoelectric effect.

A parameterization of the photoelectric cross sections in two different energy intervals, formulated as

\[ \sigma(E) = a_1 / E + a_2 / E^2 + a_3 / E^3 + a_4 / E^4 + a_5 / E^5 + a_6 / E^6 \: , \]

is used. The coefficients for this model are calculated by fitting the tabulated EPICS2014 subshell cross sections. The parameterized model applies above approximately 5 keV; below this limit (which depends on the atomic number) the tabulated cross sections are used. The angle of the emitted photoelectron is sampled from the Sauter-Gavrila distribution.

Note

This performs the same sampling routine as in Geant4’s G4LivermorePhotoElectricModel class, as documented in section 6.3.5 of the Geant4 Physics Reference (release 10.6).

class RayleighInteractor#: Apply the Livermore model of Rayleigh scattering to photons.

Note

This performs the same sampling routine as in Geant4’s G4LivermoreRayleighModel class, as documented in section 6.2.2 of the Geant4 Physics Reference (release 10.6).

class RelativisticBremInteractor#

Perform a high-energy Bremsstrahlung interaction.

This is a relativistic Bremsstrahlung model for high-energy (> 1 GeV) electrons and positrons.

Note

This performs the same sampling routine as in Geant4’s G4eBremsstrahlungRelModel class, as documented in section 10.2.2 of the Geant4 Physics Reference (release 10.7).

class SeltzerBergerInteractor#

Seltzer-Berger model for electron and positron bremsstrahlung processes.

Given an incoming electron or positron of sufficient energy (as per CutOffView), this class provides the energy loss of these particles due to radiation of photons in the field of a nucleus. This model improves accuracy using cross sections based on interpolation of published tables from Seltzer and Berger given in Nucl. Instr. and Meth. in Phys. Research B, 12(1):95–134 (1985) and Atomic Data and Nuclear Data Tables, 35():345 (1986). The cross sections are obtained from SBEnergyDistribution and are appropriately scaled in the case of positrons via SBPositronXsCorrector.

Note

This interactor performs an analogous sampling as in Geant4’s G4SeltzerBergerModel, documented in 10.2.1 of the Geant Physics Reference (release 10.6). The implementation is based on Geant4 10.4.3.

class AtomicRelaxation#

Simulate particle emission from atomic deexcitation.

The EADL radiative and non-radiative transition data is used to simulate the emission of fluorescence photons and (optionally) Auger electrons given an initial shell vacancy created by a primary process.

class EnergyLossHelper#

Precalculate energy loss fluctuation properties.

Fluctuations in the energy loss of charged particles over a given thickness of material arise from statistical variation in both the number of collisions and the energy lost in each collision. Above a given energy threshold, fluctuations are simulated through the explicit sampling of secondaries. However, the continuous energy loss below the cutoff energy also has fluctuations, and these are not taken into account in the calculation of the mean loss. For continuous energy loss, fluctuation models are used to sample the actual restricted energy loss given the mean loss.

Different models are used depending on the value of the parameter \( \kappa = \xi / T_{max} \), the ratio of the mean energy loss to the maximum allowed energy transfer in a single collision.

For large \( \kappa \), when the particle loses all or most of its energy along the step, the number of collisions is large and the straggling function can be approximated by a Gaussian distribution.

Otherwise, the Urban model for energy loss fluctuations in thin layers is used.

Note

This performs the same sampling routine as in Geant4’s G4UniversalFluctuation, as documented in section 7.3 in the Geant4 Physics Reference Manual and in PHYS332 and PHYS333 in GEANT3, CERN Program Library Long Writeup, W5013 (1993).

class UrbanMscStepLimit#

Sample a step limit for the Urban MSC model.

This distribution is to be used for tracks that have non-negligible steps and are near the boundary. Otherwise, no displacement or step limiting is needed.

Note

This code performs the same method as in ComputeTruePathLengthLimit of G4UrbanMscModel, as documented in section 8.1.6 of the Geant4 10.7 Physics Reference Manual or CERN-OPEN-2006-077 by L. Urban.

class UrbanMscScatter#: Sample angular change and lateral displacement with the Urban multiple scattering model.

Note

This code performs the same method as in G4VMultipleScattering::AlongStepDoIt and G4UrbanMscModel::SampleScattering of the Geant4 10.7 release.

template<class XSCorrector> class SBEnergyDistribution#

Sample exiting photon energy from Bremsstrahlung.

The SB energy distribution uses tabulated Seltzer-Berger differential cross section data from G4EMLOW, which stores scaled cross sections as a function of incident particle energy and exiting gamma energy (see SeltzerBergerModel for details). The sampling procedure is roughly laid out in section [PHYS341] of the GEANT3 physics reference manual, although like Geant4 we use raw tabulated SB data rather than a parameter fit. Also like Geant4 we include the extra density correction factor.

Once an element and incident energy have been selected, the exiting energy distribution is the differential cross section, which is stored as a scaled tabulated value. The reconstructed cross section gives the pdf

\[ p(\kappa) \propto \frac{\dif \sigma}{\dif k}s \propto \frac{1}{\kappa} \chi_Z(E, \kappa) \]

where \( \kappa = k / E \) is the ratio of the emitted photon energy to the incident charged particle energy, and the domain of p is nominally restricted by the allowable energy range \( \kappa_c < \kappa < 1 \), where \( \kappa_c \) is from the cutoff energy \( E_c \) for secondary gamma production. This distribution is sampled by decomposing it into an analytic sampling of \( 1/\kappa \) and a rejection sample on the scaled cross section \( \chi \). The analytic sample over the domain is from

\[ p_1(\kappa) \frac{1}{\ln 1/\kappa_c} \frac{1}{\kappa} \]

by sampling

\[ \kappa = \exp( \xi \ln \kappa_c ) \,, \]

and the rejection sample is the ratio of the scaled differential cross section at the exiting energy to the maximum across all exiting energies.

\[ p_2(\kappa) = \frac{\chi_Z(E, \kappa)}{\max_\kappa \chi_Z(E, \kappa)} \]

Since the tabulated data is bilinearly interpolated in incident and exiting energy, we can calculate a bounding maximum by precalculating (at setup time) the index of the maximum value of \( \chi \) and linearly interpolating those maximum values based on the incident energy.

The above algorithm is idealized; in practice, the minimum and maximum values are adjusted for a constant factor \(d_\rho\), which depends on the incident particle mass + energy and the material’s electron density:

\[ \kappa_\mathrm{min} = \ln (E_c^2 + d_\rho E^2) \]

and

\[ \kappa_\mathrm{max} = \ln (E^2 + d_\rho E^2) \, . \]

With this correction, the sample of the exiting gamma energy \( k = \kappa / E \) is transformed from the simple exponential above to:

\[ k = \sqrt{ \exp(\kappa_\mathrm{min} + \xi [\kappa_\mathrm{max} - \kappa_\mathrm{min}]) - d_\rho E^2} \]

Most of the mechanics of the sampling are in the template-free SBEnergyDistHelper, which is passed as a construction argument to this sampler. The separate class exists here to minimize duplication of templated code, which is required to provide for an on-the-fly correction of the cross section sampling.

class SBPositronXsCorrector#

Scale SB differential cross sections for positrons.

This correction factor is the ratio of the positron to electron cross sections. It appears in the bowels of G4SeltzerBergerModel::SampleEnergyTransfer in Geant4 and scales the SB cross section by a factor

\[ \frac{\sigma^+}{\sigma^-} = \exp( 2 \pi \alpha Z [ \beta^{-1}(E - k_c) - \beta^{-1}(E - k) ] ) \]

where the inverse positron speed is :

\[ \beta^{-1}(E) = \frac{c}{v} = \sqrt{1 - \left( \frac{m_e c^2}{E + m_e c^2} \right)^2} \,, \]

\( \alpha \) is the fine structure constant, \(E\) is the incident positron kinetic energy, \(k_c\) is the gamma production cutoff energy, and \( k \) is the provisionally sampled energy of the emitted photon.

\( \frac{\sigma^+}{\sigma^-} = 1 \) at the low end of the spectrum (where \( k / E = 0 \)) and is 0 at the tip of the spectrum where \( k / E \to 1 \).

The correction factor is described in:

Kim, Longhuan, R. H. Pratt, S. M. Seltzer, and M. J. Berger. “Ratio of Positron to Electron Bremsstrahlung Energy Loss: An Approximate Scaling Law.” Physical Review A 33, no. 5 (May 1, 1986): 3002–9. https://doi.org/10.1103/PhysRevA.33.3002.

Physics data#

Celeritas reads physics data from Geant4 (or from a ROOT file exported from data previously loaded into Geant4). Different versions of Geant4 (and Geant4 data) can be used seamlessly with any version of Celeritas, allowing differences to be isolated without respect to machine or model implementation. The following classes enumerate all the data used at runtime.

struct ImportData#

Store imported physics data from external sources.

All the data imported to Celeritas is stored in this single entity. This struct can be used in memory or recorded in a ROOT TBranch as a single TTree entry, which will be read by RootImporter to load the data into Celeritas. Currently, the TTree and TBranch names are hardcoded as geant4_data and ImportData in RootImporter .

Each entity’s id is defined by its vector position. An ImportElement with id = 3 is stored at elements[3] . Same for materials and volumes.

Seltzer-Berger, Livermore PE, and atomic relaxation data are loaded based on atomic numbers, and thus are stored in maps. To retrieve specific data use find(atomic_number) .

The “processes” field may be empty for testing applications.

Type aliases

using ZInt = int#

using ImportSBMap = std::map<ZInt, ImportSBTable>#

using ImportLivermorePEMap = std::map<ZInt, ImportLivermorePE>#

using ImportAtomicRelaxationMap = std::map<ZInt, ImportAtomicRelaxation>#

Public Members

std::vector<ImportParticle> particles#

std::vector<ImportIsotope> isotopes#

std::vector<ImportElement> elements#

std::vector<ImportMaterial> materials#

std::vector<ImportProcess> processes#

std::vector<ImportMscModel> msc_models#

std::vector<ImportVolume> volumes#

ImportEmParameters em_params#

ImportTransParameters trans_params#

ImportSBMap sb_data#

ImportLivermorePEMap livermore_pe_data#

ImportAtomicRelaxationMap atomic_relaxation_data#

Material and geometry properties#

struct ImportElement#

Store element data.

IsotopeIndex maps the isotope in the ImportData::isotopes vector.

type aliases

using IsotopeIndex = unsigned int#

using IsotopeFrac = std::pair<IsotopeIndex, double>#

using VecIsotopeFrac = std::vector<IsotopeFrac>#

Public Members

std::string name#

int atomic_number#

double atomic_mass#: [amu]

double radiation_length_tsai#: [mass/length^2]

double coulomb_factor#

VecIsotopeFrac isotopes_fractions#: Isotopic fractional abundance.

struct ImportProductionCut#

Store particle production cut.

Public Members

double energy = {}#: [MeV]

double range = {}#: [length]

struct ImportMatElemComponent#

Store elemental composition of a given material.

Public Members

unsigned int element_id = {}#: Index of element in ImportElement.

double mass_fraction = {}#: [mass/length^3]

double number_fraction = {}#

struct ImportMaterial#

Store material data.

Type aliases

using MapIntCutoff = std::map<int, ImportProductionCut>#

using VecComponent = std::vector<ImportMatElemComponent>#

Public Members

std::string name = {}#

ImportMaterialState state = {ImportMaterialState::size_}#

double temperature#: [K]

double density#: [mass/length^3]

double electron_density#: [1/length^3]

double number_density#: [1/length^3]

double radiation_length#: [length]

double nuclear_int_length#: [length]

MapIntCutoff pdg_cutoffs#: Cutoff per PDG.

VecComponent elements#

struct ImportVolume#

Store logical volume properties.

Public Functions

inline explicit operator bool() const#: Whether this represents a physical volume or is just a placeholder.

Public Members

int material_id = {-1}#

std::string name#

std::string solid_name#

struct ImportTransParameters#

Parameters related to transportation.

The looping thresholds are particle-dependent and stored in a map where the keys are the PDG number.

Type aliases

using PDGInt = int#

using ImportLoopingMap = std::unordered_map<PDGInt, ImportLoopingThreshold>#

Public Functions

inline explicit operator bool() const#: Whether parameters are assigned and valid.

Public Members

ImportLoopingMap looping#: Thresholds for killing looping tracks.

int max_substeps = {1000}#: Maximum number of substeps in the field propagator.

struct ImportLoopingThreshold#

Particle-dependent parameters for killing looping tracks.

Public Functions

inline explicit operator bool() const#: Whether parameters are assigned and valid.

Public Members

int threshold_trials = {10}#: Number of steps a higher-energy looping track takes before it’s killed.

double important_energy = {250}#: Energy below which looping tracks are immediately killed [MeV].

enum class celeritas::ImportMaterialState#

Enum for storing G4State enumerators.

Values:

enumerator other#

enumerator solid#

enumerator liquid#

enumerator gas#

enumerator size_#

Physics properties#

struct ImportParticle#

Store particle data.

Public Members

std::string name#

int pdg = {0}#

double mass = {0}#: [MeV]

double charge = {0}#: [Multiple of electron charge value]

double spin = {0}#: [Multiple of hbar]

double lifetime = {0}#: [time]

bool is_stable = {false}#

struct ImportProcess#

Store physics process data.

Note

ImportPhysicsTable is process and type (lambda, dedx, and so on) dependent, with each table type including physics vectors for all materials. Therefore, the physics vector of a given material is retrieved by finding the appropriate table_type in the tables vector and selecting the material: table.physics_vectors.at(material_id) .

Public Functions

inline explicit operator bool() const#

Public Members

int particle_pdg = {0}#

int secondary_pdg = {0}#

ImportProcessType process_type = {ImportProcessType::size_}#

ImportProcessClass process_class = {ImportProcessClass::size_}#

std::vector<ImportModel> models#

std::vector<ImportPhysicsTable> tables#

struct ImportModel#

Imported data for one model of a process.

This is always for a particular particle type since we import Processes as being for a particular particle.

The materials vector must always be assigned since we want the lower cutoff energy for each model.

Public Functions

inline explicit operator bool() const#

Public Members

ImportModelClass model_class = {ImportModelClass::size_}#

std::vector<ImportModelMaterial> materials#

struct ImportMscModel#

Store imported data for multiple scattering.

Public Functions

inline explicit operator bool() const#

Public Members

int particle_pdg = {0}#

ImportModelClass model_class = {ImportModelClass::size_}#

ImportPhysicsTable xs_table#

struct ImportModelMaterial#

Imported data for one material in a particular model.

Microscopic cross-section data are stored in the element-selector physics vector is in length^2. They will not be present for all model types, as some models only do on-the-fly calculation (e.g., photoelectric effect) or don’t depend on elemental interactions (e.g., compton scattering). The needs_micro_xs function indicates which models should store the cross section data.

The energy grid’s boundaries determine the model’s energy bounds and will always be set.

Type aliases

using MicroXs = std::vector<double>#

Public Functions

inline explicit operator bool() const#

Public Members

std::vector<double> energy#: Energy grid for the material.

std::vector<MicroXs> micro_xs#: Cross sections for each element.

Public Static Attributes

static constexpr auto energy_units = {ImportUnits::mev}#

static constexpr auto xs_units = {ImportUnits::len_sq}#

struct ImportPhysicsTable#

Imported physics table.

Each table stores physics vectors for all materials.

Public Functions

inline explicit operator bool() const#

Public Members

ImportTableType table_type = {ImportTableType::size_}#

ImportUnits x_units = {ImportUnits::none}#

ImportUnits y_units = {ImportUnits::none}#

std::vector<ImportPhysicsVector> physics_vectors#

struct ImportPhysicsVector#

Store imported physics vector data [see Geant4’s G4PhysicsVector.hh].

Each vector’s x axis is structured according to the vector_type. X is usually energy, but (as in the case of “inverse range”) can be distance or any other arbitrary value.

Public Members

ImportPhysicsVectorType vector_type#

std::vector<double> x#: Geant4 binVector.

std::vector<double> y#: Geant4 dataVector.

enum class celeritas::ImportProcessType#

Category of physics process.

See Geant4’s G4ProcessType.hh for the equivalent enum.

Values:

enumerator other#

enumerator transportation#

enumerator electromagnetic#

enumerator optical#

enumerator hadronic#

enumerator photolepton_hadron#

enumerator decay#

enumerator general#

enumerator parameterisation#

enumerator user_defined#

enumerator parallel#

enumerator phonon#

enumerator ucn#

enumerator size_#

enum class celeritas::ImportProcessClass#

Enumerator for the available physics processes.

This enum was created to safely access the many physics tables imported.

Values:

enumerator other#

enumerator ion_ioni#

enumerator msc#

enumerator h_ioni#

enumerator h_brems#

enumerator h_pair_prod#

enumerator coulomb_scat#

enumerator e_ioni#

enumerator e_brems#

enumerator photoelectric#

enumerator compton#

enumerator conversion#

enumerator rayleigh#

enumerator annihilation#

enumerator mu_ioni#

enumerator mu_brems#

enumerator mu_pair_prod#

enumerator gamma_general#

enumerator size_#

enum class celeritas::ImportModelClass#

Enumerator for the available physics models.

This enum was created to safely access the many imported physics tables.

Todo:: reorganize by physics list (major) and particle (minor) so that newly supported models are appended cleanly to the end of the list.

Values:

enumerator other#

enumerator bragg_ion#

enumerator bethe_bloch#

enumerator urban_msc#

enumerator icru_73_qo#

enumerator wentzel_vi_uni#

enumerator h_brems#

enumerator h_pair_prod#

enumerator e_coulomb_scattering#

enumerator bragg#

enumerator moller_bhabha#

enumerator e_brems_sb#

enumerator e_brems_lpm#

enumerator e_plus_to_gg#

enumerator livermore_photoelectric#

enumerator klein_nishina#

enumerator bethe_heitler#

enumerator bethe_heitler_lpm#

enumerator livermore_rayleigh#

enumerator mu_bethe_bloch#

enumerator mu_brems#

enumerator mu_pair_prod#

enumerator fluo_photoelectric#

enumerator goudsmit_saunderson#

enumerator size_#

enum class celeritas::ImportTableType#

Property being described by the physics table.

These are named based on accessors in G4VEnergyLossProcess, with one new table type, dedx_process, introduced to disambiguate the tables. In Geant4, the dedx table belonging to the ionization process is actually the sum of the de/dx for all processes that contribute to energy loss for the given particle, while the dedx tables for the remaining processes are the per-process energy loss. Here the tables are named to distinguish the summed energy loss (dedx) from the energy loss for an individual process (dedx_process). The ionization table is really just the dedx_process table for ionization, so it is redundant. The range table is calculated from the summed dedx table.

Values:

enumerator lambda#: Macroscopic cross section.

enumerator lambda_prim#: Cross section scaled by energy.

enumerator dedx#: Energy loss summed over processes.

enumerator range#: Integrated inverse energy loss.

enumerator msc_xs#: Scaled transport cross section.

enumerator size_#

enum class celeritas::ImportUnits#

Units of a physics table.

These depend on the unit selection in the main import data

Values:

enumerator none#: Unitless.

enumerator mev#: Energy [MeV].

enumerator mev_per_len#: Energy loss [MeV/len].

enumerator mev_per_cm#: Deprecated.

enumerator len#: Range [len].

enumerator cm#: Deprecated.

enumerator len_inv#: Macroscopic xs [1/len].

enumerator cm_inv#: Deprecated.

enumerator len_mev_inv#: Scaled [1/E] macroscopic xs [1/len-MeV].

enumerator cm_mev_inv#: Deprecated.

enumerator mev_sq_per_len#: Scaled [E^2] macroscopic xs [MeV^2/len].

enumerator mev_2_per_cm#: Deprecated.

enumerator len_sq#: Microscopic cross section [len^2].

enumerator cm_2#: Deprecated.

enumerator size_#

enum class celeritas::ImportPhysicsVectorType#

Geant4 equivalent enum for Physics vector types.

[See Geant4’s G4PhysicsVectorType.hh]

Values:

enumerator unknown#

enumerator linear#: Uniform and linear in x.

enumerator log#: Uniform and logarithmic in x.

enumerator free#: Nonuniform in x.

EM data#

struct ImportEmParameters#

Common electromagnetic physics parameters (see G4EmParameters.hh).

Note

Geant4 v11 removed the Spline() option from G4EmParameters.hh.

Note

The Geant4 MSC models use the values in G4EmParameters as the defaults; however, the MSC parameters can also be set directly using the model setter methods (there is no way to retrieve the values from the model in that case).

Public Functions

inline explicit operator bool() const#: Whether parameters are assigned and valid.

Public Members

bool energy_loss_fluct = {false}#: Energy loss fluctuation.

bool lpm = {true}#: LPM effect for bremsstrahlung and pair production.

bool integral_approach = {true}#: Integral cross section rejection.

double linear_loss_limit = {0.01}#: Slowing down threshold for linearity assumption.

double lowest_electron_energy = {0.001}#: Lowest e-/e+ kinetic energy [MeV].

bool auger = {false}#: Whether auger emission should be enabled (valid only for relaxation)

double msc_range_factor = {0.04}#: MSC range factor for e-/e+.

double msc_safety_factor = {0.6}#: MSC safety factor.

double msc_lambda_limit = {1 * units::millimeter}#: MSC lambda limit [length].

bool apply_cuts = {false}#: Kill secondaries below production cut.

double screening_factor = {1}#: Nuclear screening factor for single/multiple Coulomb scattering.

struct ImportAtomicTransition#

EADL transition data for atomic relaxation for a single element.

Public Members

int initial_shell = {}#: Originating shell designator.

int auger_shell = {}#: Auger shell designator.

double probability = {}#: Transition probability.

double energy = {}#: Transition energy [MeV].

struct ImportAtomicSubshell#

Public Members

int designator = {}#: Subshell designator.

std::vector<ImportAtomicTransition> fluor#: Radiative transitions.

std::vector<ImportAtomicTransition> auger#: Non-radiative transitions.

struct ImportAtomicRelaxation#

Public Members

std::vector<ImportAtomicSubshell> shells#

struct ImportLivermoreSubshell#

Livermore EPICS2014 photoelectric cross section data for a single element.

Public Members

double binding_energy#: Ionization energy [MeV].

std::vector<double> param_lo#: Low energy xs fit parameters.

std::vector<double> param_hi#: High energy xs fit parameters.

std::vector<double> xs#: Tabulated cross sections [b].

std::vector<double> energy#: Tabulated energies [MeV].

struct ImportLivermorePE#

Public Members

ImportPhysicsVector xs_lo#: Low energy range tabulated xs [b].

ImportPhysicsVector xs_hi#: High energy range tabulated xs [b].

double thresh_lo#: Threshold for low energy fit [MeV].

double thresh_hi#: Threshold for high energy fit [MeV].

std::vector<ImportLivermoreSubshell> shells#

struct ImportSBTable#

Seltzer Berger differential cross sections for a single element.

Public Members

std::vector<double> x#: Log energy of incident particle / MeV.

std::vector<double> y#: Ratio of gamma energy to incident energy.

std::vector<double> value#: Scaled DCS [mb] [ny * x + y].