diff --git a/NuRadioReco/framework/base_trace.py b/NuRadioReco/framework/base_trace.py index 020a98167..fba66e1d9 100644 --- a/NuRadioReco/framework/base_trace.py +++ b/NuRadioReco/framework/base_trace.py @@ -258,6 +258,7 @@ def apply_time_shift(self, delta_t, silent=False): """ if delta_t > .1 * self.get_number_of_samples() / self.get_sampling_rate() and not silent: logger.warning('Trace is shifted by more than 10% of its length') + spec = self.get_frequency_spectrum() spec *= np.exp(-2.j * np.pi * delta_t * self.get_frequencies()) self.set_frequency_spectrum(spec, self._sampling_rate) @@ -299,9 +300,9 @@ def deserialize(self, data_pkl): def add_to_trace(self, channel): """ - Adds the trace of another channel to the trace of this channel. The trace is only added within the + Adds the trace of another channel to the trace of this channel. The trace is only added within the time window of "this" channel. - If this channel is an empty trace with a defined _sampling_rate and _trace_start_time, and a + If this channel is an empty trace with a defined _sampling_rate and _trace_start_time, and a _time_trace containing zeros, this function can be seen as recording a channel in the specified readout window. @@ -354,7 +355,7 @@ def add_to_trace(self, channel): t_end_readout = tt_readout[i_end_readout] i_end_channel = n_samples_channel - 1 t_end_channel = t1_channel - + # Determine the remaining time between the binning of the two traces and use time shift as interpolation: residual_time_offset = t_start_channel - t_start_readout tmp_channel = copy.deepcopy(channel) diff --git a/NuRadioReco/modules/analogToDigitalConverter.py b/NuRadioReco/modules/analogToDigitalConverter.py index b7482906b..1712b79ee 100644 --- a/NuRadioReco/modules/analogToDigitalConverter.py +++ b/NuRadioReco/modules/analogToDigitalConverter.py @@ -306,12 +306,12 @@ def get_digital_trace(self, station, det, channel, trace = np.fft.irfft(trace_fft * trigger_filter) - # Random clock offset - delayed_samples = len(trace) - int(np.round(MC_sampling_frequency / adc_sampling_frequency)) - 1 - trace = delay_trace(trace, MC_sampling_frequency, adc_time_delay, delayed_samples) - - times = times + 1.0 / adc_sampling_frequency - times = times[:len(trace)] + if adc_time_delay: + # Random clock offset + trace, dt_tstart = delay_trace(trace, MC_sampling_frequency, adc_time_delay) + if dt_tstart > 0: + times = times[dt_tstart / MC_sampling_frequency:] + times = times[:len(trace)] # Upsampling to 5 GHz before downsampling using interpolation. # We cannot downsample with a Fourier method because we want to keep diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 23ea274cd..7e91d53e7 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -1,10 +1,9 @@ +from NuRadioReco.utilities import units, ice, geometryUtilities as geo_utl, fft +import NuRadioReco.framework.base_trace + import numpy as np import scipy.constants import scipy.signal -from NuRadioReco.utilities import units -from NuRadioReco.utilities import ice -from NuRadioReco.utilities import geometryUtilities as geo_utl -from NuRadioReco.utilities import fft import logging logger = logging.getLogger('NuRadioReco.trace_utilities') @@ -20,7 +19,7 @@ def get_efield_antenna_factor(station, frequencies, channels, detector, zenith, Parameters ---------- - + station: Station frequencies: array of complex frequencies of the radio signal for which the antenna response is needed @@ -195,7 +194,7 @@ def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 ** Parameters ---------- - + trace: array of floats Trace to be upsampled original_sampling_frequency: float @@ -304,49 +303,82 @@ def apply_butterworth(spectrum, frequencies, passband, order=8): return filtered_spectrum -def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): +def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True): """ Delays a trace by transforming it to frequency and multiplying by phases. - Since this method is cyclic, the trace has to be cropped. It only accepts - positive delays, so some samples from the beginning are thrown away and then - some samples from the end so that the total number of samples is equal to + + Since this method is cyclic, the delayed trace has to be cropped. A positive + delay means that the trace is shifted to the right, i.e., its delayed. A + negative delay would mean that the trace is shifted to the left. Samples + from the beginning or end are thrown away. Optionally one can crop the + trace (from the end) so that the total number of samples is equal to the argument delayed samples. Parameters ---------- - trace: array of floats + trace: array of floats or `NuRadioReco.framework.base_trace.BaseTrace` Array containing the trace sampling_frequency: float Sampling rate for the trace time_delay: float Time delay used for transforming the trace. Must be positive or 0 - delayed_samples: integer or None - Number of samples that the delayed trace must contain - if None: the trace is not cut + crop_trace: bool + If True (default), the trace is cropped to remove samples + what are unphysical after delaying (rolling) the trace. Returns ------- delayed_trace: array of floats The delayed, cropped trace + dt_start: float (optional) + The delta t of the trace start time. Only returned if crop_trace is True. """ - - if time_delay < 0: - msg = 'Time delay must be positive' - raise ValueError(msg) - - n_samples = len(trace) - - spectrum = fft.time2freq(trace, sampling_frequency) - frequencies = np.fft.rfftfreq(n_samples, 1 / sampling_frequency) + # Do nothing if time_delay is 0 + if not time_delay: + if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace): + if crop_trace: + return trace.get_trace(), 0 + else: + return trace.get_trace() + else: + if crop_trace: + return trace, 0 + else: + return trace + + if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace): + spectrum = trace.get_frequency_spectrum() + frequencies = trace.get_frequencies() + else: + n_samples = len(trace) + spectrum = fft.time2freq(trace, sampling_frequency) + frequencies = np.fft.rfftfreq(n_samples, 1 / sampling_frequency) spectrum *= np.exp(-1j * 2 * np.pi * frequencies * time_delay) delayed_trace = fft.freq2time(spectrum, sampling_frequency) - init_sample = int(time_delay * sampling_frequency) + 1 + cycled_samples = int(round(time_delay * sampling_frequency)) - if delayed_samples is not None: - delayed_trace = delayed_trace[init_sample:None] - delayed_trace = delayed_trace[:delayed_samples] + if crop_trace: + if time_delay >= 0: + delayed_trace = delayed_trace[cycled_samples:] + dt_start = cycled_samples * sampling_frequency + else: + delayed_trace = delayed_trace[:-cycled_samples] + dt_start = 0 + + return delayed_trace, dt_start + + else: + # Check if unphysical samples contain any signal and if so, throw a warning + if time_delay > 0: + if np.any(np.abs(delayed_trace[:cycled_samples]) > 0.01 * units.microvolt): + logger.warning("The delayed trace has unphysical samples that contain signal. " + "Consider cropping the trace to remove these samples.") + else: + if np.any(np.abs(delayed_trace[-cycled_samples:]) > 0.01 * units.microvolt): + logger.warning("The delayed trace has unphysical samples that contain signal. " + "Consider cropping the trace to remove these samples.") - return delayed_trace + return delayed_trace