Advanced documentation

The “advanced documentation” is the almost complete documentation of all modules.

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='')

A class to check a file integrity against provided checksum

filename

path to the file

Type:str

checksum

hex of sha256 checksum

Type:str

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)

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)

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

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)

Energy dependent Z-factor according to Lipari (1993).

get_solution(particle_name, mag=0.0, grid_idx=None, integrate=False, return_as='kinetic energy')

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 kaon k_
• conventional leptons originating neither from pion nor from kaon decay are collected in a category without any prefix, e.g. numu or mu+

Parameters:

• particle_name (str) – The name of the particle such, e.g. total_mu+ for the total flux spectrum of positive muons or pr_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)

n_e(grid_idx=None, min_energy_cutoff=0.1)

Returns the number of electrons plus positrons at a grid step above min_energy_cutoff.

Parameters:

• grid_idx (int) – Depth grid index (for profiles)
• min_energy_cutoff (float) – Energy threshold > mceq_config.e_min

n_mu(grid_idx=None, min_energy_cutoff=0.1)

Returns the number of positive and negative muons at a grid step above min_energy_cutoff.

Parameters:

• grid_idx (int) – Depth grid index (for profiles)
• min_energy_cutoff (float) – Energy threshold > mceq_config.e_min

n_particles(label, grid_idx=None, min_energy_cutoff=0.1)

Returns number of particles of type label at a grid step above an energy threshold for counting.

Parameters:

• label (str) – Particle name
• grid_idx (int) – Depth grid index (for profiles)
• min_energy_cutoff (float) – Energy threshold > mceq_config.e_min

regenerate_matrices(skip_decay_matrix=False)

Call this function after applying particle prod. modifications aka Barr parameters

set_density_model(density_config)

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 or MCEq.geometry.density_profiles.GeneralizedTarget.

Parameters:density_config (tuple of strings) – (parametrization type, arguments)

set_initial_spectrum(spectrum, pdg_id=None, append=False)

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)

Sets interaction model and/or an external charm model for calculation.

Decay and interaction matrix will be regenerated automatically after performing this call.

Parameters:

• interaction_model (str) – name of interaction model
• charm_model (str, optional) – name of charm model
• force (bool) – force loading interaction model

set_mod_pprod(prim_pdg, sec_pdg, x_func, x_func_args, delay_init=False)

Sets combination of projectile/secondary for error propagation.

The production spectrum of sec_pdg in interactions of prim_pdg is modified according to the function passed to InteractionYields.init_mod_matrix()

Parameters:

• prim_pdg (int) – interacting (primary) particle PDG ID
• sec_pdg (int) – secondary particle PDG ID
• x_func (object) – reference to function
• x_func_args (tuple) – arguments passed to x_func
• delay_init (bool) – Prevent init of mceq matrices if you are planning to add more modifications

set_primary_model(mclass, tag)

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

set_single_primary_particle(E, corsika_id=None, pdg_id=None, append=False)

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:

• E (float) – kinetic energy of a nucleus in GeV
• corsika_id (int) – ID of a nucleus (see text)
• pdg_id (int) – PDG ID of a particle
• append (bool) – If True, keep previous state and append a new particle.

set_theta_deg(theta_deg)

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)

Launches the solver.

The setting integrator in the config file decides which solver to launch.

Parameters:

• int_grid (list) – list of depths at which results are recorded
• grid_var (str) – Can be depth X or something else (currently only X supported)
• kwargs (dict) – Arguments are passed directly to the solver methods.

unset_mod_pprod(dont_fill=False)

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

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)

class MCEq.core.MatrixBuilder(particle_manager)

This class constructs the interaction and decay matrices.

construct_matrices(skip_decay_matrix=False)

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)

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)

Bundles different particle properties for simplified availability of particle properties in MCEq.core.MCEqRun.

Parameters:

• pdg_id (int) – PDG ID of the particle
• egrid (np.array, optional) – energy grid (centers)
• cs_db (object, optional) – reference to an instance of InteractionYields

add_decay_channel(child, dec_matrix, force=False)

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)

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)

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:

• kin_energy (float) – approximate interaction energy
• prim_pdg (int) – PDG ID of projectile
• sec_pdg (int) – PDG ID of secondary particle
• verbose (bool) – print out the closest energy

Returns:

,

Return type:

(numpy.array, numpy.array)

dN_dxf(energy, prim_pdg, sec_pdg, pos_only=True, verbose=True, **kwargs)

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:

• energy (float) – approximate interaction lab. energy
• prim_pdg (int) – PDG ID of projectile
• sec_pdg (int) – PDG ID of secondary particle
• verbose (bool) – print out the closest energy

Returns:

,

Return type:

(numpy.array, numpy.array)

dN_dxlab(kin_energy, sec_pdg, verbose=True, **kwargs)

Returns for interaction energy close to kin_energy for hadron-air collisions.

The function respects modifications applied via _set_mod_pprod().

Parameters:

• kin_energy (float) – approximate interaction kin_energy
• prim_pdg (int) – PDG ID of projectile
• sec_pdg (int) – PDG ID of secondary particle
• verbose (bool) – print out the closest enerkin_energygy

Returns:

,

Return type:

(numpy.array, numpy.array)

dNdec_dxlab(kin_energy, sec_pdg, verbose=True, **kwargs)

Returns for interaction energy close to kin_energy for hadron-air collisions.

The function respects modifications applied via _set_mod_pprod().

Parameters:

• kin_energy (float) – approximate interaction energy
• prim_pdg (int) – PDG ID of projectile
• sec_pdg (int) – PDG ID of secondary particle
• verbose (bool) – print out the closest energy

Returns:

,

Return type:

(numpy.array, numpy.array)

inel_cross_section(mbarn=False)

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)

Add custom particle type. (Incomplete and not debugged)

inverse_decay_length(cut=True)

Returns inverse decay length (or infinity (np.inf), if particle is stable), where the air density is factorized out.

Parameters:

• E (float) – energy in laboratory system in GeV
• cut (bool) – set to zero in ‘resonance’ regime

Returns:

in 1/cm

Return type:

(float)

inverse_interaction_length()

Returns inverse interaction length for A_target given by config.

Returns: in cm**2/g Return type:(float)

is_child(particle_ref)

True if this particle decays into particle_ref.

is_secondary(particle_ref)

True if this projectile and produces particle particle_ref.

set_cs(cs_db)

Set cross section adn recalculate the dependent variables

set_decay_channels(decay_db, pmanager)

Populates decay channel and energy distributions

set_hadronic_channels(hadronic_db, pmanager)

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)

Database for objects of MCEqParticle.

Authors:

Anatoli Fedynitch (DESY) Jonas Heinze (DESY)

add_tracking_particle(parent_list, child_pdg, alias_name, from_interactions=False)

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:

• alias (str) – Name alias under which the result is accessible in get_solution
• parents (list) – list of parent particle PDG ID’s
• child (int) – Child particle
• from_interactions (bool) – track particles from interactions

keys()

Returns pdg_ids of all particles

set_continuous_losses(contloss_db)

Set continuous losses terms to particles with ionization and radiation losses.

set_cross_sections_db(cs_db)

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)

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)

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)

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

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

class MCEq.data.Decays(mceq_hdf_db, default_decay_dset='full_decays')

Class for managing the dictionary of decay yield matrices.

Parameters:mceq_hdf_db (object) – instance of MCEq.data.HDF5Backend

get_matrix(parent, child)

Returns a DIM x DIM decay matrix.

Parameters:

• parent (int) – PDG ID of parent particle
• child (int) – PDG ID of final state child particle

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

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

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)

Returns inelastic parent-air cross-section as vector spanned over the energy grid.

Parameters:

• parent (int) – PDG ID of parent particle
• mbarn (bool,optional) – if True, the units of the cross-section will be , else

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

class MCEq.data.Interactions(mceq_hdf_db)

Class for managing the dictionary of interaction yield matrices.

Args: mceq_hdf_db (object): instance of MCEq.data.HDF5Backend

get_matrix(parent, child)

Returns a DIM x DIM yield matrix.

Parameters:

• parent (int) – PDG ID of parent particle
• child (int) – PDG ID of final state child/secondary particle

Returns:

yield matrix

Return type:

numpy.array

print_mod_pprod()

Prints the active particle production modification.

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)

This class handles the transfer between CPU and GPU memory, and the calling of GPU kernels. Initialized by MCEq.core.MCEqRun and used by solv_CUDA_sparse().

alloc_grid_sol(dim, nsols)

Allocates memory for intermediate if grid solution requested.

dump_sol()

Saves current solution to a new index in grid solution memory.

get_gridsol()

Downloads grid solution to main memory.

get_phi()

Downloads current solution from GPU memory.

set_matrices(int_m, dec_m)

Upload sparce matrices to GPU memory

set_phi(phi)

Uploads initial condition to GPU memory.

solve_step(rho_inv, dX)

Makes one solver step on GPU using cuSparse (BLAS)

MCEq.solvers.solv_CUDA_sparse(nsteps, dX, rho_inv, context, phi, grid_idcs)

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)

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)

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)

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

From https://gist.github.com/techtonik/2151727

MCEq.misc.corsikaid2pdg(corsika_id)

Conversion of CORSIKA nuclear code to PDG nuclear code

MCEq.misc.gen_xmat(energy_grid)

Generates x_lab matrix for a given energy grid

MCEq.misc.getAZN(pdg_id)

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)

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)

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)

Returns True if particle ID belongs to a heavy (charm) hadron.

MCEq.misc.pdg2corsikaid(pdg_id)

Conversion from nuclear PDG ID to CORSIKA ID.

Note:

PDG ID for nuclei is coded according to 10LZZZAAAI.
For iron-52 it is 1000260520.

MCEq.misc.print_in_rows(min_dbg_level, str_list, n_cols=5)

Prints contents of a list in rows n_cols entries per row.

MCEq.misc.theta_deg(cos_theta)

Converts to in degrees.

MCEq.misc.theta_rad(theta)

Converts from rad to degrees.