Advanced documentation¶
The “advanced documentation” is the almost complete documentation of all modules.
mceq_config
– default configuration optionsMCEq.core
– Core moduleMCEq.particlemanager
– Particle managerMCEq.data
– Data handlingMCEq.solvers
– ODE solver implementations- Miscellaneous
mceq_config
– default configuration options¶
These are all options MCEq accepts. Usually there is no need to change, except for advanced scenarios. Check out the file for a better formatted description and some advanced settings not contained in the list below.
-
class
mceq_config.
FileIntegrityCheck
(filename, checksum='')[source]¶ A class to check a file integrity against provided checksum
-
is_passed():
returns True if checksum and calculated checksum of the file are equal
-
get_file_checksum():
returns checksum of the file
-
hashlib
= <module 'hashlib' from '/home/docs/.pyenv/versions/2.7.18/lib/python2.7/hashlib.pyc'>¶
-
-
class
mceq_config.
MCEqConfigCompatibility
(namespace)[source]¶ This class provides access to the attributes of the module as a dictionary, as it was in the previous versions of MCEq
This method is deprecated and will be removed in future.
-
mceq_config.
A_target
= 14.65672¶ Average mass of target (for interaction length calculations) Change parameter only in combination with interaction model setting. By default all particle production matrices are calculated for air targets expect those for models with ‘_pp’ suffix. These are valid for hydrogen targets. <A> = 14.6568 for air as below (source https://en.wikipedia.org/wiki/Atmosphere_of_Earth)
-
mceq_config.
adv_set
= {'allowed_projectiles': [], 'disable_charm_pprod': False, 'disable_decays': [], 'disable_direct_leptons': False, 'disable_interactions_of_unstable': False, 'disable_leading_mesons': False, 'disabled_particles': [], 'exclude_from_mixing': [13], 'force_resonance': [], 'no_mixing': False}¶ Advanced settings (some options might be obsolete/not working)
-
mceq_config.
assume_nucleon_interactions_for_exotics
= True¶ Assume nucleon, pion and kaon cross sections for interactions of rare or exotic particles (mostly relevant for non-compact mode)
-
mceq_config.
average_loss_operator
= True¶ Improve (explicit solver) stability by averaging the continous loss operator
-
mceq_config.
cuda_fp_precision
= 32¶ CUDA Floating point precision (default 32-bit ‘float’)
-
mceq_config.
cuda_gpu_id
= 0¶ Select CUDA device ID if you have multiple GPUs
-
mceq_config.
dXmax
= 10.0¶ Maximal integration step dX in g/cm2. No limit necessary in most cases, use for debugging purposes when searching for stability issues.
-
mceq_config.
data_dir
= '/home/docs/checkouts/readthedocs.org/user_builds/mceq/checkouts/latest/MCEq/data'¶ Directory where the data files for the calculation are stored
-
mceq_config.
debug_level
= 1¶ Debug flag for verbose printing, 0 silences MCEq entirely
-
mceq_config.
dedx_material
= 'air'¶ Material for ionization and radiation (=continuous) loss terms Currently available choices: ‘air’, ‘water’, ‘ice’
-
mceq_config.
density_model
= ('CORSIKA', ('BK_USStd', None))¶ (model, (arguments))
Type: Atmospheric model in the format
-
mceq_config.
e_max
= 100000000000.0¶ The maximal energy is 1e12 GeV, but not all interaction models run at such high energies. If you are interested in lower energies, reduce this value for inclusive calculations to max. energy of interest + 4-5 orders of magnitude. For single primaries the maximal energy is directly limited by this value. Smaller grids speed up the initialization and integration.
-
mceq_config.
e_min
= 0.1¶ Minimal energy for grid The minimal energy (technically) is 1e-2 GeV. Currently you can run into stability problems with the integrator with such low thresholds. Use with care and check results for oscillations and feasibility.
-
mceq_config.
em_db_fname
= 'mceq_db_EM_Tsai-Max_Z7.31.h5'¶ File name of the MCEq database
-
mceq_config.
enable_default_tracking
= True¶ Enable default tracking particles, such as pi_numu, pr_mu+, etc. If only total fluxes are of interest, disable this feature to gain performance since the eqution system becomes smaller and sparser
-
mceq_config.
enable_em
= False¶ Enable electromagnetic cascade with matrices from EmCA
-
mceq_config.
enable_em_ion
= False¶ enable EM ionization loss
-
mceq_config.
enable_muon_energy_loss
= True¶ Muon energy loss according to Kokoulin et al.
-
mceq_config.
env_density
= 0.001225¶ density of default material in g/cm^3
-
mceq_config.
excpt_on_missing_particle
= False¶ Raise exception when requesting unknown particles from get_solution
-
mceq_config.
hybrid_crossover
= 0.5¶ Ratio of decay_length/interaction_length where particle interactions are neglected and the resonance approximation is used 0.5 ~ precision loss <+3% speed gain ~ factor 10 If smoothness and shape accuracy for prompt flux is crucial, use smaller values around 0.1 or 0.05
-
mceq_config.
integrator
= 'euler'¶ Selection of integrator (euler/odepack)
-
mceq_config.
kernel_config
= 'numpy'¶ euler kernel implementation (numpy/MKL/CUDA). With serious nVidia GPUs CUDA a few times faster than MKL autodetection of fastest kernel below
-
mceq_config.
leading_process
= 'auto'¶ The leading process is can be either “decays” or “interactions”. This depends on the target density and it is usually chosen automatically. For advanced applications one can force “interactions” to be the dominant process. Essentially this affects how the adaptive step size is computed. There is also the choice of “auto” that takes both processes into account
-
mceq_config.
len_target
= 1000.0¶ Default parameters for GeneralizedTarget Total length of the target [m]
-
mceq_config.
loss_step_for_average
= 0.1¶ Step size (dX) for averaging
-
mceq_config.
low_energy_extension
= {'he_le_transition': 80, 'nbins_interp': 3, 'use_unknown_cs': True}¶ This is not used in the code as before, instead the low energy extension is compiled into the HDF backend files.
-
mceq_config.
max_density
= (0.001225,)¶ Approximate value for the maximum density expected. Needed for the resonance approximation. Default value: air at the surface
-
mceq_config.
mceq_db_fname
= 'mceq_db_lext_dpm191_v12.h5'¶ File name of the MCEq database
-
mceq_config.
mkl_threads
= 8¶ Number of MKL threads (for sparse matrix multiplication the performance advantage from using more than a few threads is limited by memory bandwidth) Irrelevant for GPU integrators, but can affect initialization speed if numpy is linked to MKL.
-
mceq_config.
muon_helicity_dependence
= True¶ Helicity dependent muons decays from analytical expressions
-
mceq_config.
ode_params
= {'method': 'bdf', 'name': 'lsoda', 'nsteps': 1000, 'rtol': 0.01}¶ parameters for the odepack integrator. More details at http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode
-
mceq_config.
override_debug_fcn
= []¶ Override debug prinput for functions listed here (just give the name, “get_solution” for instance) Warning, this option slows down initialization by a lot. Use only when needed.
-
mceq_config.
override_max_level
= 10¶ Override debug printout for debug levels < value for the functions above
-
mceq_config.
pf
= 'Linux-5.15.0-1004-aws-x86_64-with-debian-buster-sid'¶ Autodetect best solver determine shared library extension and MKL path
-
mceq_config.
print_module
= False¶ Print module name in debug output
-
mceq_config.
prompt_ctau
= 0.123¶ default ctau < 0.123 cm (that of D0)
Type: Definition of prompt
-
mceq_config.
r_E
= 6391000.0¶ parameters for EarthGeometry
-
mceq_config.
return_as
= 'kinetic energy'¶ The latest versions of MCEq work in kinetic energy not total energy If you want the result to be compatible with the previous choose ‘total energy’ else ‘kinetic energy’
-
mceq_config.
stability_margin
= 0.95¶ Stability margin for the integrator. The default 0.95 means that step sizes are chosen 5% away from the stability circle. Usually no need to change, except you know what it does.
-
mceq_config.
standard_particles
= [11, 12, 13, 14, 16, 211, 321, 2212, 2112, 3122, 411, 421, 431, -11, -12, -13, -14, -16, -211, -321, -2212, -2112, -3122, -411, -421, -431, 22, 111, 130, 310]¶ Particles for compact mode
-
mceq_config.
use_isospin_sym
= True¶ When using modified particle production matrices use isospin symmetries to determine the corresponding modification for neutrons and K0L/K0S
MCEq.core
– Core module¶
This module contains the main program features. Instantiating MCEq.core.MCEqRun
will initialize the data structures and particle tables, create and fill the
interaction and decay matrix and check if all information for the calculation
of inclusive fluxes in the atmosphere is available.
-
class
MCEq.core.
MCEqRun
(interaction_model, primary_model, theta_deg, **kwargs)[source]¶ Main class for handling the calculation.
This class is the main user interface for the caclulation. It will handle initialization and various error/configuration checks. The setup has to be accomplished before invoking the integration routine is
MCeqRun.solve()
. Changes of configuration, such as:- interaction model in
MCEqRun.set_interaction_model()
, - primary flux in
MCEqRun.set_primary_model()
, - zenith angle in
MCEqRun.set_theta_deg()
, - density profile in
MCEqRun.set_density_model()
, - member particles of the special
obs_
group inMCEqRun.set_obs_particles()
,
can be made on an active instance of this class, while calling
MCEqRun.solve()
subsequently to calculate the solution corresponding to the settings.The result can be retrieved by calling
MCEqRun.get_solution()
.Parameters: - interaction_model (string) – PDG ID of the particle
- density_model (string,sting,string) – model type, location, season
- primary_model (class, param_tuple) – classes derived from
crflux.models.PrimaryFlux
and its parameters as tuple - theta_deg (float) – zenith angle in degrees, measured positively from vertical direction
- adv_set (dict) – advanced settings, see
mceq_config
- obs_ids (list) – list of particle name strings. Those lepton decay
products will be scored in the special
obs_
categories
-
closest_energy
(kin_energy)[source]¶ Convenience function to obtain the nearest grid energy to the energy argument, provided as kinetik energy in lab. frame.
-
decay_z_factor
(parent_pdg, child_pdg)[source]¶ Energy dependent Z-factor according to Lipari (1993).
-
get_solution
(particle_name, mag=0.0, grid_idx=None, integrate=False, return_as='kinetic energy')[source]¶ Retrieves solution of the calculation on the energy grid.
Some special prefixes are accepted for lepton names:
- the total flux of muons, muon neutrinos etc. from all sources/mothers
can be retrieved by the prefix
total_
, i.e.total_numu
- the conventional flux of muons, muon neutrinos etc. from all sources
can be retrieved by the prefix
conv_
, i.e.conv_numu
- correspondigly, the flux of leptons which originated from the decay
of a charged pion carries the prefix
pi_
and from a kaonk_
- conventional leptons originating neither from pion nor from kaon
decay are collected in a category without any prefix, e.g.
numu
ormu+
Parameters: - particle_name (str) – The name of the particle such, e.g.
total_mu+
for the total flux spectrum of positive muons orpr_antinumu
for the flux spectrum of prompt anti muon neutrinos - mag (float, optional) – ‘magnification factor’: the solution is
multiplied by
sol
- grid_idx (int, optional) – if the integrator has been configured to save
intermediate solutions on a depth grid, then
grid_idx
specifies the index of the depth grid for which the solution is retrieved. If not specified the flux at the surface is returned - integrate (bool, optional) – return averge particle number instead of
- flux (multiply by bin width) –
Returns: flux of particles on energy grid
e_grid
Return type: (numpy.array)
- the total flux of muons, muon neutrinos etc. from all sources/mothers
can be retrieved by the prefix
-
n_e
(grid_idx=None, min_energy_cutoff=0.1)[source]¶ Returns the number of electrons plus positrons at a grid step above min_energy_cutoff.
Parameters:
-
n_mu
(grid_idx=None, min_energy_cutoff=0.1)[source]¶ Returns the number of positive and negative muons at a grid step above min_energy_cutoff.
Parameters:
-
n_particles
(label, grid_idx=None, min_energy_cutoff=0.1)[source]¶ Returns number of particles of type label at a grid step above an energy threshold for counting.
Parameters:
-
regenerate_matrices
(skip_decay_matrix=False)[source]¶ Call this function after applying particle prod. modifications aka Barr parameters
-
set_density_model
(density_config)[source]¶ Sets model of the atmosphere.
To choose, for example, a CORSIKA parametrization for the Southpole in January, do the following:
mceq_instance.set_density_model(('CORSIKA', ('PL_SouthPole', 'January')))
More details about the choices can be found in
MCEq.geometry.density_profiles
. Calling this method will issue a recalculation of the interpolation and the integration path.From version 1.2 and above, the density_config parameter can be a reference to an instance of a density class directly. The class has to be derived either from
MCEq.geometry.density_profiles.EarthsAtmosphere
orMCEq.geometry.density_profiles.GeneralizedTarget
.Parameters: density_config (tuple of strings) – (parametrization type, arguments)
-
set_initial_spectrum
(spectrum, pdg_id=None, append=False)[source]¶ Set a user-defined spectrum for an arbitrary species as initial condition.
This function is an equivalent to
set_single_primary_particle()
. It allows to define an arbitrary spectrum for each available particle species as initial condition for the integration. Set the append argument to True for subsequent species to define initial spectra combined from different particles.The (differential) spectrum has to be distributed on the energy grid as dN/dptot, i.e. divided by the bin widths and with the total momentum units in GeV(/c).
Parameters: - spectrum (np.array) – spectrum dN/dptot
- pdg_id (int) – PDG ID in case of a particle
-
set_interaction_model
(interaction_model, particle_list=None, update_particle_list=True, force=False, build_matrices=True)[source]¶ Sets interaction model and/or an external charm model for calculation.
Decay and interaction matrix will be regenerated automatically after performing this call.
Parameters:
-
set_mod_pprod
(prim_pdg, sec_pdg, x_func, x_func_args, delay_init=False)[source]¶ Sets combination of projectile/secondary for error propagation.
The production spectrum of
sec_pdg
in interactions ofprim_pdg
is modified according to the function passed toInteractionYields.init_mod_matrix()
Parameters:
-
set_primary_model
(mclass, tag)[source]¶ Sets primary flux model.
This functions is quick and does not require re-generation of matrices.
Parameters: - interaction_model (
CRFluxModel.PrimaryFlux
) – reference - primary model class (to) –
- tag (tuple) – positional argument list for model class
- interaction_model (
-
set_single_primary_particle
(E, corsika_id=None, pdg_id=None, append=False)[source]¶ Set type and kinetic energy of a single primary nucleus to calculation of particle yields.
The functions uses the superposition theorem, where the flux of a nucleus with mass A and charge Z is modeled by using Z protons and A-Z neutrons at energy The nucleus type is defined via . For example iron has the CORSIKA ID 5226.
Single leptons or hadrons can be defined by specifiying pdg_id instead of corsika_id.
The append argument can be used to compose an initial state with multiple particles. If it is False the initial condition is reset to zero before adding the particle.
A continuous input energy range is allowed between .
Parameters:
-
set_theta_deg
(theta_deg)[source]¶ Sets zenith angle as seen from a detector.
Currently only ‘down-going’ angles (0-90 degrees) are supported.
Parameters: theta_deg (float) – zenith angle in the range 0-90 degrees
-
solve
(int_grid=None, grid_var='X', **kwargs)[source]¶ Launches the solver.
The setting integrator in the config file decides which solver to launch.
Parameters:
-
unset_mod_pprod
(dont_fill=False)[source]¶ Removes modifications from
MCEqRun.set_mod_pprod()
.Parameters: - skip_fill (bool) – If true do not regenerate matrices
- to be done at a later step by hand) ((has) –
-
z_factor
(projectile_pdg, secondary_pdg, definition='primary_e')[source]¶ Energy dependent Z-factor according to Thunman et al. (1996)
-
dim
¶ Energy grid (dimension)
-
dim_states
¶ Number of cascade particles times dimension of grid (dimension of the equation system)
-
e_bins
¶ Energy grid (bin edges)
-
e_grid
¶ Energy grid (bin centers)
-
e_widths
¶ Energy grid (bin widths)
-
pman
= None¶ Particle manager (initialized/updated in set_interaction_model)
- interaction model in
-
class
MCEq.core.
MatrixBuilder
(particle_manager)[source]¶ This class constructs the interaction and decay matrices.
-
construct_matrices
(skip_decay_matrix=False)[source]¶ Constructs the matrices for calculation.
These are:
- ,
- .
For debug_levels >= 2 some general information about matrix shape and the number of non-zero elements is printed. The intermediate matrices and are deleted afterwards to save memory.
Set the
skip_decay_matrix
flag to avoid recreating the decay matrix. This is not necessary if, for example, particle production is modified, or the interaction model is changed.Parameters: skip_decay_matrix (bool) – Omit re-creating D matrix
-
cont_loss_operator
(pdg_id)[source]¶ Returns continuous loss operator that can be summed with appropriate position in the C matrix.
-
dim
¶ Energy grid (dimension)
-
dim_states
¶ Number of cascade particles times dimension of grid (dimension of the equation system)
-
MCEq.particlemanager
– Particle manager¶
The MCEq.particlemanager.ParticleManager
handles the bookkeeping of
MCEq.particlemanager.MCEqParticle`s. It feeds the parameterizations of interactions
and decays from :mod:`MCEq.data
into the corresponding variables and validates certain relations.
The construction of the interaction and decay matrices proceeds by iterating over the particles
in MCEq.particlemanager.ParticleManager
, querying the interaction and decay yields
for child particles. Therefore, there is usually no need to directly access any of the
classes in MCEq.data
.
-
class
MCEq.particlemanager.
MCEqParticle
(pdg_id, helicity, energy_grid=None, cs_db=None, init_pdata_defaults=True)[source]¶ Bundles different particle properties for simplified availability of particle properties in
MCEq.core.MCEqRun
.Parameters: -
add_decay_channel
(child, dec_matrix, force=False)[source]¶ Add a decay channel.
The branching ratios are not renormalized and one needs to take care of this externally.
-
add_hadronic_production_channel
(child, int_matrix)[source]¶ Add a new particle that is produced in hadronic interactions.
The int_matrix is expected to be in the correct shape and scale as the other interaction (dN/dE(i,j)) matrices. Energy conservation is not checked.
-
dN_dEkin
(kin_energy, sec_pdg, verbose=True, **kwargs)[source]¶ Returns in lab frame for an interaction energy close to
kin_energy
(total) for hadron-air collisions.The function respects modifications applied via
_set_mod_pprod()
.Parameters: Returns: ,
Return type: (numpy.array, numpy.array)
-
dN_dxf
(energy, prim_pdg, sec_pdg, pos_only=True, verbose=True, **kwargs)[source]¶ Returns in c.m. for interaction energy close to
energy
(lab. not kinetic) for hadron-air collisions.The function respects modifications applied via
_set_mod_pprod()
.Parameters: Returns: ,
Return type: (numpy.array, numpy.array)
-
dN_dxlab
(kin_energy, sec_pdg, verbose=True, **kwargs)[source]¶ Returns for interaction energy close to
kin_energy
for hadron-air collisions.The function respects modifications applied via
_set_mod_pprod()
.Parameters: Returns: ,
Return type: (numpy.array, numpy.array)
-
dNdec_dxlab
(kin_energy, sec_pdg, verbose=True, **kwargs)[source]¶ Returns for interaction energy close to
kin_energy
for hadron-air collisions.The function respects modifications applied via
_set_mod_pprod()
.Parameters: Returns: ,
Return type: (numpy.array, numpy.array)
-
inel_cross_section
(mbarn=False)[source]¶ Returns inelastic cross section.
Parameters: mbarn (bool) – if True cross section in mb otherwise in cm**2 Returns: in mb or cm**2 Return type: (float)
-
init_custom_particle_data
(name, pdg_id, helicity, ctau, mass, **kwargs)[source]¶ Add custom particle type. (Incomplete and not debugged)
-
inverse_decay_length
(cut=True)[source]¶ Returns inverse decay length (or infinity (np.inf), if particle is stable), where the air density is factorized out.
Parameters: Returns: in 1/cm
Return type: (float)
-
inverse_interaction_length
()[source]¶ Returns inverse interaction length for A_target given by config.
Returns: in cm**2/g Return type: (float)
-
set_hadronic_channels
(hadronic_db, pmanager)[source]¶ Changes the hadronic interaction model.
Replaces indexing of the yield dictionary from PDG IDs with references from partilcle manager.
-
A
= None¶ Mass, charge, neutron number
-
E_crit
= None¶ (float) critical energy in air at the surface
-
N
= None¶ Mass, charge, neutron number
-
Z
= None¶ Mass, charge, neutron number
-
can_interact
= None¶ (bool) can_interact
-
ctau
= None¶ (float) ctau in cm
-
dEdX
= None¶ (np.array) continuous losses in GeV/(g/cm2)
-
decay_dists
= None¶ decay channels if any
-
hadridx
¶ Returns index range where particle behaves as hadron.
Returns: range on energy grid Return type: tuple()
(int,int)
-
has_contloss
= None¶ (bool) has continuous losses dE/dX defined
-
helicity
= None¶ (int) helicity -1, 0, 1 (0 means undefined or average)
-
is_em
= None¶ (bool) if it’s an electromagnetic particle
-
is_hadron
= None¶ (bool) particle is a hadron
-
is_lepton
= None¶ (bool) particle is a lepton
-
is_mixed
= None¶ (bool) particle has both, hadron and resonance properties
-
is_nucleus
= None¶ (bool) particle is a nucleus (not yet implemented)
-
is_projectile
= None¶ (bool) particle is interacting projectile
-
is_resonance
= None¶ (bool) if particle has just resonance behavior
-
is_stable
= None¶ (bool) particle is stable
-
is_tracking
= None¶ (bool) is a tracking particle
-
lidx
¶ Returns lower index of particle range in state vector.
Returns: lower index in state vector MCEqRun.phi
Return type: (int)
-
mass
= None¶ (float) mass in GeV
-
mceqidx
= None¶ (int) MCEq ID
-
name
= None¶ (str) species name in string representation
-
pdg_id
= None¶ (int) Particle Data Group Monte Carlo particle ID
-
residx
¶ Returns index range where particle behaves as resonance.
Returns: range on energy grid Return type: tuple()
(int,int)
-
uidx
¶ Returns upper index of particle range in state vector.
Returns: upper index in state vector MCEqRun.phi
Return type: (int)
-
unique_pdg_id
= None¶ (int) Unique PDG ID that is different for tracking particles
-
-
class
MCEq.particlemanager.
ParticleManager
(pdg_id_list, energy_grid, cs_db, mod_table=None)[source]¶ Database for objects of
MCEqParticle
.- Authors:
- Anatoli Fedynitch (DESY) Jonas Heinze (DESY)
-
add_tracking_particle
(parent_list, child_pdg, alias_name, from_interactions=False)[source]¶ Allows tracking decay and particle production chains.
Replaces previous
obs_particle
function that allowed to track only leptons from decays certain particles. This present feature removes the special PDG IDs 71XX, 72XX, etc and allows to define any channel like:$ particleManagerInstance.add_tracking_particle([211], 14, 'pi_numu')
This will store muon neutrinos from pion decays under the alias ‘pi_numu’. Multiple parents are allowed:
$ particleManagerInstance.add_tracking_particle( [411, 421, 431], 14, 'D_numu')
Parameters:
-
set_continuous_losses
(contloss_db)[source]¶ Set continuous losses terms to particles with ionization and radiation losses.
-
set_cross_sections_db
(cs_db)[source]¶ Sets the inelastic cross section to each interacting particle.
This applies to most of the hadrons and does not imply that the particle becomes a projectile. parents need in addition defined hadronic channels.
-
set_decay_channels
(decay_db)[source]¶ Attaches the references to the decay yield tables to each unstable particle
-
set_interaction_model
(cs_db, hadronic_db, updated_parent_list=None, force=False)[source]¶ Attaches the references to the hadronic yield tables to each projectile particle
-
track_leptons_from
(parent_pdg_list, prefix, exclude_em=True, from_interactions=False, use_helicities=False)[source]¶ Adds tracking particles for all leptons coming from decays of parents in parent_pdg_list.
-
mceqidx2pdg
= None¶ (dict) Converts index in state vector to PDG ID
-
mceqidx2pref
= None¶ (dict) Converts MCEq index to reference of class:particlemanager.MCEqParticle
-
nspec
= None¶ (int) Total number of species
-
pdg2mceqidx
= None¶ (dict) Converts PDG ID to index in state vector
-
pname2mceqidx
= None¶ (dict) Converts particle name to index in state vector
-
pname2pref
= None¶ (dict) Converts particle name to reference of class:particlemanager.MCEqParticle
MCEq.data
– Data handling¶
The tabulated data in MCEq is handled by HDF5Backend
. The HDF5 file densly packed
data, where matrices are stored as vectors of a sparse CSR data structure. Index dictionaries
and other metadata are stored as attributes. The other classes pf this module know how to
interact with the backend and provide an intermediate step to the ParticleManager
that propagates data further to the MCEqParticle
objects.
-
class
MCEq.data.
ContinuousLosses
(mceq_hdf_db, material='air')[source]¶ Class for managing the dictionary of hadron-air cross-sections.
Parameters: - mceq_hdf_db (object) – instance of
MCEq.data.HDF5Backend
- material (str) – name of the material (not fully implemented)
-
energy_grid
= None¶ reference to energy grid
-
index_d
= None¶ Dictionary containing the distribuiton matrices
-
mceq_db
= None¶ MCEq HDF5Backend reference
-
parents
= None¶ List of active parents
- mceq_hdf_db (object) – instance of
-
class
MCEq.data.
Decays
(mceq_hdf_db, default_decay_dset='full_decays')[source]¶ Class for managing the dictionary of decay yield matrices.
Parameters: mceq_hdf_db (object) – instance of MCEq.data.HDF5Backend
-
get_matrix
(parent, child)[source]¶ Returns a
DIM x DIM
decay matrix.Parameters: Returns: decay matrix
Return type: numpy.array
-
mceq_db
= None¶ MCEq HDF5Backend reference
-
parent_list
= None¶ (list) List of particles in the decay matrices
-
-
class
MCEq.data.
HDF5Backend
[source]¶ Provides access to tabulated data stored in an HDF5 file.
The file contains all necessary ingredients to run MCEq, i.e. no other files are required. This database is not maintained in git and it will change infrequently.
-
class
MCEq.data.
InteractionCrossSections
(mceq_hdf_db, interaction_model='SIBYLL2.3c')[source]¶ Class for managing the dictionary of hadron-air cross-sections.
Parameters: - mceq_hdf_db (object) – instance of
MCEq.data.HDF5Backend
- interaction_model (str) – name of the interaction model
-
get_cs
(parent, mbarn=False)[source]¶ Returns inelastic
parent
-air cross-section as vector spanned over the energy grid.Parameters: Returns: cross-section in or
Return type: numpy.array
-
GeV2mbarn
= 0.38937930376300284¶ unit -
-
GeVcm
= 1.9732696312541852e-14¶ unit -
-
GeVfm
= 0.19732696312541853¶ unit -
-
energy_grid
= None¶ reference to energy grid
-
iam
= None¶ (str) Interaction Model name
-
index_d
= None¶ Dictionary containing the distribuiton matrices
-
mbarn2cm2
= 9.999999999999999e-28¶ unit conversion -
-
mceq_db
= None¶ MCEq HDF5Backend reference
-
parents
= None¶ List of active parents
- mceq_hdf_db (object) – instance of
-
class
MCEq.data.
Interactions
(mceq_hdf_db)[source]¶ Class for managing the dictionary of interaction yield matrices.
Args: mceq_hdf_db (object): instance of
MCEq.data.HDF5Backend
-
get_matrix
(parent, child)[source]¶ Returns a
DIM x DIM
yield matrix.Parameters: Returns: yield matrix
Return type: numpy.array
-
description
= None¶ String containing the desciption of the model
-
energy_grid
= None¶ reference to energy grid
-
iam
= None¶ (str) Interaction Model name
-
index_d
= None¶ Dictionary containing the distribuiton matrices
-
mceq_db
= None¶ MCEq HDF5Backend reference
-
mod_pprod
= None¶ (tuple) modified particle combination for error prop.
-
parents
= None¶ List of active parents
-
particles
= None¶ List of all known particles
-
relations
= None¶ Dictionary parent/child relations
-
MCEq.solvers
– ODE solver implementations¶
The module contains functions which are called by MCEq.core.MCEqRun.solve()
method.
The implementation is a simple Forward-Euler stepper. The stability is under control since the smallest Eigenvalues are known a priory. The step size is “adaptive”, but it is deterministic and known before the integration starts.
The steps that each solver routine does are:
with
(1)¶
and
(2)¶
As one can easily see, each step can be represented by two sparse gemv calls and one vector addition. This is what happens in the MKL and CUDA functions below.
The fastest solver is using NVidia’s cuSparse library provided via the cupy matrix library. Intel MKL is recommended for Intel CPUs, in particular since MKL is using AVX instructions. The plain numpy solver is for compatibility and hacking, but not recommended for general use.
-
class
MCEq.solvers.
CUDASparseContext
(int_m, dec_m, device_id=0)[source]¶ This class handles the transfer between CPU and GPU memory, and the calling of GPU kernels. Initialized by
MCEq.core.MCEqRun
and used bysolv_CUDA_sparse()
.
-
MCEq.solvers.
solv_CUDA_sparse
(nsteps, dX, rho_inv, context, phi, grid_idcs)[source]¶ NVIDIA CUDA cuSPARSE implementation of forward-euler integration.
Function requires a working
accelerate
installation.Parameters: - nsteps (int) – number of integration steps
- dX (numpy.array[nsteps]) – vector of step-sizes in g/cm**2
- rho_inv (numpy.array[nsteps]) – vector of density values
- int_m (numpy.array) – interaction matrix (1) in dense or sparse representation
- dec_m (numpy.array) – decay matrix (2) in dense or sparse representation
- phi (numpy.array) – initial state vector
- mu_loss_handler (object) – object of type
SemiLagrangianEnergyLosses
Returns: state vector after integration
Return type: numpy.array
-
MCEq.solvers.
solv_MKL_sparse
(nsteps, dX, rho_inv, int_m, dec_m, phi, grid_idcs)[source]¶ Intel MKL sparse BLAS implementation of forward-euler integration.
Function requires that the path to the MKL runtime library
libmkl_rt.[so/dylib]
defined in the config file.Parameters: - nsteps (int) – number of integration steps
- dX (numpy.array[nsteps]) – vector of step-sizes in g/cm**2
- rho_inv (numpy.array[nsteps]) – vector of density values
- int_m (numpy.array) – interaction matrix (1) in dense or sparse representation
- dec_m (numpy.array) – decay matrix (2) in dense or sparse representation
- phi (numpy.array) – initial state vector
- grid_idcs (list) – indices at which longitudinal solutions have to be saved.
Returns: state vector after integration
Return type: numpy.array
-
MCEq.solvers.
solv_numpy
(nsteps, dX, rho_inv, int_m, dec_m, phi, grid_idcs)[source]¶ numpy
implementation of forward-euler integration.Parameters: - nsteps (int) – number of integration steps
- dX (numpy.array[nsteps]) – vector of step-sizes in g/cm**2
- rho_inv (numpy.array[nsteps]) – vector of density values
- int_m (numpy.array) – interaction matrix (1) in dense or sparse representation
- dec_m (numpy.array) – decay matrix (2) in dense or sparse representation
- phi (numpy.array) – initial state vector
Returns: state vector after integration
Return type: numpy.array
Miscellaneous¶
Different helper functions.
-
class
MCEq.misc.
energy_grid
(c, b, w, d)¶ Energy grid (centers, bind widths, dimension)
-
b
¶ Alias for field number 1
-
c
¶ Alias for field number 0
-
d
¶ Alias for field number 3
-
w
¶ Alias for field number 2
-
-
MCEq.misc.
caller_name
(skip=2)[source]¶ Get a name of a caller in the format module.class.method
skip specifies how many levels of stack to skip while getting caller name. skip=1 means “who calls me”, skip=2 “who calls my caller” etc. An empty string is returned if skipped levels exceed stack height.abs
-
MCEq.misc.
getAZN
(pdg_id)[source]¶ Returns mass number , charge and neutron number of
pdg_id
.Note:
PDG ID for nuclei is coded according to 10LZZZAAAI. For iron-52 it is 1000260520.
Parameters: pdgid (int) – PDG ID of nucleus/mass group Returns: (Z,A) tuple Return type: (int,int,int)
-
MCEq.misc.
getAZN_corsika
(corsikaid)[source]¶ Returns mass number , charge and neutron number of
corsikaid
.Parameters: corsikaid (int) – corsika id of nucleus/mass group Returns: (Z,A) tuple Return type: (int,int,int)
-
MCEq.misc.
info
(min_dbg_level, *message, **kwargs)[source]¶ Print to console if min_debug_level <= config.debug_level
The fuction determines automatically the name of caller and appends the message to it. Message can be a tuple of strings or objects which can be converted to string using str().
Parameters: - min_dbg_level (int) – Minimum debug level in config for printing
- message (tuple) – Any argument or list of arguments that casts to str
- condition (bool) – Print only if condition is True
- blank_caller (bool) – blank the caller name (for multiline output)
- no_caller (bool) – don’t print the name of the caller
- Authors:
- Anatoli Fedynitch (DESY) Jonas Heinze (DESY)
-
MCEq.misc.
is_charm_pdgid
(pdgid)[source]¶ Returns True if particle ID belongs to a heavy (charm) hadron.
-
MCEq.misc.
pdg2corsikaid
(pdg_id)[source]¶ Conversion from nuclear PDG ID to CORSIKA ID.
Note:
PDG ID for nuclei is coded according to 10LZZZAAAI. For iron-52 it is 1000260520.