Geometry package¶

In MCEq, geometry is everything related to the medium in which the particle cascade develops. The very basic geometrical functions for the polar coordinate system of the Earth - no it’s not flat, but just azimuth symmetric - are located in MCEq.geometry.geometry. The density parameterizations and interfaces are in MCEq.geometry.density_profiles

`MCEq.geometry.density_profiles`¶

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 EarthsAtmosphere, which provides the functions for other parts of the program. In particular the function EarthsAtmosphere.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::

EarthsAtmosphere.set_theta(),
EarthsAtmosphere.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.geometry.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.

Interpolates table at requested value for previously set year and day of year (doy).

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

get_temperature(h_cm)[source]¶

Returns the temperature in K.

Interpolates table at requested value for previously set year and day of year (doy).

Parameters:	h_cm (float) – height in cm
Returns:	temperature $T(h_{cm})$ in K
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.geometry.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¶

(5x5) Stores 5 atmospheric parameters _aatm, _batm, _catm, _thickl, _hlay for each of the 5 layers

Type:	numpy.array

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. Parameters from ANTARES/KM3NET are based on the work of T. Heid (see this issue)

location	CORSIKA Table	Description/season
“USStd”	23	US Standard atmosphere
“BK_USStd”	37	Bianca Keilhauer’s USStd
“Karlsruhe”	24	AT115 / Karlsruhe
“SouthPole”	26 and 28	MSIS-90-E for Dec and June
“PL_SouthPole”	29 and 30	Lipari’s Jan and Aug
“ANTARES/KM3NeT-ORCA”	NA	PhD T. Heid
“KM3NeT-ARCA”	NA	PhD T. Heid

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.geometry.density_profiles.EarthsAtmosphere(*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¶

current zenith angle $\theta$ in radiants

Type:	float

theta_deg¶

current zenith angle $\theta$ in degrees

Type:	float

max_X¶

Slant depth at the surface according to the geometry defined in the MCEq.geometry

Type:	float

geometry¶

Can be a custom instance of EarthGeometry

Type:	object

X2h(X)[source]¶

Returns the height above surface as a function of slant depth for currently selected zenith angle.

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

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

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

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

max_X: Depth at altitude 0.

max_den¶: Density at altitude 0.

s_X2rho¶: Spline for conversion from depth to density.

s_h2X¶: Spline for conversion from altitude to depth.

s_lX2h¶: Spline for conversion from depth to altitude.

class MCEq.geometry.density_profiles.GeneralizedTarget(len_target=100000.0, env_density=0.001225, env_name='air')[source]¶

This class provides a way to run MCEq on piece-wise constant one-dimenional density profiles.

The default values for the average density are taken from config file variables len_target, env_density and env_name. The density profile has to be built by calling subsequently add_material(). The current composition of the target can be checked with draw_materials() or print_table().

Note

If the target is not air or hydrogen, the result is approximate, since seconray particle yields are provided for nucleon-air or proton-proton collisions. Depending on this choice one has to adjust the nuclear mass in mceq_config.

Parameters:	len_target (float) – total length of the target in meters env_density (float) – density of the default material in g/cm3 env_name** (str) – title for this environment

add_material(start_position_cm, density, name)[source]¶

Adds one additional material to a composite target.

Parameters:	start_position_cm (float) – position where the material starts counted from target origin l\|X = 0 in cm density (float) – density of material in g/cm3 name** (str) – any user defined name
Raises:	`Exception` – If requested start_position_cm is not properly defined.

draw_materials(axes=None, logx=False)[source]¶

Makes a plot of depth and density profile as a function of the target length. The list of materials is printed out, too.

Parameters:	axes (plt.axes, optional) – handle for matplotlib axes

get_density(l_cm)[source]¶

Returns the density in g/cm**3 as a function of position l in cm.

Parameters:	l (float) – position in target in cm
Returns:	density in g/cm**3
Return type:	float
Raises:	`Exception` – If requested position exceeds target length.

get_density_X(X)[source]¶

Returns the density in g/cm**3 as a function of depth X.

Parameters:	X (float) – depth in g/cm**2
Returns:	density in g/cm**3
Return type:	float
Raises:	`Exception` – If requested depth exceeds target.

print_table(min_dbg_lev=0)[source]¶: Prints table of materials to standard output.

r_X2rho(X)[source]¶

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

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

reset()[source]¶: Resets material list to defaults.

set_length(new_length_cm)[source]¶

Updates the total length of the target.

Usually the length is set

set_theta(*args)[source]¶

This method is not defined for the generalized target. The purpose is to catch usage errors.

Raises:	`NotImplementedError` – always

max_X¶: Maximal depth of target.

s_X2h¶: Spline for depth at distance.

s_h2X¶: Spline for distance at depth.

class MCEq.geometry.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.geometry.density_profiles.MSIS00Atmosphere(location, season=None, doy=None, use_loc_altitudes=False)[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

get_temperature(h_cm)[source]¶

Returns the temperature of air in K.

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

Parameters:	h_cm (float) – height in cm
Returns:	density $T(h_{cm})$ in K
Return type:	float

init_parameters(location, season, doy, use_loc_altitudes)[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. use_loc_altitudes (bool) – If to use default altitudes from location

set_doy(day_of_year)[source]¶

Changes MSIS season by day of year.

Parameters:	day_of_year (int) – Jan.=0, 1.Feb=32

set_location(location)[source]¶

Changes MSIS location by strings defined in _msis_wrapper.

Parameters:	location (str) – location as defined in `NRLMSISE-00.`

set_location_coord(longitude, latitude)[source]¶

Changes MSIS location by longitude, latitude in _msis_wrapper

Parameters:	longitude (float) – longitude of the location with abs(longitude) <= 180 latitude (float) – latitude of the location with abs(latitude) <= 90

set_season(month)[source]¶

Changes MSIS location by month strings defined in _msis_wrapper.

Parameters:	location (str) – month as defined in `NRLMSISE-00.`

update_parameters(**kwargs)[source]¶

Updates parameters of the density model

Parameters:	location_coord (tuple of str) – (longitude, latitude) season (str) – months of the year: January, February, etc. doy (int) – day of the year. ‘doy’ takes precedence over ‘season’ if both are set

class MCEq.geometry.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

set_theta(theta_deg)[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

`MCEq.geometry.geometry`¶

The module contains the geometry for an azimuth symmetric Earth.

class MCEq.geometry.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)

h_obs¶

observation level height [cm]

Type:	float

h_atm¶

top of the atmosphere [cm]

Type:	float

r_E¶

radius Earth [cm]

Type:	float

r_top¶

radius at top of the atmosphere [cm]

Type:	float

r_obs¶

radius at observation level [cm]

Type:	float

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.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.geometry.nrlmsise00`¶

CTypes interface to the C-version of the NRLMSISE-00 code, originally developed by Picone et al.. The C-translation is by Dominik Brodowski <https://www.brodo.de/space/nrlmsise/index.html>_.

`MCEq.geometry.corsikaatm`¶

This set of functions are C implementations of the piecewise defined exponential profiles as used in CORSIKA. An efficient implementation is difficult in plain numpy.

Geometry package¶

MCEq.geometry.density_profiles¶

MCEq.geometry.geometry¶

MCEq.geometry.nrlmsise00¶

MCEq.geometry.corsikaatm¶

`MCEq.geometry.density_profiles`¶

`MCEq.geometry.geometry`¶

`MCEq.geometry.nrlmsise00`¶

`MCEq.geometry.corsikaatm`¶