Add sanity check to subtract_time_delay_from_response

[fschlueter]

Hi all, i noticed that the subtraction of a time delay can go wrong if the response from which the time delay is subracted is to corsely sampled, i.e., the frequency binning is insufficient such that time_delay * Delta frequencies * 2 * pi > pi which caused issues.

I used this little program to convince myself of it.

import numpy as np
import matplotlib.pyplot as plt
from NuRadioReco.framework import electric_field
from NuRadioReco.utilities import analytic_pulse
from NuRadioReco.detector.response import subtract_time_delay_from_response, Response
import signal
signal.signal(signal.SIGINT, signal.SIG_DFL)

ef = electric_field.ElectricField(channel_ids=[0])
ef.set_frequency_spectrum(
    analytic_pulse.get_analytic_pulse_freq(
        1, -0.5, 1 * np.pi, 2048, 1, 2.4, bandpass=(0.1, 0.7)), sampling_rate=2.4)

ef.set_trace(np.roll(ef.get_trace(), 480), "same")

fig, ax = plt.subplots()
ax.grid()

dt = []
df = []
for n in np.arange(80, 150, 1):
    resp = np.ones(n) + 1j * np.zeros(n)
    freq = np.linspace(0, 1.2, len(resp))
    df.append(np.diff(2  * 50 * freq)[0])

    resp_shifted = subtract_time_delay_from_response(freq, resp, time_delay=-50)
    fake_response = Response(freq, [np.abs(resp_shifted), np.angle(resp_shifted)], ["mag", "rad"])
    dt.append(fake_response._calculate_time_delay())

ax.plot(df, dt)
ax.axhline(50, color="k", lw=1, ls="--")
ax.legend()
plt.show()

[fschlueter]

If you are fine with that we should also implement it in the function in trace_utilities

[sjoerd-bouma]

If you are fine with that we should also implement it in the function in trace_utilities

Have you checked this is actually a problem for the function in trace_utilities? I thought this was an issue that only arises because of your interpolation; it might therefore even be preferable to implement some clever fix in the interpolation (if that's possible).

[fschlueter]

More or less yes - the subtraction comes before any interpolation and I think there is a physical reason for that. However, I have not the time to make a code example and give a better answer atm. For trace_utilities is should not be a serious issue because it shifts traces which are typically recorded with a much higher sampling rate.

[sjoerd-bouma]

What am I missing?

from NuRadioReco.utilities import trace_utilities
import numpy as np
import matplotlib.pyplot as plt

trace = np.zeros(2048)
trace[1024] = 1
plt.plot(trace, color='k', label='original')

prop_cycle = plt.rcParams['axes.prop_cycle']
colors = list(prop_cycle.by_key()['color'])

shift_in_ns = 500
for i, sampling_rate in enumerate([2.4, 1.2, 0.6, 0.3, 0.1]): # GHz
    delayed_trace = trace_utilities.delay_trace(trace, sampling_rate, shift_in_ns)

    plt.plot(delayed_trace, color=colors[i]) # shifted trace
    plt.axvline((1024 + shift_in_ns * sampling_rate) % 2048, ls=':', color=colors[i])

plt.ylim(0, 1.5)
plt.legend()
plt.xlabel('sample')