Simulations

The simulation module is a tool to setup arrival simulations in a few lines of code. It is a wrapper for the core functions and is based on the data container provided by the cosmic_rays module. The simulation module is a parametrized simulation based on healpy, a tool to pixelize the sphere into pixels with equally solid angle (https://healpy.readthedocs.io/en/latest/index.html). For examples have a look at the tutorial!

Module to setup a parametrized simulation, that work on probability distributions

class simulations.BaseSimulation(nsets, ncrs)

Bases: object

Base class where the classes ObservedBound() and SourceBound() inherit from.

apply_uncertainties(err_e=None, err_a=None, err_xmax=None)

Apply measurement uncertainties.

Parameters

• err_e – relative uncertainty on the energy (typical 0.14)

• err_a – angular uncertainty in degree on the arrival directions (typical 0.5 degree)

• err_xmax – absolute uncertainty on the shower depth xmax (typical 15 g/cm^2)

set_xmax(z2a='double', model='EPOS-LHC')

Calculate Xmax bei gumbel distribution for the simulated energies and charges.

Parameters

• z2a – How the charge is converted to mass number [‘double’, ‘empiric’, ‘stable’, ‘abundance’]

• model – Hadronic interaction for gumbel distribution

Returns

no return

shuffle_events()

Independently shuffle the cosmic rays of each set.

class simulations.CompositionModel(shape, log10e=None)

Bases: object

Predefined composition models

auger(smoothed=True, model='EPOS-LHC')

Simple estimate from AUGER Xmax measurements

auger_exponential()

Simple exponential estimate from AUGER Xmax measurements

equal()

Assumes a equal distribution in (H, He, N, Fe) groups.

mixed()

Simple estimate of the composition above ~20 EeV by M. Erdmann (2017)

mixed_clipped()

mixed from above, but CNO group only Z=6 because of no lenses at low rigidities

class simulations.EnergySpectrum(shape, log10e_min, log10e_max=20.5)

Bases: object

Predefined energy spectra

auger_fit()

Energies following the AUGER spectrum above log10e_min 17.5.

power_law(gamma=- 3)

Power law spectrum, with spectral index corresponding to non differential spectrum, where gamma=-3.7 corresponds to the AUGER fit at intermediate energies.

Parameters

gamma – non-differential spectral index (E ~ E^(gamma))

Returns

energies in shape self.shape

class simulations.ObservedBound(nside, nsets, ncrs)

Bases: simulations.BaseSimulation

Class to simulate cosmic ray arrival scenario by sources located at the sky, including energies, charges, smearing and galactic magnetic field effects. This is an observed bound simulation, thus energies and composition is set on Earth and might differ at sources.

apply_exposure(a0=- 35.25, zmax=60)

Apply the exposure (coverage) of any experiment (default: AUGER) to the observed maps.

Parameters

• a0 (float, int) – equatorial declination [deg] of the experiment (default: AUGER, a0=-35.25 deg)

• zmax – maximum zenith angle [deg] for the events

Returns

no return

apply_uncertainties(err_e=None, err_a=None, err_xmax=None)

Apply measurement uncertainties (see BaseSimulation.apply_uncertainties()).

arrival_setup(fsig, seed=1)

Creates the realizations of the arrival maps.

Parameters

fsig (float) – signal fraction of cosmic rays per set (signal = originate from sources)

Returns

no return

convert_pixel(keyword='vecs', convert_all=False)

Converts pixelized information stored under key ‘pixel’ to vectors (keyword=’vecs’) or angles (keyword=’angles’), accessible via ‘lon’, ‘lat’. When convert_all is True, both are saved. This can be used at a later stage, if convert_all was set to False originally.

get_data(convert_all=False, shuffle=False)

Returns the data in the form of the cosmic_rays.CosmicRaysSets() container.

Parameters

convert_all (bool) – if True, also vectors and lon/lat of the cosmic ray sets are saved (more memory expensive)

Returns

Instance of cosmic_rays.CosmicRaysSets() classes

Example: sim = ObservedBound() … crs = sim.get_data(convert_all=True) pixel = crs[‘pixel’] lon = crs[‘lon’] lat = crs[‘lat’] log10e = crs[‘log10e’] charge = crs[‘charge’]

lensing_map(lens, cache=None)

Apply a galactic magnetic field to the extragalactic map.

Parameters

• lens – Instance of astrotools.gamale.Lens class, for the galactic magnetic field

• cache – Caches all the loaded lens parts (increases speed, but may consume a lot of memory!)

Returns

no return

set_charges(charge, **kwargs)

Setting the charges of the simulated cosmic ray set.

Parameters

charge – Either charge number that is used or numpy.array of charges in shape (nsets, ncrs) or keyword

Type

charge: Union[np.ndarray, str, float]

Returns

charges

set_energy(log10e_min, log10e_max=None, energy_spectrum=None, **kwargs)

Setting the energies of the simulated cosmic ray set.

Parameters

• log10e_min (Union[np.ndarray, float]) – Either minimum energy (in log10e) for AUGER setup or power-law or numpy.array of energies in shape (nsets, ncrs)

• log10e_max – Maximum energy for AUGER setup

• energy_spectrum – model that is defined in below class EnergySpectrum

• kwargs – additional named keywords which will be passed to class EnergySpectrum

Returns

energies in log10e

set_rigidity_bins(lens_or_bins, cover_rigidity=True)

Defines the bin centers of the rigidities.

Parameters

lens_or_bins – Either the binning of the lens (left bin borders) or the lens itself

Returns

no return

set_sources(sources, fluxes=None)

Define source position and optional weights (cosmic ray luminosity).

Parameters

• sources – array of shape (3, n_sources) that point towards the center of the sources or integer for number of random sources or keyword [‘sbg’]

• fluxes – corresponding cosmic ray fluxes of the sources of shape (n_sources).

Returns

no return

set_xmax(z2a='double', model='EPOS-LHC')

Calculate Xmax bei gumbel distribution for the simulated energies and charges. For more information refer to BaseSimulation.set_xmax().

smear_sources(delta, dynamic=None)

Smears the source positions with a fisher distribution of width delta (optional dynamic smearing).

Parameters

• delta – either the constant width of the fisher distribution or dynamic (delta / R[10 EV]), in radians

• dynamic – if True, function applies dynamic smearing (delta / R[EV]) with value delta at 10 EV rigidity

Returns

no return

class simulations.SourceBound(nsets, ncrs)

Bases: simulations.BaseSimulation

Class to simulate cosmic ray arrival scenario by sources located at the sky, including energies, charges, smearing. This is a source bound simulation, thus energies and composition is set at the sources and might differ at Earth.

attenuate(library_path=None, inside_fraction=None, evolution=0)

Apply attenuation for far away sources based on CRPropa simulations

Parameters

• library_path – Input library file to use.

• evolution – m for source evolution of the homogeneous background (right now: not for foreground sources)

• cosmological_expansion – if the cosmological expansion of the universe is taken into account. If this is True, all distances (e.g. rmax) have to be given in comoving distance.

get_data(shuffle=False)

Returns the data in the form of the cosmic_rays.CosmicRaysSets() container.

Returns

Instance of cosmic_rays.CosmicRaysSets() classes

plot_arrivals(idx=None, opath=None, emin=None)

Plot arrival map

plot_composition(idx=None, opath=None)

Plots composition (comparison measurements/simulation)

plot_composition_skymap(idx=None, opath=None, emin=None)

Plot arrival map

plot_distance(sets='all', opath=None, emin=None, sig='both')

Plot histogram of travel distances (either mean of all sets or one specific)

plot_spectrum(opath=None)

Plot energy spectrum

set_charges(charges)

Define fraction of charge groups in form of dictionary (e.g. {‘h’:0.5, ‘fe’:0.5}) at source or as keyword ‘first_minimum’/’second_minimum’ from Auger’s combined fit paper (arXiv:1612.07155) If string is given, gamma and Rcut are also set to the respective best fit values.

Parameters

charges – dictionary hosting the fractions of injected elements (H, He, N, Si, Fe possible) or string (‘first_minimum’/’second_minimum’)

set_energy(log10e_min, gamma=- 2, log10_cut=20.0, rig_cut=True)

Define energy spectrum and cut off energy of sources.

Parameters

• log10e_min – Minimum threshold energy for observation

• gamma – Spectral index of energy spectrum at sources

• log10_cut – Maximum cut-off energy or rigidity for sources

• rig_cut – if True, log10_cut refers to a rigidity cut

set_sources(source_density=None, sources=None, fluxes=None, n_src=None, rmax=None)

Define source density or directly positions and optional weights (cosmic ray luminosity).

Parameters

• source_density – source density (in 1 / Mpc^3) or array of shape (3, n_sources) that point towards the center of the sources or keyword ‘sbg’

• fluxes – corresponding cosmic ray fluxes of the sources of shape (n_sources).

• n_src – Number of point sources to be considered.

Returns

no return

set_xmax(z2a='double', model='EPOS-LHC')

Calculate Xmax by Gumbel distribution for the simulated energies and charges. For more information refer to BaseSimulation.set_xmax().

smear_sources(delta, dynamic=None)

Smears the source positions with a fisher distribution of width delta (optional dynamic smearing).

Parameters

• delta – either the constant width of the fisher distribution or dynamic (delta / R[10 EV]), in radians

• dynamic – if True, function applies dynamic smearing (delta / R[EV]) with value delta at 10 EV rigidity

Returns

no return

class simulations.SourceGeometry(nsets, cosmo=True)

Bases: object

Class to set up a 3d Universe out of isotropically distributed sources.

distance_indices(distance_bins)

Get indices of given distance bins in the shape of the sources (nsets, n_src).

Parameters

distance_bins – Distance bins which refer to indices

Return indices

indices of distance_bins in shape (nsets, n_src)

has_sources()

Check if geometry features sources

set_sources(source_density=None, sources=None, fluxes=None, n_src=None, rmax=None)

Define source density or directly positions of sources and optional weights (cosmic ray luminosity).

Parameters

• source_density – sets source density (in 1 / Mpc^3) as float or array with shape (self.nsets, ) can be keyword of source catalog like ‘sbg’

• sources – explicitly set source coordinates of shape (3, None) Here, distances are always lookback distances, never comoving!

• fluxes – corresponding cosmic ray fluxes of the sources of shape (n_sources).

• n_src – Number of point sources to be considered.

• rmax – distance where background sources start.

Returns

no return

class simulations.SourceScenario

Bases: object

Predefined source scenarios

static agn_3fhl()

catalog of gamma-AGN emitters from 3FHL catalog, downloaded 15.12.2020 from https://www.auger.unam.mx/AugerWiki/ArrivalDirections/HEfd

static cena()

Only Cen A as source https://www.auger.unam.mx/AugerWiki/ArrivalDirections/HEfd

static gamma_agn()

AGN scenario used in GAP note 2017_007

static sbg()

Starburst Galaxy Scenario used in GAP note 2017_007

static sbg_NGC1068()

Lunardini SBG catalog + Circinus (44 sources), without NGC 253 for cross checks

static sbg_NGC253()

Lunardini SBG catalog + Circinus (44 sources), without NGC 253 for cross checks

static sbg_lunardini()

Lunardini SBG catalog + Circinus (44 sources), from J.Biteau, received May 2019 per email

simulations.sample_from_m_distributions(weight_matrix, n)

Sample n events from m distributions each having k bins and defined in weight_matrix.

Parameters

• weight_matrix – shape (m, k), hosting m different distributions each with k bins.

• n – Number of events drawn for each of the m distributions.

Return indices

Indice array with values 0..(k-1), in shape (m, n)

simulations.set_fisher_smeared_sources(nside, sources, delta, source_fluxes=None)

Smears the source positions (optional fluxes) with a fisher distribution of width delta.

Parameters

• nside (int) – nside of the HEALPix pixelization (default: 64)

• sources – array of shape (3, n_sources) that point towards the center of the sources

• delta – float or array with same length as sources: width of the fisher distribution (in radians)

• source_fluxes – corresponding cosmic ray fluxes of the sources of shape (n_sources).

Returns

healpy map (with npix entries) for the smeared sources normalized to sum 1