Trigger Readout Window

[cg-laser]

the trace_start_time will include the (channel-specific) trigger time delays at the end of a NuRadioMC simulation.
For a clean implementation and to avoid bugs in old code, a new module channelReadoutWindowCutter was written. The previous triggerTimeAdjuster module will be deprecated in a future release.

[fschlueter]

Do we want to move the deprecated modules into a folder deprecated (still allowing to import them from NuRadioReco.modules ofc?

[cg-laser]

Do we want to move the deprecated modules into a folder deprecated (still allowing to import them from NuRadioReco.modules ofc?

no, because then the usercode will fail and they won't get the deprecation message ;-)
for the next release, I would just delete them.

[sjoerd-bouma]

This looks fine to me - as far as I can tell it really only changes the trace start times back to how they were pre-2.2.0.

I just have some nitpicky comments about multiline strings and a request to use the warnings module for DeprecationWarnings as this will make it easier to catch this stuff in the tests in the future.

[sjoerd-bouma]

Do we want to move the deprecated modules into a folder deprecated (still allowing to import them from NuRadioReco.modules ofc?

no, because then the usercode will fail and they won't get the deprecation message ;-) for the next release, I would just delete them.

I think Felix's suggestion was to move the modules to a deprecated folder, and change the original module to from deprecated import *. You keep the functionality but move deprecated code away. I don't have strong feelings about this either way (though maybe we would want to name the module _deprecated to avoid it showing up as a suggestion for users).

[fschlueter]

Do we want to move the deprecated modules into a folder deprecated (still allowing to import them from NuRadioReco.modules ofc?

no, because then the usercode will fail and they won't get the deprecation message ;-) for the next release, I would just delete them.

I think Felix's suggestion was to move the modules to a deprecated folder, and change the original module to from deprecated import *. You keep the functionality but move deprecated code away. I don't have strong feelings about this either way (though maybe we would want to name the module _deprecated to avoid it showing up as a suggestion for users).

I would have added to the __init__.py something like:

from deprecated import *

which should avoid that user code fails.

[cg-laser]

thanks for the comments! I fixed them. We still need to update the tests once we agree on the implementation

[fschlueter]

What is with the channelLengthAdjuster do we still need this module?

[fschlueter]

What is with the channelLengthAdjuster do we still need this module?

Probably we can keep it because it is a simple module which one can use in reconstruction pipelines?

[cg-laser]

yes, this module might be useful for some studies where running a trigger module would be overkill

[fschlueter]

Do we need this warning at all? We ensure that we have the correct number of samples with:

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

Maybe we can only check that detector_sampling_rate <= channel_sampling_rate?

[sjoerd-bouma]

I think when I implemented this I thought about it sufficiently to realize that if sampling_rate < detector_sampling_rate we cannot guarantee the correct number of samples (with the way resampling is currently implemented); I am not actually sure if the correct number of samples is guaranteed in all cases if sampling_rate > detector_sampling_rate; my feeling was that this is only true in most cases which is why I left the warning as it is. If you can work out the maths and figure out that this is wrong and we don't need the warning we can change this : )

[sjoerd-bouma]

I'm a bit confused, and think potentially more changes are needed. The triggerTimeAdjuster (now channelReadoutWindowCutter) is called in line 1522; later non-trigger channels are added still, with (in the new convention) wrong trace start times; I am also not 100% sure that the slightly hacky looking adjustments to pre_trigger_time by making a deep copy of the sim channels were correct, but in any case I guess with the new convention for trace_start_time they aren't correct now.

Is there any reason why we do not just call the triggerTimeAdjuster (now channelReadoutWindowCutter) at the end of the simulation, after downsampling to the detector sampling rate has already occurred? I feel like this would be much less error-prone. I guess there will be a small performance impact because we have to generate more noise, but probably this only makes a small difference in the grand scheme of things?

[sjoerd-bouma]

Also also (seeing as I looked at this today) - please update the documentation page documentation/NuRadioReco/pages/times.rst before merging, bad documentation is worse than no documentation : )

[fschlueter]

I also think this is not yet correct. Let me elaborate a bit. The SimChannels from the non-trigger channels do not all have the same trace_start_time since they were simulated with apply_det_response_sim (maybe an obvious fact but I wanted to put it here). For those non-trigger channels, the Channel object is created with self._add_empty_channel(station, channel_id) which does:

        # 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())

So it creates a channel with the same trace_start_time as the first channel (which at this point has to be a trigger channel I think/hope). For this trigger channel the trace_start_time is correct and already takes the pre_trigger_time into account. However, a new channel can have for the same trigger a different pre_trigger_time (we want to simulate different readout delays with that right?). So if we add the traces of all the SimChannel to this Channel we have to correct the trace_start_time for this channel and not for the SimChannels as currently done:

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

Hence I am suggesting to change the function _add_empty_channel to:

    def _add_empty_channel(self, station, channel_id, primary_trigger):
        """ 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'])
        
        trigger_channel_id = station.get_channel_ids()[0]
        trigger_channel_trace_start_time = station.get_channel(trigger_channel_id).get_trace_start_time()
        channel_trace_start_time = trigger_channel_trace_start_time + (
            primary_trigger.get_pre_trigger_time_channel(channel_id) - primary_trigger.get_pre_trigger_time_channel(trigger_channel_id))
        
        channel.set_trace_start_time(channel_trace_start_time)
        station.add_channel(channel)

And the code block in the run method to:

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.
            self._add_empty_channel(station, channel_id)

        # Add the sim_channel to the station channel:
        channel = station.get_channel(sim_channel.get_id())
        channel.add_to_trace(sim_channel_copy)

[fschlueter]

FYI, I noticed that the function _add_empty_channel just exists in develop yet ... My code snippet should still work I hope.

[fschlueter]

I started implementing my ideas into a new branch: https://github.com/nu-radio/NuRadioMC/tree/refactor-triggertimeadjuster-fs. Alongside the changes described above, I removed everything from the deprecated module + raising an Error now, and drastically reduced the code in the new module which now also uses the add_to_trace function.

[fschlueter]

Running the test which fails on #820 on this branch gives us a lot of:

WARNING - 2025-02-04 14:16:36,095 - NuRadioReco.channelReadoutWindowCutter - trigger time is before 'trigger offset window' (requested samples before trigger = 275,trigger time sample = 166), the trace needs to be rolled by 109 samples first = 21.80ns

This could be the reason why the unit test fails? @cg-laser you updated at one point the unit tests with values which should match this simulations?

[fschlueter]

What is also confusing me. Why can we roll the traces without adjusting trace_start_time? This module should not effect the signal time, correct? So when you roll a trace you would need to adjust for that by correcting the trace_start_time? However, the module and also the processor seemed to have had this feature since always. So is that a bug? It does of course not matter to much if you always roll all channels by the same amount but this is or was never guaranteed.

[fschlueter]

What is also confusing me. Why can we roll the traces without adjusting trace_start_time? This module should not effect the signal time, correct? So when you roll a trace you would need to adjust for that by correcting the trace_start_time? However, the module and also the processor seemed to have had this feature since always. So is that a bug? It does of course not matter to much if you always roll all channels by the same amount but this is or was never guaranteed.

Okay, apologies. I confused myself again. The rolling is just applied to bring the pulse to the correct position in the trace which matches the trace_start_time.

[fschlueter]

Okay, I merged a branch with PR #820 into this branch to fix the issue with this branch. Originally I tried to implement the channelReadoutWindowCutter in a completely different way in PR #820. This failed so far but I felt the code commented in the module (only 5 lines). I think we should discuss about it on Thursday.

Completely unrelated question @cg-laser, this PR also edits two scripts. Is that on purpose?

[fschlueter]

@sjoerd-bouma I updated the documentation as you wished. However, I think now we are getting errors because we can not reference NuRadioReco.modules._deprecated.triggerTimeAdjuster. Any suggestion?

[cg-laser]

Thanks a lot Felix, I had some discussion with Martin this morning and it would be good to further improve the add_to_trace function. The function should get a precision parameter which determines if an additional residual time correction is performed. Now this is done inexplicitely by rounding to 5 digits. We can also allow more tolerance by default, e.g. 1ps should be sufficiently precise.

[fschlueter]

Fine, but that can be done in a different PR?

[fschlueter]

I just added a flag to add_to_trace(channel) that by default an error is raised if the trace of channel does not cover the entire readout window. The reason behind this is simple: If we have simulations with noise and this happens the resulting trace of the readout window will contain a bunches of zeros. Now the error actually is raised for non-trigger channels in one test. I actually implemented the flag exactly for this scenario. Background is that for non-trigger channels we do not ensure that traces are long enough. However, what I forgot is that for those channels the noise is not yet added. Which means we should be good. If everyone is okay with that I will just set the flag to False for this particular case of non-trigger channels.

[fschlueter]

Okay, the trace_start_time failed for one test for 55 ns which sounds like it is the default pre_trigger_time. Shall we just update the reference values than?

[MartinRavn]

I have two remaining concerns about the channelReadoutWindowCuttter:

1. It may not be obvious to a user that the function should not be used when the readout window is not almost completely inside the simulated trace. Now, the function will fail if the window is completely outside the trace, but maybe we should raise an error when the overlap is less than e.g. half the number_of_samples.
2. In channelReadoutWindowCuttter we round both the trigger_time_sample and samples_before_trigger to integers, but in in simulation.py when we treat the non-trigger channels, continuous time shifts are allowed. This means that the relative timing between antennas can be up to 1.5 samples wrong if the trigger time and/or pre-trigger times are not multiples of 1/sampling_rate. I believe it can be solved by calculating the residual time offset in channelReadoutWindowCuttter and adding it to the final trace start time.

Are these issues big enough that we should solve them, or can we leave it as it is?

[fschlueter]

I have two remaining concerns about the channelReadoutWindowCuttter:

1. It may not be obvious to a user that the function should not be used when the readout window is not almost completely inside the simulated trace. Now, the function will fail if the window is completely outside the trace, but maybe we should raise an error when the overlap is less than e.g. half the number_of_samples.
2. In channelReadoutWindowCuttter we round both the trigger_time_sample and samples_before_trigger to integers, but in in simulation.py when we treat the non-trigger channels, continuous time shifts are allowed. This means that the relative timing between antennas can be up to half a sample wrong if the trigger time and/or pre-trigger times are not multiples of 1/sampling_rate. I believe it can be solved by calculating the residual time offset in channelReadoutWindowCuttter and adding it to the final trace start time.

Are these issues big enough that we should solve them, or can we leave it as it is?

@MartinRavn, to you first point: I think we discussed such a check and I am agree we can implement it. And to your second point, will instead of using trigger_time - pre_trigger_time using trigger_time - int(pre_trigger_time * sampling_rate) / sampling_rate for the trace_start_time solve the problem?

[MartinRavn]

We also need to take into account that trigger_time_channel has been rounded, so we could set the trace_start_time to:

int(np.round((trigger_time - channel.get_trace_start_time()) * sampling_rate)) / sampling_rate - int(pre_trigger_time * sampling_rate) / sampling_rate

What I came up with was to calculate the residual time offset:

residual_time_offset = channel.get_times()[trigger_time_sample - samples_before_trigger] - (trigger_time - pre_trigger_time)

and add it to the trace_start_time at the end.

[fschlueter]

Are you sure this is necessary? It is only important that t_trigger_in_channel_trace = trig.get_trigger_time() - channel.get_trace_start_time() gives us the correct "location", it is not important that this location coincides with a sample.

[MartinRavn]

I would say it is necessary. We set:
channel.set_trace_start_time(trigger_time - pre_trigger_time)
Which will not be the true trace start time since we took the bin nearest to trigger_time - pre_trigger_time

[MartinRavn]

The validate_separate_trigger_channels test fails if I run it with a pre-trigger time which is not a multiple of 1/sampling_rate. I did it by adding pre_trigger_times = 55 * units.ns + 1.17341 * units.ns to this line. However, none of the solutions I suggested solves the problem, so the test maybe fails for unrelated reasons.

[fschlueter]

What is the error?

[MartinRavn]

I get this error: