Fix ara trigger

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

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

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

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.
@@ -142,7 +188,7 @@ def has_triggered(self, channel):
                                             max_freq=1000 * units.MHz,
                                             n_samples=10000,
                                             sampling_rate=channel.get_sampling_rate(),
-                                            amplitude=20.*units.mV,
+                                            amplitude=20 * units.mV,
                                             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,49 +247,59 @@ 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_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:
-                has_triggered = True
-                trigger_time_sample = min(trigger_times)
-                break
+        trigger_time = None
+
+        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,
                                          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_time = 0
+            trigger.set_trigger_time(trigger_time)
         else:
             trigger.set_triggered(True)
-            trigger.set_trigger_time(trigger_time_sample / sampling_rate)
+            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)
 

NuRadioReco/modules/phasedarray/triggerSimulator.py

@@ -231,7 +231,6 @@ 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:

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