From 001d2de97862f9cc48da7c7d71eb7d3ca18329bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Fri, 22 Nov 2024 12:02:19 +0100 Subject: [PATCH 01/10] Allow negative delays. Always crop the trace by the cycled samples --- NuRadioReco/utilities/trace_utilities.py | 25 ++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 57bc82fb4..8d43ffcf2 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -20,7 +20,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 @@ -120,7 +120,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 @@ -232,9 +232,12 @@ def apply_butterworth(spectrum, frequencies, passband, order=8): def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): """ 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 @@ -254,11 +257,6 @@ def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): delayed_trace: array of floats The delayed, cropped trace """ - - if time_delay < 0: - msg = 'Time delay must be positive' - raise ValueError(msg) - n_samples = len(trace) spectrum = fft.time2freq(trace, sampling_frequency) @@ -268,10 +266,13 @@ def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): delayed_trace = fft.freq2time(spectrum, sampling_frequency) - init_sample = int(time_delay * sampling_frequency) + 1 + cycled_samples = int(time_delay * sampling_frequency) + 1 + if time_delay > 0: + delayed_trace = delayed_trace[cycled_samples:] + else: + delayed_trace = delayed_trace[:cycled_samples] if delayed_samples is not None: - delayed_trace = delayed_trace[init_sample:None] delayed_trace = delayed_trace[:delayed_samples] return delayed_trace From 6ab4e138ab733f4b455cdd5189c478969d44018a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Fri, 22 Nov 2024 12:15:41 +0100 Subject: [PATCH 02/10] Add a return argument to correctly addjust the trace start time after cropping. Also change calculation of the number of samples to be cropped. The original code was wron when time_delay / sampling_frequency was a round number. --- NuRadioReco/utilities/trace_utilities.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 8d43ffcf2..5e7c60f01 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -256,6 +256,8 @@ def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): ------- delayed_trace: array of floats The delayed, cropped trace + dt_start: float + The time change of the trace start time. """ n_samples = len(trace) @@ -266,13 +268,15 @@ def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): delayed_trace = fft.freq2time(spectrum, sampling_frequency) - cycled_samples = int(time_delay * sampling_frequency) + 1 + cycled_samples = int(round(time_delay * sampling_frequency)) 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 if delayed_samples is not None: delayed_trace = delayed_trace[:delayed_samples] - return delayed_trace + return delayed_trace, dt_start From bea7db3ba18a61e45ffaf56eca482953fb04c0a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Fri, 29 Nov 2024 17:07:09 +0100 Subject: [PATCH 03/10] Add option to not crop delayed traces. Add warning if non-physical samples contain a signal --- .../modules/analogToDigitalConverter.py | 6 +-- NuRadioReco/utilities/trace_utilities.py | 42 ++++++++++++------- 2 files changed, 30 insertions(+), 18 deletions(-) diff --git a/NuRadioReco/modules/analogToDigitalConverter.py b/NuRadioReco/modules/analogToDigitalConverter.py index b7482906b..8a46b6483 100644 --- a/NuRadioReco/modules/analogToDigitalConverter.py +++ b/NuRadioReco/modules/analogToDigitalConverter.py @@ -307,10 +307,10 @@ 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) + trace, dt_tstart = delay_trace(trace, MC_sampling_frequency, adc_time_delay) - times = times + 1.0 / adc_sampling_frequency + if dt_tstart > 0: + times = times[dt_tstart / MC_sampling_frequency:] times = times[:len(trace)] # Upsampling to 5 GHz before downsampling using interpolation. diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 5e7c60f01..9b4e26063 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -229,7 +229,7 @@ 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. @@ -248,16 +248,16 @@ def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): 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 - The time change of the trace start time. + dt_start: float (optional) + The delta t of the trace start time. Only returned if crop_trace is True. """ n_samples = len(trace) @@ -269,14 +269,26 @@ def delay_trace(trace, sampling_frequency, time_delay, delayed_samples=None): delayed_trace = fft.freq2time(spectrum, sampling_frequency) cycled_samples = int(round(time_delay * sampling_frequency)) - 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 - if delayed_samples is not 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, dt_start + return delayed_trace From 867a0390d44cc4540e2c573a9b9550e32acb0ba3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Fri, 29 Nov 2024 17:23:18 +0100 Subject: [PATCH 04/10] Allow pass a base_trace object --- NuRadioReco/utilities/trace_utilities.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 9b4e26063..e52163b59 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') @@ -242,7 +241,7 @@ def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True): 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 @@ -259,10 +258,13 @@ def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True): dt_start: float (optional) The delta t of the trace start time. Only returned if crop_trace is True. """ - n_samples = len(trace) - - spectrum = fft.time2freq(trace, sampling_frequency) - frequencies = np.fft.rfftfreq(n_samples, 1 / sampling_frequency) + 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) From 8def4f98b6f346b526efea18c92f2a008817686b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Fri, 29 Nov 2024 17:25:45 +0100 Subject: [PATCH 05/10] Use trace_utility in apply_time_shift --- NuRadioReco/framework/base_trace.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/NuRadioReco/framework/base_trace.py b/NuRadioReco/framework/base_trace.py index 020a98167..2b4570252 100644 --- a/NuRadioReco/framework/base_trace.py +++ b/NuRadioReco/framework/base_trace.py @@ -5,7 +5,7 @@ import decimal import numbers import functools -from NuRadioReco.utilities import fft, bandpass_filter +from NuRadioReco.utilities import fft, bandpass_filter, trace_utilities import NuRadioReco.detector.response import scipy.signal import copy @@ -258,9 +258,9 @@ 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) + + delayed_trace = trace_utilities.apply_time_shift(self, self._sampling_rate, delta_t, crop_trace=False) + self.set_trace(delayed_trace, self._sampling_rate) def resample(self, sampling_rate): if sampling_rate == self.get_sampling_rate(): @@ -299,9 +299,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 +354,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) From 20aedf3de8029258939b073266ba0fd3ac08f525 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Mon, 2 Dec 2024 08:24:15 +0100 Subject: [PATCH 06/10] Fix function call --- NuRadioReco/framework/base_trace.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NuRadioReco/framework/base_trace.py b/NuRadioReco/framework/base_trace.py index 2b4570252..385312b44 100644 --- a/NuRadioReco/framework/base_trace.py +++ b/NuRadioReco/framework/base_trace.py @@ -259,7 +259,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') - delayed_trace = trace_utilities.apply_time_shift(self, self._sampling_rate, delta_t, crop_trace=False) + delayed_trace = trace_utilities.delay_trace(self, self._sampling_rate, delta_t, crop_trace=False) self.set_trace(delayed_trace, self._sampling_rate) def resample(self, sampling_rate): From c010d0f50b044a2688dfddf968995bf5d4fcca00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Mon, 2 Dec 2024 11:55:15 +0100 Subject: [PATCH 07/10] Do nothing when delay is 0 --- NuRadioReco/utilities/trace_utilities.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index e52163b59..970651858 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -258,6 +258,13 @@ def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True): dt_start: float (optional) The delta t of the trace start time. Only returned if crop_trace is True. """ + if not time_delay: + # If time delay is 0 + if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace): + return trace.get_trace(), 0 + else: + return trace, 0 + if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace): spectrum = trace.get_frequency_spectrum() frequencies = trace.get_frequencies() @@ -273,7 +280,7 @@ def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True): cycled_samples = int(round(time_delay * sampling_frequency)) if crop_trace: - if time_delay > 0: + if time_delay >= 0: delayed_trace = delayed_trace[cycled_samples:] dt_start = cycled_samples * sampling_frequency else: From 89a791a3c5fa2e08a0823d9aee59b94615bf8fb2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Mon, 2 Dec 2024 12:18:58 +0100 Subject: [PATCH 08/10] return correct number of objects --- NuRadioReco/utilities/trace_utilities.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 970651858..59a973f32 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -258,12 +258,18 @@ def delay_trace(trace, sampling_frequency, time_delay, crop_trace=True): dt_start: float (optional) The delta t of the trace start time. Only returned if crop_trace is True. """ + # Do nothing if time_delay is 0 if not time_delay: - # If time delay is 0 if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace): - return trace.get_trace(), 0 + if crop_trace: + return trace.get_trace(), 0 + else: + return trace.get_trace() else: - return trace, 0 + if crop_trace: + return trace, 0 + else: + return trace if isinstance(trace, NuRadioReco.framework.base_trace.BaseTrace): spectrum = trace.get_frequency_spectrum() From b3499fee7bad699ac4afda711bcab25c98e685f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Mon, 2 Dec 2024 14:32:00 +0100 Subject: [PATCH 09/10] Make it explicit in the module --- NuRadioReco/modules/analogToDigitalConverter.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/NuRadioReco/modules/analogToDigitalConverter.py b/NuRadioReco/modules/analogToDigitalConverter.py index 8a46b6483..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 - 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)] + 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 From f1a04cb7d201867a3319cb27e627de4afe653ed1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= <30903175+fschlueter@users.noreply.github.com> Date: Thu, 12 Dec 2024 17:28:56 +0100 Subject: [PATCH 10/10] Update base_trace.py: Revert changes introduced with this PR. It should be rather happening in an dedicated PR --- NuRadioReco/framework/base_trace.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/NuRadioReco/framework/base_trace.py b/NuRadioReco/framework/base_trace.py index 385312b44..fba66e1d9 100644 --- a/NuRadioReco/framework/base_trace.py +++ b/NuRadioReco/framework/base_trace.py @@ -5,7 +5,7 @@ import decimal import numbers import functools -from NuRadioReco.utilities import fft, bandpass_filter, trace_utilities +from NuRadioReco.utilities import fft, bandpass_filter import NuRadioReco.detector.response import scipy.signal import copy @@ -259,8 +259,9 @@ 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') - delayed_trace = trace_utilities.delay_trace(self, self._sampling_rate, delta_t, crop_trace=False) - self.set_trace(delayed_trace, self._sampling_rate) + 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) def resample(self, sampling_rate): if sampling_rate == self.get_sampling_rate():