Trigger Readout Window

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 = []

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 = []

NuRadioMC/simulation/simulation.py

@@ -29,7 +29,7 @@
 import NuRadioReco.modules.channelSignalReconstructor
 import NuRadioReco.modules.channelResampler
 import NuRadioReco.modules.channelGenericNoiseAdder
-import NuRadioReco.modules.triggerTimeAdjuster
+import NuRadioReco.modules.channelReadoutWindowCutter
 
 from NuRadioReco.detector import detector, antennapattern
 import NuRadioReco.framework.sim_station
@@ -59,7 +59,7 @@
 channelResampler = NuRadioReco.modules.channelResampler.channelResampler()
 channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder()
 eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter()
-triggerTimeAdjuster = NuRadioReco.modules.triggerTimeAdjuster.triggerTimeAdjuster()
+channelReadoutWindowCutter = NuRadioReco.modules.channelReadoutWindowCutter.channelReadoutWindowCutter()
 
 
 def merge_config(user, default):
@@ -1578,7 +1578,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
@@ -1628,30 +1628,32 @@ def run(self):
 
                         logger.debug(f"Adding sim_station to station {station_id} for event group {event_group.get_run_number()}, channel {channel_id}")
                         station.add_sim_station(sim_station)  # this will add the channels and efields to the existing sim_station object
-                        for evt in output_buffer[station_id].values():
-                            # determine the trigger that was used to determine the readout window
 
+                        # The non-triggered channels were simulated using the efieldToVoltageConverterPerEfield
+                        # (notice the "PerEfield" in the name). This means that each electric field was converted to
+                        # a sim channel. Now we still have to add together all sim channels associated with one "physical"
+                        # channel. Furthermore we have to cut out the correct readout window. For the trigger channels
+                        # this is done with the channelReadoutWindowCutter, here we have to do it manually.
+                        for evt in output_buffer[station_id].values():
                             for sim_channel in sim_station.get_channels_by_channel_id(channel_id):
                                 if not station.has_channel(sim_channel.get_id()):
-                                    # add empty channel with the correct length and time if it doesn't exist yet.
+                                    # For each physical channel we first create a "empty" trace (all zeros)
+                                    # with the start time and length ....
                                     self._add_empty_channel(station, channel_id)
 
-                                # Add the sim_channel to the station channel:
                                 channel = station.get_channel(sim_channel.get_id())
-                                # we need to account for the pre trigger time of the trigger that was used to determine the readout window
-                                pre_trigger_time = station.get_primary_trigger().get_pre_trigger_time_channel(channel_id)
-                                sim_channel_copy = copy.deepcopy(sim_channel)
-                                sim_channel_copy.set_trace_start_time(sim_channel.get_trace_start_time() + pre_trigger_time)
-                                channel.add_to_trace(sim_channel_copy)
+                                # ... and now add the sim channel to the correct window defined by the "empty trace"
+                                # At this point the traces are noiseless, hence, we do not have to raise an error.
+                                channel.add_to_trace(sim_channel, raise_error=False)
 
                 for evt in output_buffer[station_id].values():
-                    # we might not have a channel object in case there was no ray tracing solution to this channel, or if the timing did not match
-                    # the readout window. In this case we need to create a channel object and add it to the station
+                    # We might not have a channel object in case there was no ray tracing solution to this channel,
+                    # i.e. no sim_channel was associated to the channel. Hence, we add an empty channel object.
                     for channel_id in non_trigger_channels:
                         if not station.has_channel(channel_id):
                             self._add_empty_channel(station, channel_id)
 
-                    # the only thing left is to add noise to the non-trigger traces
+                    # The only thing left is to add noise to the non-trigger traces
                     # we need to do it a bit differently than for the trigger traces,
                     # because we need to add noise to traces where the amplifier response
                     # was already applied to.
@@ -1723,13 +1725,10 @@ def run(self):
 
     def _add_empty_channel(self, station, channel_id):
         """ Adds a channel with an empty trace (all zeros) to the station with the correct length and trace_start_time """
-        channel = NuRadioReco.framework.channel.Channel(channel_id)
-        n_samples = int(round(self._det.get_number_of_samples(station.get_id(), channel_id))
-                        * self._config['sampling_rate'] / self._det.get_sampling_frequency(station.get_id(), channel_id))
-        channel.set_trace(np.zeros(n_samples), self._config['sampling_rate'])
-        # we need to use any other channel to get the correct trace_start_time. All channels have the same start time at the end
-        # of the simulation.
-        channel.set_trace_start_time(station.get_channel(station.get_channel_ids()[0]).get_trace_start_time())
+        trigger = station.get_primary_trigger()
+        channel = NuRadioReco.modules.channelReadoutWindowCutter.get_empty_channel(
+            station.get_id(), channel_id, self._det, trigger, self._config['sampling_rate'])
+
         station.add_channel(channel)
 
     def _is_in_fiducial_volume(self, pos):

NuRadioMC/test/SingleEvents/T05validate_nur_file.py

@@ -20,10 +20,11 @@
 
 print("Testing the files {} and {} for equality".format(file1, file2))
 
-def all_traces(file):
+def all_traces(file, return_trace_start_times=True):
     eventReader1 = NuRadioReco.modules.io.eventReader.eventReader()
     eventReader1.begin(file)
     i = 0
+    trace_start_times = []
     for iE1, event1 in enumerate(eventReader1.run()):
         for st1, station1 in enumerate(event1.get_stations()):
             # print(f"eventid {event1.get_run_number()} station id: {station1.get_id()}")
@@ -35,11 +36,14 @@ def all_traces(file):
                 else:
                     all_traces = np.append(all_traces, trace1)
                     # print(f"apending trace {len(trace1)} to all_traces {len(all_traces)}")
+                trace_start_times += [channel1.get_trace_start_time()]
                 i += 1
+    if return_trace_start_times:
+        return all_traces, np.array(trace_start_times)
     return all_traces
 
-all_traces_1 = all_traces(file1)
-all_traces_2 = all_traces(file2)
+all_traces_1, trace_start_times_1 = all_traces(file1)
+all_traces_2, trace_start_times_2 = all_traces(file2)
 
 diff = all_traces_1 - all_traces_2
 
@@ -50,6 +54,11 @@ def all_traces(file):
 
 testing.assert_almost_equal(all_traces_1, all_traces_2,decimal=precision)
 
+# check that the trace_start_times are all equal
+testing.assert_almost_equal(
+    trace_start_times_1, trace_start_times_2, decimal=precision,
+    err_msg=f"Trace start times are not equal (maximum difference: {max(np.abs(trace_start_times_1-trace_start_times_2))})")
+
 try:
     testing.assert_equal(all_traces_1, all_traces_2)
 except:

NuRadioReco/framework/base_trace.py

@@ -299,20 +299,26 @@ def deserialize(self, data_pkl):
         if 'trace_start_time' in data.keys():
             self.set_trace_start_time(data['trace_start_time'])
 
-    def add_to_trace(self, channel):
+    def add_to_trace(self, channel, min_residual_time_offset=1e-5, raise_error=True):
         """
-        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
-        _time_trace containing zeros, this function can be seen as recording a channel in the specified
-        readout window.
+        Adds the trace of another channel to the trace of this channel.
+
+        The trace of "channel" is only added within the time window of "this". If "this" is an empty
+        trace (i.e., a trace containing zeros) with a defined trace_start_time, this function can be
+        seen as recording a channel in the specified readout window. Hence, "this" is referred
+        to as the "readout" in the comments of this function.
 
         Parameters
         ----------
         channel: BaseTrace
             The channel whose trace is to be added to the trace of this channel.
+        min_residual_time_offset: float (default: 1e-5)
+            Minimum risdual time between the target bin of this channel and the target bin of the channel
+            to be added. Below this threshold the residual time shift is not applied to increase performance
+            and minimize numerical artifacts from Fourier transforms.
+        raise_error: bool (default: True)
+            If True, an error is raised if "self" is not fully contained in "channel".
         """
-
         assert self.get_number_of_samples() is not None, "No trace is set for this channel"
         assert self.get_sampling_rate() == channel.get_sampling_rate(), "Sampling rates of the two channels do not match"
 
@@ -331,37 +337,59 @@ def add_to_trace(self, channel):
         # We handle 1+2x2 cases:
         # 1. Channel is completely outside readout window:
         if t1_channel < t0_readout or t1_readout < t0_channel:
+            if raise_error:
+                logger.error("The channel is completely outside the readout window")
+                raise ValueError('The channel is completely outside the readout window')
             return
+
+        def floor(x):
+            return int(np.floor(round(x, 5)))
+
+        def ceil(x):
+            return int(np.ceil(round(x, 5)))
+
         # 2. Channel starts before readout window:
         if t0_channel < t0_readout:
             i_start_readout = 0
             t_start_readout = t0_readout
-            i_start_channel = int((t0_readout-t0_channel) * sampling_rate_channel) + 1 # The first bin of channel inside readout
+            i_start_channel = ceil((t0_readout - t0_channel) * sampling_rate_channel) # The first bin of channel inside readout
             t_start_channel = tt_channel[i_start_channel]
         # 3. Channel starts after readout window:
         elif t0_channel >= t0_readout:
-            i_start_readout = int((t0_channel-t0_readout) * sampling_rate_readout) # The bin of readout right before channel starts
+            if t0_channel > t0_readout and raise_error:
+                logger.error("The channel starts after the readout window")
+                raise ValueError('The channel starts after the readout window')
+
+            i_start_readout = floor((t0_channel - t0_readout) * sampling_rate_readout) # The bin of readout right before channel starts
             t_start_readout = tt_readout[i_start_readout]
             i_start_channel = 0
             t_start_channel = t0_channel
+
         # 4. Channel ends after readout window:
         if t1_channel >= t1_readout:
-            i_end_readout = n_samples_readout - 1
-            t_end_readout = t1_readout
-            i_end_channel = int((t1_readout - t0_channel) * sampling_rate_channel) + 1 # The bin of channel right after readout ends
-            t_end_channel = tt_channel[i_end_channel]
+            i_end_readout = n_samples_readout
+            i_end_channel = ceil((t1_readout - t0_channel) * sampling_rate_channel) + 1 # The bin of channel right after readout ends
         # 5. Channel ends before readout window:
         elif t1_channel < t1_readout:
-            i_end_readout = int((t1_channel - t0_readout) * sampling_rate_readout) # The bin of readout right before channel ends
-            t_end_readout = tt_readout[i_end_readout]
-            i_end_channel = n_samples_channel - 1
-            t_end_channel = t1_channel
+            if raise_error:
+                logger.error("The channel ends before the readout window")
+                raise ValueError('The channel ends before the readout window')
+
+            i_end_readout = floor((t1_channel - t0_readout) * sampling_rate_readout) + 1 # The bin of readout right before channel ends
+            i_end_channel = n_samples_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)
-        tmp_channel.apply_time_shift(residual_time_offset)
-        trace_to_add = tmp_channel.get_trace()[i_start_channel:i_end_channel]
+        if np.abs(residual_time_offset) >= min_residual_time_offset:
+            tmp_channel = copy.deepcopy(channel)
+            tmp_channel.apply_time_shift(residual_time_offset)
+            trace_to_add = tmp_channel.get_trace()[i_start_channel:i_end_channel]
+        else:
+            trace_to_add = channel.get_trace()[i_start_channel:i_end_channel]
+
+        if i_end_readout - i_start_readout != i_end_channel - i_start_channel:
+            logger.error("The traces do not have the same length. This should not happen.")
+            raise ValueError('The traces do not have the same length. This should not happen.')
 
         # Add the trace to the original trace:
         original_trace = self.get_trace()

NuRadioReco/modules/__init__.py

@@ -0,0 +1 @@
+from NuRadioReco.modules._deprecated import triggerTimeAdjuster

NuRadioReco/modules/_deprecated/triggerTimeAdjuster.py

@@ -0,0 +1,9 @@
+import warnings
+
+class triggerTimeAdjuster:
+
+    def __init__(self, *args, **kwargs):
+        warnings.error("The module `triggerTimeAdjuster `is deprecated. 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)
+        raise NotImplementedError

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):