diff --git a/NuRadioMC/simulation/examples/A01calculate_sim_efield.py b/NuRadioMC/simulation/examples/A01calculate_sim_efield.py index 4b691c9e18..5d1b1a2c18 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 18d71b87e3..9c372593b0 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/NuRadioMC/simulation/simulation.py b/NuRadioMC/simulation/simulation.py index a69bd4ccfe..ffc7ddb5b5 100644 --- a/NuRadioMC/simulation/simulation.py +++ b/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): diff --git a/NuRadioMC/test/SingleEvents/T05validate_nur_file.py b/NuRadioMC/test/SingleEvents/T05validate_nur_file.py index 5a799727cf..7da35f1c43 100755 --- a/NuRadioMC/test/SingleEvents/T05validate_nur_file.py +++ b/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: diff --git a/NuRadioReco/framework/base_trace.py b/NuRadioReco/framework/base_trace.py index 3894a5d465..3c16ad1eaf 100644 --- a/NuRadioReco/framework/base_trace.py +++ b/NuRadioReco/framework/base_trace.py @@ -2,6 +2,7 @@ from NuRadioReco.utilities import fft, bandpass_filter import NuRadioReco.detector.response +from NuRadioReco.utilities import units import numpy as np import logging @@ -243,7 +244,7 @@ def get_number_of_samples(self): length = (self._frequency_spectrum.shape[-1] - 1) * 2 return length - def apply_time_shift(self, delta_t, silent=False): + def apply_time_shift(self, delta_t, silent=False, fourier_shift_threshold = 1e-5 * units.ns): """ Uses the fourier shift theorem to apply a time shift to the trace Note that this is a cyclic shift, which means the trace will wrap @@ -257,12 +258,23 @@ def apply_time_shift(self, delta_t, silent=False): Turn off warnings if time shift is larger than 10% of trace length Only use this option if you are sure that your trace is long enough to acommodate the time shift + fourier_shift_threshold: float (default: 1e-5 * units.ns) + Threshold for the Fourier shift. If the shift is closer to a multiple of + 1 / sampling_rate than this, the trace is rolled instead of using the Fourier + shift theorem to save time and avoid numerical errors in the Fourier transforms. """ - if delta_t > .1 * self.get_number_of_samples() / self.get_sampling_rate() and not silent: + if np.abs(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) + + if np.abs(np.round(delta_t * self.get_sampling_rate()) - delta_t * self.get_sampling_rate()) < fourier_shift_threshold: + roll_by = int(np.round(delta_t * self.get_sampling_rate())) + trace = self.get_trace() + trace = np.roll(trace, roll_by, axis=-1) + self.set_trace(trace, self.get_sampling_rate()) + else: + spec = self.get_frequency_spectrum() + spec *= np.exp(-2.j * np.pi * delta_t * self.get_frequencies()) + self.set_frequency_spectrum(spec, self.get_sampling_rate()) def resample(self, sampling_rate): if sampling_rate == self.get_sampling_rate(): @@ -299,20 +311,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 +349,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() diff --git a/NuRadioReco/modules/__init__.py b/NuRadioReco/modules/__init__.py index e69de29bb2..01b401f89d 100644 --- a/NuRadioReco/modules/__init__.py +++ b/NuRadioReco/modules/__init__.py @@ -0,0 +1 @@ +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 0000000000..e69de29bb2 diff --git a/NuRadioReco/modules/_deprecated/triggerTimeAdjuster.py b/NuRadioReco/modules/_deprecated/triggerTimeAdjuster.py new file mode 100644 index 0000000000..2cab7d8ce5 --- /dev/null +++ b/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 diff --git a/NuRadioReco/modules/channelLengthAdjuster.py b/NuRadioReco/modules/channelLengthAdjuster.py index 0252829187..d570fdebb0 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 0000000000..5103919021 --- /dev/null +++ b/NuRadioReco/modules/channelReadoutWindowCutter.py @@ -0,0 +1,171 @@ +from NuRadioReco.modules.base.module import register_run +import NuRadioReco.framework.channel +from NuRadioReco.utilities import units + +import numpy as np +import logging +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': + counter += 1 + if counter > 1: + logger.warning('channelReadoutWindowCutter was called twice. ' + '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 or not trigger.has_triggered(): + logger.info('No trigger found (which triggered)! Channel timings will not be changed.') + return + + trigger_time = trigger.get_trigger_time() + for channel in station.iter_channels(): + + detector_sampling_rate = detector.get_sampling_frequency(station.get_id(), channel.get_id()) + sampling_rate = channel.get_sampling_rate() + n_samples = detector.get_number_of_samples(station.get_id(), channel.get_id()) + self.__check_sampling_rates(channel, detector_sampling_rate, sampling_rate, n_samples) + + + # 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(n_samples / 2 * sampling_rate / detector_sampling_rate)) + + trace = channel.get_trace() + 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 + + # windowed_channel = get_empty_channel( + # station.get_id(), channel.get_id(), detector, trigger, sampling_rate) + # windowed_channel.add_to_trace(channel) + + # channel.set_trace(windowed_channel.get_trace(), windowed_channel.get_sampling_rate()) + # channel.set_trace_start_time(windowed_channel.get_trace_start_time()) + + channel_id = channel.get_id() + pre_trigger_time = trigger.get_pre_trigger_time_channel(channel_id) + + pre_trigger_time_channel = trigger_time - pre_trigger_time - channel.get_trace_start_time() + + # throw warinings or errors depending on if the readout window is outside or too far outside the trace + trace_length = len(trace) + if pre_trigger_time_channel < -pre_trigger_time: + msg = ("Start of the readout window is more than pre_trigger_time before the start of the " + "trace (trigger time = {}ns, pre-trigger time = {}ns, start of trace {}ns, requested " + "time before trace = {}ns), this would result in rolling over the edge of the trace " + "and is not the intended use of this function").format( + trigger_time, pre_trigger_time, channel.get_trace_start_time(), pre_trigger_time_channel) + logger.error(msg) + raise AttributeError(msg) + elif pre_trigger_time_channel < 0: + msg = ("Start of the readout window is before the start of the trace (trigger time = {}ns, " + "pre-trigger time = {}ns, start of trace {}ns, requested time before trace = {}ns), " + "the trace will be rolled over the edge to fit in the readout window").format( + trigger_time, pre_trigger_time, channel.get_trace_start_time(), pre_trigger_time_channel) + logger.warning(msg) + elif trigger_time > channel.get_trace_start_time() + trace_length / sampling_rate: + msg = ("Trigger time is after the end of the trace (trigger time = {}ns, end of trace {}ns), " + "this would result in rolling over the edge of the trace and is not the intended use " + "of this function").format( + trigger_time, channel.get_trace_start_time() + trace_length / sampling_rate, + pre_trigger_time_channel + number_of_samples / sampling_rate - trace_length / sampling_rate) + logger.error(msg) + raise AttributeError(msg) + elif pre_trigger_time_channel + number_of_samples / sampling_rate > trace_length / sampling_rate: + msg = ("End of the readout window is outside the end of the trace (trigger time = {}ns, " + "pre-trigger time = {}ns, length of readout window {}ns, end of trace {}, requested time after trace = {}ns), " + "the trace will be rolled over the edge to fit in the readout window").format( + trigger_time, pre_trigger_time, number_of_samples/sampling_rate, channel.get_trace_start_time() + trace_length/sampling_rate, + pre_trigger_time_channel + number_of_samples / sampling_rate - trace_length / sampling_rate) + logger.warning(msg) + + # "roll" the start of the readout window to the start of the trace + channel.apply_time_shift(-pre_trigger_time_channel, silent=True) + + # cut the trace + trace = channel.get_trace() + trace = trace[:number_of_samples] + + channel.set_trace(trace, channel.get_sampling_rate()) + channel.set_trace_start_time(trigger_time - pre_trigger_time) + + + def __check_sampling_rates(self, channel, detector_sampling_rate, channel_sampling_rate, n_samples): + 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( + 'channelReadoutWindowCutter 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 + + +def get_empty_channel(station_id, channel_id, detector, trigger, sampling_rate): + channel = NuRadioReco.framework.channel.Channel(channel_id) + + # Get the correct number of sample for the final sampling rate + sampling_rate_ratio = sampling_rate / detector.get_sampling_frequency(station_id, channel_id) + n_samples = int(round(detector.get_number_of_samples(station_id, channel_id) * sampling_rate_ratio)) + + # get the correct trace start time taking into account different `pre_trigger_times` + channel_trace_start_time = trigger.get_trigger_time() - trigger.get_pre_trigger_time_channel(channel_id) + + channel.set_trace(np.zeros(n_samples), sampling_rate) + channel.set_trace_start_time(channel_trace_start_time) + + return channel diff --git a/NuRadioReco/modules/io/LOFAR/readLOFARData.py b/NuRadioReco/modules/io/LOFAR/readLOFARData.py index d2d9ff3b95..a70d90bf60 100644 --- a/NuRadioReco/modules/io/LOFAR/readLOFARData.py +++ b/NuRadioReco/modules/io/LOFAR/readLOFARData.py @@ -222,7 +222,7 @@ def nrrID_to_tbbID(channel_id): fourth element of the channelID with a "0". Parameters - ------------ + ---------- channel_id: str or int Channel ID @@ -550,7 +550,7 @@ def begin(self, event_id, logger_level=logging.NOTSET): # Read in energy estimate from LORA energy = lora_dict["LORA"]["energy_GeV"] - + # The LORA coordinate system has x pointing East -> set this through magnetic field vector (values from 2015) self.__hybrid_shower.set_parameter(showerParameters.magnetic_field_vector, np.array([0.004675, 0.186270, -0.456412])) @@ -560,7 +560,7 @@ def begin(self, event_id, logger_level=logging.NOTSET): # Add LORA core and energy to parameters. The z-Position of the core is always at 7.6m for LOFAR self.__hybrid_shower.set_parameter(showerParameters.core, np.array([core_pos_x * units.m, core_pos_y * units.m, 7.6 * units.m])) self.__hybrid_shower.set_parameter(showerParameters.energy, energy * units.GeV) - + # Go through TBB directory and identify all files for this event tbb_filename_pattern = tbb_filetag_from_unix(self.__lora_timestamp) @@ -720,7 +720,7 @@ def run(self, detector, trace_length=65536): evt.set_station(station) lofar_trace_access.close_file() - + # Add general event radio shower to event to store reconstruction values later radio_shower = NuRadioReco.framework.radio_shower.RadioShower( shower_id=evt.get_id(), station_ids=evt.get_station_ids() diff --git a/NuRadioReco/modules/triggerTimeAdjuster.py b/NuRadioReco/modules/triggerTimeAdjuster.py deleted file mode 100644 index d9df0767ef..0000000000 --- a/NuRadioReco/modules/triggerTimeAdjuster.py +++ /dev/null @@ -1,180 +0,0 @@ -from NuRadioReco.modules.base.module import register_run -import numpy as np -import logging -from NuRadioReco.utilities import units - -logger = logging.getLogger('NuRadioReco.triggerTimeAdjuster') - - -class triggerTimeAdjuster: - """ - 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 - - - def get_pre_trigger_time(self, trigger_name, channel_id): - """ Get the pre_trigger_time for a given trigger_name and channel_id """ - if isinstance(self.__pre_trigger_time, float): - return self.__pre_trigger_time - - pre_trigger_time = self.__pre_trigger_time - while isinstance(pre_trigger_time, dict): - if trigger_name in pre_trigger_time: # keys are different triggers - pre_trigger_time = pre_trigger_time[trigger_name] - elif channel_id in pre_trigger_time: # keys are channel_ids - pre_trigger_time = pre_trigger_time[channel_id] - else: - logger.error( - 'pre_trigger_time was specified as a dictionary, ' - f'but neither the trigger_name {trigger_name} ' - f'nor the channel id {channel_id} are present as keys' - ) - raise KeyError - - return pre_trigger_time - - @register_run() - def run(self, event, station, detector, mode='sim_to_data'): - """ - Run the trigger time adjuster. - - This module can be used either to 'cut' the simulated traces into - the appropriate readout windows, or to adjust the trace start times - of simulated / real data to account for the different trigger readout - delays. - - 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` - - mode: 'sim_to_data' (default) | 'data_to_sim' - If 'sim_to_data', cuts the (arbitrary-length) simulated traces - to the appropriate readout windows. - If 'data_to_sim', looks through all triggers in the station and adjusts the - trace_start_time according to the different readout delays - - If the ``trigger_name`` was specified in the ``begin`` function, - only this trigger is considered. - """ - counter = 0 - for i, (name, instance, kwargs) in enumerate(event.iter_modules(station.get_id())): - if name == 'triggerTimeAdjuster': - if(kwargs['mode'] == mode): - counter += 1 - if counter > 1: - logger.warning('triggerTimeAdjuster 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 mode == 'sim_to_data': - 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) - # 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}") - - elif mode == 'data_to_sim': - for channel in station.iter_channels(): - pre_trigger_time = trigger.get_pre_trigger_time_channel(channel.get_id()) - channel.set_trace_start_time(channel.get_trace_start_time()-pre_trigger_time) - else: - raise ValueError(f"Argument '{mode}' for mode is not valid. Options are 'sim_to_data' or 'data_to_sim'.") - - 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/documentation/source/NuRadioReco/pages/times.rst b/documentation/source/NuRadioReco/pages/times.rst index 779db32a79..40bb016f7c 100644 --- a/documentation/source/NuRadioReco/pages/times.rst +++ b/documentation/source/NuRadioReco/pages/times.rst @@ -93,14 +93,18 @@ We list all relevant modules that are used for a MC simulation and reconstructio * unfolds amplifier -> also implies a time delay in the channel trace * cable delay is subtracted from the trace start time (due to the limited trace length, the trace is not rolled to account for cable delays) -* `NuRadioReco.modules.triggerTimeAdjuster` - * ``sim_to_data`` mode: This modules cuts the trace to the correct length (as specified in the detector description) around the trigger time with a pre-trigger time as defined by the respective trigger module. In the case of multiple triggers it used the primary trigger. If no primary trigger is defined, it uses the trigger with the earliest trigger time. In the end, the `trace_start_time ` is set to the trigger time. This is done because this reflects what raw experimental data looks like. - * ``data_to_sim`` mode: The module determines the trigger that was used to cut the trace to its current length (the 'sim_to_data' step above in case of simulations) and adjusts the `trace_start_time ` according to the different readout delays. The "primary trigger" defines the readout delays. **After** applying this module in the "data_to_sim" direction, the position in the trace that caused the trigger can be found via the `trigger_time `. - * `NuRadioReco.modules.channelStopFilter`: this module prepends and appends all channels with a fixed length (128ns by default). The 'prepend' time is subtracted from the trace start time (because all channels get the same time delay). It additionally applies a tukey window to taper off the start and end (by default, the first and last 5%) of the trace. * `NuRadioReco.modules.voltageToEfieldConverter`: * the traces from all used channels are cut to the overlapping region (including delays due to geometry and differences in delays due to different group delays in hardware, e.g. different antenna/amplifier responses) - * the E-field `trace_start_time ` is set accordingly \ No newline at end of file + * the E-field `trace_start_time ` is set accordingly + +* `NuRadioReco.modules.channelReadoutWindowCutter`: + * Cuts out the readout windows from simulated traces according to the `trigger_time`, `pre_trigger_time`, and `number_of_samples` (i.e., length of trace) parameters. + +* NuRadioReco.modules._deprecated.triggerTimeAdjuster (deprecated): + * This module is now deprecated. It was replaced by the `NuRadioReco.modules.channelReadoutWindowCutter` module for simulation (doing essential what the mode ``sim_to_data`` did with the notable difference that it sets the trace start time to `trigger_time - pre_trigger_time`. The mode ``data_to_sim`` is not needed anymore, i.e., the experiment specific IO modules have to make sure that the trace start time is set correctly. + * ``sim_to_data`` mode: This modules cuts the trace to the correct length (as specified in the detector description) around the trigger time with a pre-trigger time as defined by the respective trigger module. In the case of multiple triggers it used the primary trigger. If no primary trigger is defined, it uses the trigger with the earliest trigger time. In the end, the `trace_start_time ` is set to the trigger time. This is done because this reflects what raw experimental data looks like. + * ``data_to_sim`` mode: The module determines the trigger that was used to cut the trace to its current length (the 'sim_to_data' step above in case of simulations) and adjusts the `trace_start_time ` according to the different readout delays. The "primary trigger" defines the readout delays. **After** applying this module in the "data_to_sim" direction, the position in the trace that caused the trigger can be found via the `trigger_time `.