Source code for NuRadioReco.modules.ARIANNA.hardwareResponseIncorporator

from NuRadioReco.modules.base.module import register_run
from NuRadioReco.detector.ARIANNA import analog_components
import numpy as np
from NuRadioReco.utilities import units, fft
import time
import logging


[docs]class hardwareResponseIncorporator: """ Incorporates the gain and phase induced by the ARIANNA hardware. """ def __init__(self): self.logger = logging.getLogger("NuRadioReco.ARIANNA.hardwareResponseIncorporator") self.__debug = False self.__time_delays = {} self.__t = 0 self.__mingainlin = None self.begin()
[docs] def begin(self, debug=False): self.__debug = debug
[docs] def get_filter(self, frequencies, station_id, channel_id, det, sim_to_data=False, phase_only=False, mode=None, mingainlin=None): """ helper function to return the filter that the module applies. Parameters ---------- frequencies: array of floats the frequency array for which the filter should be returned station_id: int the station id channel_id: int the channel id det: detector instance the detector sim_to_data: bool (default False) if False, deconvolve the hardware response if True, convolve with the hardware response phase_only: bool (default False) if True, only the phases response is applied but not the amplitude response mode: string or None, default None Options: * 'phase_only': only the phases response is applied but not the amplitude response (identical to phase_only=True ) * 'relativ': gain of amp is divided by maximum of the gain, i.e. at the maximum of the filter response is 1 (before applying cable response). This makes it easier to compare the filtered to unfiltered signal * None : default, gain and phase effects are applied 'normally' mingainlin: float In frequency ranges where the gain gets very small, the reconstruction of the original signal (obtained by dividing the measured signal by the gain) leads to excessively high values, due to the effect of post-amplifier noise. In order to mitigate this effect, a minimum gain (linear scale!) as fraction of the maximum gain can be defined. If specified, any gain value smaller than mingainlin will be replaced by mingainlin. Note: The adjustment to the minimal gain is NOT visible when getting the amp response from analog_components.get_amplifier_response() Returns ------- array of complex floats the complex filter amplitudes """ amp_type = det.get_amplifier_type(station_id, channel_id) amp_measurement = det.get_amplifier_measurement(station_id, channel_id) amp_response = analog_components.get_amplifier_response(frequencies, amp_type=amp_type, amp_measurement=amp_measurement) if mingainlin is not None: mingainlin = float(mingainlin) ampmax = np.max(np.abs(amp_response)) iamp_gain_low = np.where(np.abs(amp_response) < (mingainlin * ampmax)) amp_response[iamp_gain_low] = (mingainlin * ampmax) * np.exp(1j * np.angle(amp_response[iamp_gain_low])) cable_response = analog_components.get_cable_response_parametrized(frequencies, * det.get_cable_type_and_length(station_id, channel_id)) if(mode == 'phase_only'): cable_response = np.ones_like(cable_response) * np.exp(1j * np.angle(cable_response)) amp_response = np.ones_like(amp_response) * np.angle(amp_response) elif(mode == 'relative'): ampmax = np.max(np.abs(amp_response)) amp_response /= ampmax if sim_to_data: return amp_response * cable_response else: return 1. / (amp_response * cable_response)
[docs] @register_run() def run(self, evt, station, det, sim_to_data=False, phase_only=False, mode=None, mingainlin=None): """ Switch sim_to_data to go from simulation to data or otherwise. The option zero_noise can be used to zero the noise around the pulse. It is unclear, how useful this is. Parameters ---------- evt: Event the event on which to run the module station: Station The station on which to run the module det: Detector or GenericDetector The detector description sim_to_data: bool (default False) if False, deconvolve the hardware response if True, convolve with the hardware response phase_only: bool (default False) if True, only the phases response is applied but not the amplitude response mode: string or None, default None Options: * 'phase_only': only the phases response is applied but not the amplitude response (identical to phase_only=True ) * 'relativ': gain of amp is divided by maximum of the gain, i.e. at the maximum of the filter response is 1 (before applying cable response). This makes it easier to compare the filtered to unfiltered signal * None : default, gain and phase effects are applied 'normally' mingainlin: float In frequency ranges where the gain gets very small, the reconstruction of the original signal (obtained by dividing the measured signal by the gain) leads to excessively high values, due to the effect of post-amplifier noise. In order to mitigate this effect, a minimum gain (linear scale!) as fraction of the maximum gain can be defined. If specified, any gain value smaller than mingainlin will be replaced by mingainlin. Note: The adjustment to the minimal gain is NOT visible when getting the amp response from analog_components.get_amplifier_response() """ self.__mingainlin = mingainlin if (phase_only): mode = 'phase_only' self.logger.warning('Please use option mode=''phase_only'' in the future, use of option phase_only will be phased out') t = time.time() for channel in station.iter_channels(): frequencies = channel.get_frequencies() trace_fft = channel.get_frequency_spectrum() trace_fft *= self.get_filter(frequencies, station.get_id(), channel.get_id(), det, sim_to_data, phase_only, mode, mingainlin) # zero first bins to avoid DC offset trace_fft[0] = 0 # hardwareResponse incorporator should always be used in conjunction with bandpassfilter # otherwise, noise will be blown up channel.set_frequency_spectrum(trace_fft, channel.get_sampling_rate()) if not sim_to_data: # Include cable delays cable_delay = det.get_cable_delay(station.get_id(), channel.get_id()) self.logger.debug("cable delay of channel {} is {}ns".format(channel.get_id(), cable_delay / units.ns)) channel.add_trace_start_time(-cable_delay) self.__t += time.time() - t
[docs] def end(self): from datetime import timedelta self.logger.setLevel(logging.INFO) dt = timedelta(seconds=self.__t) self.logger.info("total time used by this module is {}".format(dt)) return dt
def __calculate_time_delays_cable(self): """ helper function to calculate the time delay of the amp for a delta pulse """ sampling_rate = 10 * units.GHz # assume a huge sampling rate to have a good time resolution n = 2 ** 12 trace = np.zeros(n) trace[n // 2] = 1 max_time = trace.argmax() / sampling_rate spec = fft.time2freq(trace, sampling_rate) ff = np.fft.rfftfreq(n, 1. / sampling_rate) response = analog_components.get_cable_response(ff) response_gain = response['gain'] response_phase = response['phase'] trace2 = fft.freq2time(spec * response_gain * response_phase, sampling_rate) # import matplotlib.pyplot as plt # fig, ax = plt.subplots(1, 1) # ax.plot(trace) # ax.plot(trace2) # plt.show() max_time2 = np.abs(trace2).argmax() / sampling_rate return max_time2 - max_time def __calculate_time_delays_amp(self, amp_type): """ helper function to calculate the time delay of the amp for a delta pulse """ amp_response_f = analog_components.amplifier_response[amp_type] sampling_rate = 10 * units.GHz # assume a huge sampling rate to have a good time resolution n = 2 ** 12 trace = np.zeros(n) trace[n // 2] = 1 max_time = trace.argmax() / sampling_rate spec = fft.time2freq(trace, sampling_rate) ff = np.fft.rfftfreq(n, 1. / sampling_rate) amp_response_gain = amp_response_f['gain'](ff) amp_response_phase = amp_response_f['phase'](ff) mask = (ff < 70 * units.MHz) & (ff > 40 * units.MHz) spec[~mask] = 0 trace2 = fft.freq2time(spec * amp_response_gain * amp_response_phase, sampling_rate) # import matplotlib.pyplot as plt # fig, ax = plt.subplots(1, 1) # ax.plot(trace) # ax.plot(trace2) # plt.show() max_time2 = np.abs(trace2).argmax() / sampling_rate return max_time2 - max_time
[docs] def get_time_delay(self, amp_type): if(amp_type not in self.__time_delays.keys()): # not yet calculated -> calculate the time delay self.__time_delays[amp_type] = self.__calculate_time_delays_amp(amp_type) self.logger.info("time delays of amp {} have not yet been calculated -> calculating -> time delay is {:.2f} ns".format(amp_type, self.__time_delays[amp_type] / units.ns)) return self.__time_delays[amp_type]
[docs] def get_mingainlin(self): return self.__mingainlin