Source code for NuRadioReco.modules.templateDirectionFitter

from NuRadioReco.modules.base.module import register_run
import numpy as np
from scipy import constants
from radiotools import helper as hp
from NuRadioReco.utilities import units, ice
from NuRadioReco.framework.parameters import stationParameters as stnp
from NuRadioReco.framework.parameters import channelParameters as chp

import logging
logger = logging.getLogger('NuRadioReco.templateDirectionFitter')


[docs]class templateDirectionFitter: """ Calculates basic signal parameters. """ def __init__(self): pass
[docs] def begin(self): pass
[docs] @register_run() def run(self, evt, station, det, channels_to_use=None, cosmic_ray=False): """ Fits the direction using templates Parameters ---------- evt: event station: station det: detector channels_to_use: list (default: [0, 1, 2, 3] antenna to use for fit cosmic_ray: bool type to set correlation template """ if channels_to_use is None: channels_to_use = [0, 1, 2, 3] if(cosmic_ray): type_str = 'cr' xcorrelations = chp.cr_xcorrelations else: type_str = 'nu' xcorrelations = chp.nu_xcorrelations station_id = station.get_id() channels = station.iter_channels(channels_to_use) times = [] positions = [] for iCh, channel in enumerate(channels): channel_id = channel.get_id() times.append(channel[xcorrelations]['{}_ref_xcorr_time'.format(type_str)] + channel.get_trace_start_time()) positions.append(det.get_relative_position(station_id, channel_id)) times = np.array(times) positions = np.array(positions) site = det.get_site(station_id) n_ice = ice.get_refractive_index(-0.01, site) from scipy import optimize as opt def obj_plane(params, positions, t_measured): zenith, azimuth = params if cosmic_ray: if((zenith < 0) or (zenith > 0.5 * np.pi)): return np.inf else: if((zenith < 0.5 * np.pi) or (zenith > np.pi)): return np.inf v = hp.spherical_to_cartesian(zenith, azimuth) c = constants.c * units.m / units.s if not cosmic_ray: c = c / n_ice logger.debug("using speed of light = {:.4g}".format(c)) t_expected = -(np.dot(v, positions.T) / c) sigma = 1 * units.ns chi2 = np.sum(((t_expected - t_expected.mean()) - (t_measured - t_measured.mean())) ** 2 / sigma ** 2) logger.debug("texp = {texp}, tm = {tmeas}, {chi2}".format(texp=t_expected, tmeas=t_measured, chi2=chi2)) return chi2 method = "Nelder-Mead" options = {'maxiter': 1000, 'disp': False} zenith_start = 135 * units.deg if cosmic_ray: zenith_start = 45 * units.deg starting_chi2 = {} for starting_az in np.array([0, 90, 180, 270]) * units.degree: starting_chi2[starting_az] = obj_plane((zenith_start, starting_az), positions, times) azimuth_start = min(starting_chi2, key=starting_chi2.get) res = opt.minimize(obj_plane, x0=[zenith_start, azimuth_start], args=(positions, times), method=method, options=options) output_str = "reconstucted angles theta = {:.1f}, phi = {:.1f}".format(res.x[0] / units.deg, hp.get_normalized_angle(res.x[1]) / units.deg) if station.has_sim_station(): sim_zen = station.get_sim_station()[stnp.zenith] sim_az = station.get_sim_station()[stnp.azimuth] dOmega = hp.get_angle(hp.spherical_to_cartesian(sim_zen, sim_az), hp.spherical_to_cartesian(res.x[0], res.x[1])) output_str += " MC theta = {:.1f}, phi = {:.1f}, dOmega = {:.2f}".format(sim_zen / units.deg, sim_az / units.deg, dOmega / units.deg) logger.info(output_str) station[stnp.zenith] = res.x[0] station[stnp.azimuth] = hp.get_normalized_angle(res.x[1]) if(cosmic_ray): station[stnp.cr_zenith] = res.x[0] station[stnp.cr_azimuth] = hp.get_normalized_angle(res.x[1]) else: station[stnp.nu_zenith] = res.x[0] station[stnp.nu_azimuth] = hp.get_normalized_angle(res.x[1])
[docs] def end(self): pass