Source code for NuRadioReco.modules.cosmicRayEnergyReconstructor

import numpy as np
from NuRadioReco.modules.base.module import register_run
from NuRadioReco.framework.parameters import stationParameters as stnp
from NuRadioReco.framework.parameters import electricFieldParameters as efp
from NuRadioReco.utilities import units
import NuRadioReco.utilities.trace_utilities
import logging
import radiotools.helper as hp
import radiotools.coordinatesystems
import radiotools.atmosphere.models


[docs]class cosmicRayEnergyReconstructor: """ Reconstructs the energy of an air shower from its radio signal Requires the following modules to be run beforehand: * a 10th order Butterworth bandpass filter with passband 80-300 MHz * a direction reconstruction * the voltageToAnalyticEfieldConverter """ def __init__(self): self.logger = logging.getLogger('NuRadioReco.cosmicRayEnergyReconstructor') self.__atmosphere = radiotools.atmosphere.models.Atmosphere() self.__parametrizations = { 'mooresbay': { 'scale': np.array([(442.46, -281.75, 324.96), (394.08, -308.36, 436.30)]), 'falloff': np.array([(-.1584, -.07943), (.8070, -1.4098)]) }, 'southpole': { 'scale': np.array([(976.30, -1213.43, 626.98), (643.39, -667.08, 478.06)]), 'falloff': np.array([(-.2273, .05627), (1.3372, -2.1653)]) }, 'auger': { 'scale': np.array([(229.96, -123.75, 110.51), (214.46, -111.01, 119.18)]), 'falloff': np.array([(-.1445, -.09820), (.5936, -1.1763)]) } } self.__elevations = { # TODO: This should be changed once we have implemented a proper coordinate system 'mooresbay': 30., 'southpole': 2800., 'auger': 1560. } self.__site = None self.__parametrization_for_site = None self.__elevation = None
[docs] def begin(self, site='mooresbay'): self.__site = site if site not in self.__parametrizations.keys(): self.logger.error('Unsupported site. Please select one of the following: {}'.format(self.__parametrizations.keys())) raise ValueError self.__parametrization_for_site = self.__parametrizations[site] self.__elevation = self.__elevations[site]
[docs] @register_run() def run(self, event, station, detector): if not station.is_cosmic_ray(): self.logger.warning('Event is not a cosmic ray!') if not station.has_parameter(stnp.zenith) or not station.has_parameter(stnp.azimuth): self.logger.error('No incoming direction available. Energy can not be reconstructed!') return zenith = station.get_parameter(stnp.zenith) azimuth = station.get_parameter(stnp.azimuth) if zenith < 30. * units.deg: self.logger.warning('Zenith angle is smaller than 30deg. Energy reconstruction is likely to be inaccurate!') n_efields = len(station.get_electric_fields()) if n_efields == 0: self.logger.error('No E-field found. Please run the voltageToAnalyticEfieldConverter beforehand!') return if n_efields > 1: self.logger.warning('Multiple E-fields were found. Only the first E-field will be used.') electric_field = station.get_electric_fields()[0] spectrum_slope = electric_field.get_parameter(efp.cr_spectrum_slope) alpha = hp.get_angle_to_magnetic_field_vector(zenith, azimuth, self.__site) cs = radiotools.coordinatesystems.cstrafo(zenith, azimuth, site=self.__site) efield_trace_vxB_vxvxB = cs.transform_to_vxB_vxvxB(cs.transform_from_onsky_to_ground(electric_field.get_trace())) efield_trace_vxB_vxvxB[0] /= np.sin(alpha) # correct energy fluence for effect of angle to magnetic field energy_fluence = NuRadioReco.utilities.trace_utilities.get_electric_field_energy_fluence(efield_trace_vxB_vxvxB, electric_field.get_times()) energy_fluence = np.abs(energy_fluence[0]) + np.abs(energy_fluence[1]) xmax_distance = self.__atmosphere.get_distance_xmax_geometric(zenith, 750., self.__elevation) # parametrization is for Xmax of 750g/cm^2 # find out if we are inside or outside of the Cherenkov ring second_order_spectrum_parameter = electric_field.get_parameter(efp.cr_spectrum_quadratic_term) if second_order_spectrum_parameter <= spectrum_slope * .1: scale_parameter = self.__parametrization_for_site['scale'][0][0] * zenith ** 2 + self.__parametrization_for_site['scale'][0][1] * zenith + self.__parametrization_for_site['scale'][0][0] falloff_parameter = self.__parametrization_for_site['falloff'][0][0] * zenith + self.__parametrization_for_site['falloff'][0][1] else: scale_parameter = self.__parametrization_for_site['scale'][1][0] * zenith ** 2 + self.__parametrization_for_site['scale'][1][1] * zenith + self.__parametrization_for_site['scale'][1][0] falloff_parameter = self.__parametrization_for_site['falloff'][1][0] * zenith + self.__parametrization_for_site['falloff'][1][1] rec_energy = 1.e18 * np.sqrt(energy_fluence) * (xmax_distance / units.km) / (scale_parameter * np.exp(falloff_parameter * np.abs(spectrum_slope) ** 0.8)) station.set_parameter(stnp.cr_energy_em, rec_energy)