From e2a65ea5efa503b5caf259ab6284abca5c0f3b71 Mon Sep 17 00:00:00 2001 From: Christian Glaser Date: Thu, 21 Nov 2024 14:25:48 +0000 Subject: [PATCH 1/5] include the readout time delays in the trace start time --- NuRadioMC/simulation/simulation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/NuRadioMC/simulation/simulation.py b/NuRadioMC/simulation/simulation.py index 8fb7a7a8c..b1a65b3b7 100644 --- a/NuRadioMC/simulation/simulation.py +++ b/NuRadioMC/simulation/simulation.py @@ -28,7 +28,7 @@ import NuRadioReco.modules.efieldToVoltageConverter import NuRadioReco.modules.channelAddCableDelay import NuRadioReco.modules.channelResampler -import NuRadioReco.modules.triggerTimeAdjuster +import NuRadioReco.modules.channelReadoutWindowCutter from NuRadioReco.detector import detector import NuRadioReco.framework.sim_station import NuRadioReco.framework.electric_field @@ -59,7 +59,7 @@ channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() channelResampler = NuRadioReco.modules.channelResampler.channelResampler() eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() -triggerTimeAdjuster = NuRadioReco.modules.triggerTimeAdjuster.triggerTimeAdjuster() +channelReadoutWindowCutter = NuRadioReco.modules.channelReadoutWindowCutter.channelReadoutWindowCutter() def merge_config(user, default): """ @@ -1519,7 +1519,7 @@ def run(self): if not evt.get_station().has_triggered(): continue - triggerTimeAdjuster.run(evt, station, self._det) + channelReadoutWindowCutter.run(evt, station, self._det) evt_group_triggered = True output_buffer[station_id][evt.get_id()] = evt # end event loop From 58d23938e778ed276ea369e13cb922399e4a0b0e Mon Sep 17 00:00:00 2001 From: Christian Glaser Date: Thu, 21 Nov 2024 14:51:36 +0000 Subject: [PATCH 2/5] commit modules --- NuRadioReco/modules/channelLengthAdjuster.py | 2 + .../modules/channelReadoutWindowCutter.py | 138 ++++++++++++++++++ NuRadioReco/modules/triggerTimeAdjuster.py | 6 + 3 files changed, 146 insertions(+) create mode 100644 NuRadioReco/modules/channelReadoutWindowCutter.py diff --git a/NuRadioReco/modules/channelLengthAdjuster.py b/NuRadioReco/modules/channelLengthAdjuster.py index 025282918..8e37b40f8 100644 --- a/NuRadioReco/modules/channelLengthAdjuster.py +++ b/NuRadioReco/modules/channelLengthAdjuster.py @@ -11,6 +11,8 @@ class channelLengthAdjuster: def __init__(self): self.number_of_samples = None self.offset = None + logger.warning("In most cases it is advisable to run a trigger module and use the channelReadoutWindowCutter module to cut the traces to the readout window \ + instead of this simple module.") self.begin(()) def begin(self, number_of_samples=256, offset=50): diff --git a/NuRadioReco/modules/channelReadoutWindowCutter.py b/NuRadioReco/modules/channelReadoutWindowCutter.py new file mode 100644 index 000000000..a483e7b91 --- /dev/null +++ b/NuRadioReco/modules/channelReadoutWindowCutter.py @@ -0,0 +1,138 @@ +from NuRadioReco.modules.base.module import register_run +import numpy as np +import logging +from NuRadioReco.utilities import units + +logger = logging.getLogger('NuRadioReco.channelReadoutWindowCutter') + + +class channelReadoutWindowCutter: + """ + Modifies channel traces to simulate the effects of the trigger + + The trace is cut to the length defined in the detector description relative to the trigger time. + If no trigger exists, nothing is done. + """ + + def __init__(self, log_level=logging.NOTSET): + logger.setLevel(log_level) + self.__sampling_rate_warning_issued = False + self.begin() + + def begin(self): + pass + + @register_run() + def run(self, event, station, detector): + """ + Cuts the traces to the readout window defined in the trigger. + + If multiple triggers exist, the primary trigger is used. If multiple + primary triggers exist, an error is raised. + If no primary trigger exists, the trigger with the earliest trigger time + is defined as the primary trigger and used to set the readout windows. + + Parameters + ---------- + event: `NuRadioReco.framework.event.Event` + + station: `NuRadioReco.framework.base_station.Station` + + detector: `NuRadioReco.detector.detector.Detector` + """ + counter = 0 + for i, (name, instance, kwargs) in enumerate(event.iter_modules(station.get_id())): + if name == 'channelReadoutWindowCutter': + if(kwargs['mode'] == mode): + counter += 1 + if counter > 1: + logger.warning('channelReadoutWindowCutter was called twice with the same mode. ' + 'This is likely a mistake. The module will not be applied again.') + return 0 + + + # determine which trigger to use + # if no primary trigger exists, use the trigger with the earliest trigger time + trigger = station.get_primary_trigger() + if trigger is None: # no primary trigger found + logger.debug('No primary trigger found. Using the trigger with the earliest trigger time.') + trigger = station.get_first_trigger() + if trigger is not None: + logger.info(f"setting trigger {trigger.get_name()} primary because it triggered first") + trigger.set_primary(True) + + if trigger is None: + logger.info('No trigger found! Channel timings will not be changed.') + return + + if trigger.has_triggered(): + trigger_time = trigger.get_trigger_time() + for channel in station.iter_channels(): + trigger_time_channel = trigger_time - channel.get_trace_start_time() + # if trigger_time_channel == 0: + # logger.warning(f"the trigger time is equal to the trace start time for channel {channel.get_id()}. This is likely because this module was already run on this station. The trace will not be changed.") + # continue + + trace = channel.get_trace() + trace_length = len(trace) + detector_sampling_rate = detector.get_sampling_frequency(station.get_id(), channel.get_id()) + sampling_rate = channel.get_sampling_rate() + self.__check_sampling_rates(detector_sampling_rate, sampling_rate) + + # this should ensure that 1) the number of samples is even and + # 2) resampling to the detector sampling rate results in the correct number of samples + # (note that 2) can only be guaranteed if the detector sampling rate is lower than the + # current sampling rate) + number_of_samples = int( + 2 * np.ceil( + detector.get_number_of_samples(station.get_id(), channel.get_id()) / 2 + * sampling_rate / detector_sampling_rate + )) + + if number_of_samples > trace.shape[0]: + logger.error("Input has fewer samples than desired output. Channels has only {} samples but {} samples are requested.".format( + trace.shape[0], number_of_samples)) + raise AttributeError + else: + trigger_time_sample = int(np.round(trigger_time_channel * sampling_rate)) + # logger.debug(f"channel {channel.get_id()}: trace_start_time = {channel.get_trace_start_time():.1f}ns, trigger time channel {trigger_time_channel/units.ns:.1f}ns, trigger time sample = {trigger_time_sample}") + channel_id = channel.get_id() + pre_trigger_time = trigger.get_pre_trigger_time_channel(channel_id) + samples_before_trigger = int(pre_trigger_time * sampling_rate) + cut_samples_beginning = 0 + if(samples_before_trigger <= trigger_time_sample): + cut_samples_beginning = trigger_time_sample - samples_before_trigger + roll_by = 0 + if(cut_samples_beginning + number_of_samples > trace_length): + logger.warning("trigger time is sample {} but total trace length is only {} samples (requested trace length is {} with an offest of {} before trigger). To achieve desired configuration, trace will be rolled".format( + trigger_time_sample, trace_length, number_of_samples, samples_before_trigger)) + roll_by = cut_samples_beginning + number_of_samples - trace_length # roll_by is positive + trace = np.roll(trace, -1 * roll_by) + cut_samples_beginning -= roll_by + rel_station_time_samples = cut_samples_beginning + roll_by + elif(samples_before_trigger > trigger_time_sample): + roll_by = -trigger_time_sample + samples_before_trigger + logger.warning(f"trigger time is before 'trigger offset window' (requested samples before trigger = {samples_before_trigger}," \ + f"trigger time sample = {trigger_time_sample}), the trace needs to be rolled by {roll_by} samples first" \ + f" = {roll_by / sampling_rate/units.ns:.2f}ns") + trace = np.roll(trace, roll_by) + + # shift trace to be in the correct location for cutting + trace = trace[cut_samples_beginning:(number_of_samples + cut_samples_beginning)] + channel.set_trace(trace, channel.get_sampling_rate()) + channel.set_trace_start_time(trigger_time - pre_trigger_time) + # channel.set_trace_start_time(channel.get_trace_start_time() + rel_station_time_samples / channel.get_sampling_rate()) + # logger.debug(f"setting trace start time to {channel.get_trace_start_time() + rel_station_time_samples / channel.get_sampling_rate():.0f} = {channel.get_trace_start_time():.0f} + {rel_station_time_samples / channel.get_sampling_rate():.0f}") + + + + def __check_sampling_rates(self, detector_sampling_rate, channel_sampling_rate): + if not self.__sampling_rate_warning_issued: # we only issue this warning once + if not np.isclose(detector_sampling_rate, channel_sampling_rate): + logger.warning( + 'triggerTimeAdjuster was called, but the channel sampling rate ' + f'({channel_sampling_rate/units.GHz:.3f} GHz) is not equal to ' + f'the target detector sampling rate ({detector_sampling_rate/units.GHz:.3f} GHz). ' + 'Traces may not have the correct trace length after resampling.' + ) + self.__sampling_rate_warning_issued = True diff --git a/NuRadioReco/modules/triggerTimeAdjuster.py b/NuRadioReco/modules/triggerTimeAdjuster.py index d9df0767e..c071487a1 100644 --- a/NuRadioReco/modules/triggerTimeAdjuster.py +++ b/NuRadioReco/modules/triggerTimeAdjuster.py @@ -18,6 +18,9 @@ def __init__(self, log_level=logging.NOTSET): logger.setLevel(log_level) self.__sampling_rate_warning_issued = False self.begin() + logger.warning(f"triggerTimeAdjuster is deprecated and will be removed soon. In most cased you can safely delete the application\ + of this module as it is automatically applied in NuRadioMC simulations. If you really need to use this module, \ + please use the channelReadoutWindowCutter module instead.") def begin(self): pass @@ -76,6 +79,9 @@ def run(self, event, station, detector, mode='sim_to_data'): If the ``trigger_name`` was specified in the ``begin`` function, only this trigger is considered. """ + logger.warning(f"triggerTimeAdjuster is deprecated and will be removed soon. In most cased you can safely delete the application\ + of this module as it is automatically applied in NuRadioMC simulations. If you really need to use this module, \ + please use the channelReadoutWindowCutter module instead.") counter = 0 for i, (name, instance, kwargs) in enumerate(event.iter_modules(station.get_id())): if name == 'triggerTimeAdjuster': From 60964b83272015e09a8f2343132a2483d5060a74 Mon Sep 17 00:00:00 2001 From: Christian Glaser Date: Thu, 21 Nov 2024 15:15:03 +0000 Subject: [PATCH 3/5] fix --- NuRadioReco/modules/channelReadoutWindowCutter.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/NuRadioReco/modules/channelReadoutWindowCutter.py b/NuRadioReco/modules/channelReadoutWindowCutter.py index a483e7b91..7c0df2e14 100644 --- a/NuRadioReco/modules/channelReadoutWindowCutter.py +++ b/NuRadioReco/modules/channelReadoutWindowCutter.py @@ -43,10 +43,9 @@ def run(self, event, station, detector): counter = 0 for i, (name, instance, kwargs) in enumerate(event.iter_modules(station.get_id())): if name == 'channelReadoutWindowCutter': - if(kwargs['mode'] == mode): - counter += 1 + counter += 1 if counter > 1: - logger.warning('channelReadoutWindowCutter was called twice with the same mode. ' + logger.warning('channelReadoutWindowCutter was called twice. ' 'This is likely a mistake. The module will not be applied again.') return 0 From 1efa01a84e3109acbc29617c3d641262a2c0c2ba Mon Sep 17 00:00:00 2001 From: Christian Glaser Date: Fri, 22 Nov 2024 12:23:44 +0000 Subject: [PATCH 4/5] address comments --- NuRadioReco/modules/__init__.py | 1 + .../{ => _deprecated}/triggerTimeAdjuster.py | 13 +++++++------ NuRadioReco/modules/channelLengthAdjuster.py | 4 ++-- 3 files changed, 10 insertions(+), 8 deletions(-) rename NuRadioReco/modules/{ => _deprecated}/triggerTimeAdjuster.py (92%) diff --git a/NuRadioReco/modules/__init__.py b/NuRadioReco/modules/__init__.py index e69de29bb..487c88c90 100644 --- a/NuRadioReco/modules/__init__.py +++ b/NuRadioReco/modules/__init__.py @@ -0,0 +1 @@ +from _deprecated import * \ No newline at end of file diff --git a/NuRadioReco/modules/triggerTimeAdjuster.py b/NuRadioReco/modules/_deprecated/triggerTimeAdjuster.py similarity index 92% rename from NuRadioReco/modules/triggerTimeAdjuster.py rename to NuRadioReco/modules/_deprecated/triggerTimeAdjuster.py index c071487a1..6628e3a4b 100644 --- a/NuRadioReco/modules/triggerTimeAdjuster.py +++ b/NuRadioReco/modules/_deprecated/triggerTimeAdjuster.py @@ -1,6 +1,7 @@ from NuRadioReco.modules.base.module import register_run import numpy as np import logging +import warnings from NuRadioReco.utilities import units logger = logging.getLogger('NuRadioReco.triggerTimeAdjuster') @@ -18,9 +19,9 @@ def __init__(self, log_level=logging.NOTSET): logger.setLevel(log_level) self.__sampling_rate_warning_issued = False self.begin() - logger.warning(f"triggerTimeAdjuster is deprecated and will be removed soon. In most cased you can safely delete the application\ - of this module as it is automatically applied in NuRadioMC simulations. If you really need to use this module, \ - please use the channelReadoutWindowCutter module instead.") + warnings.warn("triggerTimeAdjuster is deprecated and will be removed soon. In most cased you can safely delete the application " + "of this module as it is automatically applied in NuRadioMC simulations. If you really need to use this module, " + "please use the channelReadoutWindowCutter module instead.", DeprecationWarning) def begin(self): pass @@ -79,9 +80,9 @@ def run(self, event, station, detector, mode='sim_to_data'): If the ``trigger_name`` was specified in the ``begin`` function, only this trigger is considered. """ - logger.warning(f"triggerTimeAdjuster is deprecated and will be removed soon. In most cased you can safely delete the application\ - of this module as it is automatically applied in NuRadioMC simulations. If you really need to use this module, \ - please use the channelReadoutWindowCutter module instead.") + warnings.warn("triggerTimeAdjuster is deprecated and will be removed soon. In most cased you can safely delete the application " + "of this module as it is automatically applied in NuRadioMC simulations. If you really need to use this module, " + "please use the channelReadoutWindowCutter module instead.", DeprecationWarning) counter = 0 for i, (name, instance, kwargs) in enumerate(event.iter_modules(station.get_id())): if name == 'triggerTimeAdjuster': diff --git a/NuRadioReco/modules/channelLengthAdjuster.py b/NuRadioReco/modules/channelLengthAdjuster.py index 8e37b40f8..d570fdebb 100644 --- a/NuRadioReco/modules/channelLengthAdjuster.py +++ b/NuRadioReco/modules/channelLengthAdjuster.py @@ -11,8 +11,8 @@ class channelLengthAdjuster: def __init__(self): self.number_of_samples = None self.offset = None - logger.warning("In most cases it is advisable to run a trigger module and use the channelReadoutWindowCutter module to cut the traces to the readout window \ - instead of this simple module.") + logger.warning("In most cases it is advisable to run a trigger module and use the channelReadoutWindowCutter module to cut the traces to the readout window " + "instead of this simple module.") self.begin(()) def begin(self, number_of_samples=256, offset=50): From 309bc222f5325b3b2733c84e9205ed5aed099a3a Mon Sep 17 00:00:00 2001 From: Christian Glaser Date: Wed, 27 Nov 2024 15:54:03 +0000 Subject: [PATCH 5/5] properly initialize propagator --- .../examples/A01calculate_sim_efield.py | 21 ++++++++++++------- .../examples/A02calculate_channel.py | 20 +++++++++++------- NuRadioReco/modules/__init__.py | 2 +- NuRadioReco/modules/_deprecated/__init__.py | 0 4 files changed, 27 insertions(+), 16 deletions(-) create mode 100644 NuRadioReco/modules/_deprecated/__init__.py diff --git a/NuRadioMC/simulation/examples/A01calculate_sim_efield.py b/NuRadioMC/simulation/examples/A01calculate_sim_efield.py index 4b691c9e1..5d1b1a2c1 100644 --- a/NuRadioMC/simulation/examples/A01calculate_sim_efield.py +++ b/NuRadioMC/simulation/examples/A01calculate_sim_efield.py @@ -26,22 +26,27 @@ time_logger = NuRadioMC.simulation.time_logger.timeLogger(logger) +# initialize the detector description (from the json file) +kwargs = dict(json_filename="surface_station_1GHz.json", antenna_by_depth=False) +det = detector.Detector(**kwargs) +det.update(datetime.now()) + +# get the general config settings +cfg = sim.get_config("config.yaml") + # set the ice model ice = medium.get_ice_model('southpole_simple') # set the propagation module -propagator = propagation.get_propagation_module("analytic")(ice) +# it is important to pass the detector object to the propagator for an accurate calculation of the attenuation length +# (if the detector is availble, the sampling rate is used to determine the maximum frequency for an efficient interpolation +# of the attenuation length, otherwise the internal sampling rate of the simulation is used which is much higher, this would +# lead to inaccuracies at low frequencies) +propagator = propagation.get_propagation_module("analytic")(ice, detector=det, config=cfg) # set the station id and channel id sid = 101 cid = 1 -# get the general config settings -cfg = sim.get_config("config.yaml") - -# initialize the detector description (from the json file) -kwargs = dict(json_filename="surface_station_1GHz.json", antenna_by_depth=False) -det = detector.Detector(**kwargs) -det.update(datetime.now()) # define the showers that should be simulated showers = [] diff --git a/NuRadioMC/simulation/examples/A02calculate_channel.py b/NuRadioMC/simulation/examples/A02calculate_channel.py index 18d71b87e..9c372593b 100644 --- a/NuRadioMC/simulation/examples/A02calculate_channel.py +++ b/NuRadioMC/simulation/examples/A02calculate_channel.py @@ -24,22 +24,28 @@ propagation module to use (e.g. the analytic ray tracer). """ time_logger = NuRadioMC.simulation.time_logger.timeLogger(logger) +# initialize the detector description (from the json file) +kwargs = dict(json_filename="surface_station_1GHz.json", antenna_by_depth=False) +det = detector.Detector(**kwargs) +det.update(datetime.now()) + +# get the general config settings +cfg = sim.get_config("config.yaml") + # set the ice model ice = medium.get_ice_model('southpole_simple') # set the propagation module -propagator = propagation.get_propagation_module("analytic")(ice) +# it is important to pass the detector object to the propagator for an accurate calculation of the attenuation length +# (if the detector is availble, the sampling rate is used to determine the maximum frequency for an efficient interpolation +# of the attenuation length, otherwise the internal sampling rate of the simulation is used which is much higher, this would +# lead to inaccuracies at low frequencies) +propagator = propagation.get_propagation_module("analytic")(ice, detector=det, config=cfg) # set the station id and channel id sid = 101 cid = 1 -# get the general config settings -cfg = sim.get_config("config.yaml") -# initialize the detector description (from the json file) -kwargs = dict(json_filename="surface_station_1GHz.json", antenna_by_depth=False) -det = detector.Detector(**kwargs) -det.update(datetime.now()) # define the showers that should be simulated showers = [] diff --git a/NuRadioReco/modules/__init__.py b/NuRadioReco/modules/__init__.py index 487c88c90..01b401f89 100644 --- a/NuRadioReco/modules/__init__.py +++ b/NuRadioReco/modules/__init__.py @@ -1 +1 @@ -from _deprecated import * \ No newline at end of file +from NuRadioReco.modules._deprecated import triggerTimeAdjuster \ No newline at end of file diff --git a/NuRadioReco/modules/_deprecated/__init__.py b/NuRadioReco/modules/_deprecated/__init__.py new file mode 100644 index 000000000..e69de29bb