From 046c072d8e6de6710c045956fcf097adb74e9e48 Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Mon, 4 Nov 2019 14:32:24 +0100 Subject: [PATCH 01/11] Changing the preferences for the secondary channels. None indicates no secondary channels from now on. Corrected a bug when checking for the positions of secondary channels. --- NuRadioReco/examples/PhasedArray/SNR_curves/T02RunSNR.py | 3 ++- NuRadioReco/modules/phasedarray/triggerSimulator.py | 2 -- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/NuRadioReco/examples/PhasedArray/SNR_curves/T02RunSNR.py b/NuRadioReco/examples/PhasedArray/SNR_curves/T02RunSNR.py index 5829b28e..582881e0 100644 --- a/NuRadioReco/examples/PhasedArray/SNR_curves/T02RunSNR.py +++ b/NuRadioReco/examples/PhasedArray/SNR_curves/T02RunSNR.py @@ -160,7 +160,8 @@ def _detector_simulation(self): triggered_channels=None, # run trigger on all channels secondary_channels=[0,1,3,4,6,7], # secondary channels trigger_name='primary_and_secondary_phasing', - set_not_triggered=(not self._station.has_triggered("simple_threshold"))) + set_not_triggered=(not self._station.has_triggered("simple_threshold")), + ref_index=1.78) if (has_triggered and not reject_event): print('Trigger for SNR', SNRs[iSNR]) diff --git a/NuRadioReco/modules/phasedarray/triggerSimulator.py b/NuRadioReco/modules/phasedarray/triggerSimulator.py index 857e3c4a..61dd79b1 100644 --- a/NuRadioReco/modules/phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/phasedarray/triggerSimulator.py @@ -231,11 +231,9 @@ def run(self, evt, station, det, beam_rolls = self.get_beam_rolls(station, det, triggered_channels, phasing_angles, ref_index=ref_index) empty_rolls = [ {} for direction in range(len(phasing_angles)) ] logger.debug("secondary_channels:", secondary_channels) - if (len(secondary_channels) == 0): only_primary = True else: - only_primary = False secondary_beam_rolls = self.get_beam_rolls(station, det, secondary_channels, secondary_phasing_angles, ref_index=ref_index) if only_primary: From 03e8b83efb65777113f221da095fe6a6d5d332ae Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Tue, 12 Nov 2019 17:39:10 +0100 Subject: [PATCH 02/11] Fixing the ARA diode trigger. New method for calculating the threshold noise parameters --- NuRadioReco/modules/ARA/triggerSimulator.py | 88 +++++++++++++++------ 1 file changed, 65 insertions(+), 23 deletions(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index 38205744..49f8f055 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -113,6 +113,52 @@ def tunnel_diode(self, channel): return trace_after_tunnel_diode + def calculate_noise_parameters(self, + sampling_rate = 1 * units.GHz, + min_freq = 50 * units.MHz, + max_freq = 1 * units.GHz, + amplitude = 10 * units.microvolt, + type='rayleigh'): + """ + Calculates the mean and the standard deviation for the diode-filtered noise. + + Parameters + ---------- + sampling_rate: float + Sampling rate + min_freq: float + Minimum frequency of the bandwidth + max_freq: float + Maximum frequency of the bandwidth + amplitude: float + Voltage amplitude (RMS) for the noise + type: string + Noise type + + Returns + ------- + power_mean: float + Mean of the diode-filtered noise + power_std: float + Standard deviation of the diode-filtered noise + """ + noise = NuRadioReco.framework.channel.Channel(0) + + long_noise = channelGenericNoiseAdder().bandlimited_noise(min_freq=min_freq, + max_freq=max_freq, + n_samples=10000, + sampling_rate=sampling_rate, + amplitude=amplitude, + type=type) + + noise.set_trace(long_noise, sampling_rate) + power_noise = self.tunnel_diode(noise) + + power_mean = np.mean(power_noise) + power_std = np.std(power_noise) + + return power_mean, power_std + def has_triggered(self, channel): """ Check if the detector system triggers on a given channel. @@ -138,11 +184,11 @@ def has_triggered(self, channel): # than taken the full ARA signal chain noise = NuRadioReco.framework.channel.Channel(0) - long_noise = channelGenericNoiseAdder().bandlimited_noise(min_freq=50 * units.MHz, - max_freq=1000 * units.MHz, + long_noise = channelGenericNoiseAdder().bandlimited_noise(min_freq=130 * units.MHz, + max_freq=1500 * units.MHz, n_samples=10000, sampling_rate=channel.get_sampling_rate(), - amplitude=20.*units.mV, + amplitude=16.70*units.microvolt, type='perfect_white') noise.set_trace(long_noise, channel.get_sampling_rate()) @@ -156,6 +202,7 @@ def has_triggered(self, channel): after_tunnel_diode = self.tunnel_diode(channel) low_trigger = (self._power_mean - self._power_std * np.abs(self.power_threshold)) + return np.min(after_tunnel_diode) < low_trigger @register_run() @@ -200,47 +247,42 @@ def run(self, evt, station, det, # No coincidence requirement yet trigger = {} + trigger_times = [] for channel in station.iter_channels(): channel_id = channel.get_id() if channel_id not in triggered_channels: continue trigger[channel_id] = self.has_triggered(channel) + if trigger[channel_id]: + times = channel.get_times() + trace_after_diode = self.tunnel_diode(channel) + arg_trigger = np.argwhere( trace_after_diode == np.min(trace_after_diode) )[0,0] + trigger_times.append(times[arg_trigger]) has_triggered = False trigger_time_sample = None - # loop over the trace with a sliding window of "coinc_window" - coinc_window_samples = np.int(np.round(coinc_window * channel.get_sampling_rate())) - trace_length = len(station.get_channel(0).get_trace()) - sampling_rate = station.get_channel(0).get_sampling_rate() - - for i in range(0, trace_length - coinc_window_samples): - istop = i + coinc_window_samples - coinc = 0 - trigger_times = [] - for iCh, tr in trigger.items(): # loops through triggers of each channel - tr = np.array(tr) - mask_trigger_in_coind_window = (tr >= i) & (tr < istop) - if(np.sum(mask_trigger_in_coind_window)): - coinc += 1 - trigger_times.append(tr[mask_trigger_in_coind_window][0]) # save time/sample of first trigger in - - if coinc >= number_concidences: + + trace_times = station.get_channel(0).get_times() + trigger_times = np.array(trigger_times) + for trace_time in trace_times: + if ( np.sum( np.abs(trace_time-trigger_times) <= coinc_window/2 ) >= number_concidences ): has_triggered = True - trigger_time_sample = min(trigger_times) + trigger_time_sample = np.min(trigger_times) break trigger = IntegratedPowerTrigger(trigger_name, power_threshold, coinc_window, channels=triggered_channels, number_of_coincidences=number_concidences, power_mean=self._power_mean, power_std=self._power_std) + if not has_triggered: trigger.set_triggered(False) logger.info("Station has NOT passed trigger") trigger_time_sample = 0 - station.get_trigger().set_trigger_time(trigger_time_sample / sampling_rate) + trigger.set_trigger_time(trigger_time_sample) else: trigger.set_triggered(True) - trigger.set_trigger_time(trigger_time_sample / sampling_rate) + trigger.set_trigger_time(trigger_time_sample) logger.info("Station has passed trigger, trigger time is {:.1f} ns (sample {})".format( trigger.get_trigger_time() / units.ns, trigger_time_sample)) From 31cd994fad15e4e1b38005fe57f0e4aac197e118 Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Tue, 12 Nov 2019 17:44:00 +0100 Subject: [PATCH 03/11] Change log --- changelog.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/changelog.txt b/changelog.txt index 62beeca5..1bfa7a35 100644 --- a/changelog.txt +++ b/changelog.txt @@ -20,5 +20,6 @@ Detector description can be stored in .nur files bugfixes: -Fixes increase in filesize caused by switch to python3 -Fixed bug when using no secondary channels on the phased array +-Fixed bug when using ARA trigger version 1.0.0 - 2019/08/30 - first python 3 release From 0c20978415c4f75f533418152c9664dfb00aeb5f Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Tue, 12 Nov 2019 17:57:24 +0100 Subject: [PATCH 04/11] Skip trace bins --- NuRadioReco/modules/ARA/triggerSimulator.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index 49f8f055..748dd4dd 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -264,7 +264,9 @@ def run(self, evt, station, det, trace_times = station.get_channel(0).get_times() trigger_times = np.array(trigger_times) - for trace_time in trace_times: + slice_left = int(coinc_window/2/(trace_times[1]-trace_times[0])) + slice_right = len(trace_times)-slice_left + for trace_time in trace_times[slice_left:slice_right]: if ( np.sum( np.abs(trace_time-trigger_times) <= coinc_window/2 ) >= number_concidences ): has_triggered = True trigger_time_sample = np.min(trigger_times) From 2e9a9ddb8f719b5918901883037e0550cc93ad50 Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Tue, 12 Nov 2019 18:02:53 +0100 Subject: [PATCH 05/11] Back to default parameters --- NuRadioReco/modules/ARA/triggerSimulator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index 748dd4dd..4630f243 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -184,11 +184,11 @@ def has_triggered(self, channel): # than taken the full ARA signal chain noise = NuRadioReco.framework.channel.Channel(0) - long_noise = channelGenericNoiseAdder().bandlimited_noise(min_freq=130 * units.MHz, - max_freq=1500 * units.MHz, + long_noise = channelGenericNoiseAdder().bandlimited_noise(min_freq=50 * units.MHz, + max_freq=1000 * units.MHz, n_samples=10000, sampling_rate=channel.get_sampling_rate(), - amplitude=16.70*units.microvolt, + amplitude=20 * units.mV, type='perfect_white') noise.set_trace(long_noise, channel.get_sampling_rate()) From 19aec33efff3a837ea87e4ff47caabed0d2d2d20 Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Tue, 12 Nov 2019 18:09:11 +0100 Subject: [PATCH 06/11] Fixed bug when accessing only_primary --- NuRadioReco/modules/phasedarray/triggerSimulator.py | 1 + 1 file changed, 1 insertion(+) diff --git a/NuRadioReco/modules/phasedarray/triggerSimulator.py b/NuRadioReco/modules/phasedarray/triggerSimulator.py index 61dd79b1..1e9ebeee 100644 --- a/NuRadioReco/modules/phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/phasedarray/triggerSimulator.py @@ -234,6 +234,7 @@ def run(self, evt, station, det, if (len(secondary_channels) == 0): only_primary = True else: + only_primary = False secondary_beam_rolls = self.get_beam_rolls(station, det, secondary_channels, secondary_phasing_angles, ref_index=ref_index) if only_primary: From f8c1d8a682ae8cead20335b3afad680fad228b16 Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Tue, 12 Nov 2019 18:19:08 +0100 Subject: [PATCH 07/11] Removed old NuRadioMC units dependences --- .../PhasedArray/Effective_volume/T01generate_event_list.py | 2 +- .../examples/PhasedArray/SNR_curves/T01generate_event_list.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NuRadioReco/examples/PhasedArray/Effective_volume/T01generate_event_list.py b/NuRadioReco/examples/PhasedArray/Effective_volume/T01generate_event_list.py index bf534599..f366f3a5 100644 --- a/NuRadioReco/examples/PhasedArray/Effective_volume/T01generate_event_list.py +++ b/NuRadioReco/examples/PhasedArray/Effective_volume/T01generate_event_list.py @@ -13,7 +13,7 @@ """ from __future__ import absolute_import, division, print_function -from NuRadioMC.utilities import units +from NuRadioReco.utilities import units from NuRadioMC.EvtGen.generator import generate_eventlist_cylinder import numpy as np import os diff --git a/NuRadioReco/examples/PhasedArray/SNR_curves/T01generate_event_list.py b/NuRadioReco/examples/PhasedArray/SNR_curves/T01generate_event_list.py index b9bc945b..247d56fc 100644 --- a/NuRadioReco/examples/PhasedArray/SNR_curves/T01generate_event_list.py +++ b/NuRadioReco/examples/PhasedArray/SNR_curves/T01generate_event_list.py @@ -11,7 +11,7 @@ """ from __future__ import absolute_import, division, print_function -from NuRadioMC.utilities import units +from NuRadioReco.utilities import units from NuRadioMC.EvtGen.generator import generate_eventlist_cylinder import numpy as np import os From 0dd7f4061a8fda860dca6c6e2b99c641e9cf66b4 Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Wed, 13 Nov 2019 10:34:54 +0100 Subject: [PATCH 08/11] Changed trigger_time_sample to trigger_time --- NuRadioReco/modules/ARA/triggerSimulator.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index 4630f243..4dba10d1 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -260,7 +260,7 @@ def run(self, evt, station, det, trigger_times.append(times[arg_trigger]) has_triggered = False - trigger_time_sample = None + trigger_time = None trace_times = station.get_channel(0).get_times() trigger_times = np.array(trigger_times) @@ -269,7 +269,7 @@ def run(self, evt, station, det, for trace_time in trace_times[slice_left:slice_right]: if ( np.sum( np.abs(trace_time-trigger_times) <= coinc_window/2 ) >= number_concidences ): has_triggered = True - trigger_time_sample = np.min(trigger_times) + trigger_time = np.min(trigger_times) break trigger = IntegratedPowerTrigger(trigger_name, power_threshold, @@ -280,13 +280,13 @@ def run(self, evt, station, det, if not has_triggered: trigger.set_triggered(False) logger.info("Station has NOT passed trigger") - trigger_time_sample = 0 - trigger.set_trigger_time(trigger_time_sample) + trigger_time = 0 + trigger.set_trigger_time(trigger_time) else: trigger.set_triggered(True) - trigger.set_trigger_time(trigger_time_sample) + trigger.set_trigger_time(trigger_time) logger.info("Station has passed trigger, trigger time is {:.1f} ns (sample {})".format( - trigger.get_trigger_time() / units.ns, trigger_time_sample)) + trigger.get_trigger_time() / units.ns, trigger_time)) station.set_trigger(trigger) From c6f932ee81f4df451dda8fcdef8e8be0454a610c Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Wed, 13 Nov 2019 17:25:55 +0100 Subject: [PATCH 09/11] Using argmin --- NuRadioReco/modules/ARA/triggerSimulator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index 4dba10d1..e13d09d3 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -256,7 +256,7 @@ def run(self, evt, station, det, if trigger[channel_id]: times = channel.get_times() trace_after_diode = self.tunnel_diode(channel) - arg_trigger = np.argwhere( trace_after_diode == np.min(trace_after_diode) )[0,0] + arg_trigger = np.argmin(trace_after_diode) trigger_times.append(times[arg_trigger]) has_triggered = False From 7065ed8168a60ceb4243bb583a7f7f36b01570fe Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Wed, 13 Nov 2019 17:52:04 +0100 Subject: [PATCH 10/11] Calculating the correct total time trace. Now we take the minimum and the maximum times from all the channels. --- NuRadioReco/modules/ARA/triggerSimulator.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index e13d09d3..a4ceed3c 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -262,7 +262,17 @@ def run(self, evt, station, det, has_triggered = False trigger_time = None - trace_times = station.get_channel(0).get_times() + times_min = [] + times_max = [] + sampling_rates = [] + for channel in station.iter_channels(): + times_min.append(np.min(channel.get_times())) + times_max.append(np.max(channel.get_times())) + sampling_rates.append(channel.get_sampling_rate()) + + trace_times = np.arange(np.min(times_min), np.max(times_max), + 1/np.min(sampling_rates)) + trigger_times = np.array(trigger_times) slice_left = int(coinc_window/2/(trace_times[1]-trace_times[0])) slice_right = len(trace_times)-slice_left From 3b0565451a4e3e4829a7900e021a20c50008be8b Mon Sep 17 00:00:00 2001 From: daniel-zid Date: Wed, 13 Nov 2019 18:09:50 +0100 Subject: [PATCH 11/11] Moving loop into another loop. Implementing a new condition to save time. If the number of triggers is less than the required coincidences, the coincidence window is not calculated --- NuRadioReco/modules/ARA/triggerSimulator.py | 41 +++++++++++---------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/NuRadioReco/modules/ARA/triggerSimulator.py b/NuRadioReco/modules/ARA/triggerSimulator.py index a4ceed3c..31ce0f8b 100644 --- a/NuRadioReco/modules/ARA/triggerSimulator.py +++ b/NuRadioReco/modules/ARA/triggerSimulator.py @@ -248,39 +248,42 @@ def run(self, evt, station, det, # No coincidence requirement yet trigger = {} trigger_times = [] + times_min = [] + times_max = [] + sampling_rates = [] + number_triggered_channels = 0 + for channel in station.iter_channels(): channel_id = channel.get_id() if channel_id not in triggered_channels: continue trigger[channel_id] = self.has_triggered(channel) if trigger[channel_id]: + number_triggered_channels += 1 times = channel.get_times() trace_after_diode = self.tunnel_diode(channel) arg_trigger = np.argmin(trace_after_diode) trigger_times.append(times[arg_trigger]) + times_min.append(np.min(times)) + times_max.append(np.max(times)) + sampling_rates.append(channel.get_sampling_rate()) has_triggered = False trigger_time = None - times_min = [] - times_max = [] - sampling_rates = [] - for channel in station.iter_channels(): - times_min.append(np.min(channel.get_times())) - times_max.append(np.max(channel.get_times())) - sampling_rates.append(channel.get_sampling_rate()) - - trace_times = np.arange(np.min(times_min), np.max(times_max), - 1/np.min(sampling_rates)) - - trigger_times = np.array(trigger_times) - slice_left = int(coinc_window/2/(trace_times[1]-trace_times[0])) - slice_right = len(trace_times)-slice_left - for trace_time in trace_times[slice_left:slice_right]: - if ( np.sum( np.abs(trace_time-trigger_times) <= coinc_window/2 ) >= number_concidences ): - has_triggered = True - trigger_time = np.min(trigger_times) - break + if (number_triggered_channels >= number_concidences): + + trace_times = np.arange(np.min(times_min), np.max(times_max), + 1/np.min(sampling_rates)) + + trigger_times = np.array(trigger_times) + slice_left = int(coinc_window/2/(trace_times[1]-trace_times[0])) + slice_right = len(trace_times)-slice_left + for trace_time in trace_times[slice_left:slice_right]: + if ( np.sum( np.abs(trace_time-trigger_times) <= coinc_window/2 ) >= number_concidences ): + has_triggered = True + trigger_time = np.min(trigger_times) + break trigger = IntegratedPowerTrigger(trigger_name, power_threshold, coinc_window, channels=triggered_channels,