From a1b9b24ee32009675fc90aabd7db570ad7124c48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Thu, 10 Oct 2024 15:59:09 +0200 Subject: [PATCH 1/8] refactor get_time_delay function in rno-g detector --- NuRadioReco/detector/RNO_G/rnog_detector.py | 36 +++++++++++---------- NuRadioReco/detector/response.py | 6 +++- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/NuRadioReco/detector/RNO_G/rnog_detector.py b/NuRadioReco/detector/RNO_G/rnog_detector.py index ffca9c09d..3026c4c0c 100644 --- a/NuRadioReco/detector/RNO_G/rnog_detector.py +++ b/NuRadioReco/detector/RNO_G/rnog_detector.py @@ -1203,31 +1203,33 @@ def get_time_delay(self, station_id, channel_id, cable_only=False, use_stored=Tr signal_chain_dict = self.get_channel_signal_chain( station_id, channel_id) - time_delay = 0 - for key, value in signal_chain_dict["response_chain"].items(): - - if re.search("cable", key) is None and re.search("fiber", key) and cable_only: - continue - - weight = value.get("weight", 1) - if use_stored: - if "time_delay" not in value and "cable_delay" not in value: - self.logger.warning( - f"The signal chain component \"{key}\" of station.channel " - f"{station_id}.{channel_id} has no cable/time delay stored... (hence return time delay without it)") + if use_stored and not cable_only: + resp = self.get_signal_chain_response(station_id, channel_id) + return resp.get_time_delay() + elif use_stored and cable_only: + time_delays = resp.get_time_delays() + names = resp.get_names() + time_delay = 0 + for name, dt in zip(names, time_delays): + if re.search("cable", name) is None and re.search("fiber", name): continue + time_delay += dt + return time_delay + else: + time_delay = 0 + for key, value in signal_chain_dict["response_chain"].items(): - try: - time_delay += weight * value["time_delay"] - except KeyError: - time_delay += weight * value["cable_delay"] + if re.search("cable", key) is None and re.search("fiber", key) and cable_only: + continue - else: ydata = [value["mag"], value["phase"]] + # This is different from within `get_signal_chain_response` because we do set the time delay here + # and thus we do not remove it from the response. response = Response(value["frequencies"], ydata, value["y-axis_units"], name=key, station_id=station_id, channel_id=channel_id, log_level=self.__log_level) + weight = value.get("weight", 1) time_delay += weight * response._calculate_time_delay() return time_delay diff --git a/NuRadioReco/detector/response.py b/NuRadioReco/detector/response.py index f1efb6de7..14f86ce4f 100644 --- a/NuRadioReco/detector/response.py +++ b/NuRadioReco/detector/response.py @@ -380,10 +380,14 @@ def plot(self, ax1=None, show=False, in_dB=True, plt_kwargs={}): else: return fig, ax - def get_time_delay(self): + def get_time_delay(self, ): """ Get time delay from DB """ return np.sum(self.__time_delays) + def get_time_delays(self, ): + """ Get time delay from DB """ + return self.__time_delays + def _calculate_time_delay(self): """ Calculate time delay from phase of the stored complex response function. From 0ea3613cf54b461601ec6d58600bf2019a4813ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= <30903175+fschlueter@users.noreply.github.com> Date: Wed, 16 Oct 2024 11:07:53 +0200 Subject: [PATCH 2/8] Update rnog_detector.py: FIx get_time_delay Addressing Steffens comments --- NuRadioReco/detector/RNO_G/rnog_detector.py | 25 +++++++++++---------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/NuRadioReco/detector/RNO_G/rnog_detector.py b/NuRadioReco/detector/RNO_G/rnog_detector.py index 3026c4c0c..8f44b2447 100644 --- a/NuRadioReco/detector/RNO_G/rnog_detector.py +++ b/NuRadioReco/detector/RNO_G/rnog_detector.py @@ -1203,23 +1203,24 @@ def get_time_delay(self, station_id, channel_id, cable_only=False, use_stored=Tr signal_chain_dict = self.get_channel_signal_chain( station_id, channel_id) - if use_stored and not cable_only: + if use_stored: resp = self.get_signal_chain_response(station_id, channel_id) - return resp.get_time_delay() - elif use_stored and cable_only: - time_delays = resp.get_time_delays() - names = resp.get_names() - time_delay = 0 - for name, dt in zip(names, time_delays): - if re.search("cable", name) is None and re.search("fiber", name): - continue - time_delay += dt - return time_delay + if not cable_only: + return resp.get_time_delay() + else: + time_delays = resp.get_time_delays() + names = resp.get_names() + time_delay = 0 + for name, dt in zip(names, time_delays): + if re.search("cable", name) is None and re.search("fiber", name) is None: + continue + time_delay += dt + return time_delay else: time_delay = 0 for key, value in signal_chain_dict["response_chain"].items(): - if re.search("cable", key) is None and re.search("fiber", key) and cable_only: + if re.search("cable", key) is None and re.search("fiber", key) is None and cable_only: continue ydata = [value["mag"], value["phase"]] From c3f2129f313997a0b35f9a0567df29b284b1841c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= <30903175+fschlueter@users.noreply.github.com> Date: Wed, 16 Oct 2024 11:08:22 +0200 Subject: [PATCH 3/8] Update response.py --- NuRadioReco/detector/response.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NuRadioReco/detector/response.py b/NuRadioReco/detector/response.py index 14f86ce4f..ac13432ff 100644 --- a/NuRadioReco/detector/response.py +++ b/NuRadioReco/detector/response.py @@ -384,7 +384,7 @@ def get_time_delay(self, ): """ Get time delay from DB """ return np.sum(self.__time_delays) - def get_time_delays(self, ): + def get_time_delays(self): """ Get time delay from DB """ return self.__time_delays From 68fe0c5f1e45da44896c3738fbf8ca41387ba530 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Tue, 26 Nov 2024 13:26:11 +0100 Subject: [PATCH 4/8] Improve binning to calculate time delay from response --- NuRadioReco/detector/response.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/NuRadioReco/detector/response.py b/NuRadioReco/detector/response.py index ac13432ff..0d8d62b59 100644 --- a/NuRadioReco/detector/response.py +++ b/NuRadioReco/detector/response.py @@ -98,6 +98,8 @@ def __init__(self, frequency, y, y_unit, time_delay=0, weight=1, if y[0] is None or y[1] is None: raise ValueError("Data for response incomplete, detected \"None\"") + # print(name, y_unit[1], y[1]) + y_ampl, y_phase = np.array(y) if y_unit[0] == "dB": gain = 10 ** (y_ampl / 20) @@ -109,10 +111,14 @@ def __init__(self, frequency, y, y_unit, time_delay=0, weight=1, if y_unit[1].lower() == "deg": if np.max(np.abs(y_phase)) < 2 * np.pi: self.logger.warning("Is the phase really in deg? Does not look like it... " - f"Do not convert to rad. Phase: {y_phase}") + f"Do not convert {name} to rad: {y_phase}") else: y_phase = np.deg2rad(y_phase) elif y_unit[1].lower() == "rad": + # We can not make this test because the phase might be already unwrapped + # if np.amax(y_phase) - np.amin(y_phase) > 2 * np.pi: + # self.logger.warning("Is the phase really in rad? Does not look like it... " + # f"Do convert {name} to rad: {y_phase}") y_phase = y_phase else: raise KeyError @@ -404,15 +410,13 @@ def _calculate_time_delay(self): The time delay at ~ 200 MHz """ - freqs = np.arange(0.05, 1.2, 0.001) * units.GHz + freqs = np.arange(50, 1200, 0.5) * units.MHz response = self(freqs) - delta_freq = np.diff(freqs) - phase = np.angle(response) - time_delay = -np.diff(np.unwrap(phase)) / delta_freq / 2 / np.pi + time_delay = -np.diff(np.unwrap(phase)) / delta_freq / 2 / np.pi mask = np.all([195 * units.MHz < freqs, freqs < 250 * units.MHz], axis=0)[:-1] time_delay1 = np.mean(time_delay[mask]) From 1ba0e0862e50e3011cbf55a645b94ed84ac629aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Tue, 26 Nov 2024 13:28:01 +0100 Subject: [PATCH 5/8] Fix get_time_delay function --- NuRadioReco/detector/RNO_G/rnog_detector.py | 27 ++++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/NuRadioReco/detector/RNO_G/rnog_detector.py b/NuRadioReco/detector/RNO_G/rnog_detector.py index 8f44b2447..52e60f4f1 100644 --- a/NuRadioReco/detector/RNO_G/rnog_detector.py +++ b/NuRadioReco/detector/RNO_G/rnog_detector.py @@ -1199,7 +1199,6 @@ def get_time_delay(self, station_id, channel_id, cable_only=False, use_stored=Tr time_delay: float Sum of the time delays of all components in the signal chain for one channel """ - signal_chain_dict = self.get_channel_signal_chain( station_id, channel_id) @@ -1254,20 +1253,24 @@ def get_site(self, station_id): from NuRadioReco.detector import detector det = detector.Detector(source="rnog_mongo", log_level=logging.DEBUG, always_query_entire_description=False, - database_connection='RNOG_public', select_stations=24) + database_connection='RNOG_public', select_stations=13) - det.update(datetime.datetime(2023, 8, 2, 0, 0)) + det.update(datetime.datetime(2023, 7, 2, 0, 0)) + print(det.get_time_delay(13, 0)) + print(det.get_time_delay(13, 0, cable_only=True)) + # print(det.get_time_delay(13, 0, cable_only=False, use_stored=False)) + print(det.get_time_delay(13, 0, cable_only=True, use_stored=False)) - response = det.get_signal_chain_response(station_id=24, channel_id=0) + # response = det.get_signal_chain_response(station_id=24, channel_id=0) - from NuRadioReco.framework import electric_field - ef = electric_field.ElectricField(channel_ids=[0]) - ef.set_frequency_spectrum(np.ones(1025, dtype=complex), sampling_rate=2.4) + # from NuRadioReco.framework import electric_field + # ef = electric_field.ElectricField(channel_ids=[0]) + # ef.set_frequency_spectrum(np.ones(1025, dtype=complex), sampling_rate=2.4) - # Multipy the response to a trace. The multiply operator takes care of everything - trace_at_readout = ef * response + # # Multipy the response to a trace. The multiply operator takes care of everything + # trace_at_readout = ef * response - # getting the complex response as array - freq = np.arange(50, 1000) * units.MHz - complex_resp = response(freq) + # # getting the complex response as array + # freq = np.arange(50, 1000) * units.MHz + # complex_resp = response(freq) From c608e2b84f10bb5527b004816f5272bd63c182fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Tue, 26 Nov 2024 13:35:17 +0100 Subject: [PATCH 6/8] Remove option to get only cable delay from rno-g detector as this is ill-defined for the DB detector. --- NuRadioReco/detector/RNO_G/rnog_detector.py | 30 ++++----------------- 1 file changed, 5 insertions(+), 25 deletions(-) diff --git a/NuRadioReco/detector/RNO_G/rnog_detector.py b/NuRadioReco/detector/RNO_G/rnog_detector.py index 52e60f4f1..2e615cce1 100644 --- a/NuRadioReco/detector/RNO_G/rnog_detector.py +++ b/NuRadioReco/detector/RNO_G/rnog_detector.py @@ -1165,18 +1165,14 @@ def is_channel_noiseless(self, station_id, channel_id): def get_cable_delay(self, station_id, channel_id, use_stored=True): - """ - Return the cable delay of a signal chain as stored in the detector description. - This interface is required by simulation.py. See get_time_delay for description of - arguments. - """ + """ Same as `get_time_delay`. Only here to keep the same interface as the other detector classes. """ # FS: For the RNO-G detector description it is not easy to determine the cable delay alone # because it is not clear which reference components may need to be substraced. # However, having the cable delay without amplifiers is anyway weird. - return self.get_time_delay(station_id, channel_id, cable_only=False, use_stored=use_stored) + return self.get_time_delay(station_id, channel_id, use_stored=use_stored) - def get_time_delay(self, station_id, channel_id, cable_only=False, use_stored=True): - """ Return the sum of the time delay of all components in the signal chain calculated from the phase + def get_time_delay(self, station_id, channel_id, use_stored=True): + """ Return the sum of the time delay of all components in the signal chain calculated from the phase. Parameters ---------- @@ -1187,9 +1183,6 @@ def get_time_delay(self, station_id, channel_id, cable_only=False, use_stored=Tr channel_id: int The channel id - cable_only: bool - If True: Consider only cables to calculate delay. (Default: False) - use_stored: bool If True, take time delay as stored in DB rather than calculated from response. (Default: True) @@ -1204,24 +1197,11 @@ def get_time_delay(self, station_id, channel_id, cable_only=False, use_stored=Tr if use_stored: resp = self.get_signal_chain_response(station_id, channel_id) - if not cable_only: - return resp.get_time_delay() - else: - time_delays = resp.get_time_delays() - names = resp.get_names() - time_delay = 0 - for name, dt in zip(names, time_delays): - if re.search("cable", name) is None and re.search("fiber", name) is None: - continue - time_delay += dt - return time_delay + return resp.get_time_delay() else: time_delay = 0 for key, value in signal_chain_dict["response_chain"].items(): - if re.search("cable", key) is None and re.search("fiber", key) is None and cable_only: - continue - ydata = [value["mag"], value["phase"]] # This is different from within `get_signal_chain_response` because we do set the time delay here # and thus we do not remove it from the response. From 8d241ef563d2a309b8ed37cc4cfcc4bbd1cb1d56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= <30903175+fschlueter@users.noreply.github.com> Date: Tue, 3 Dec 2024 09:37:40 +0100 Subject: [PATCH 7/8] Update rnog_detector.py: Remove test code --- NuRadioReco/detector/RNO_G/rnog_detector.py | 23 ++++++++------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/NuRadioReco/detector/RNO_G/rnog_detector.py b/NuRadioReco/detector/RNO_G/rnog_detector.py index 9a9cb2a25..953a634c5 100644 --- a/NuRadioReco/detector/RNO_G/rnog_detector.py +++ b/NuRadioReco/detector/RNO_G/rnog_detector.py @@ -1264,20 +1264,15 @@ def get_site_coordinates(self, station_id=None): det.update(datetime.datetime(2023, 7, 2, 0, 0)) - print(det.get_time_delay(13, 0)) - print(det.get_time_delay(13, 0, cable_only=True)) - # print(det.get_time_delay(13, 0, cable_only=False, use_stored=False)) - print(det.get_time_delay(13, 0, cable_only=True, use_stored=False)) + response = det.get_signal_chain_response(station_id=13, channel_id=0) - # response = det.get_signal_chain_response(station_id=24, channel_id=0) + from NuRadioReco.framework import electric_field + ef = electric_field.ElectricField(channel_ids=[0]) + ef.set_frequency_spectrum(np.ones(1025, dtype=complex), sampling_rate=2.4) - # from NuRadioReco.framework import electric_field - # ef = electric_field.ElectricField(channel_ids=[0]) - # ef.set_frequency_spectrum(np.ones(1025, dtype=complex), sampling_rate=2.4) + # Multipy the response to a trace. The multiply operator takes care of everything + trace_at_readout = ef * response - # # Multipy the response to a trace. The multiply operator takes care of everything - # trace_at_readout = ef * response - - # # getting the complex response as array - # freq = np.arange(50, 1000) * units.MHz - # complex_resp = response(freq) + # getting the complex response as array + freq = np.arange(50, 1000) * units.MHz + complex_resp = response(freq) From d4e29548cd701c398637617624362d20aff8206d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= <30903175+fschlueter@users.noreply.github.com> Date: Tue, 3 Dec 2024 09:42:10 +0100 Subject: [PATCH 8/8] Update response.py: Clean up --- NuRadioReco/detector/response.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/NuRadioReco/detector/response.py b/NuRadioReco/detector/response.py index 2cf57c53e..3867b506d 100644 --- a/NuRadioReco/detector/response.py +++ b/NuRadioReco/detector/response.py @@ -98,8 +98,6 @@ def __init__(self, frequency, y, y_unit, time_delay=0, weight=1, if y[0] is None or y[1] is None: raise ValueError("Data for response incomplete, detected \"None\"") - # print(name, y_unit[1], y[1]) - y_ampl, y_phase = np.array(y) if y_unit[0] == "dB": gain = 10 ** (y_ampl / 20) @@ -386,7 +384,7 @@ def plot(self, ax1=None, show=False, in_dB=True, plt_kwargs={}): else: return fig, ax - def get_time_delay(self, ): + def get_time_delay(self): """ Get time delay from DB """ return np.sum(self.__time_delays)