Physical models¶

The cascade equation can be applied whenever particles interact and decay. The main purpose of the current program version is to calculate inclusive fluxes of leptons ( $\mu,\ \nu_{\mu},\ \nu_e\ \text{and}\ \nu_{\tau}$ ) in the Earth’s atmosphere. The calculation requires a spherical model of the Earth’s atmosphere which on hand is based on its MCEq.geometry and on the other hand parameterizations or numerical models of the MCEq.density_profiles.

Alternatiely, this code could be used for calculations of the high energy hadron and lepton flux in astrophysical environments, where the gas density and the (re-)interaction probability are very low. Prediction of detailed photon spectra is possible but additional extensions.

At very high energies, i.e. beyond 100 TeV, the decays of very short-lived particles become an important contribution to the flux. However, heavy-flavor production at these energies is not well known and it can not be consitently predicted within theoretical frameworks. Some of the default interaction models (SIBYLL-2.3 or SIBYLL-2.3c) contain models of charmed particle production but they represent only hypotheses based on experimental data and phenomenology. Other MCEq.charm_models exist, such as the MCEq.charm_models.MRS_charm, which can be coupled with any other interaction model for normal hadron production. In practice any kind of model which predicts a $x$ distribution can be employed in this code as extension of the MCEq.charm_models module.

`MCEq.density_profiles` - models of the Earth’s atmosphere¶

This module includes classes and functions modeling the Earth’s atmosphere. Currently, two different types models are supported:

Linsley-type/CORSIKA-style parameterization
Numerical atmosphere via external routine (NRLMSISE-00)

Both implementations have to inherit from the abstract class EarthAtmosphere, which provides the functions for other parts of the program. In particular the function EarthAtmosphere.get_density()

Typical interaction:

$ atm_object = CorsikaAtmosphere("BK_USStd")
$ atm_object.set_theta(90)
$ print 'density at X=100', atm_object.X2rho(100.)

The class MCEqRun will only the following routines::

EarthAtmosphere.set_theta(),
EarthAtmosphere.r_X2rho().

If you are extending this module make sure to provide these functions without breaking compatibility.

Example

An example can be run by executing the module:

$ python MCEq/atmospheres.py

class MCEq.density_profiles.AIRSAtmosphere(location, season, extrapolate=True, *args, **kwargs)[source]¶

Interpolation class for tabulated atmospheres.

This class is intended to read preprocessed AIRS Satellite data.

Parameters:	location (str) – see `init_parameters()` season (str,optional) – see `init_parameters()`

get_density(h_cm)[source]¶

Returns the density of air in g/cm**3.

Wraps around ctypes calls to the NRLMSISE-00 C library.

Parameters:	h_cm (float) – height in cm
Returns:	density $\rho(h_{cm})$ in g/cm**3
Return type:	float

init_parameters(location, **kwargs)[source]¶

Loads tables and prepares interpolation.

Parameters:	location (str) – supported is only “SouthPole” doy (int) – Day Of Year

class MCEq.density_profiles.CorsikaAtmosphere(location, season=None)[source]¶

Class, holding the parameters of a Linsley type parameterization similar to the Air-Shower Monte Carlo CORSIKA.

The parameters pre-defined parameters are taken from the CORSIKA manual. If new sets of parameters are added to init_parameters(), the array _thickl can be calculated using calc_thickl() .

_atm_param¶: numpy.array – (5x5) Stores 5 atmospheric parameters _aatm, _batm, _catm, _thickl, _hlay for each of the 5 layers

Parameters:	location (str) – see `init_parameters()` season (str,optional) – see `init_parameters()`

calc_thickl()[source]¶

Calculates thickness layers for depth2height()

The analytical inversion of the CORSIKA parameterization relies on the knowledge about the depth $X$ , where trasitions between layers/exponentials occur.

Example

Create a new set of parameters in init_parameters() inserting arbitrary values in the _thikl array:

$ cor_atm = CorsikaAtmosphere(new_location, new_season)
$ cor_atm.calc_thickl()

Replace _thickl values with printout.

depth2height(x_v)[source]¶

Converts column/vertical depth to height.

Parameters:	x_v (float) – column depth $X_v$ in g/cm**2
Returns:	height in cm
Return type:	float

get_density(h_cm)[source]¶

Returns the density of air in g/cm**3.

Uses the optimized module function corsika_get_density_jit().

Parameters:	h_cm (float) – height in cm
Returns:	density $\rho(h_{cm})$ in g/cm**3
Return type:	float

get_mass_overburden(h_cm)[source]¶

Returns the mass overburden in atmosphere in g/cm**2.

Uses the optimized module function corsika_get_m_overburden_jit().

Parameters:	h_cm (float) – height in cm
Returns:	column depth $T(h_{cm})$ in g/cm**2
Return type:	float

init_parameters(location, season)[source]¶

Initializes _atm_param.

location	CORSIKA Table	Description/season
“USStd”	1	US Standard atmosphere
“BK_USStd”	31	Bianca Keilhauer’s USStd
“Karlsruhe”	18	AT115 / Karlsruhe
“SouthPole”	26 and 28	MSIS-90-E for Dec and June
“PL_SouthPole”	29 and 30	Lipari’s Jan and Aug

Parameters:	location (str) – see table season (str, optional) – choice of season for supported locations
Raises:	`Exception` – if parameter set not available

rho_inv(X, cos_theta)[source]¶

Returns reciprocal density in cm**3/g using planar approximation.

This function uses the optimized function planar_rho_inv_jit()

Parameters:	h_cm (float) – height in cm
Returns:	$\frac{1}{\rho}(X,\cos{\theta})$ cm**3/g
Return type:	float

class MCEq.density_profiles.EarthAtmosphere(*args, **kwargs)[source]¶

Abstract class containing common methods on atmosphere. You have to inherit from this class and implement the virtual method get_density().

Note

Do not instantiate this class directly.

thrad¶: float – current zenith angle $\theta$ in radiants

theta_deg¶: float – current zenith angle $\theta$ in degrees

max_X¶: float – Slant depth at the surface according to the geometry defined in the MCEq.geometry

X2rho(X)[source]¶

Returns the density $\rho(X)$ .

The spline s_X2rho is used, which was calculated or retrieved from cache during the set_theta() call.

Parameters:	X (float) – slant depth in g/cm**2
Returns:	$\rho$ in cm**3/g
Return type:	float

calculate_density_spline(n_steps=2000)[source]¶

Calculates and stores a spline of $\rho(X)$ .

Parameters:	n_steps (int, optional) – number of $X$ values to use for interpolation
Raises:	`Exception` – if `set_theta()` was not called before.

gamma_cherenkov_air(h_cm)[source]¶: Returns the Lorentz factor gamma of Cherenkov threshold in air (MeV).

get_density(h_cm)[source]¶

Abstract method which implementation should return the density in g/cm**3.

Parameters:	h_cm (float) – height in cm
Returns:	density in g/cm**3
Return type:	float
Raises:	`NotImplementedError`

h2X(h)[source]¶

Returns the depth along path as function of height above surface.

The spline s_X2rho is used, which was calculated or retrieved from cache during the set_theta() call.

Parameters:	h (float) – vertical height above surface in cm
Returns:	X slant depth in g/cm**2
Return type:	float

moliere_air(h_cm)[source]¶: Returns the Moliere unit of air for US standard atmosphere.

nref_rel_air(h_cm)[source]¶: Returns the refractive index - 1 in air (density parametrization as in CORSIKA).

r_X2rho(X)[source]¶

Returns the inverse density $\frac{1}{\rho}(X)$ .

The spline s_X2rho is used, which was calculated or retrieved from cache during the set_theta() call.

Parameters:	X (float) – slant depth in g/cm**2
Returns:	$1/\rho$ in cm**3/g
Return type:	float

set_theta(theta_deg, force_spline_calc=False)[source]¶

Configures geometry and initiates spline calculation for $\rho(X)$ .

If the option ‘use_atm_cache’ is enabled in the config, the function will check, if a corresponding spline is available in the cache and use it. Otherwise it will call calculate_density_spline(), make the function r_X2rho() available to the core code and store the spline in the cache.

Parameters:	theta_deg (float) – zenith angle $\theta$ at detector force_spline_calc (bool) – forces (re-)calculation of the spline for each call

theta_cherenkov_air(h_cm)[source]¶: Returns the Cherenkov angle in air (degrees).

class MCEq.density_profiles.IsothermalAtmosphere(location, season, hiso_km=6.3, X0=1300.0)[source]¶

Isothermal model of the atmosphere.

This model is widely used in semi-analytical calculations. The isothermal approximation is valid in a certain range of altitudes and usually one adjust the parameters to match a more realistic density profile at altitudes between 10 - 30 km, where the high energy muon production rate peaks. Such parametrizations are given in the book “Cosmic Rays and Particle Physics”, Gaisser, Engel and Resconi (2016). The default values are from M. Thunman, G. Ingelman, and P. Gondolo, Astropart. Physics 5, 309 (1996).

Parameters:	location (str) – no effect season (str) – no effect hiso_km (float) – isothermal scale height in km X0 (float) – Ground level overburden

get_density(h_cm)[source]¶

Returns the density of air in g/cm**3.

Parameters:	h_cm (float) – height in cm
Returns:	density $\rho(h_{cm})$ in g/cm**3
Return type:	float

get_mass_overburden(h_cm)[source]¶

Returns the mass overburden in atmosphere in g/cm**2.

Parameters:	h_cm (float) – height in cm
Returns:	column depth $T(h_{cm})$ in g/cm**2
Return type:	float

class MCEq.density_profiles.MSIS00Atmosphere(location, season)[source]¶

Wrapper class for a python interface to the NRLMSISE-00 model.

NRLMSISE-00 is an empirical model of the Earth’s atmosphere. It is available as a FORTRAN 77 code or as a verson traslated into C by Dominik Borodowski. Here a PYTHON wrapper has been used.

_msis¶: NRLMSISE-00 python wrapper object handler

Parameters:	location (str) – see `init_parameters()` season (str,optional) – see `init_parameters()`

get_density(h_cm)[source]¶

Returns the density of air in g/cm**3.

Wraps around ctypes calls to the NRLMSISE-00 C library.

Parameters:	h_cm (float) – height in cm
Returns:	density $\rho(h_{cm})$ in g/cm**3
Return type:	float

init_parameters(location, season)[source]¶

Sets location and season in NRLMSISE-00.

Translates location and season into day of year and geo coordinates.

Parameters:	location (str) – Supported are “SouthPole” and “Karlsruhe” season (str) – months of the year: January, February, etc.

class MCEq.density_profiles.MSIS00IceCubeCentered(location, season)[source]¶

Extension of MSIS00Atmosphere which couples the latitude setting with the zenith angle of the detector.

Parameters:	location (str) – see `init_parameters()` season (str,optional) – see `init_parameters()`

latitude(det_zenith_deg)[source]¶

Returns the geographic latitude of the shower impact point.

Assumes a spherical earth. The detector is 1948m under the surface.

Credits: geometry fomulae by Jakob van Santen, DESY Zeuthen.

Parameters:	det_zenith_deg (float) – zenith angle at detector in degrees
Returns:	latitude of the impact point in degrees
Return type:	float

`MCEq.geometry` — Extensive-Air-Shower geometry¶

This module includes classes and functions modeling the geometry of the cascade trajectory.

class MCEq.geometry.EarthGeometry[source]¶

A model of the Earth’s geometry, approximating it: by a sphere. The figure below illustrates the meaning of the parameters.

Curved geometry as it is used in the code (not to scale!).

Example

The plots below will be produced by executing the module:

$ python geometry.py

(Source code)

(png, hires.png, pdf)

(png, hires.png, pdf)

(png, hires.png, pdf)

(png, hires.png, pdf)

The module currently contains only module level functions. To implement a different geometry, e.g. in an astrophysical content, you can create a new module providing similar functions or place the current content into a class.

h_obs¶: float – observation level height [cm]

h_atm¶: float – top of the atmosphere [cm]

r_E¶: float – radius Earth [cm]

r_top¶: float – radius at top of the atmosphere [cm]

r_obs¶: float – radius at observation level [cm]

cos_th_star(theta)[source]¶: Returns the zenith angle at atmospheric boarder $\cos(\theta^*)$ in [rad] as a function of zenith at detector.

delta_l(h, theta)[source]¶: Distance $dl$ covered along path $l(\theta)$ as a function of current height. Inverse to h().

h(dl, theta)[source]¶: Height above surface at distance $dl$ counted from the beginning of path $l(\theta)$ in cm.

l(theta)[source]¶: Returns path length in [cm] for given zenith angle $\theta$ [rad].

MCEq.geometry.chirkin_cos_theta_star(costheta)[source]¶

$\cos(\theta^*)$ parameterization.

This function returns the equivalent zenith angle for for very inclined showers. It is based on a CORSIKA study by D. Chirkin, hep-ph/0407078v1, 2004.

Parameters:	costheta (float) – $\cos(\theta)$ in [rad]
Returns:	$\cos(\theta*)$ in [rad]
Return type:	float

`MCEq.charm_models` — charmed particle production¶

This module includes classes for custom charmed particle production. Currently only the MRS model is implemented as the class MRS_charm. The abstract class CharmModel guides the implementation of custom classes.

The Yields instantiates derived classes of CharmModel and calls CharmModel.get_yield_matrix() when overwriting a model yield file in Yields.set_custom_charm_model().

class MCEq.charm_models.CharmModel[source]¶

Abstract class, from which implemeted charm models can inherit.

Note

Do not instantiate this class directly.

get_yield_matrix(proj, sec)[source]¶

The implementation of this abstract method returns the yield matrix spanned over the energy grid of the calculation.

Parameters:	proj (int) – PDG ID of the interacting particle (projectile) sec (int) – PDG ID of the final state charmed meson (secondary)
Returns:	yield matrix
Return type:	np.array
Raises:	`NotImplementedError`

class MCEq.charm_models.MRS_charm(e_grid, csm)[source]¶

Martin-Ryskin-Stasto charm model.

The model is described in A. D. Martin, M. G. Ryskin, and A. M. Stasto, Acta Physica Polonica B 34, 3273 (2003). The parameterization of the inclusive $c\bar{c}$ cross-section is given in the appendix of the paper. This formula provides the behavior of the cross-section, while fragmentation functions and certain scales are needed to obtain meson and baryon fluxes as a function of the kinematic variable $x_F$ . At high energies and $x_F > 0.05$ , where this model is valid, $x_F \approx x=E_c/E_{proj}$ . Here, these fragmentation functions are used:

$D$ -mesons $\frac{4}{3} x$

$\Lambda$ -baryons $\frac{1}{1.47} x$

The production ratios between the different types of $D$ -mesons are stored in the attribute cs_scales and D0_scale, where D0_scale is the $c\bar{c}$ to $D^0$ ratio and cs_scales stores the production ratios of $D^\pm/D^0$ , $D_s/D^0$ and $\Lambda_c/D^0$ .

Since the model employs only perturbartive production of charm, the charge conjugates are symmetric, i.e. $\sigma_{D^+} = \sigma_{D^-}$ etc.

Parameters:	e_grid (np.array) – energy grid as it is defined in `MCEqRun`. csm (np.array) – inelastic cross-sections as used in `MCEqRun`.

D_dist(x, E, mes)[source]¶

Returns the Feynman- $x_F$ distribution of $\sigma_{D-mesons}$ in mb

Parameters:	x (float or np.array) – $x_F$ E (float) – center-of-mass energy in GeV mes (int) – PDG ID of D-meson: $\pm421, \pm431, \pm411$
Returns:	$\sigma_{D-mesons}$ in mb
Return type:	float

LambdaC_dist(x, E)[source]¶

Returns the Feynman- $x_F$ distribution of $\sigma_{\Lambda_C}$ in mb

Parameters:	x (float or np.array) – $x_F$ E (float) – center-of-mass energy in GeV mes (int) – PDG ID of D-meson: $\pm421, \pm431, \pm411$
Returns:	$\sigma_{D-mesons}$ in mb
Return type:	float

dsig_dx(x, E)[source]¶

Returns the Feynman- $x_F$ distribution of $\sigma_{c\bar{c}}$ in mb

Parameters:	x (float or np.array) – $x_F$ E (float) – center-of-mass energy in GeV
Returns:	$\sigma_{c\bar{c}}$ in mb
Return type:	float

get_yield_matrix(proj, sec)[source]¶

Returns the yield matrix in proper format for MCEqRun.

Parameters:

proj (int) – projectile PDG ID $\pm$ [2212, 211, 321]
sec (int) – charmed particle PDG ID $\pm$ [411, 421, 431, 4122]

Returns:

yield matrix if (proj,sec) combination allowed,: else zero matrix

Return type:

np.array

sigma_cc(E)[source]¶: Returns the integrated ccbar cross-section in mb.

Note

Integration is not going over complete phase-space due to limitations of the parameterization.

test()[source]¶: Plots the meson, baryon and charm quark distribution as shown in the plot below.

D0_scale = 0.47619047619047616¶: D0 cross-section wrt to the ccbar cross-section

allowed_proj = [2212, -2212, 2112, -2112, 211, -211, 321, -321]¶: hadron projectiles, which are allowed to produce charm

allowed_sec = [411, 421, 431, 4122]¶: charm secondaries, which are predicted by this model

cs_scales = {4122: 0.45, 411: 0.5, 421: 1.0, 431: 0.15}¶: fractions of cross-section wrt to D0 cross-section

class MCEq.charm_models.WHR_charm(e_grid, csm)[source]¶

Logan Wille, Francis Halzen, Hall Reno.

The approach is the same as in A. D. Martin, M. G. Ryskin, and A. M. Stasto, Acta Physica Polonica B 34, 3273 (2003). The parameterization of the inclusive $c\bar{c}$ cross-section is replaced by interpolated tables from the calculation. Fragmentation functions and certain scales are needed to obtain meson and baryon fluxes as a function of the kinematic variable $x_F$ . At high energies and $x_F > 0.05$ , where this model is valid, $x_F \approx x=E_c/E_{proj}$ . Here, these fragmentation functions are used:

$D$ -mesons $\frac{4}{3} x$

$\Lambda$ -baryons $\frac{1}{1.47} x$

The production ratios between the different types of $D$ -mesons are stored in the attribute cs_scales and D0_scale, where D0_scale is the $c\bar{c}$ to $D^0$ ratio and cs_scales stores the production ratios of $D^\pm/D^0$ , $D_s/D^0$ and $\Lambda_c/D^0$ .

Since the model employs only perturbartive production of charm, the charge conjugates are symmetric, i.e. $\sigma_{D^+} = \sigma_{D^-}$ etc.

Parameters:	e_grid (np.array) – energy grid as it is defined in `MCEqRun`. csm (np.array) – inelastic cross-sections as used in `MCEqRun`.

dsig_dx(x, E)[source]¶

Returns the Feynman- $x_F$ distribution of $\sigma_{c\bar{c}}$ in mb

Parameters:	x (float or np.array) – $x_F$ E (float) – center-of-mass energy in GeV
Returns:	$\sigma_{c\bar{c}}$ in mb
Return type:	float

Physical models¶

MCEq.density_profiles - models of the Earth’s atmosphere¶

MCEq.geometry — Extensive-Air-Shower geometry¶

MCEq.charm_models — charmed particle production¶

`MCEq.density_profiles` - models of the Earth’s atmosphere¶

`MCEq.geometry` — Extensive-Air-Shower geometry¶

`MCEq.charm_models` — charmed particle production¶