Source code for NuRadioReco.utilities.cr_flux

import numpy as np
from NuRadioReco.utilities import units
import scipy.interpolate as interpolate
from scipy.integrate import quad
import os.path


# References
# data_auger_19: https://www.auger.org/document-centre-public?task=download.send&id=5045&catid=78&m=0
# data_TA_21: data provided by van Vliet
# analytic spectrum TA_19: https://pos.sissa.it/358/298/pdf
# analytic spectrum auger: https://git.rwth-aachen.de/astro/astrotools/-/blob/master/astrotools

[docs]def get_measured_data(type='data_auger_19'): """ Helper function get data points in NuRadio base units Parameters ---------- type: string choose between data_auger_19 and data_TA_21 Returns ------- E: array of floats energie in eV J: array of floats flux of cosmic ray Err_up: array of floats upper error bar on flux Err_low: array of floats lower error bar on flux """ if type == 'data_auger_19': abspath = os.path.dirname(os.path.abspath(__file__)) data = np.loadtxt(os.path.join(abspath, 'data/Auger_combined_spectrum_ICRC_2019.txt'), skiprows=3) # from PRL paper 2020 based on ICRC 2019 E = 10 ** (data[:, 0]) * units.eV E_J = data[:, 1] * units.m ** -2 * units.second ** -1 * units.steradian ** -1 J = E_J / E J_err_low = data[:, 3] * units.m ** -2 * units.second ** -1 * units.steradian ** -1 / E J_err_up = data[:, 2] * units.m ** -2 * units.second ** -1 * units.steradian ** -1 / E elif type == 'data_TA_21': abspath = os.path.dirname(os.path.abspath(__file__)) data = np.loadtxt(os.path.join(abspath, 'data/TA_combined_spectrum_ICRC_2021.txt'), skiprows=10) E = 10 ** (data[:, 0]) * units.eV J = data[:, 2] * units.eV ** -1 * units.m ** -2 * units.second ** -1 * units.steradian ** -1 J_band_low = data[:, 3] * units.eV ** -1 * units.m ** -2 * units.second ** -1 * units.steradian ** -1 J_band_up = data[:, 4] * units.eV ** -1 * units.m ** -2 * units.second ** -1 * units.steradian ** -1 J_err_low = J - J_band_low J_err_up = J_band_up - J else: raise NotImplementedError(f"Type {type} is not supported") return E, J, J_err_low, J_err_up
[docs]def get_interpolated_cr_flux(log10_energy, type='data_auger_19'): """ Evaluates the scipy interpolation of the measured data at a given energy Parameters ---------- log10_energy: float energies (in log10(E / eV))) type: string choose between data_auger_19 and data_TA_21 Return ------- scipy interpolation of data in NuRadio base units (which are 1/(eV m^2 sr ns) ) """ E, J, J_err_low, J_err_up = get_measured_data(type) log10_e = np.log10(E) flux_interpolation = interpolate.interp1d(log10_e, J, fill_value=0, bounds_error=True) return flux_interpolation(log10_energy)
[docs]def get_analytic_cr_flux(log10_energy, type="auger_19"): """ Parameters ---------- log10_energy: float energies (in log10(E / eV))) type: string choose between auger_17, auger_19 and TA_19 Returns ------- analytic parametrization of spectrum in NuRadio base units (which are 1/(eV m^2 sr ns) """ energy = 10 ** log10_energy if type == "auger_17": p = np.array([2.8e-19, 5.08e18, 39e18, 3.293, 2.53, 2.5]) spectrum = np.where(energy < p[1], p[0] * (energy / p[1]) ** (-p[3]), p[0] * (energy / p[1]) ** (-p[4]) * (1 + (p[1] / p[2]) ** p[5]) * (1 + (energy / p[2]) ** p[5]) ** -1) spectrum = spectrum * (units.eV * units.km ** 2 * units.sr * units.year) ** -1 elif type == "auger_19": p = np.array([3.46e12, 1.5e17, 6.2e18, 12e18, 50e18, 2.92, 3.27, 2.2, 3.2, 5.4]) spectrum = (energy / p[0]) ** (-p[5]) * \ (1 + (energy / p[1]) ** p[5]) / (1 + (energy / p[1]) ** p[6]) * \ (1 + (energy / p[2]) ** p[6]) / (1 + (energy / p[2]) ** p[7]) * \ (1 + (energy / p[3]) ** p[7]) / (1 + (energy / p[3]) ** p[8]) * \ (1 + (energy / p[4]) ** p[8]) / (1 + (energy / p[4]) ** p[9]) spectrum = spectrum * (units.eV * units.km ** 2 * units.sr * units.year) ** -1 elif type == "TA_19": p1 = -3.28 p2 = -2.68 p3 = -4.84 E1 = 10 ** 18.69 E2 = 10 ** 19.81 c = 2.24e-30 c1 = c * (E1 / 1e18) ** p1 c2 = c1 * (E2 / E1) ** p2 spectrum = np.where(energy < E1, c * (energy / 1e18) ** p1, np.where(energy < E2, c1 * (energy / E1) ** p2, c2 * (energy / E2) ** p3)) spectrum = spectrum * (units.eV * units.m ** 2 * units.sr * units.s) ** -1 else: raise NotImplementedError(f"Type {type} is not supported") return spectrum
[docs]def get_flux_per_energy_bin(log10e_min, log10e_max, type='auger_19'): """ Returns an scipy integration of the measured data or the analytic spectrum over given interval. Parameters ---------- log10e_min: float Input energy lower bound (in log10(E / eV))) log10e_max: float Input energy upper bound (in log10(E / eV))) type: string choose between data_auger_19, data_TA_21, auger_17, auger_19, TA_19 Returns ------- scipy integration of data in NuRadio base units (which are 1/(eV m^2 sr ns)) """ if type in ['auger_17', 'auger_19', 'TA_19']: def flux(x): """ Bring parametrized energy spectrum in right shape for quad() function """ return get_analytic_cr_flux(np.log10(np.array([x])), type)[0] integrated_flux = quad(flux, 10 ** log10e_min, 10 ** log10e_max) else: raise NotImplementedError(f"Type {type} is not supported") return integrated_flux[0]
[docs]def get_cr_event_rate(log10energy=18, zenith=50*units.deg, a_eff=1, type="auger_19"): """ Cosmic ray event rate at a specific energy and zenith angle assuming a detector with effective area 'A_eff'. The detector projection and solid range are taken into account. The flux is calculated with the analytic spectrum. Parameters ---------- log10energy: float energy in units log10(energy/eV) zenith: float zenith angle a_eff: effective area of detector type: string choose between auger_17, auger_19 and TA_19 Returns ------- eventrate for an isotropic flux with given energy and zenith angle. The differential flux is returned in d/d zenith_angle, not d/d solid_angle. """ # The projected area, is the area visible to CR. # None of the horizontal CR reach the flat detector, all CR from above reach the detector projected_area = np.cos(zenith) solid_angle = 2 * np.pi * np.sin(zenith) flux = get_analytic_cr_flux(log10energy, type=type) return flux * projected_area * solid_angle * a_eff
[docs]def plot_measured_spectrum(ax=None, scale=2.7, type='data_auger_19', base_units=False): """ Plot measured spectrum. Attention: time unit is year instead of ns. Parameters ---------- ax: axis on which the data is plotted scale: scale factor for energy, default = 2.7 type: string choose between data_auger_19 and data_TA_21 base_units: bool if False spectrum will be plotted in [m^{-2} yr^{-1} sr^{-1} eV^{%scale-1}] instead of NuRadio base units Returns ------- plot of data without plt.show() """ import matplotlib.pyplot as plt if ax is None: fig = plt.figure() ax = fig.add_subplot(111) E, J, J_err_low, J_err_up = get_measured_data(type) if not base_units: E_J_scale_units = units.m ** -2 * units.year ** -1 * units.steradian ** -1 * units.eV ** (scale - 1) yl = r'$J(E)$ [m$^{-2}$ yr$^{-1}$ sr$^{-1}$ eV$^{%g}$]' % (scale - 1) else: E_J_scale_units = units.m ** -2 * units.ns ** -1 * units.steradian ** -1 * units.eV ** (scale - 1) yl = r'$J(E)$ [m$^{-2}$ ns$^{-1}$ sr$^{-1}$ eV$^{%g}$]' % (scale - 1) ax.errorbar(E, E ** scale * J / E_J_scale_units, yerr=[(E ** scale * J_err_low) / E_J_scale_units, (E ** scale * J_err_up) / E_J_scale_units], marker='x', linewidth=1, markersize=8, ls='None', label=type) if scale != 0: yl = r'$E^{%g}\,$' % scale + yl ax.set_ylabel(yl) ax.set_xlabel(r'$E$ [eV]') ax.set_xscale('log') ax.set_yscale('log') plt.legend()