MCEqRun#

class MCEq.core.MCEqRun(interaction_model, primary_model, theta_deg, **kwargs)[source]#

Bases: object

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) – interaction model name, e.g. SIBYLL2.3E
primary_model (class, param_tuple) – classes derived from crflux.models.PrimaryFlux and its parameters as tuple
theta_deg (float) – zenith angle \(\theta\) in degrees, measured positively from vertical direction
medium (string, optional) – “air”, “water”, “rock”, “co2”, “hydrogen”, “iron”
density_model (instance or tuple) – Instance of initialized density model or tuple of strings, such as (‘CORSIKA’, (‘BK_USStd’, None))
particle_list (list, optional) – Construct a system for only these partices including their decay products.

Attributes Summary

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

Methods Summary

`closest_energy`(kin_energy)	Convenience function to obtain the nearest grid energy to the energy argument, provided as kinetik energy in lab.
`decay_z_factor`(parent_pdg, child_pdg)	Energy dependent Z-factor according to Lipari (1993).
`etot_grid`(particle_name[, return_bins])	Computes and returns the total energy grid.
`get_solution`(particle_name[, mag, grid_idx, ...])	Retrieves solution of the calculation on the energy grid.
`inject_ddm`(ddm)	Set a DDM object as interaction model.
`n_e`([grid_idx, min_energy_cutoff])	Returns the number of electrons plus positrons at a grid step above min_energy_cutoff.
`n_mu`([grid_idx, min_energy_cutoff])	Returns the number of positive and negative muons at a grid step above min_energy_cutoff.
`n_particles`(label[, grid_idx, min_energy_cutoff])	Returns number of particles of type label at a grid step above an energy threshold for counting.
`ptot_grid`(particle_name[, return_bins])	Computes and returns the total momentum grid.
`regenerate_matrices`([skip_decay_matrix])	Call this function after applying particle prod.
`set_density_model`(density_model_or_config)	Sets model of the atmosphere.
`set_initial_spectrum`(spectrum, pdg_id[, append])	Set a user-defined spectrum for an arbitrary species as initial condition.
`set_interaction_model`(interaction_model[, ...])	Sets interaction model and/or an external charm model for calculation.
`set_mod_pprod`(prim_pdg, sec_pdg, x_func, ...)	Sets combination of projectile/secondary for error propagation.
`set_primary_model`(model_class_or_object[, tag])	Sets primary flux model.
`set_single_primary_particle`(E[, corsika_id, ...])	Set type and kinetic energy of a single primary nucleus to calculation of particle yields.
`set_theta_deg`(theta_deg)	Sets zenith angle \(\theta\) as seen from a detector.
`solve`([int_grid, grid_var])	Launches the solver.
`solve_from_integration_path`(nsteps, dX, ...)	Launches the solver directly for parameters of the integration path.
`unset_mod_pprod`([dont_fill])	Removes modifications from `MCEqRun.set_mod_pprod()`.
`xgrid`(particle_name, return_as[, return_bins])	Uniform access to the spectrum variable, depending on the same return_as argument as in get_solution.
`z_factor`(projectile_pdg, secondary_pdg[, ...])	Energy dependent Z-factor according to Thunman et al. (1996).

Attributes Documentation

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)

Methods Documentation

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

etot_grid(particle_name, return_bins=False)[source]#

Computes and returns the total energy grid.

If return_bins = True return bins and centers, otherwise just the bin centers.

get_solution(particle_name, mag=0.0, grid_idx=None, integrate=False, return_as='kinetic energy', dont_sum_helicities=False)[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 without a prefix mu+ or with the prefix total_mu+, total_numu
the conventional flux of muons, muon neutrinos etc. from all sources can be retrieved by the prefix conv_, i.e. conv_numu
the prompt flux of muons, muon neutrinos etc. from all sources can be retrieved by the prefix pr_, i.e. pr_numu
correspondigly, the flux of leptons which originated from the decay of a charged pion carries the prefix pi_ and from a kaon k_

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 \(= \Phi \cdot E^{mag}\)
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)
return_as (str, optional) – the flux can be returned as total energy, kinetic energy, or total momentum flux. This defaults to kinetic energy and is in general taken from MCEq.config.return_as
dont_sum_helicities (bool, optional) – Per default the lepton flux is summed over the available helicities, e.g. total_mu+ is the muon flux from (-1, 0, +1) helicity for mu+.

Returns:

func: numpy.array): flux of particles on energy grid e_grid

Return type:

(

inject_ddm(ddm)[source]#

Set a DDM object as interaction model.

The argument requires a DDM model object. Calling set_interaction_model overwrites DDM with a different model.

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:

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)[source]#

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)[source]#

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

ptot_grid(particle_name, return_bins=False)[source]#

Computes and returns the total momentum grid.

If return_bins True, return bins, centers, otherwise just the bin centers.

regenerate_matrices(skip_decay_matrix=False)[source]#: Call this function after applying particle prod. modifications aka Barr parameters

set_density_model(density_model_or_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_model_or_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_model_or_config (obj or tuple of strings) – (parametrization type, arguments)

set_initial_spectrum(spectrum, pdg_id, 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:

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)[source]#

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(model_class_or_object, tag=None)[source]#

Sets primary flux model.

This functions is quick and does not require re-generation of matrices.

Parameters:

interaction_model (CRFluxModel.PrimaryFlux) – reference
class (to primary model)
tag (tuple) – positional argument list for model class

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 \(E_{nucleon}= E_{nucleus} / A\) The nucleus type is defined via \(\text{CORSIKA ID} = A*100 + Z\). 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 \(50*A~ \text{GeV} < E_\text{nucleus} < 10^{10}*A \text{GeV}\).

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)[source]#

Sets zenith angle \(\theta\) 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:

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.

solve_from_integration_path(nsteps, dX, rho_inv, grid_idcs)[source]#

Launches the solver directly for parameters of the integration path.

The helper function is useful if you want to skip the calculation of the integration path every time. This function is intended for expert use and is not required for normal operation.

The parameters can be obtained after calling _calculate_integration_path with correct settings for density and angle parameters:

nsteps, dX, rho_inv, grid_idcs = self.integration_path

Parameters:

phi0 (np.array) – initial condition
nsteps (int) – number of integration steps
dX (list) – the delta_X’s
rho_inv (list) – the inverse of the density at each step
grid_idcs (list) – list of steps at which the solution
grid_sol (is dumped into)

unset_mod_pprod(dont_fill=False)[source]#

Removes modifications from MCEqRun.set_mod_pprod().

Parameters:

skip_fill (bool) – If true do not regenerate matrices
hand) ((has to be done at a later step by)

xgrid(particle_name, return_as, return_bins=False)[source]#: Uniform access to the spectrum variable, depending on the same return_as argument as in get_solution.

z_factor(projectile_pdg, secondary_pdg, definition='primary_e', min_energy=0.3, use_cs_scaling=True)[source]#: Energy dependent Z-factor according to Thunman et al. (1996)