From 697c0af0ac25cc72c9665c2c7ac5b72d65fe8478 Mon Sep 17 00:00:00 2001 From: Ryan Krebs Date: Fri, 20 Dec 2024 21:30:43 -0500 Subject: [PATCH 1/7] update phased trigger and digitizer bug fixes --- NuRadioMC/simulation/simulation.py | 9 - NuRadioMC/utilities/runner.py | 2 + .../RNO_G/RNO_single_station_only_PA.json | 8 +- NuRadioReco/detector/detector_base.py | 16 + .../RNO_G/hardwareResponseIncorporator.py | 19 +- .../modules/RNO_G/triggerBoardResponse.py | 124 ++++++- .../modules/analogToDigitalConverter.py | 6 +- .../modules/phasedarray/triggerSimulator.py | 339 +++++++++++++----- NuRadioReco/utilities/signal_processing.py | 13 +- 9 files changed, 399 insertions(+), 137 deletions(-) diff --git a/NuRadioMC/simulation/simulation.py b/NuRadioMC/simulation/simulation.py index d0f69bac2e..99e674f864 100644 --- a/NuRadioMC/simulation/simulation.py +++ b/NuRadioMC/simulation/simulation.py @@ -1152,15 +1152,6 @@ def __init__(self, inputfilename, self.detector_simulation_part2 = getattr(self, "_detector_simulation_part2", None) self.detector_simulation_filter_amp = getattr(self, "_detector_simulation_filter_amp", None) self.detector_simulation_trigger = getattr(self, "_detector_simulation_trigger", None) - if((self.detector_simulation_filter_amp is None) ^ (self.detector_simulation_trigger is None)): - logger.error("Please provide both detector_simulation_filter_amp and detector_simulation_trigger or detector_simulation_part1 and detector_simulation_part2") - raise ValueError("Please provide both detector_simulation_filter_amp and detector_simulation_trigger or detector_simulation_part1 and detector_simulation_part2") - if((self.detector_simulation_part1 is None) ^ (self.detector_simulation_part2 is None)): - logger.error("Please provide both detector_simulation_filter_amp and detector_simulation_trigger or detector_simulation_part1 and detector_simulation_part2") - raise ValueError("Please provide both detector_simulation_filter_amp and detector_simulation_trigger or detector_simulation_part1 and detector_simulation_part2") - if((self.detector_simulation_part1 is not None) and (self.detector_simulation_filter_amp is not None)): - logger.error("Please provide either detector_simulation_filter_amp and detector_simulation_trigger or detector_simulation_part1 and detector_simulation_part2") - raise ValueError("Please provide either detector_simulation_filter_amp and detector_simulation_trigger or detector_simulation_part1 and detector_simulation_part2") # Initialize detector self._antenna_pattern_provider = antennapattern.AntennaPatternProvider() diff --git a/NuRadioMC/utilities/runner.py b/NuRadioMC/utilities/runner.py index 2182c99bef..124c244cf9 100644 --- a/NuRadioMC/utilities/runner.py +++ b/NuRadioMC/utilities/runner.py @@ -70,6 +70,8 @@ def run(self): else: if(not self.q.empty()): n_trig = self.q.get_nowait() + if n_trig is None: + n_trig=0 self.n_triggers += n_trig print(f"{iN} has finished with {n_trig} events, total number of triggered events is {self.n_triggers}", flush=True) outputfilename = self.get_outputfilename() diff --git a/NuRadioReco/detector/RNO_G/RNO_single_station_only_PA.json b/NuRadioReco/detector/RNO_G/RNO_single_station_only_PA.json index 4d4f7a3a0b..0c39da45f1 100644 --- a/NuRadioReco/detector/RNO_G/RNO_single_station_only_PA.json +++ b/NuRadioReco/detector/RNO_G/RNO_single_station_only_PA.json @@ -15,7 +15,7 @@ "ant_type": "RNOG_vpol_4inch_center_n1.73", "trigger_adc_sampling_frequency" : 0.472, "trigger_adc_nbits": 8, - "trigger_adc_noise_nbits": 3.321, + "trigger_adc_noise_nbits": 2.585, "trigger_amp_type": "ULP_216", "cab_time_delay": 0.0, "channel_id": 0, @@ -34,7 +34,7 @@ "ant_type": "RNOG_vpol_4inch_center_n1.73", "trigger_adc_sampling_frequency" : 0.472, "trigger_adc_nbits": 8, - "trigger_adc_noise_nbits": 3.321, + "trigger_adc_noise_nbits": 2.585, "trigger_amp_type": "ULP_216", "reference_channel": 0, "channel_id": 1, @@ -53,7 +53,7 @@ "ant_type": "RNOG_vpol_4inch_center_n1.73", "trigger_adc_sampling_frequency" : 0.472, "trigger_adc_nbits": 8, - "trigger_adc_noise_nbits": 3.321, + "trigger_adc_noise_nbits": 2.585, "trigger_amp_type": "ULP_216", "reference_channel": 0, "channel_id": 2, @@ -72,7 +72,7 @@ "ant_type": "RNOG_vpol_4inch_center_n1.73", "trigger_adc_sampling_frequency" : 0.472, "trigger_adc_nbits": 8, - "trigger_adc_noise_nbits": 3.321, + "trigger_adc_noise_nbits": 2.585, "trigger_amp_type": "ULP_216", "reference_channel": 0, "channel_id": 3, diff --git a/NuRadioReco/detector/detector_base.py b/NuRadioReco/detector/detector_base.py index dcc93488ae..8aabc1cf37 100644 --- a/NuRadioReco/detector/detector_base.py +++ b/NuRadioReco/detector/detector_base.py @@ -875,6 +875,22 @@ def get_amplifier_response(self, station_id, channel_id, frequencies): amp_phase = amp_response_functions['phase'](frequencies) return amp_gain * amp_phase + def get_trigger_amplifier_type(self, station_id, channel_id): + """ + returns the type of the amplifier + + Parameters + ---------- + station_id: int + the station id + channel_id: int + the channel id + + Returns string + """ + res = self.__get_channel(station_id, channel_id) + return res['trigger_amp_type'] + def get_sampling_frequency(self, station_id, channel_id): """ returns the sampling frequency diff --git a/NuRadioReco/modules/RNO_G/hardwareResponseIncorporator.py b/NuRadioReco/modules/RNO_G/hardwareResponseIncorporator.py index 25edada741..95d16e1b9b 100644 --- a/NuRadioReco/modules/RNO_G/hardwareResponseIncorporator.py +++ b/NuRadioReco/modules/RNO_G/hardwareResponseIncorporator.py @@ -101,6 +101,12 @@ def get_filter(self, frequencies, station_id, channel_id, det, amp_response = analog_components.load_amp_response(amp_type) amp_response = amp_response['gain']( frequencies, temp) * amp_response['phase'](frequencies) + if is_trigger: + trig_amp_type = det.get_trigger_amplifier_type(station_id, channel_id) + trig_amp_response = analog_components.load_amp_response(trig_amp_type) + trig_amp_response= trig_amp_response['gain']( + frequencies, temp) * trig_amp_response['phase'](frequencies) + amp_response=amp_response*trig_amp_response else: raise NotImplementedError("Detector type not implemented") @@ -184,9 +190,6 @@ def run(self, evt, station, det, temp=293.15, sim_to_data=False, phase_only=Fals t = time.time() - if self.trigger_channels is not None and not isinstance(det, detector.rnog_detector.Detector): - raise ValueError("Simulating extra trigger channels is only possible with the `rnog_detector.Detector` class.") - has_trigger_channels = False for channel in station.iter_channels(): frequencies = channel.get_frequencies() @@ -280,22 +283,22 @@ def get_mingainlin(self): file_dir = os.path.dirname(__file__) detectorfile = os.path.join( - file_dir, "../../detector/RNO_G/RNO_single_station.json") + file_dir, "../../detector/RNO_G/RNO_single_station_only_PA.json") det_old = detector.generic_detector.GenericDetector( json_filename=detectorfile, default_station=11, antenna_by_depth=False) det = detector.rnog_detector.Detector(log_level=logging.DEBUG, over_write_handset_values={ - "sampling_frequency": 2.4 * units.GHz}, always_query_entire_description=True) + "sampling_frequency": 2.4 * units.GHz}, always_query_entire_description=False) det.update(datetime.datetime(2022, 8, 2, 0, 0)) hri = hardwareResponseIncorporator() - + hri.begin(trigger_channels=[0,1,2,3]) frequencies = np.linspace(0, 1) * units.GHz filter_old = hri.get_filter( - frequencies, station_id=11, channel_id=0, det=det_old, sim_to_data=True) + frequencies, station_id=11, channel_id=0, det=det_old, sim_to_data=True,is_trigger=True) filter = hri.get_filter(frequencies, station_id=11, - channel_id=0, det=det, sim_to_data=True) + channel_id=0, det=det, sim_to_data=True,is_trigger=True) fig, ax = plt.subplots() diff --git a/NuRadioReco/modules/RNO_G/triggerBoardResponse.py b/NuRadioReco/modules/RNO_G/triggerBoardResponse.py index 6e2eea31df..62c07c3f29 100644 --- a/NuRadioReco/modules/RNO_G/triggerBoardResponse.py +++ b/NuRadioReco/modules/RNO_G/triggerBoardResponse.py @@ -3,7 +3,8 @@ from NuRadioReco.modules.base.module import register_run from NuRadioReco.modules.analogToDigitalConverter import analogToDigitalConverter -from NuRadioReco.utilities import units +from NuRadioReco.utilities import units, fft +#from scipy.signal import resample, firwin, hilbert logger = logging.getLogger("NuRadioReco.triggerBoardResponse") @@ -74,7 +75,7 @@ def apply_trigger_filter(self, station, trigger_channels, trigger_filter): filt = trigger_filter(freqs) channel.set_frequency_spectrum(channel.get_frequency_spectrum() * filt, channel.get_sampling_rate()) - def get_avg_vrms(self, station, trigger_channels, trace_split=20): + def get_vrms(self, station, trigger_channels, trace_split=20): """ Estimates the RMS voltage of the triggering antennas by splitting the waveforms into chunks and taking the median of standard deviation of the chunks @@ -96,18 +97,17 @@ def get_avg_vrms(self, station, trigger_channels, trace_split=20): """ - avg_vrms = 0 + vrms=[] for channel_id in trigger_channels: channel = station.get_trigger_channel(channel_id) trace = np.array(channel.get_trace()) trace = trace[: int(trace_split * int(len(trace) / trace_split))].reshape((trace_split, -1)) approx_vrms = np.median(np.std(trace, axis=1)) logger.debug(f" Ch: {channel_id}\tObser Vrms: {approx_vrms / units.mV:0.3f} mV") - avg_vrms += approx_vrms + vrms.append(approx_vrms) - avg_vrms /= len(trigger_channels) - self.logger.debug(f"Average Vrms: {avg_vrms / units.mV:0.3f} mV") - return approx_vrms + self.logger.debug(vrms) + return vrms def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): """ @@ -138,7 +138,7 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): """ if avg_vrms is None: - avg_vrms = self.get_avg_vrms(station, trigger_channels) + avg_vrms = self.get_vrms(station, trigger_channels) self.logger.debug("Applying gain at ADC level") @@ -151,8 +151,8 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): noise_bits = det_channel["trigger_adc_noise_nbits"] total_bits = det_channel["trigger_adc_nbits"] - volts_per_adc = self._adc_input_range / (2 ** total_bits - 1) - ideal_vrms = volts_per_adc * (2 ** (noise_bits - 1)) + volts_per_adc = self._adc_input_range / (2 ** total_bits) + ideal_vrms = volts_per_adc * (2 ** (noise_bits) - 1) msg = f"\t Ch: {channel_id}\t Target Vrms: {ideal_vrms / units.mV:0.3f} mV" msg += f"\t V/ADC: {volts_per_adc / units.mV:0.3f} mV" @@ -232,16 +232,17 @@ def run(self, evt, station, det, trigger_channels, vrms=None, apply_adc_gain=Tru self.logger.debug("Applying the RNO-G trigger board response") if vrms is None: - vrms = self.get_avg_vrms(station, trigger_channels) + vrms = self.get_vrms(station, trigger_channels) if apply_adc_gain: trigger_board_vrms, ideal_vrms = self.apply_adc_gain(station, det, trigger_channels, vrms) else: trigger_board_vrms = vrms - ideal_vrms = vrms + ideal_vrms = np.mean(vrms) if digitize_trace: self.digitize_trace(station, det, trigger_channels, ideal_vrms) + trigger_board_vrms=self.get_vrms(station, trigger_channels) return trigger_board_vrms @@ -252,3 +253,102 @@ def end(self): dt = timedelta(seconds=self.__t) self.logger.info("total time used by this module is {}".format(dt)) return dt + +if __name__=='__main__': + + from NuRadioReco.detector import detector + import NuRadioReco.modules.channelGenericNoiseAdder + import NuRadioReco.modules.channelBandPassFilter + from NuRadioReco.modules.RNO_G import hardwareResponseIncorporator + import matplotlib.pyplot as plt + from scipy import constants + + rnogHarwareResponse = hardwareResponseIncorporator.hardwareResponseIncorporator() + rnogHarwareResponse.begin(trigger_channels=[0,1,2,3]) + rnogADCResponse = triggerBoardResponse(log_level=logging.ERROR) + rnogADCResponse.begin(adc_input_range=2 * units.volt, clock_offset=0.0, adc_output="counts") + channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() + + det_file='RNO_G/RNO_single_station_only_PA.json' + det = detector.Detector(source='json',json_filename=det_file) + + station_id=11 + channel_ids = np.arange(4) + + n_samples=1024 + sampling_rate=472*units.MHz + dt = 1 / sampling_rate + ff = np.fft.rfftfreq(n_samples, dt) + max_freq = ff[-1] + min_freq = 0 + fff = np.linspace(min_freq, max_freq, 10000) + + four_filters={} + rf_filter=rnogHarwareResponse.get_filter(ff,station_id,channel_ids[0],det,sim_to_data=True,is_trigger=True) + chain_filter=rf_filter + + for i in range(4): + four_filters[i]=chain_filter + + four_filters_highres={} + rf_filter_highres=rnogHarwareResponse.get_filter(fff,station_id,channel_ids[0],det,sim_to_data=True,is_trigger=True) + chain_filter_highres=rf_filter_highres + + for i in channel_ids: + four_filters_highres[i]=chain_filter_highres + + Vrms = 1 + noise_temp=300 + bandwidth={} + Vrms_ratio = {} + amplitude = {} + per_channel_vrms=[] + for i in channel_ids: + integrated_channel_response = np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) + rel_channel_response=np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) + bandwidth[i]=integrated_channel_response + Vrms_ratio[i] = np.sqrt(rel_channel_response / (max_freq - min_freq)) + chan_vrms=(noise_temp * 50 * constants.k * integrated_channel_response / units.Hz) ** 0.5 + per_channel_vrms.append(chan_vrms) + amplitude[i] = chan_vrms / Vrms_ratio[i] + + station = NuRadioReco.framework.station.Station(station_id) + evt = NuRadioReco.framework.event.Event(0, 0) + + channelGenericNoiseAdder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() + channelGenericNoiseAdder.begin() + + for channel_id in channel_ids: + + spectrum = channelGenericNoiseAdder.bandlimited_noise(min_freq, max_freq, n_samples, sampling_rate, amplitude[channel_id], + type="rayleigh", time_domain=False) + + trace=fft.freq2time(spectrum*four_filters[channel_id], sampling_rate) + + channel = NuRadioReco.framework.channel.Channel(channel_id) + channel.set_trace(trace, sampling_rate) + station.add_trigger_channel(channel) + + fig,ax=plt.subplots(1,2,sharex=True,figsize=(11,7)) + + for channel_id in channel_ids: + ch=station.get_channel(channel_id) + ax[0].plot(ch.get_times(),ch.get_trace(),label='ch %i'%channel_id) + + chan_rms=rnogADCResponse.run(evt,station,det,requested_channels=channel_ids, + digitize_trace=True,apply_adc_gain=True) + + for channel_id in channel_ids: + ch=station.get_channel(channel_id) + ax[1].plot(ch.get_times(),ch.get_trace(),label='ch %i'%channel_id) + + ax[0].set_title('Raw Voltage Trace') + ax[0].set_ylabel('Voltage [V]') + ax[0].set_xlabel('Samples') + ax[0].legend(loc='upper right') + ax[1].set_xlabel('Samples') + ax[1].set_ylabel('Voltage [ADC]') + ax[1].legend(loc='upper right') + ax[1].set_title('Gain Eq. and Digitized Trace') + fig.tight_layout() + plt.show() diff --git a/NuRadioReco/modules/analogToDigitalConverter.py b/NuRadioReco/modules/analogToDigitalConverter.py index b7482906be..911a375ad9 100644 --- a/NuRadioReco/modules/analogToDigitalConverter.py +++ b/NuRadioReco/modules/analogToDigitalConverter.py @@ -36,9 +36,7 @@ def perfect_comparator(trace, adc_n_bits, adc_ref_voltage, mode='floor', output= digital_trace: array of floats Digitised voltage trace in volts or ADC counts """ - - lsb_voltage = adc_ref_voltage / (2 ** (adc_n_bits - 1) - 1) - + lsb_voltage = adc_ref_voltage / (2 ** (adc_n_bits) - 1) if (mode == 'floor'): digital_trace = np.floor(trace / lsb_voltage) elif (mode == 'ceiling'): @@ -290,7 +288,7 @@ def get_digital_trace(self, station, det, channel, adc_ref_voltage = det_channel[adc_ref_voltage_label] * units.V else: - adc_ref_voltage = Vrms * (2 ** (adc_n_bits - 1) - 1) / (2 ** (adc_noise_n_bits - 1) - 1) + adc_ref_voltage = Vrms * (2 ** (adc_n_bits) - 1) / (2 ** (adc_noise_n_bits ) - 1) if(adc_sampling_frequency > channel.get_sampling_rate()): error_msg = 'The ADC sampling rate is greater than ' diff --git a/NuRadioReco/modules/phasedarray/triggerSimulator.py b/NuRadioReco/modules/phasedarray/triggerSimulator.py index 8d2b1d6acb..fc5dc4fcdc 100644 --- a/NuRadioReco/modules/phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/phasedarray/triggerSimulator.py @@ -7,8 +7,9 @@ import scipy import numpy as np from scipy import constants +from scipy.signal import hilbert -logger = logging.getLogger('NuRadioReco.phasedTriggerSimulator') +logger = logging.getLogger('phasedTriggerSimulator') cspeed = constants.c * units.m / units.s @@ -20,14 +21,14 @@ class triggerSimulator: """ Calculates the trigger for a phased array with a primary beam. - + The channels that participate in both beams and the pointing angle for each subbeam can be specified. See https://arxiv.org/pdf/1809.04573.pdf """ - def __init__(self, log_level=logging.NOTSET): + def __init__(self, log_level=logging.WARNING): self.__t = 0 self.__pre_trigger_time = None self.__debug = None @@ -98,7 +99,7 @@ def calculate_time_delays(self, station, det, """ if(triggered_channels is None): - triggered_channels = [channel.get_id() for channel in station.iter_channels()] + triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] time_step = 1. / sampling_frequency @@ -109,7 +110,7 @@ def calculate_time_delays(self, station, det, # Need to add in delay for trigger delay cable_delays = {} - for channel in station.iter_channels(use_channels=triggered_channels): + for channel in station.iter_trigger_channels(use_channels=triggered_channels): cable_delays[channel.get_id()] = det.get_cable_delay(station.get_id(), channel.get_id()) beam_rolls = [] @@ -152,7 +153,7 @@ def get_channel_trace_start_time(self, station, triggered_channels): """ channel_trace_start_time = None - for channel in station.iter_channels(use_channels=triggered_channels): + for channel in station.iter_trigger_channels(use_channels=triggered_channels): if channel_trace_start_time is None: channel_trace_start_time = channel.get_trace_start_time() elif channel_trace_start_time != channel.get_trace_start_time(): @@ -186,7 +187,7 @@ def check_vertical_string(self, station, det, triggered_channels): if (sum(diff_x) > cut or sum(diff_y) > cut): raise NotImplementedError('The phased triggering array should lie on a vertical line') - def power_sum(self, coh_sum, window, step, adc_output='voltage'): + def power_sum(self, coh_sum, window, step, adc_output='voltage',rnog_like=False): """ Calculate power summed over a length defined by 'window', overlapping at intervals defined by 'step' @@ -223,14 +224,53 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage'): if(adc_output == 'voltage'): coh_sum_squared = (coh_sum * coh_sum).astype(float) - elif(adc_output == 'counts'): + else: #(adc_output == 'counts'): + if rnog_like: + coh_sum[coh_sum>2**6-1]=2**6-1 + coh_sum[coh_sum<-2**6]=-2**6 coh_sum_squared = (coh_sum * coh_sum).astype(int) coh_sum_windowed = np.lib.stride_tricks.as_strided(coh_sum_squared, (num_frames, window), (coh_sum_squared.strides[0] * step, coh_sum_squared.strides[0])) + power = np.sum(coh_sum_windowed, axis=1) + return_power=power.astype(float) / window + + if adc_output=='counts': return_power=np.floor(return_power) + + return return_power, num_frames + + def hilbert_envelope(self,coh_sum,adc_output='voltage',coeff_gain=1,rnog_like=False): - return power.astype(float) / window, num_frames + if rnog_like: + coeff_gain=128 + + #31 sample fir transformer + #hil=[ -0.0424413 , 0. , -0.0489708 , 0. , -0.0578745 , 0. , -0.0707355 , 0. , -0.0909457 , 0. , -0.127324 , 0. + # , -0.2122066 , 0. , -0.6366198 , 0., 0.6366198 , 0. , 0.2122066 , 0. , 0.127324 ,0. , 0.0909457 , 0. , .0707355 + # , 0. , 0.0578745 , 0. , 0.0489708 , 0. , 0.0424413 ] + + #middle 15 coefficients ^ + hil=np.array([ -0.0909457 , 0. , -0.127324 , 0. , -0.2122066 , 0. , -0.6366198 , 0. , 0.6366198 , 0. , 0.2122066 , + 0. , 0.127324 , 0. , 0.0909457 ]) + + if coeff_gain!=1: + hil=np.round(hil*coeff_gain)/coeff_gain + + imag_an=np.convolve(coh_sum,hil,mode='full')[len(hil)//2:len(coh_sum)+len(hil)//2] + + if adc_output: + imag_an=np.round(imag_an) + + if rnog_like: + envelope=np.max(np.array((coh_sum,imag_an)),axis=0)+3/8*min(np.array((coh_sum,imag_an)),axis=0) + else: + envelope=np.sqrt(coh_sum**2+imag_an**2) + + if adc_output=='counts': + envelope=np.round(envelope) + + return envelope def phase_signals(self, traces, beam_rolls): """ @@ -253,8 +293,8 @@ def phase_signals(self, traces, beam_rolls): running_i = 0 for subbeam_rolls in beam_rolls: - phased_trace = np.zeros(len(list(traces.values())[0])) + for channel_id in traces: trace = traces[channel_id] @@ -265,22 +305,25 @@ def phase_signals(self, traces, beam_rolls): return phased_traces - def phased_trigger( - self, station, det, - Vrms=None, - threshold=60 * units.mV, - triggered_channels=None, - phasing_angles=default_angles, - ref_index=1.75, - trigger_adc=False, # by default, assumes the trigger ADC is the same as the channels ADC - clock_offset=0, - adc_output='voltage', - trigger_filter=None, - upsampling_factor=1, - window=32, - step=16, - apply_digitization=True, - ): + def phased_trigger(self, station, det, + Vrms=None, + threshold=60 * units.mV, + triggered_channels=None, + phasing_angles=default_angles, + ref_index=1.75, + trigger_adc=False, # by default, assumes the trigger ADC is the same as the channels ADC + clock_offset=0, + adc_output='voltage', + trigger_filter=None, + upsampling_factor=1, + window=32, + step=16, + apply_digitization=False, + upsampling_method='fft', + coeff_gain=128, + rnog_like=False, + trig_type='power_integration' + ): """ simulates phased array trigger for each event @@ -315,9 +358,8 @@ def phased_trigger( clock_offset: float (default 0) Overall clock offset, for adc clock jitter reasons (see `apply_digitization`) adc_output: string (default 'voltage') - - 'voltage' to store the ADC output as discretised voltage trace - - 'counts' to store the ADC output in ADC counts + - 'counts' to store the ADC output in ADC counts and apply integer based math trigger_filter: array floats (default None) Freq. domain of the response to be applied to post-ADC traces Must be length for "MC freq" @@ -335,6 +377,19 @@ def phased_trigger( apply_digitization: bool (default True) Perform the quantization of the ADC. If set to true, should also set options `trigger_adc`, `adc_output`, `clock_offset` + upsampling_method: str (default 'fft') + Choose between FFT, FIR, or Linear Interpolaion based upsampling methods + coeff_gain: int (default 1) + If using the FIR upsampling, this will convert the floating point output of the + scipy filter to a fixed point value by multiplying by this factor and rounding to an + int. + rnog_like: bool (default False) + If true, this will apply the RNO-G FLOWER based math/rounding done in firmware. + trig_type: str (default "power_integration") + - "power_integration" do the power integration for the given window size and + step length + - "envelope" perform a hilbrt envelope threshold trigger on the beamformed + traces Returns ------- @@ -349,10 +404,14 @@ def phased_trigger( all time bins that fulfil the trigger condition per beam. The key is the beam number. Time with respect to first interaction time. maximum_amps: list of floats (length equal to that of `phasing_angles`) the maximum value of all the integration windows for each of the phased waveforms + n_trigs: int + total number of triggers that happened for all beams across the full traces + triggered_beams: list + list of bools for which beams triggered """ if(triggered_channels is None): - triggered_channels = [channel.get_id() for channel in station.iter_channels()] + triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] if(adc_output != 'voltage' and adc_output != 'counts'): error_msg = 'ADC output type must be "counts" or "voltage". Currently set to:' + str(adc_output) @@ -366,7 +425,7 @@ def phased_trigger( logger.debug(f"trigger channels: {triggered_channels}") traces = {} - for channel in station.iter_channels(use_channels=triggered_channels): + for channel in station.iter_trigger_channels(use_channels=triggered_channels): channel_id = channel.get_id() trace = np.array(channel.get_trace()) @@ -383,7 +442,6 @@ def phased_trigger( else: adc_sampling_frequency = channel.get_sampling_rate() - # Upsampling here, linear interpolate to mimic an FPGA internal upsampling if not isinstance(upsampling_factor, int): try: upsampling_factor = int(upsampling_factor) @@ -391,9 +449,41 @@ def phased_trigger( raise ValueError("Could not convert upsampling_factor to integer. Exiting.") if(upsampling_factor >= 2): - # FFT upsampling - new_len = len(trace) * upsampling_factor - upsampled_trace = scipy.signal.resample(trace, new_len) + + new_len = len(trace) * upsampling_factor + cur_t=np.arange(0,1/adc_sampling_frequency*len(trace),1/adc_sampling_frequency) + new_t=np.arange(0,1/adc_sampling_frequency*len(trace),1/adc_sampling_frequency/upsampling_factor) + + if upsampling_method=='fft': + upsampled_trace = scipy.signal.resample(trace, new_len) + + elif upsampling_method=='lin': + upsampled_trace = np.interp(new_t,cur_t,trace) + + elif upsampling_method=='fir': + if rnog_like: + up_filt=np.array([ 0.0, 0.0 , 0.0 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , -0.0 , -0.015625 , + -0.0390625 , -0.03125 , 0.0 , 0.0703125 , 0.15625 , 0.2265625 , 0.25 , 0.2265625 , 0.15625 , + 0.0703125 , 0.0 , -0.03125 , -0.0390625 , -0.015625 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , + 0.0 , 0.0 , 0.0 , 0.0 ]) + else: + cutoff=.5 + base_freq_filter_length=8 + filter_length=base_freq_filter_length*upsampling_factor-1 + up_filt=scipy.signal.firwin(filter_length,adc_sampling_frequency*cutoff,pass_zero='lowpass', + fs=adc_sampling_frequency*upsampling_factor) + if coeff_gain!=1: + up_filt=np.round(up_filt*coeff_gain)/coeff_gain + + zero_pad=np.zeros(len(trace)*upsampling_factor) + zero_pad[::upsampling_factor]=trace[:] + upsampled_trace=np.convolve(zero_pad,up_filt,mode='full')[len(up_filt)//2:len(zero_pad)+len(up_filt)//2]*upsampling_factor + + else: + error_msg = 'Interpolation method must be lin, fft, fir, ...' + raise ValueError(error_msg) + + if adc_output=='counts' and rnog_like==True: upsampled_trace=np.trunc(upsampled_trace) # If upsampled is performed, the final sampling frequency changes trace = upsampled_trace[:] @@ -412,7 +502,6 @@ def phased_trigger( phasing_angles, ref_index=ref_index, sampling_frequency=adc_sampling_frequency) - phased_traces = self.phase_signals(traces, beam_rolls) trigger_time = None @@ -420,28 +509,51 @@ def phased_trigger( channel_trace_start_time = self.get_channel_trace_start_time(station, triggered_channels) trigger_delays = {} - maximum_amps = np.zeros_like(phased_traces) + maximum_amps = np.zeros(len(phased_traces)) + n_trigs=0 + triggered_beams=[] for iTrace, phased_trace in enumerate(phased_traces): + is_triggered=False + + if trig_type=='power_integration': + squared_mean, num_frames = self.power_sum(coh_sum=phased_trace, window=window, step=step, adc_output=adc_output, rnog_like=rnog_like) + maximum_amps[iTrace] = np.max(squared_mean) + + if True in (squared_mean > threshold): + n_trigs+=len(np.where((squared_mean > threshold)==True)[0]) + trigger_delays[iTrace] = {} + + for channel_id in beam_rolls[iTrace]: + trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step + + triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(squared_mean > threshold))) + logger.debug(f"Station has triggered, at bins {triggered_bins}") + logger.debug(trigger_delays) + logger.debug(f"trigger_delays {trigger_delays[iTrace][triggered_channels[0]]}") + is_triggered = True + trigger_times[iTrace] = trigger_delays[iTrace][triggered_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time + logger.debug(f"trigger times = {trigger_times[iTrace]}") + + elif trig_type=='envelope': + hilbert_env = self.hilbert_envelope(coh_sum=phased_trace,adc_output=adc_output,rnog_like=rnog_like) + maximum_amps[iTrace] = np.max(hilbert_env) + + if True in (hilbert_env>threshold): + n_trigs+=len(np.where((hilbert_env > threshold)==True)[0]) + trigger_delays[iTrace] = {} + for channel_id in beam_rolls[iTrace]: + trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step + triggered_bins=np.atleast_1d(np.squeeze(np.argwhere(hilbert_env > threshold))) + is_triggered=True + trigger_times[iTrace] = trigger_delays[iTrace][triggered_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time - # Create a sliding window - squared_mean, num_frames = self.power_sum(coh_sum=phased_trace, window=window, step=step, adc_output=adc_output) - - maximum_amps[iTrace] = np.max(squared_mean) - - if True in (squared_mean > threshold): - trigger_delays[iTrace] = {} + else: + raise NotImplementedError("not a good trigger type: options (power_integration, envelope)") - for channel_id in beam_rolls[iTrace]: - trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step + triggered_beams.append(is_triggered) - triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(squared_mean > threshold))) - logger.debug(f"Station has triggered, at bins {triggered_bins}") - logger.debug(trigger_delays) - logger.debug(f"trigger_delays {trigger_delays[iTrace][triggered_channels[0]]}") - is_triggered = True - trigger_times[iTrace] = trigger_delays[iTrace][triggered_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time - logger.debug(f"trigger times = {trigger_times[iTrace]}") + is_triggered=np.any(triggered_beams) if is_triggered: logger.debug("Trigger condition satisfied!") @@ -449,7 +561,7 @@ def phased_trigger( trigger_time = min([x.min() for x in trigger_times.values()]) logger.debug(f"minimum trigger time is {trigger_time:.0f}ns") - return is_triggered, trigger_delays, trigger_time, trigger_times, maximum_amps + return is_triggered, trigger_delays, trigger_time, trigger_times, maximum_amps, n_trigs, triggered_beams @register_run() def run(self, evt, station, det, @@ -468,10 +580,15 @@ def run(self, evt, station, det, window=32, step=16, apply_digitization=True, + upsampling_method='fft', + coeff_gain=128, + rnog_like=False, + trig_type='power_integration', + return_n_triggers=False ): """ - Simulates phased array trigger for each event. + simulates phased array trigger for each event Several channels are phased by delaying their signals by an amount given by a pointing angle. Several pointing angles are possible in order to cover @@ -500,7 +617,7 @@ def run(self, evt, station, det, phasing_angles: array of float pointing angles for the primary beam set_not_triggered: bool (default False) - If True, no trigger simulation will be performed and this trigger will be set to not_triggered + if True not trigger simulation will be performed and this trigger will be set to not_triggered ref_index: float (default 1.75) refractive index for beam forming trigger_adc: bool, (default True) @@ -530,73 +647,109 @@ def run(self, evt, station, det, apply_digitization: bool (default True) Perform the quantization of the ADC. If set to true, should also set options `trigger_adc`, `adc_output`, `clock_offset` + upsampling_method: str (default 'fft') + Choose between FFT, FIR, or Linear Interpolaion based upsampling methods + coeff_gain: int (default 1) + If using the FIR upsampling, this will convert the floating point output of the + scipy filter to a fixed point value by multiplying by this factor and rounding to an + int. + rnog_like: bool (default False) + If true, this will apply the RNO-G FLOWER based math/rounding done in firmware. + trig_type: str (default "power_integration") + - "power_integration" do the power integration for the given window size and + step length + - "envelope" perform a hilbrt envelope threshold trigger on the beamformed + traces + return_n_triggers: bool (default False) + To better estimate simulated thresholds one should count the total triggers + in the entire trace for each beam. If true, this return the total trigger number. Returns ------- is_triggered: bool True if the triggering condition is met + n_triggers: int (Optional) + Count of the total number of triggers in all beamformed traces """ - if triggered_channels is None: - triggered_channels = [channel.get_id() for channel in station.iter_channels()] + if(triggered_channels is None): + triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] - if adc_output != 'voltage' and adc_output != 'counts': + if(adc_output != 'voltage' and adc_output != 'counts'): error_msg = 'ADC output type must be "counts" or "voltage". Currently set to:' + str(adc_output) raise ValueError(error_msg) is_triggered = False trigger_delays = {} - if set_not_triggered: + if(set_not_triggered): is_triggered = False trigger_delays = {} - maximum_amps = np.zeros_like(phasing_angles) - + triggered_beams = [] else: - is_triggered, trigger_delays, trigger_time, trigger_times, maximum_amps = self.phased_trigger( - station=station, - det=det, - Vrms=Vrms, - threshold=threshold, - triggered_channels=triggered_channels, - phasing_angles=phasing_angles, - ref_index=ref_index, - trigger_adc=trigger_adc, - clock_offset=clock_offset, - adc_output=adc_output, - trigger_filter=trigger_filter, - upsampling_factor=upsampling_factor, - window=window, - step=step, - apply_digitization=apply_digitization, - ) + is_triggered, trigger_delays, trigger_time, trigger_times,\ + maximum_amps, n_triggers, triggered_beams = self.phased_trigger(station=station, + det=det, + Vrms=Vrms, + threshold=threshold, + triggered_channels=triggered_channels, + phasing_angles=phasing_angles, + ref_index=ref_index, + trigger_adc=trigger_adc, + clock_offset=clock_offset, + adc_output=adc_output, + trigger_filter=trigger_filter, + upsampling_factor=upsampling_factor, + window=window, + step=step, + apply_digitization=apply_digitization, + upsampling_method=upsampling_method, + coeff_gain=coeff_gain, + rnog_like=rnog_like, + trig_type=trig_type + ) # Create a trigger object to be returned to the station - trigger = SimplePhasedTrigger( - trigger_name, - threshold, - channels=triggered_channels, - primary_angles=phasing_angles, - trigger_delays=trigger_delays, - window_size=window, - step_size=step, - maximum_amps=maximum_amps, - ) + if trig_type=='power_integration': + trigger = SimplePhasedTrigger( + trigger_name, + threshold, + channels=triggered_channels, + primary_angles=phasing_angles, + trigger_delays=trigger_delays, + window_size=window, + step_size=step, + maximum_amps=maximum_amps + ) + + elif trig_type=='envelope': + trigger = SimplePhasedTrigger( + trigger_name, + threshold, + channels=triggered_channels, + primary_angles=phasing_angles, + trigger_delays=trigger_delays, + maximum_amps=maximum_amps + ) + + else: + raise NotImplementedError("invalid phased trigger type") trigger.set_triggered(is_triggered) if is_triggered: - # trigger_time(s)= time(s) from start of trace + start time of trace with respect to moment of first - # interaction = trigger time from moment of first interaction; time offset to interaction time - # (channel_trace_start_time) already recognized in self.phased_trigger - trigger.set_trigger_time(trigger_time) + #trigger_time(s)= time(s) from start of trace + start time of trace with respect to moment of first interaction = trigger time from moment of first interaction; time offset to interaction time (channel_trace_start_time) already recognized in self.phased_trigger + trigger.set_trigger_time(trigger_time)# trigger.set_trigger_times(trigger_times) else: trigger.set_trigger_time(None) station.set_trigger(trigger) - return is_triggered + if return_n_triggers: + return is_triggered,n_triggers + else: + return is_triggered def end(self): pass diff --git a/NuRadioReco/utilities/signal_processing.py b/NuRadioReco/utilities/signal_processing.py index 77c56e73be..1fb091ca53 100644 --- a/NuRadioReco/utilities/signal_processing.py +++ b/NuRadioReco/utilities/signal_processing.py @@ -54,20 +54,20 @@ def add_cable_delay(station, det, sim_to_data, trigger=False, logger=None): If set, use `logger.debug(..)` to log the cable delay. """ - if trigger and not isinstance(det, detector.rnog_detector.Detector): - raise ValueError("Simulating extra trigger channels is only possible with the `rnog_detector.Detector` class.") - if sim_to_data: new_trace_start_times = [] for channel in station.iter_channels(): if trigger: if not channel.has_extra_trigger_channel(): continue - cable_delay = det.get_cable_delay(station.get_id(), channel.get_id(), trigger=True) + + if isinstance(det,detector.detector_base.DetectorBase): + cable_delay = det.get_cable_delay(station.get_id(), channel.get_id()) + else: + cable_delay = det.get_cable_delay(station.get_id(), channel.get_id(), trigger=True) new_trace_start_times.append(channel.get_trigger_channel().get_trace_start_time() + cable_delay) else: - # Only the RNOG detector has the argument `trigger`. Default is false cable_delay = det.get_cable_delay(station.get_id(), channel.get_id()) new_trace_start_times.append(channel.get_trace_start_time() + cable_delay) @@ -88,11 +88,10 @@ def add_cable_delay(station, det, sim_to_data, trigger=False, logger=None): if delta_time: # tmp code block to try to get the same results as before - roll_samples = int(delta_time / channel.get_sampling_rate()) + roll_samples = int(delta_time * channel.get_sampling_rate()) delta_time = delta_time - roll_samples * channel.get_sampling_rate() trace = np.roll(channel.get_trace(), roll_samples) channel.set_trace(trace, channel.get_sampling_rate()) - channel.apply_time_shift(delta_time) channel.set_trace_start_time(new_trace_start_time) From a7dcc799e268ca3de68343247f827ac34a6c27c0 Mon Sep 17 00:00:00 2001 From: Ryan Krebs Date: Sat, 21 Dec 2024 02:47:44 -0500 Subject: [PATCH 2/7] track gain to help efficiency sims and fix formatting --- .../modules/RNO_G/triggerBoardResponse.py | 149 ++++++++---------- .../modules/phasedarray/triggerSimulator.py | 79 +++++----- 2 files changed, 109 insertions(+), 119 deletions(-) diff --git a/NuRadioReco/modules/RNO_G/triggerBoardResponse.py b/NuRadioReco/modules/RNO_G/triggerBoardResponse.py index 62c07c3f29..6aa81b9900 100644 --- a/NuRadioReco/modules/RNO_G/triggerBoardResponse.py +++ b/NuRadioReco/modules/RNO_G/triggerBoardResponse.py @@ -51,30 +51,6 @@ def begin(self, adc_input_range=2 * units.volt, clock_offset=0.0, adc_output="vo self._triggerBoardAmplifications = np.array([1, 1.25, 2, 2.5, 4, 5, 8, 10, 12.5, 16, 20, 25, 32, 50]) self._adc_input_range = adc_input_range - def apply_trigger_filter(self, station, trigger_channels, trigger_filter): - """ - Applies the requested trigger filter to the `trigger_channels` - - Parameters - ---------- - station : Station - Station to use - trigger_channels : list - Channels that this function should be applied to - trigger_filter : function - set of interpolations describing the `gain` and `phase` of the filter - (see function `load_amp_response` in file `./detector/RNO_G/analog_components.py`) - - """ - - for channel_id in trigger_channels: - channel = station.get_trigger_channel(channel_id) - - # calculate and apply trigger filters - freqs = channel.get_frequencies() - filt = trigger_filter(freqs) - channel.set_frequency_spectrum(channel.get_frequency_spectrum() * filt, channel.get_sampling_rate()) - def get_vrms(self, station, trigger_channels, trace_split=20): """ Estimates the RMS voltage of the triggering antennas by splitting the waveforms @@ -109,7 +85,7 @@ def get_vrms(self, station, trigger_channels, trace_split=20): self.logger.debug(vrms) return vrms - def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): + def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None, gain_values=None): """ Calculates and applies the gain adjustment such that the correct number of "noise bits" are realized. The ADC has fixed possible gain values and @@ -126,15 +102,19 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): avg_rms : float (default: None) The Vrms of the trigger channels including the trigger board filters If set to `None`, this will be estimated using the waveforms - + gain_values : list (default: None) + If set these will be applied to the channel. Otherwise gains are recalculated. Returns ------- - vrms_after_gain : float + vrms_after_gain : list the RMS voltage of the waveforms after the gain has been applied ideal_vrms: float the ideal vrms, as measured on the ADC capacitors + ret_gain_values: list + gain values applied to each channel + """ if avg_vrms is None: @@ -146,6 +126,8 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): avg_vrms = np.full_like(trigger_channels, avg_vrms, dtype=float) vrms_after_gain = [] + ret_gain_values = [] + for channel_id, vrms in zip(trigger_channels, avg_vrms): det_channel = det.get_channel(station.get_id(), channel_id) @@ -154,30 +136,39 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None): volts_per_adc = self._adc_input_range / (2 ** total_bits) ideal_vrms = volts_per_adc * (2 ** (noise_bits) - 1) - msg = f"\t Ch: {channel_id}\t Target Vrms: {ideal_vrms / units.mV:0.3f} mV" - msg += f"\t V/ADC: {volts_per_adc / units.mV:0.3f} mV" - self.logger.debug(msg) - - # find the ADC gain from the possible values that makes the realized - # vrms closest-to-but-greater-than the ideal value - amplified_vrms_values = vrms * self._triggerBoardAmplifications - mask = amplified_vrms_values > ideal_vrms - if np.any(mask): - gain_to_use = self._triggerBoardAmplifications[mask][0] - vrms_after_gain.append(amplified_vrms_values[mask][0]) + if gain_values is not None: + vrms_after_gain.append(vrms * gain_values[channel_id]) + channel = station.get_trigger_channel(channel_id) + channel.set_trace(channel.get_trace() * gain_values[channel_id], channel.get_sampling_rate()) + else: - gain_to_use = self._triggerBoardAmplifications[-1] - vrms_after_gain.append(amplified_vrms_values[-1]) - channel = station.get_trigger_channel(channel_id) - self.logger.debug(f"\t Ch: {channel_id}\t Actuall Vrms: {np.std(channel.get_trace() * gain_to_use) / units.mV:0.3f} mV") - channel.set_trace(channel.get_trace() * gain_to_use, channel.get_sampling_rate()) - self.logger.debug(f"\t Used Vrms: {vrms_after_gain[-1] / units.mV:0.3f} mV" + f"\tADC Gain {gain_to_use}") - eff_noise_bits = np.log2(vrms_after_gain[-1] / volts_per_adc) + 1 - self.logger.debug(f"\t Eff noise bits: {eff_noise_bits:0.2f}\tRequested: {noise_bits}") + msg = f"\t Ch: {channel_id}\t Target Vrms: {ideal_vrms / units.mV:0.3f} mV" + msg += f"\t V/ADC: {volts_per_adc / units.mV:0.3f} mV" + self.logger.debug(msg) + + # find the ADC gain from the possible values that makes the realized + # vrms closest-to-but-greater-than the ideal value + amplified_vrms_values = vrms * self._triggerBoardAmplifications + mask = amplified_vrms_values > ideal_vrms - return np.array(vrms_after_gain), ideal_vrms + if np.any(mask): + gain_to_use = self._triggerBoardAmplifications[mask][0] + vrms_after_gain.append(amplified_vrms_values[mask][0]) + else: + gain_to_use = self._triggerBoardAmplifications[-1] + vrms_after_gain.append(amplified_vrms_values[-1]) + ret_gain_values.append(gain_to_use) + channel = station.get_trigger_channel(channel_id) + channel.set_trace(channel.get_trace() * gain_to_use, channel.get_sampling_rate()) + eff_noise_bits = np.log2(vrms_after_gain[-1] / volts_per_adc) + 1 + + self.logger.debug(f"\t Ch: {channel_id}\t Actuall Vrms: {np.std(channel.get_trace() * gain_to_use) / units.mV:0.3f} mV") + self.logger.debug(f"\t Used Vrms: {vrms_after_gain[-1] / units.mV:0.3f} mV" + f"\tADC Gain {gain_to_use}") + self.logger.debug(f"\t Eff noise bits: {eff_noise_bits:0.2f}\tRequested: {noise_bits}") + + return np.array(vrms_after_gain), ideal_vrms, ret_gain_values def digitize_trace(self, station, det, trigger_channels, vrms): for channel_id in trigger_channels: @@ -201,7 +192,7 @@ def digitize_trace(self, station, det, trigger_channels, vrms): @register_run() def run(self, evt, station, det, trigger_channels, vrms=None, apply_adc_gain=True, - digitize_trace=True): + digitize_trace=True, gain_values=None): """ Applies the additional filters on the trigger board and performs a gain amplification to get the correct number of trigger bits. @@ -228,6 +219,8 @@ def run(self, evt, station, det, trigger_channels, vrms=None, apply_adc_gain=Tru ------- trigger_board_vrms : float the RMS voltage of the waveforms on the trigger board after applying the ADC gain + ret_gain_values : list + the gain values applied to each channel """ self.logger.debug("Applying the RNO-G trigger board response") @@ -235,16 +228,17 @@ def run(self, evt, station, det, trigger_channels, vrms=None, apply_adc_gain=Tru vrms = self.get_vrms(station, trigger_channels) if apply_adc_gain: - trigger_board_vrms, ideal_vrms = self.apply_adc_gain(station, det, trigger_channels, vrms) + trigger_board_vrms, ideal_vrms, ret_gain_values = self.apply_adc_gain(station, det, trigger_channels, vrms, gain_values) else: trigger_board_vrms = vrms ideal_vrms = np.mean(vrms) + ret_gain_values = None if digitize_trace: self.digitize_trace(station, det, trigger_channels, ideal_vrms) - trigger_board_vrms=self.get_vrms(station, trigger_channels) + trigger_board_vrms = self.get_vrms(station, trigger_channels) - return trigger_board_vrms + return trigger_board_vrms, ret_gain_values def end(self): from datetime import timedelta @@ -269,46 +263,39 @@ def end(self): rnogADCResponse.begin(adc_input_range=2 * units.volt, clock_offset=0.0, adc_output="counts") channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() - det_file='RNO_G/RNO_single_station_only_PA.json' + det_file = 'RNO_G/RNO_single_station_only_PA.json' det = detector.Detector(source='json',json_filename=det_file) - station_id=11 + station_id = 11 channel_ids = np.arange(4) - n_samples=1024 - sampling_rate=472*units.MHz + n_samples = 1024 + sampling_rate = 472 * units.MHz dt = 1 / sampling_rate ff = np.fft.rfftfreq(n_samples, dt) max_freq = ff[-1] min_freq = 0 fff = np.linspace(min_freq, max_freq, 10000) - four_filters={} - rf_filter=rnogHarwareResponse.get_filter(ff,station_id,channel_ids[0],det,sim_to_data=True,is_trigger=True) - chain_filter=rf_filter - - for i in range(4): - four_filters[i]=chain_filter - - four_filters_highres={} - rf_filter_highres=rnogHarwareResponse.get_filter(fff,station_id,channel_ids[0],det,sim_to_data=True,is_trigger=True) - chain_filter_highres=rf_filter_highres + four_filters_highres = {} + rf_filter_highres = rnogHarwareResponse.get_filter(fff, station_id, channel_ids[0], det,sim_to_data=True, is_trigger=True) + chain_filter_highres = rf_filter_highres for i in channel_ids: - four_filters_highres[i]=chain_filter_highres - + four_filters_highres[i] = chain_filter_highres Vrms = 1 - noise_temp=300 - bandwidth={} + noise_temp = 300 + bandwidth = {} Vrms_ratio = {} amplitude = {} - per_channel_vrms=[] + per_channel_vrms = [] + for i in channel_ids: integrated_channel_response = np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) - rel_channel_response=np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) - bandwidth[i]=integrated_channel_response + rel_channel_response = np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) + bandwidth[i] = integrated_channel_response Vrms_ratio[i] = np.sqrt(rel_channel_response / (max_freq - min_freq)) - chan_vrms=(noise_temp * 50 * constants.k * integrated_channel_response / units.Hz) ** 0.5 + chan_vrms = (noise_temp * 50 * constants.k * integrated_channel_response / units.Hz) ** 0.5 per_channel_vrms.append(chan_vrms) amplitude[i] = chan_vrms / Vrms_ratio[i] @@ -323,24 +310,24 @@ def end(self): spectrum = channelGenericNoiseAdder.bandlimited_noise(min_freq, max_freq, n_samples, sampling_rate, amplitude[channel_id], type="rayleigh", time_domain=False) - trace=fft.freq2time(spectrum*four_filters[channel_id], sampling_rate) + trace = fft.freq2time(spectrum * four_filters_highres[channel_id], sampling_rate) channel = NuRadioReco.framework.channel.Channel(channel_id) channel.set_trace(trace, sampling_rate) station.add_trigger_channel(channel) - fig,ax=plt.subplots(1,2,sharex=True,figsize=(11,7)) + fig , ax = plt.subplots(1, 2, sharex=True, figsize=(11,7)) for channel_id in channel_ids: - ch=station.get_channel(channel_id) - ax[0].plot(ch.get_times(),ch.get_trace(),label='ch %i'%channel_id) + ch = station.get_channel(channel_id) + ax[0].plot(ch.get_times(), ch.get_trace(), label='ch %i' % channel_id) - chan_rms=rnogADCResponse.run(evt,station,det,requested_channels=channel_ids, - digitize_trace=True,apply_adc_gain=True) + chan_rms, gain_values = rnogADCResponse.run(evt, station, det, requested_channels=channel_ids, + digitize_trace=True, apply_adc_gain=True) for channel_id in channel_ids: - ch=station.get_channel(channel_id) - ax[1].plot(ch.get_times(),ch.get_trace(),label='ch %i'%channel_id) + ch = station.get_channel(channel_id) + ax[1].plot(ch.get_times(), ch.get_trace(), label='ch %i' % channel_id) ax[0].set_title('Raw Voltage Trace') ax[0].set_ylabel('Voltage [V]') diff --git a/NuRadioReco/modules/phasedarray/triggerSimulator.py b/NuRadioReco/modules/phasedarray/triggerSimulator.py index fc5dc4fcdc..a43549e54e 100644 --- a/NuRadioReco/modules/phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/phasedarray/triggerSimulator.py @@ -121,9 +121,7 @@ def calculate_time_delays(self, station, det, delays += [-(ant_z[key] - ref_z) / cspeed * ref_index * np.sin(angle) - cable_delays[key]] delays -= np.max(delays) - roll = np.array(np.round(np.array(delays) / time_step)).astype(int) - subbeam_rolls = dict(zip(triggered_channels, roll)) # logger.debug("angle:", angle / units.deg) @@ -154,8 +152,10 @@ def get_channel_trace_start_time(self, station, triggered_channels): channel_trace_start_time = None for channel in station.iter_trigger_channels(use_channels=triggered_channels): + if channel_trace_start_time is None: channel_trace_start_time = channel.get_trace_start_time() + elif channel_trace_start_time != channel.get_trace_start_time(): error_msg = 'Phased array channels do not have matching trace start times. ' error_msg += 'This module is not prepared for this case.' @@ -184,6 +184,7 @@ def check_vertical_string(self, station, det, triggered_channels): diff_x = np.abs(ant_x - ant_x[0]) ant_y = np.fromiter(self.get_antenna_positions(station, det, triggered_channels, 1).values(), dtype=float) diff_y = np.abs(ant_y - ant_y[0]) + if (sum(diff_x) > cut or sum(diff_y) > cut): raise NotImplementedError('The phased triggering array should lie on a vertical line') @@ -226,8 +227,8 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage',rnog_like=False) coh_sum_squared = (coh_sum * coh_sum).astype(float) else: #(adc_output == 'counts'): if rnog_like: - coh_sum[coh_sum>2**6-1]=2**6-1 - coh_sum[coh_sum<-2**6]=-2**6 + coh_sum[coh_sum>2**6-1] = 2**6 - 1 + coh_sum[coh_sum<-2**6] = -2**6 coh_sum_squared = (coh_sum * coh_sum).astype(int) coh_sum_windowed = np.lib.stride_tricks.as_strided(coh_sum_squared, (num_frames, window), @@ -236,14 +237,14 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage',rnog_like=False) power = np.sum(coh_sum_windowed, axis=1) return_power=power.astype(float) / window - if adc_output=='counts': return_power=np.floor(return_power) + if adc_output=='counts': return_power = np.floor(return_power) return return_power, num_frames - def hilbert_envelope(self,coh_sum,adc_output='voltage',coeff_gain=1,rnog_like=False): + def hilbert_envelope(self, coh_sum, adc_output='voltage', coeff_gain=1, rnog_like=False): if rnog_like: - coeff_gain=128 + coeff_gain = 128 #31 sample fir transformer #hil=[ -0.0424413 , 0. , -0.0489708 , 0. , -0.0578745 , 0. , -0.0707355 , 0. , -0.0909457 , 0. , -0.127324 , 0. @@ -251,24 +252,24 @@ def hilbert_envelope(self,coh_sum,adc_output='voltage',coeff_gain=1,rnog_like=Fa # , 0. , 0.0578745 , 0. , 0.0489708 , 0. , 0.0424413 ] #middle 15 coefficients ^ - hil=np.array([ -0.0909457 , 0. , -0.127324 , 0. , -0.2122066 , 0. , -0.6366198 , 0. , 0.6366198 , 0. , 0.2122066 , - 0. , 0.127324 , 0. , 0.0909457 ]) + hil = np.array([ -0.0909457 , 0. , -0.127324 , 0. , -0.2122066 , 0. , -0.6366198 , 0. , 0.6366198 , 0. , 0.2122066 , + 0. , 0.127324 , 0. , 0.0909457 ]) if coeff_gain!=1: - hil=np.round(hil*coeff_gain)/coeff_gain + hil = np.round(hil * coeff_gain) / coeff_gain - imag_an=np.convolve(coh_sum,hil,mode='full')[len(hil)//2:len(coh_sum)+len(hil)//2] + imag_an = np.convolve(coh_sum, hil, mode='full')[len(hil)//2 : len(coh_sum) + len(hil)//2] if adc_output: - imag_an=np.round(imag_an) + imag_an = np.round(imag_an) if rnog_like: - envelope=np.max(np.array((coh_sum,imag_an)),axis=0)+3/8*min(np.array((coh_sum,imag_an)),axis=0) + envelope = np.max(np.array((coh_sum,imag_an)), axis=0) + (3 / 8) * min(np.array((coh_sum,imag_an)), axis=0) else: - envelope=np.sqrt(coh_sum**2+imag_an**2) + envelope = np.sqrt(coh_sum**2 + imag_an**2) if adc_output=='counts': - envelope=np.round(envelope) + envelope = np.round(envelope) return envelope @@ -451,8 +452,8 @@ def phased_trigger(self, station, det, if(upsampling_factor >= 2): new_len = len(trace) * upsampling_factor - cur_t=np.arange(0,1/adc_sampling_frequency*len(trace),1/adc_sampling_frequency) - new_t=np.arange(0,1/adc_sampling_frequency*len(trace),1/adc_sampling_frequency/upsampling_factor) + cur_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/adc_sampling_frequency) + new_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/adc_sampling_frequency/upsampling_factor) if upsampling_method=='fft': upsampled_trace = scipy.signal.resample(trace, new_len) @@ -462,26 +463,27 @@ def phased_trigger(self, station, det, elif upsampling_method=='fir': if rnog_like: - up_filt=np.array([ 0.0, 0.0 , 0.0 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , -0.0 , -0.015625 , - -0.0390625 , -0.03125 , 0.0 , 0.0703125 , 0.15625 , 0.2265625 , 0.25 , 0.2265625 , 0.15625 , - 0.0703125 , 0.0 , -0.03125 , -0.0390625 , -0.015625 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , - 0.0 , 0.0 , 0.0 , 0.0 ]) + #I think I can do better than this in firmware + up_filt = np.array([ 0.0, 0.0 , 0.0 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , -0.0 , -0.015625 , + -0.0390625 , -0.03125 , 0.0 , 0.0703125 , 0.15625 , 0.2265625 , 0.25 , 0.2265625 , 0.15625 , + 0.0703125 , 0.0 , -0.03125 , -0.0390625 , -0.015625 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , + 0.0 , 0.0 , 0.0 , 0.0 ]) else: - cutoff=.5 - base_freq_filter_length=8 - filter_length=base_freq_filter_length*upsampling_factor-1 - up_filt=scipy.signal.firwin(filter_length,adc_sampling_frequency*cutoff,pass_zero='lowpass', + cutoff = .5 + base_freq_filter_length = 8 + filter_length = base_freq_filter_length * upsampling_factor - 1 + up_filt = scipy.signal.firwin(filter_length, adc_sampling_frequency*cutoff, pass_zero='lowpass', fs=adc_sampling_frequency*upsampling_factor) if coeff_gain!=1: - up_filt=np.round(up_filt*coeff_gain)/coeff_gain + up_filt = np.round(up_filt*coeff_gain)/coeff_gain - zero_pad=np.zeros(len(trace)*upsampling_factor) - zero_pad[::upsampling_factor]=trace[:] - upsampled_trace=np.convolve(zero_pad,up_filt,mode='full')[len(up_filt)//2:len(zero_pad)+len(up_filt)//2]*upsampling_factor + zero_pad = np.zeros(len(trace)*upsampling_factor) + zero_pad[::upsampling_factor] = trace[:] + upsampled_trace = np.convolve(zero_pad, up_filt, mode='full')[len(up_filt)//2 : len(zero_pad) + len(up_filt)//2] * upsampling_factor else: - error_msg = 'Interpolation method must be lin, fft, fir, ...' - raise ValueError(error_msg) + error_msg = 'Interpolation method must be lin, fft, or fir' + raise NotImplementedError(error_msg) if adc_output=='counts' and rnog_like==True: upsampled_trace=np.trunc(upsampled_trace) @@ -502,6 +504,7 @@ def phased_trigger(self, station, det, phasing_angles, ref_index=ref_index, sampling_frequency=adc_sampling_frequency) + phased_traces = self.phase_signals(traces, beam_rolls) trigger_time = None @@ -510,8 +513,8 @@ def phased_trigger(self, station, det, trigger_delays = {} maximum_amps = np.zeros(len(phased_traces)) - n_trigs=0 - triggered_beams=[] + n_trigs = 0 + triggered_beams = [] for iTrace, phased_trace in enumerate(phased_traces): is_triggered=False @@ -521,7 +524,7 @@ def phased_trigger(self, station, det, maximum_amps[iTrace] = np.max(squared_mean) if True in (squared_mean > threshold): - n_trigs+=len(np.where((squared_mean > threshold)==True)[0]) + n_trigs += len(np.where((squared_mean > threshold)==True)[0]) trigger_delays[iTrace] = {} for channel_id in beam_rolls[iTrace]: @@ -536,16 +539,16 @@ def phased_trigger(self, station, det, logger.debug(f"trigger times = {trigger_times[iTrace]}") elif trig_type=='envelope': - hilbert_env = self.hilbert_envelope(coh_sum=phased_trace,adc_output=adc_output,rnog_like=rnog_like) + hilbert_env = self.hilbert_envelope(coh_sum=phased_trace, adc_output=adc_output, rnog_like=rnog_like) maximum_amps[iTrace] = np.max(hilbert_env) if True in (hilbert_env>threshold): - n_trigs+=len(np.where((hilbert_env > threshold)==True)[0]) + n_trigs += len(np.where((hilbert_env > threshold)==True)[0]) trigger_delays[iTrace] = {} for channel_id in beam_rolls[iTrace]: trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step - triggered_bins=np.atleast_1d(np.squeeze(np.argwhere(hilbert_env > threshold))) - is_triggered=True + triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(hilbert_env > threshold))) + is_triggered = True trigger_times[iTrace] = trigger_delays[iTrace][triggered_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time else: From 52c786ffb14f0fd46bcd663049f34b710b29f9da Mon Sep 17 00:00:00 2001 From: Ryan Krebs Date: Tue, 7 Jan 2025 13:32:03 -0500 Subject: [PATCH 3/7] refactor phased array, add things for eff sims, add upsampling to trace utils --- NuRadioMC/simulation/output_writer_hdf5.py | 4 +- NuRadioMC/simulation/simulation.py | 4 +- NuRadioReco/framework/trigger.py | 61 ++- .../modules/phasedarray/phasedArray.py | 211 +++++++++ .../analog_beamformed_envelope_trigger.py} | 149 ++++--- .../beamformed_power_integration.py} | 416 +++--------------- .../digital_beamformed_envelope_trigger.py | 411 +++++++++++++++++ .../modules/trigger/highLowThreshold.py | 2 + NuRadioReco/utilities/trace_utilities.py | 104 ++++- 9 files changed, 915 insertions(+), 447 deletions(-) create mode 100644 NuRadioReco/modules/phasedarray/phasedArray.py rename NuRadioReco/modules/{envelope_phasedarray/triggerSimulator.py => trigger/analog_beamformed_envelope_trigger.py} (66%) rename NuRadioReco/modules/{phasedarray/triggerSimulator.py => trigger/beamformed_power_integration.py} (55%) create mode 100644 NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py diff --git a/NuRadioMC/simulation/output_writer_hdf5.py b/NuRadioMC/simulation/output_writer_hdf5.py index 4acca361fa..91bc1a2fbf 100644 --- a/NuRadioMC/simulation/output_writer_hdf5.py +++ b/NuRadioMC/simulation/output_writer_hdf5.py @@ -542,6 +542,6 @@ def calculate_Veff(self): V = self._mout_attributes['volume'] Veff = V * n_triggered_weighted / n_events logger.status(f"Veff = {Veff / units.km ** 3:.4g} km^3, Veffsr = {Veff * 4 * np.pi/units.km**3:.4g} km^3 sr") - return Veff + return Veff, n_triggered except: - return None + return 0, 0 diff --git a/NuRadioMC/simulation/simulation.py b/NuRadioMC/simulation/simulation.py index 99e674f864..84e84c9eb1 100644 --- a/NuRadioMC/simulation/simulation.py +++ b/NuRadioMC/simulation/simulation.py @@ -1653,11 +1653,13 @@ def run(self): eventWriter.end() logger.debug("closing nur file") - self._output_writer_hdf5.calculate_Veff() + Veff, n_triggered = self._output_writer_hdf5.calculate_Veff() if not self._output_writer_hdf5.write_output_file(): logger.warning("No events were triggered. Writing empty HDF5 output file.") self._output_writer_hdf5.write_empty_output_file(self._fin_attrs) + return n_triggered + 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) diff --git a/NuRadioReco/framework/trigger.py b/NuRadioReco/framework/trigger.py index fb5b801fa6..f870c3d966 100644 --- a/NuRadioReco/framework/trigger.py +++ b/NuRadioReco/framework/trigger.py @@ -17,14 +17,16 @@ def deserialize(triggers_pkl): trigger = SimpleThresholdTrigger(None, None) elif(trigger_type == 'high_low'): trigger = HighLowTrigger(None, None, None, None, None) - elif(trigger_type == 'simple_phased'): - trigger = SimplePhasedTrigger(None, None) + elif(trigger_type == 'power_integration_phased'): + trigger = PowerIntegrationPhasedTrigger(None, None) elif(trigger_type == 'envelope_trigger'): trigger = EnvelopeTrigger(None, None, None, None) elif trigger_type == 'int_power': trigger = IntegratedPowerTrigger(None, None, None) - elif trigger_type == 'envelope_phased': - trigger = EnvelopePhasedTrigger(None, None, None, None) + elif trigger_type == 'digital_envelope_phased': + trigger = AnalogEnvelopePhasedTrigger(None, None, None, None) + elif trigger_type == 'analog_envelope_phased': + trigger = DigitalEnvelopePhasedTrigger(None, None, None) elif(trigger_type == 'rnog_surface_trigger'): trigger = RNOGSurfaceTrigger(None, None, None, None) else: @@ -238,7 +240,7 @@ def __init__(self, name, threshold, channels=None, number_of_coincidences=1, self._coinc_window = channel_coincidence_window -class EnvelopePhasedTrigger(Trigger): +class AnalogEnvelopePhasedTrigger(Trigger): def __init__(self, name, threshold_factor, power_mean, power_std, triggered_channels=None, phasing_angles=None, trigger_delays=None, @@ -270,7 +272,7 @@ def __init__(self, name, threshold_factor, power_mean, power_std, Frequencies for a 6th-order Butterworth filter to be applied after the diode filtering. """ - Trigger.__init__(self, name, triggered_channels, 'envelope_phased', pre_trigger_times=pre_trigger_times) + Trigger.__init__(self, name, triggered_channels, 'analog_envelope_phased', pre_trigger_times=pre_trigger_times) self._triggered_channels = triggered_channels self._phasing_angles = phasing_angles self._threshold_factor = threshold_factor @@ -279,10 +281,49 @@ def __init__(self, name, threshold_factor, power_mean, power_std, self._trigger_delays = trigger_delays self._output_passband = output_passband +class DigitalEnvelopePhasedTrigger(Trigger): -class SimplePhasedTrigger(Trigger): + def __init__(self, name, threshold, trigger_channels=None, + phasing_angles=None, trigger_delays=None, + maximum_amps=None, pre_trigger_times=55 * units.ns): + """ + initialize trigger class + + Parameters + ---------- + name: string + unique name of the trigger + threshold_factor: float + the threshold factor + power_mean: float + mean of the noise trace after being filtered with the diode + power_std: float + standard deviation of the noise trace after being filtered with the + diode. power_mean and power_std can be calculated with the function + calculate_noise_parameters from utilities.diodeSimulator + triggered_channels: array of ints or None + the channels that are involved in the main phased beam + default: None, i.e. all channels + phasing_angles: array of floats or None + the angles for each beam + trigger_delays: dictionary + the delays for the channels that have caused a trigger. + If there is no trigger, it's an empty dictionary + output_passband: (float, float) tuple + Frequencies for a 6th-order Butterworth filter to be applied after + the diode filtering. + """ + Trigger.__init__(self, name, trigger_channels, 'digital_envelope_phased', pre_trigger_times=pre_trigger_times) + self._trigger_channels = trigger_channels + self._phasing_angles = phasing_angles + self._threshold = threshold + self._trigger_delays = trigger_delays + self._maximum_amps = maximum_amps + + +class PowerIntegrationPhasedTrigger(Trigger): - def __init__(self, name, threshold, channels=None, secondary_channels=None, + def __init__(self, name, threshold, trigger_channels=None, secondary_channels=None, primary_angles=None, secondary_angles=None, trigger_delays=None, sec_trigger_delays=None, window_size=None, step_size=None, @@ -319,8 +360,8 @@ def __init__(self, name, threshold, channels=None, secondary_channels=None, maximum_amps: list of floats (length equal to that of `phasing_angles`) the maximum value of all the integration windows for each of the phased waveforms """ - Trigger.__init__(self, name, channels, 'simple_phased', pre_trigger_times=pre_trigger_times) - self._primary_channels = channels + Trigger.__init__(self, name, trigger_channels, 'power_integration_phased', pre_trigger_times=pre_trigger_times) + self._primary_channels = trigger_channels self._primary_angles = primary_angles self._secondary_channels = secondary_channels self._secondary_angles = secondary_angles diff --git a/NuRadioReco/modules/phasedarray/phasedArray.py b/NuRadioReco/modules/phasedarray/phasedArray.py new file mode 100644 index 0000000000..83c8cacd38 --- /dev/null +++ b/NuRadioReco/modules/phasedarray/phasedArray.py @@ -0,0 +1,211 @@ +from NuRadioReco.modules.base.module import register_run +from NuRadioReco.utilities import units +import logging +import scipy +import numpy as np +from scipy import constants + +logger = logging.getLogger('phasedTriggerBase') + +cspeed = constants.c * units.m / units.s + +class phasedArray(): + """ + Base class for all phased array trigger modules. + """ + + def __init__(self, log_level=logging.WARNING): + self.__t = 0 + self.__pre_trigger_time = None + self.__debug = None + logger.setLevel(log_level) + self.begin() + + def begin(self, debug=False, pre_trigger_time=100 * units.ns): + self.__pre_trigger_time = pre_trigger_time + self.__debug = debug + + def get_antenna_positions(self, station, det, triggered_channels=None, component=2): + """ + Calculates the vertical coordinates of the antennas of the detector + + Parameters + ---------- + station: Station object + Description of the current station + det: Detector object + Description of the current detector + triggered_channels: array of ints + channels ids of the channels that form the primary phasing array + if None, all channels are taken + component: int + Which cartesian coordinate to return + + Returns + ------- + ant_pos: array of floatss + Desired antenna position in requested coordinate + + """ + + ant_pos = {} + for channel in station.iter_channels(use_channels=triggered_channels): + ant_pos[channel.get_id()] = det.get_relative_position(station.get_id(), channel.get_id())[component] + + return ant_pos + + def calculate_time_delays(self, station, det, + triggered_channels, + phasing_angles=None, + ref_index=1.75, + sampling_frequency=None): + + """ + Calculates the delays needed for phasing the array. + + Parameters + ---------- + station: Station object + Description of the current station + det: Detector object + Description of the current detector + triggered_channels: array of ints + channels ids of the channels that form the primary phasing array + if None, all channels are taken + phasing_angles: array of float + pointing angles for the primary beam + ref_index: float + refractive index for beam forming + sampling_frequency: float + Rate of the ADC used + + Returns + ------- + beam_rolls: array of dicts of keys=antenna and content=delay + """ + + if(triggered_channels is None): + triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] + + time_step = 1. / sampling_frequency + + ant_z = self.get_antenna_positions(station, det, triggered_channels, 2) + + self.check_vertical_string(station, det, triggered_channels) + ref_z = np.max(np.fromiter(ant_z.values(), dtype=float)) + + # Need to add in delay for trigger delay + cable_delays = {} + for channel in station.iter_trigger_channels(use_channels=triggered_channels): + cable_delays[channel.get_id()] = det.get_cable_delay(station.get_id(), channel.get_id()) + + beam_rolls = [] + for angle in phasing_angles: + + delays = [] + for key in ant_z: + delays += [-(ant_z[key] - ref_z) / cspeed * ref_index * np.sin(angle) - cable_delays[key]] + + delays -= np.max(delays) + roll = np.array(np.round(np.array(delays) / time_step)).astype(int) + subbeam_rolls = dict(zip(triggered_channels, roll)) + + # logger.debug("angle:", angle / units.deg) + # logger.debug(subbeam_rolls) + + beam_rolls.append(subbeam_rolls) + + return beam_rolls + + def get_channel_trace_start_time(self, station, triggered_channels): + """ + Finds the start time of the desired traces. + Throws an error if all the channels dont have the same start time. + + Parameters + ---------- + station: Station object + Description of the current station + triggered_channels: array of ints + channels ids of the channels that form the primary phasing array + if None, all channels are taken + + Returns + ------- + channel_trace_start_time: float + Channel start time + """ + + channel_trace_start_time = None + for channel in station.iter_trigger_channels(use_channels=triggered_channels): + + if channel_trace_start_time is None: + channel_trace_start_time = channel.get_trace_start_time() + + elif channel_trace_start_time != channel.get_trace_start_time(): + error_msg = 'Phased array channels do not have matching trace start times. ' + error_msg += 'This module is not prepared for this case.' + raise ValueError(error_msg) + + return channel_trace_start_time + + def check_vertical_string(self, station, det, triggered_channels): + """ + Checks if the triggering antennas lie in a straight vertical line + Throws error if not. + + Parameters + ---------- + station: Station object + Description of the current station + det: Detector object + Description of the current detector + triggered_channels: array of ints + channels ids of the channels that form the primary phasing array + if None, all channels are taken + """ + + cut = 1.e-3 * units.m + ant_x = np.fromiter(self.get_antenna_positions(station, det, triggered_channels, 0).values(), dtype=float) + diff_x = np.abs(ant_x - ant_x[0]) + ant_y = np.fromiter(self.get_antenna_positions(station, det, triggered_channels, 1).values(), dtype=float) + diff_y = np.abs(ant_y - ant_y[0]) + + if (sum(diff_x) > cut or sum(diff_y) > cut): + raise NotImplementedError('The phased triggering array should lie on a vertical line') + + def phase_signals(self, traces, beam_rolls): + """ + Phase signals together given the rolls + + Parameters + ---------- + traces: 2D array of floats + Signals from the antennas to be phased together. + beam_rolls: 2D array of floats + The amount to shift each signal before phasing the + traces together + + Returns + ------- + phased_traces: array of arrays + """ + + phased_traces = [[] for i in range(len(beam_rolls))] + + running_i = 0 + for subbeam_rolls in beam_rolls: + phased_trace = np.zeros(len(list(traces.values())[0])) + + for channel_id in traces: + + trace = traces[channel_id] + phased_trace += np.roll(trace, subbeam_rolls[channel_id]) + + phased_traces[running_i] = phased_trace + running_i += 1 + + return phased_traces + + def end(self): + pass diff --git a/NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py b/NuRadioReco/modules/trigger/analog_beamformed_envelope_trigger.py similarity index 66% rename from NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py rename to NuRadioReco/modules/trigger/analog_beamformed_envelope_trigger.py index b6808a307e..4fd8b7cbae 100644 --- a/NuRadioReco/modules/envelope_phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/trigger/analog_beamformed_envelope_trigger.py @@ -1,8 +1,7 @@ from NuRadioReco.modules.base.module import register_run from NuRadioReco.utilities import units -from NuRadioReco.framework.trigger import EnvelopePhasedTrigger -from NuRadioReco.modules.phasedarray.triggerSimulator import triggerSimulator as phasedTrigger -from NuRadioReco.modules.phasedarray.triggerSimulator import get_beam_rolls, get_channel_trace_start_time +from NuRadioReco.framework.trigger import AnalogEnvelopePhasedTrigger +from NuRadioReco.modules.phasedarray.phasedArray import phasedArray from NuRadioReco.utilities.diodeSimulator import diodeSimulator from NuRadioReco.modules.analogToDigitalConverter import analogToDigitalConverter import numpy as np @@ -17,7 +16,7 @@ default_angles = np.arcsin(np.linspace(np.sin(main_low_angle), np.sin(main_high_angle), 30)) -class triggerSimulator(phasedTrigger): +class triggerSimulator(phasedArray): """ Calculates the trigger for a envelope phased array. The channels that participate in the beam forming and the pointing angle for each @@ -27,18 +26,23 @@ class triggerSimulator(phasedTrigger): See https://arxiv.org/pdf/1903.11043.pdf and https://elog.phys.hawaii.edu/elog/anita_notes/080827_041639/powertrigger.pdf """ - + def envelope_trigger(self, station, det, - beam_rolls, + phasing_angles, + ref_index, triggered_channels, - threshold_factor, - power_mean, - power_std, + envelope_type="diode", + threshold_factor=None, + power_mean=None, + power_std=None, output_passband=(None, 200 * units.MHz), - cut_times=(None, None), - trigger_adc=False): + threshold=None, + trigger_adc=False, + apply_digitization=False, + adc_output="voltage" + ): """ Calculates the envelope trigger for a certain phasing configuration. Beams are formed. Then, each channel to be phased is filtered with a @@ -87,20 +91,20 @@ def envelope_trigger(self, station_id = station.get_id() diode = diodeSimulator(output_passband) + ADC = analogToDigitalConverter() traces = {} for channel in station.iter_channels(use_channels=triggered_channels): channel_id = channel.get_id() + adc_sampling_frequency = channel.get_sampling_rate() time_step = 1 / channel.get_sampling_rate() if trigger_adc: - - ADC = analogToDigitalConverter() trace = ADC.get_digital_trace(station, det, channel, - trigger_adc=trigger_adc, - random_clock_offset=True, - adc_type='perfect_floor_comparator', - diode=diode) + trigger_adc=trigger_adc, + random_clock_offset=True, + adc_type='perfect_floor_comparator', + diode=diode) time_step = 1 / det.get_channel(station_id, channel_id)['trigger_adc_sampling_frequency'] times = np.arange(len(trace), dtype=float) * time_step times += channel.get_trace_start_time() @@ -110,28 +114,32 @@ def envelope_trigger(self, trace = diode.tunnel_diode(channel) # get the enveloped trace times = np.copy(channel.get_times()) # get the corresponding time bins - if cut_times != (None, None): - left_bin = np.argmin(np.abs(times - cut_times[0])) - right_bin = np.argmin(np.abs(times - cut_times[1])) - trace[0:left_bin] = 0 - trace[right_bin:None] = 0 - traces[channel_id] = trace[:] - for subbeam_rolls in beam_rolls: + beam_rolls = self.calculate_time_delays(station, det, + triggered_channels, + phasing_angles, + ref_index=ref_index, + sampling_frequency=adc_sampling_frequency) + + phased_traces = self.phase_signals(traces, beam_rolls) - phased_trace = None - # Number of antennas: primary beam antennas - Nant = len(beam_rolls[0]) - for channel_id in traces: + trigger_time = None + trigger_times = {} + channel_trace_start_time = self.get_channel_trace_start_time(station, triggered_channels) + + trigger_delays = {} + maximum_amps = np.zeros(len(phased_traces)) + n_trigs = 0 + triggered_beams = [] + + for iTrace, phased_trace in enumerate(phased_traces): + is_triggered = False - trace = traces[channel_id] + # Number of antennas: primary beam antennas + Nant = len(beam_rolls[iTrace]) - if(phased_trace is None): - phased_trace = np.roll(trace, subbeam_rolls[channel_id]) - else: - phased_trace += np.roll(trace, subbeam_rolls[channel_id]) low_trigger = power_mean - power_std * np.abs(threshold_factor) low_trigger *= Nant @@ -139,27 +147,41 @@ def envelope_trigger(self, threshold_passed = np.min(phased_trace) < low_trigger if threshold_passed: - trigger_delays = {} - for channel_id in subbeam_rolls: - trigger_delays[channel_id] = subbeam_rolls[channel_id] * time_step + is_triggered = True + for channel_id in beam_rolls: + trigger_delays[channel_id] = beam_rolls[channel_id] * time_step logger.debug("Station has triggered") - return True, trigger_delays - return False, {} + triggered_beams.append(is_triggered) + + is_triggered=np.any(triggered_beams) + + if is_triggered: + logger.debug("Trigger condition satisfied!") + logger.debug("all trigger times", trigger_times) + trigger_time = min([x.min() for x in trigger_times.values()]) + logger.debug(f"minimum trigger time is {trigger_time:.0f}ns") + + return is_triggered, trigger_delays, trigger_time, trigger_times, maximum_amps, n_trigs, triggered_beams @register_run() def run(self, evt, station, det, - threshold_factor=6.5, - power_mean=None, - power_std=None, - triggered_channels=None, trigger_name='envelope_phased_threshold', + triggered_channels=None, phasing_angles=default_angles, set_not_triggered=False, ref_index=1.75, + envelope_type="diode", + threshold_factor=None, + power_mean=None, + power_std=None, output_passband=(None, 200 * units.MHz), - cut_times=(None, None), - trigger_adc=False): + threshold=None, + trigger_adc=False, + apply_digitization=False, + adc_output="voltage", + return_n_triggers=False + ): """ simulates phased array trigger for each event @@ -231,28 +253,39 @@ def run(self, evt, station, det, else: - channel_trace_start_time = get_channel_trace_start_time(station, triggered_channels) - logger.debug("primary channels:", triggered_channels) - beam_rolls = get_beam_rolls(station, det, triggered_channels, - phasing_angles, ref_index=ref_index) - - is_triggered, trigger_delays = self.envelope_trigger(station, det, beam_rolls, - triggered_channels, threshold_factor, - power_mean, power_std, output_passband, cut_times, - trigger_adc) - trigger = EnvelopePhasedTrigger(trigger_name, threshold_factor, power_mean, power_std, + is_triggered, trigger_delays, trigger_time, \ + trigger_times, n_triggers = self.envelope_trigger(station, + det, + phasing_angles, + ref_index, + triggered_channels=triggered_channels, + threshold_factor=threshold_factor, + power_mean=power_mean, + power_std=power_std, + output_passband=output_passband, + threshold=threshold, + trigger_adc=trigger_adc, + apply_digitization=apply_digitization, + adc_output=adc_output) + + trigger = AnalogEnvelopePhasedTrigger(trigger_name, threshold_factor, power_mean, power_std, triggered_channels, phasing_angles, trigger_delays, output_passband) trigger.set_triggered(is_triggered) + if is_triggered: - trigger.set_trigger_time(channel_trace_start_time) + #trigger_time(s)= time(s) from start of trace + start time of trace with respect to moment of first interaction = trigger time from moment of first interaction; time offset to interaction time (channel_trace_start_time) already recognized in self.phased_trigger + trigger.set_trigger_time(trigger_time) + trigger.set_trigger_times(trigger_times) else: trigger.set_trigger_time(None) + station.set_trigger(trigger) - return is_triggered + if return_n_triggers: + return is_triggered, n_triggers + else: + return is_triggered - def end(self): - pass diff --git a/NuRadioReco/modules/phasedarray/triggerSimulator.py b/NuRadioReco/modules/trigger/beamformed_power_integration.py similarity index 55% rename from NuRadioReco/modules/phasedarray/triggerSimulator.py rename to NuRadioReco/modules/trigger/beamformed_power_integration.py index a43549e54e..6e42895111 100644 --- a/NuRadioReco/modules/phasedarray/triggerSimulator.py +++ b/NuRadioReco/modules/trigger/beamformed_power_integration.py @@ -1,8 +1,9 @@ from NuRadioReco.modules.base.module import register_run from NuRadioReco.utilities import units - -from NuRadioReco.framework.trigger import SimplePhasedTrigger +from NuRadioReco.utilities import trace_utilities +from NuRadioReco.framework.trigger import PowerIntegrationPhasedTrigger from NuRadioReco.modules.analogToDigitalConverter import analogToDigitalConverter +from NuRadioReco.modules.phasedarray.phasedArray import phasedArray import logging import scipy import numpy as np @@ -18,7 +19,7 @@ default_angles = np.arcsin(np.linspace(np.sin(main_low_angle), np.sin(main_high_angle), 11)) -class triggerSimulator: +class triggerSimulator(phasedArray): """ Calculates the trigger for a phased array with a primary beam. @@ -28,166 +29,6 @@ class triggerSimulator: See https://arxiv.org/pdf/1809.04573.pdf """ - def __init__(self, log_level=logging.WARNING): - self.__t = 0 - self.__pre_trigger_time = None - self.__debug = None - logger.setLevel(log_level) - self.begin() - - def begin(self, debug=False, pre_trigger_time=100 * units.ns): - self.__pre_trigger_time = pre_trigger_time - self.__debug = debug - - def get_antenna_positions(self, station, det, triggered_channels=None, component=2): - """ - Calculates the vertical coordinates of the antennas of the detector - - Parameters - ---------- - station: Station object - Description of the current station - det: Detector object - Description of the current detector - triggered_channels: array of ints - channels ids of the channels that form the primary phasing array - if None, all channels are taken - component: int - Which cartesian coordinate to return - - Returns - ------- - ant_pos: array of floatss - Desired antenna position in requested coordinate - - """ - - ant_pos = {} - for channel in station.iter_channels(use_channels=triggered_channels): - ant_pos[channel.get_id()] = det.get_relative_position(station.get_id(), channel.get_id())[component] - - return ant_pos - - def calculate_time_delays(self, station, det, - triggered_channels, - phasing_angles=default_angles, - ref_index=1.75, - sampling_frequency=None): - - """ - Calculates the delays needed for phasing the array. - - Parameters - ---------- - station: Station object - Description of the current station - det: Detector object - Description of the current detector - triggered_channels: array of ints - channels ids of the channels that form the primary phasing array - if None, all channels are taken - phasing_angles: array of float - pointing angles for the primary beam - ref_index: float - refractive index for beam forming - sampling_frequency: float - Rate of the ADC used - - Returns - ------- - beam_rolls: array of dicts of keys=antenna and content=delay - """ - - if(triggered_channels is None): - triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] - - time_step = 1. / sampling_frequency - - ant_z = self.get_antenna_positions(station, det, triggered_channels, 2) - - self.check_vertical_string(station, det, triggered_channels) - ref_z = np.max(np.fromiter(ant_z.values(), dtype=float)) - - # Need to add in delay for trigger delay - cable_delays = {} - for channel in station.iter_trigger_channels(use_channels=triggered_channels): - cable_delays[channel.get_id()] = det.get_cable_delay(station.get_id(), channel.get_id()) - - beam_rolls = [] - for angle in phasing_angles: - - delays = [] - for key in ant_z: - delays += [-(ant_z[key] - ref_z) / cspeed * ref_index * np.sin(angle) - cable_delays[key]] - - delays -= np.max(delays) - roll = np.array(np.round(np.array(delays) / time_step)).astype(int) - subbeam_rolls = dict(zip(triggered_channels, roll)) - - # logger.debug("angle:", angle / units.deg) - # logger.debug(subbeam_rolls) - - beam_rolls.append(subbeam_rolls) - - return beam_rolls - - def get_channel_trace_start_time(self, station, triggered_channels): - """ - Finds the start time of the desired traces. - Throws an error if all the channels dont have the same start time. - - Parameters - ---------- - station: Station object - Description of the current station - triggered_channels: array of ints - channels ids of the channels that form the primary phasing array - if None, all channels are taken - - Returns - ------- - channel_trace_start_time: float - Channel start time - """ - - channel_trace_start_time = None - for channel in station.iter_trigger_channels(use_channels=triggered_channels): - - if channel_trace_start_time is None: - channel_trace_start_time = channel.get_trace_start_time() - - elif channel_trace_start_time != channel.get_trace_start_time(): - error_msg = 'Phased array channels do not have matching trace start times. ' - error_msg += 'This module is not prepared for this case.' - raise ValueError(error_msg) - - return channel_trace_start_time - - def check_vertical_string(self, station, det, triggered_channels): - """ - Checks if the triggering antennas lie in a straight vertical line - Throws error if not. - - Parameters - ---------- - station: Station object - Description of the current station - det: Detector object - Description of the current detector - triggered_channels: array of ints - channels ids of the channels that form the primary phasing array - if None, all channels are taken - """ - - cut = 1.e-3 * units.m - ant_x = np.fromiter(self.get_antenna_positions(station, det, triggered_channels, 0).values(), dtype=float) - diff_x = np.abs(ant_x - ant_x[0]) - ant_y = np.fromiter(self.get_antenna_positions(station, det, triggered_channels, 1).values(), dtype=float) - diff_y = np.abs(ant_y - ant_y[0]) - - if (sum(diff_x) > cut or sum(diff_y) > cut): - raise NotImplementedError('The phased triggering array should lie on a vertical line') - def power_sum(self, coh_sum, window, step, adc_output='voltage',rnog_like=False): """ Calculate power summed over a length defined by 'window', overlapping at intervals defined by 'step' @@ -227,9 +68,10 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage',rnog_like=False) coh_sum_squared = (coh_sum * coh_sum).astype(float) else: #(adc_output == 'counts'): if rnog_like: - coh_sum[coh_sum>2**6-1] = 2**6 - 1 - coh_sum[coh_sum<-2**6] = -2**6 - coh_sum_squared = (coh_sum * coh_sum).astype(int) + sat_bit = 7 + coh_sum[coh_sum>2**sat_bit-1] = 2**sat_bit - 1 + coh_sum[coh_sum<-2**sat_bit] = -2**sat_bit + coh_sum_squared = (coh_sum * coh_sum) coh_sum_windowed = np.lib.stride_tricks.as_strided(coh_sum_squared, (num_frames, window), (coh_sum_squared.strides[0] * step, coh_sum_squared.strides[0])) @@ -237,79 +79,14 @@ def power_sum(self, coh_sum, window, step, adc_output='voltage',rnog_like=False) power = np.sum(coh_sum_windowed, axis=1) return_power=power.astype(float) / window - if adc_output=='counts': return_power = np.floor(return_power) + if adc_output=='counts': return_power = np.round(return_power) return return_power, num_frames - def hilbert_envelope(self, coh_sum, adc_output='voltage', coeff_gain=1, rnog_like=False): - - if rnog_like: - coeff_gain = 128 - - #31 sample fir transformer - #hil=[ -0.0424413 , 0. , -0.0489708 , 0. , -0.0578745 , 0. , -0.0707355 , 0. , -0.0909457 , 0. , -0.127324 , 0. - # , -0.2122066 , 0. , -0.6366198 , 0., 0.6366198 , 0. , 0.2122066 , 0. , 0.127324 ,0. , 0.0909457 , 0. , .0707355 - # , 0. , 0.0578745 , 0. , 0.0489708 , 0. , 0.0424413 ] - - #middle 15 coefficients ^ - hil = np.array([ -0.0909457 , 0. , -0.127324 , 0. , -0.2122066 , 0. , -0.6366198 , 0. , 0.6366198 , 0. , 0.2122066 , - 0. , 0.127324 , 0. , 0.0909457 ]) - - if coeff_gain!=1: - hil = np.round(hil * coeff_gain) / coeff_gain - - imag_an = np.convolve(coh_sum, hil, mode='full')[len(hil)//2 : len(coh_sum) + len(hil)//2] - - if adc_output: - imag_an = np.round(imag_an) - - if rnog_like: - envelope = np.max(np.array((coh_sum,imag_an)), axis=0) + (3 / 8) * min(np.array((coh_sum,imag_an)), axis=0) - else: - envelope = np.sqrt(coh_sum**2 + imag_an**2) - - if adc_output=='counts': - envelope = np.round(envelope) - - return envelope - - def phase_signals(self, traces, beam_rolls): - """ - Phase signals together given the rolls - - Parameters - ---------- - traces: 2D array of floats - Signals from the antennas to be phased together. - beam_rolls: 2D array of floats - The amount to shift each signal before phasing the - traces together - - Returns - ------- - phased_traces: array of arrays - """ - - phased_traces = [[] for i in range(len(beam_rolls))] - - running_i = 0 - for subbeam_rolls in beam_rolls: - phased_trace = np.zeros(len(list(traces.values())[0])) - - for channel_id in traces: - - trace = traces[channel_id] - phased_trace += np.roll(trace, subbeam_rolls[channel_id]) - - phased_traces[running_i] = phased_trace - running_i += 1 - - return phased_traces - def phased_trigger(self, station, det, Vrms=None, threshold=60 * units.mV, - triggered_channels=None, + trigger_channels=None, phasing_angles=default_angles, ref_index=1.75, trigger_adc=False, # by default, assumes the trigger ADC is the same as the channels ADC @@ -330,7 +107,7 @@ def phased_trigger(self, station, det, Several channels are phased by delaying their signals by an amount given by a pointing angle. Several pointing angles are possible in order to cover - the sky. The array triggered_channels controls the channels that are phased, + the sky. The array trigger_channels controls the channels that are phased, according to the angles phasing_angles. Parameters @@ -345,7 +122,7 @@ def phased_trigger(self, station, det, int the detector description file. threshold: float threshold above (or below) a trigger is issued, absolute amplitude - triggered_channels: array of ints + trigger_channels: array of ints channels ids of the channels that form the primary phasing array if None, all channels are taken phasing_angles: array of float @@ -411,8 +188,8 @@ def phased_trigger(self, station, det, list of bools for which beams triggered """ - if(triggered_channels is None): - triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] + if(trigger_channels is None): + trigger_channels = [channel.get_id() for channel in station.iter_trigger_channels()] if(adc_output != 'voltage' and adc_output != 'counts'): error_msg = 'ADC output type must be "counts" or "voltage". Currently set to:' + str(adc_output) @@ -423,10 +200,10 @@ def phased_trigger(self, station, det, is_triggered = False trigger_delays = {} - logger.debug(f"trigger channels: {triggered_channels}") + logger.debug(f"trigger channels: {trigger_channels}") traces = {} - for channel in station.iter_trigger_channels(use_channels=triggered_channels): + for channel in station.iter_trigger_channels(use_channels=trigger_channels): channel_id = channel.get_id() trace = np.array(channel.get_trace()) @@ -450,42 +227,9 @@ def phased_trigger(self, station, det, raise ValueError("Could not convert upsampling_factor to integer. Exiting.") if(upsampling_factor >= 2): - - new_len = len(trace) * upsampling_factor - cur_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/adc_sampling_frequency) - new_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/adc_sampling_frequency/upsampling_factor) - - if upsampling_method=='fft': - upsampled_trace = scipy.signal.resample(trace, new_len) - - elif upsampling_method=='lin': - upsampled_trace = np.interp(new_t,cur_t,trace) - - elif upsampling_method=='fir': - if rnog_like: - #I think I can do better than this in firmware - up_filt = np.array([ 0.0, 0.0 , 0.0 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , -0.0 , -0.015625 , - -0.0390625 , -0.03125 , 0.0 , 0.0703125 , 0.15625 , 0.2265625 , 0.25 , 0.2265625 , 0.15625 , - 0.0703125 , 0.0 , -0.03125 , -0.0390625 , -0.015625 , 0.0 , 0.0078125 , 0.0078125 , 0.0078125 , - 0.0 , 0.0 , 0.0 , 0.0 ]) - else: - cutoff = .5 - base_freq_filter_length = 8 - filter_length = base_freq_filter_length * upsampling_factor - 1 - up_filt = scipy.signal.firwin(filter_length, adc_sampling_frequency*cutoff, pass_zero='lowpass', - fs=adc_sampling_frequency*upsampling_factor) - if coeff_gain!=1: - up_filt = np.round(up_filt*coeff_gain)/coeff_gain - - zero_pad = np.zeros(len(trace)*upsampling_factor) - zero_pad[::upsampling_factor] = trace[:] - upsampled_trace = np.convolve(zero_pad, up_filt, mode='full')[len(up_filt)//2 : len(zero_pad) + len(up_filt)//2] * upsampling_factor - - else: - error_msg = 'Interpolation method must be lin, fft, or fir' - raise NotImplementedError(error_msg) - - if adc_output=='counts' and rnog_like==True: upsampled_trace=np.trunc(upsampled_trace) + upsampled_trace, new_sampling_frequency = trace_utilities.digital_upsampling(trace, adc_sampling_frequency, upsampling_method=upsampling_method, + upsampling_factor=upsampling_factor, rnog_like=rnog_like, coeff_gain=coeff_gain, + adc_output='voltage', filter_taps=31) # If upsampled is performed, the final sampling frequency changes trace = upsampled_trace[:] @@ -493,14 +237,13 @@ def phased_trigger(self, station, det, if(len(trace) % 2 == 1): trace = trace[:-1] - adc_sampling_frequency *= upsampling_factor - - time_step = 1.0 / adc_sampling_frequency - traces[channel_id] = trace[:] + + adc_sampling_frequency *= upsampling_factor + time_step = 1.0 / adc_sampling_frequency beam_rolls = self.calculate_time_delays(station, det, - triggered_channels, + trigger_channels, phasing_angles, ref_index=ref_index, sampling_frequency=adc_sampling_frequency) @@ -509,7 +252,7 @@ def phased_trigger(self, station, det, trigger_time = None trigger_times = {} - channel_trace_start_time = self.get_channel_trace_start_time(station, triggered_channels) + channel_trace_start_time = self.get_channel_trace_start_time(station, trigger_channels) trigger_delays = {} maximum_amps = np.zeros(len(phased_traces)) @@ -519,40 +262,24 @@ def phased_trigger(self, station, det, for iTrace, phased_trace in enumerate(phased_traces): is_triggered=False - if trig_type=='power_integration': - squared_mean, num_frames = self.power_sum(coh_sum=phased_trace, window=window, step=step, adc_output=adc_output, rnog_like=rnog_like) - maximum_amps[iTrace] = np.max(squared_mean) - - if True in (squared_mean > threshold): - n_trigs += len(np.where((squared_mean > threshold)==True)[0]) - trigger_delays[iTrace] = {} - - for channel_id in beam_rolls[iTrace]: - trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step - - triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(squared_mean > threshold))) - logger.debug(f"Station has triggered, at bins {triggered_bins}") - logger.debug(trigger_delays) - logger.debug(f"trigger_delays {trigger_delays[iTrace][triggered_channels[0]]}") - is_triggered = True - trigger_times[iTrace] = trigger_delays[iTrace][triggered_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time - logger.debug(f"trigger times = {trigger_times[iTrace]}") - - elif trig_type=='envelope': - hilbert_env = self.hilbert_envelope(coh_sum=phased_trace, adc_output=adc_output, rnog_like=rnog_like) - maximum_amps[iTrace] = np.max(hilbert_env) - - if True in (hilbert_env>threshold): - n_trigs += len(np.where((hilbert_env > threshold)==True)[0]) - trigger_delays[iTrace] = {} - for channel_id in beam_rolls[iTrace]: - trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step - triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(hilbert_env > threshold))) - is_triggered = True - trigger_times[iTrace] = trigger_delays[iTrace][triggered_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time - else: - raise NotImplementedError("not a good trigger type: options (power_integration, envelope)") + squared_mean, num_frames = self.power_sum(coh_sum=phased_trace, window=window, step=step, adc_output=adc_output, rnog_like=rnog_like) + maximum_amps[iTrace] = np.max(squared_mean) + + if True in (squared_mean > threshold): + n_trigs += len(np.where((squared_mean > threshold)==True)[0]) + trigger_delays[iTrace] = {} + + for channel_id in beam_rolls[iTrace]: + trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step + + triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(squared_mean > threshold))) + logger.debug(f"Station has triggered, at bins {triggered_bins}") + logger.debug(trigger_delays) + logger.debug(f"trigger_delays {trigger_delays[iTrace][trigger_channels[0]]}") + is_triggered = True + trigger_times[iTrace] = trigger_delays[iTrace][trigger_channels[0]] + triggered_bins * step * time_step + channel_trace_start_time + logger.debug(f"trigger times = {trigger_times[iTrace]}") triggered_beams.append(is_triggered) @@ -570,7 +297,7 @@ def phased_trigger(self, station, det, def run(self, evt, station, det, Vrms=None, threshold=60 * units.mV, - triggered_channels=None, + trigger_channels=None, trigger_name='simple_phased_threshold', phasing_angles=default_angles, set_not_triggered=False, @@ -586,7 +313,6 @@ def run(self, evt, station, det, upsampling_method='fft', coeff_gain=128, rnog_like=False, - trig_type='power_integration', return_n_triggers=False ): @@ -595,7 +321,7 @@ def run(self, evt, station, det, Several channels are phased by delaying their signals by an amount given by a pointing angle. Several pointing angles are possible in order to cover - the sky. The array triggered_channels controls the channels that are phased, + the sky. The array trigger_channels controls the channels that are phased, according to the angles phasing_angles. Parameters @@ -612,7 +338,7 @@ def run(self, evt, station, det, int the detector description file. threshold: float threshold above (or below) a trigger is issued, absolute amplitude - triggered_channels: array of ints + trigger_channels: array of ints channels ids of the channels that form the primary phasing array if None, all channels are taken trigger_name: string @@ -658,11 +384,6 @@ def run(self, evt, station, det, int. rnog_like: bool (default False) If true, this will apply the RNO-G FLOWER based math/rounding done in firmware. - trig_type: str (default "power_integration") - - "power_integration" do the power integration for the given window size and - step length - - "envelope" perform a hilbrt envelope threshold trigger on the beamformed - traces return_n_triggers: bool (default False) To better estimate simulated thresholds one should count the total triggers in the entire trace for each beam. If true, this return the total trigger number. @@ -675,8 +396,8 @@ def run(self, evt, station, det, Count of the total number of triggers in all beamformed traces """ - if(triggered_channels is None): - triggered_channels = [channel.get_id() for channel in station.iter_trigger_channels()] + if(trigger_channels is None): + trigger_channels = [channel.get_id() for channel in station.iter_trigger_channels()] if(adc_output != 'voltage' and adc_output != 'counts'): error_msg = 'ADC output type must be "counts" or "voltage". Currently set to:' + str(adc_output) @@ -695,7 +416,7 @@ def run(self, evt, station, det, det=det, Vrms=Vrms, threshold=threshold, - triggered_channels=triggered_channels, + trigger_channels=trigger_channels, phasing_angles=phasing_angles, ref_index=ref_index, trigger_adc=trigger_adc, @@ -708,41 +429,28 @@ def run(self, evt, station, det, apply_digitization=apply_digitization, upsampling_method=upsampling_method, coeff_gain=coeff_gain, - rnog_like=rnog_like, - trig_type=trig_type + rnog_like=rnog_like ) # Create a trigger object to be returned to the station - if trig_type=='power_integration': - trigger = SimplePhasedTrigger( - trigger_name, - threshold, - channels=triggered_channels, - primary_angles=phasing_angles, - trigger_delays=trigger_delays, - window_size=window, - step_size=step, - maximum_amps=maximum_amps - ) - - elif trig_type=='envelope': - trigger = SimplePhasedTrigger( - trigger_name, - threshold, - channels=triggered_channels, - primary_angles=phasing_angles, - trigger_delays=trigger_delays, - maximum_amps=maximum_amps - ) - else: - raise NotImplementedError("invalid phased trigger type") + trigger = PowerIntegrationPhasedTrigger( + trigger_name, + threshold, + trigger_channels=trigger_channels, + primary_angles=phasing_angles, + trigger_delays=trigger_delays, + window_size=window, + step_size=step, + maximum_amps=maximum_amps + ) + trigger.set_triggered(is_triggered) if is_triggered: #trigger_time(s)= time(s) from start of trace + start time of trace with respect to moment of first interaction = trigger time from moment of first interaction; time offset to interaction time (channel_trace_start_time) already recognized in self.phased_trigger - trigger.set_trigger_time(trigger_time)# + trigger.set_trigger_time(trigger_time) trigger.set_trigger_times(trigger_times) else: trigger.set_trigger_time(None) @@ -750,9 +458,7 @@ def run(self, evt, station, det, station.set_trigger(trigger) if return_n_triggers: - return is_triggered,n_triggers + return is_triggered, n_triggers else: return is_triggered - def end(self): - pass diff --git a/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py b/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py new file mode 100644 index 0000000000..4154956074 --- /dev/null +++ b/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py @@ -0,0 +1,411 @@ +from NuRadioReco.modules.base.module import register_run +from NuRadioReco.utilities import units +from NuRadioReco.utilities import trace_utilities +from NuRadioReco.framework.trigger import DigitalEnvelopePhasedTrigger +from NuRadioReco.modules.phasedarray.phasedArray import phasedArray +from NuRadioReco.modules.analogToDigitalConverter import analogToDigitalConverter +import logging +import scipy +import numpy as np +from scipy import constants +from scipy.signal import hilbert + +logger = logging.getLogger('phasedEnvelopeTriggerSimulator') + +cspeed = constants.c * units.m / units.s + +main_low_angle = np.deg2rad(-55.0) +main_high_angle = -1.0 * main_low_angle +default_angles = np.arcsin(np.linspace(np.sin(main_low_angle), np.sin(main_high_angle), 11)) + + +class triggerSimulator(phasedArray): + """ + Calculates the trigger for a phased array with primary beams. + + The channels that participate in both beams and the pointing angle for each + subbeam can be specified. + + Calculates Hilbert Envelope using digital FIR filters. + """ + + def hilbert_envelope(self, coh_sum, adc_output='voltage', coeff_gain=1, rnog_like=False): + + if rnog_like: + coeff_gain = 128 + + #31 sample fir transformer + #hil=[ -0.0424413 , 0. , -0.0489708 , 0. , -0.0578745 , 0. , -0.0707355 , 0. , -0.0909457 , 0. , -0.127324 , 0. + # , -0.2122066 , 0. , -0.6366198 , 0., 0.6366198 , 0. , 0.2122066 , 0. , 0.127324 ,0. , 0.0909457 , 0. , .0707355 + # , 0. , 0.0578745 , 0. , 0.0489708 , 0. , 0.0424413 ] + + #middle 15 coefficients ^ + hil = np.array([ -0.0909457 , 0. , -0.127324 , 0. , -0.2122066 , 0. , -0.6366198 , 0. , 0.6366198 , 0. , 0.2122066 , + 0. , 0.127324 , 0. , 0.0909457 ]) + + if coeff_gain!=1: + hil = np.round(hil * coeff_gain) / coeff_gain + + imag_an = np.convolve(coh_sum, hil, mode='full')[len(hil)//2 : len(coh_sum) + len(hil)//2] + #imag_an=np.imag(hilbert(coh_sum)) + if adc_output=='counts': + imag_an = np.round(imag_an) + + if rnog_like: + envelope = np.max(np.array((coh_sum,imag_an)), axis=0) + (3 / 8) * np.min(np.array((coh_sum,imag_an)), axis=0) + else: + envelope = np.sqrt(coh_sum**2 + imag_an**2) + + if adc_output=='counts': + envelope = np.round(envelope) + + return envelope + + def envelope_trigger(self, station, det, + Vrms=None, + threshold=60 * units.mV, + trigger_channels=None, + phasing_angles=default_angles, + ref_index=1.75, + trigger_adc=False, # by default, assumes the trigger ADC is the same as the channels ADC + clock_offset=0, + adc_output='voltage', + trigger_filter=None, + upsampling_factor=1, + apply_digitization=False, + upsampling_method='fft', + coeff_gain=128, + rnog_like=False, + trig_type='power_integration' + ): + """ + simulates phased array trigger for each event + + Several channels are phased by delaying their signals by an amount given + by a pointing angle. Several pointing angles are possible in order to cover + the sky. The array trigger_channels controls the channels that are phased, + according to the angles phasing_angles. + + Parameters + ---------- + station: Station object + Description of the current station + det: Detector object + Description of the current detector + Vrms: float + RMS of the noise on a channel, used to automatically create the digitizer + reference voltage. If set to None, tries to use reference voltage as defined + int the detector description file. + threshold: float + threshold above (or below) a trigger is issued, absolute amplitude + trigger_channels: array of ints + channels ids of the channels that form the primary phasing array + if None, all channels are taken + phasing_angles: array of float + pointing angles for the primary beam + ref_index: float (default 1.75) + refractive index for beam forming + trigger_adc: bool (default True) + If True, uses the ADC settings from the trigger. It must be specified in the + detector file. See analogToDigitalConverter module for information + (see option `apply_digitization`) + clock_offset: float (default 0) + Overall clock offset, for adc clock jitter reasons (see `apply_digitization`) + adc_output: string (default 'voltage') + - 'voltage' to store the ADC output as discretised voltage trace + - 'counts' to store the ADC output in ADC counts and apply integer based math + trigger_filter: array floats (default None) + Freq. domain of the response to be applied to post-ADC traces + Must be length for "MC freq" + upsampling_factor: integer (default 1) + Upsampling factor. The trace will be a upsampled to a + sampling frequency int_factor times higher than the original one + after conversion to digital + apply_digitization: bool (default True) + Perform the quantization of the ADC. If set to true, should also set options + `trigger_adc`, `adc_output`, `clock_offset` + upsampling_method: str (default 'fft') + Choose between FFT, FIR, or Linear Interpolaion based upsampling methods + coeff_gain: int (default 1) + If using the FIR upsampling, this will convert the floating point output of the + scipy filter to a fixed point value by multiplying by this factor and rounding to an + int. + rnog_like: bool (default False) + If true, this will apply the RNO-G FLOWER based math/rounding done in firmware. + + Returns + ------- + is_triggered: bool + True if the triggering condition is met + trigger_delays: dictionary + the delays for the primary channels that have caused a trigger. + If there is no trigger, it's an empty dictionary + trigger_time: float + the earliest trigger time with respect to first interaction time. + trigger_times: dictionary + all time bins that fulfil the trigger condition per beam. The key is the beam number. Time with respect to first interaction time. + maximum_amps: list of floats (length equal to that of `phasing_angles`) + the maximum value of all the integration windows for each of the phased waveforms + n_trigs: int + total number of triggers that happened for all beams across the full traces + triggered_beams: list + list of bools for which beams triggered + """ + + if(trigger_channels is None): + trigger_channels = [channel.get_id() for channel in station.iter_trigger_channels()] + + if(adc_output != 'voltage' and adc_output != 'counts'): + error_msg = 'ADC output type must be "counts" or "voltage". Currently set to:' + str(adc_output) + raise ValueError(error_msg) + + ADC = analogToDigitalConverter() + + is_triggered = False + trigger_delays = {} + + logger.debug(f"trigger channels: {trigger_channels}") + + traces = {} + for channel in station.iter_trigger_channels(use_channels=trigger_channels): + channel_id = channel.get_id() + + trace = np.array(channel.get_trace()) + + if apply_digitization: + trace, adc_sampling_frequency = ADC.get_digital_trace(station, det, channel, + Vrms=Vrms, + trigger_adc=trigger_adc, + clock_offset=clock_offset, + return_sampling_frequency=True, + adc_type='perfect_floor_comparator', + adc_output=adc_output, + trigger_filter=None) + else: + adc_sampling_frequency = channel.get_sampling_rate() + + if not isinstance(upsampling_factor, int): + try: + upsampling_factor = int(upsampling_factor) + except: + raise ValueError("Could not convert upsampling_factor to integer. Exiting.") + + if(upsampling_factor >= 2): + upsampled_trace, new_sampling_frequency = trace_utilities.digital_upsampling(trace, adc_sampling_frequency, upsampling_method=upsampling_method, + upsampling_factor=upsampling_factor, rnog_like=rnog_like, coeff_gain=coeff_gain, + adc_output='voltage', filter_taps=31) + + # If upsampled is performed, the final sampling frequency changes + trace = upsampled_trace[:] + + if(len(trace) % 2 == 1): + trace = trace[:-1] + + traces[channel_id] = trace[:] + + adc_sampling_frequency *= upsampling_factor + + time_step = 1.0 / adc_sampling_frequency + beam_rolls = self.calculate_time_delays(station, det, + trigger_channels, + phasing_angles, + ref_index=ref_index, + sampling_frequency=adc_sampling_frequency) + + phased_traces = self.phase_signals(traces, beam_rolls) + + trigger_time = None + trigger_times = {} + channel_trace_start_time = self.get_channel_trace_start_time(station, trigger_channels) + + trigger_delays = {} + maximum_amps = np.zeros(len(phased_traces)) + n_trigs = 0 + triggered_beams = [] + + for iTrace, phased_trace in enumerate(phased_traces): + is_triggered=False + hilbert_env = self.hilbert_envelope(coh_sum=phased_trace, adc_output=adc_output, rnog_like=rnog_like) + maximum_amps[iTrace] = np.max(hilbert_env) + + if True in (hilbert_env>threshold): + n_trigs += len(np.where((hilbert_env > threshold)==True)[0]) + trigger_delays[iTrace] = {} + for channel_id in beam_rolls[iTrace]: + trigger_delays[iTrace][channel_id] = beam_rolls[iTrace][channel_id] * time_step + triggered_bins = np.atleast_1d(np.squeeze(np.argwhere(hilbert_env > threshold))) + is_triggered = True + trigger_times[iTrace] = trigger_delays[iTrace][trigger_channels[0]] + triggered_bins * time_step + channel_trace_start_time + + triggered_beams.append(is_triggered) + + is_triggered=np.any(triggered_beams) + + if is_triggered: + logger.debug("Trigger condition satisfied!") + logger.debug("all trigger times", trigger_times) + trigger_time = min([x.min() for x in trigger_times.values()]) + logger.debug(f"minimum trigger time is {trigger_time:.0f}ns") + + return is_triggered, trigger_delays, trigger_time, trigger_times, maximum_amps, n_trigs, triggered_beams + + @register_run() + def run(self, evt, station, det, + Vrms=None, + threshold=60 * units.mV, + trigger_channels=None, + trigger_name='digital_envelope_phased_threshold', + phasing_angles=default_angles, + set_not_triggered=False, + ref_index=1.75, + trigger_adc=False, # by default, assumes the trigger ADC is the same as the channels ADC + clock_offset=0, + adc_output='voltage', + trigger_filter=None, + upsampling_factor=1, + apply_digitization=True, + upsampling_method='fft', + coeff_gain=128, + rnog_like=False, + return_n_triggers=False + ): + + """ + simulates phased array trigger for each event + + Several channels are phased by delaying their signals by an amount given + by a pointing angle. Several pointing angles are possible in order to cover + the sky. The array triggered_channels controls the channels that are phased, + according to the angles phasing_angles. + + Parameters + ---------- + evt: Event object + Description of the current event + station: Station object + Description of the current station + det: Detector object + Description of the current detector + Vrms: float + RMS of the noise on a channel, used to automatically create the digitizer + reference voltage. If set to None, tries to use reference voltage as defined + int the detector description file. + threshold: float + threshold above (or below) a trigger is issued, absolute amplitude + trigger_channels: array of ints + channels ids of the channels that form the primary phasing array + if None, all channels are taken + trigger_name: string + name for the trigger + phasing_angles: array of float + pointing angles for the primary beam + set_not_triggered: bool (default False) + if True not trigger simulation will be performed and this trigger will be set to not_triggered + ref_index: float (default 1.75) + refractive index for beam forming + trigger_adc: bool, (default True) + If True, uses the ADC settings from the trigger. It must be specified in the + detector file. See analogToDigitalConverter module for information + clock_offset: float (default 0) + Overall clock offset, for adc clock jitter reasons + adc_output: string (default 'voltage') + + - 'voltage' to store the ADC output as discretised voltage trace + - 'counts' to store the ADC output in ADC counts + + trigger_filter: array floats (default None) + Freq. domain of the response to be applied to post-ADC traces + Must be length for "MC freq" + upsampling_factor: integer (default 1) + Upsampling factor. The trace will be a upsampled to a + sampling frequency int_factor times higher than the original one + after conversion to digital + window: int (default 32) + Power integral window + Units of ADC time ticks + step: int (default 16) + Time step in power integral. If equal to window, there is no time overlap + in between neighboring integration windows. + Units of ADC time ticks + apply_digitization: bool (default True) + Perform the quantization of the ADC. If set to true, should also set options + `trigger_adc`, `adc_output`, `clock_offset` + upsampling_method: str (default 'fft') + Choose between FFT, FIR, or Linear Interpolaion based upsampling methods + coeff_gain: int (default 1) + If using the FIR upsampling, this will convert the floating point output of the + scipy filter to a fixed point value by multiplying by this factor and rounding to an + int. + rnog_like: bool (default False) + If true, this will apply the RNO-G FLOWER based math/rounding done in firmware. + return_n_triggers: bool (default False) + To better estimate simulated thresholds one should count the total triggers + in the entire trace for each beam. If true, this return the total trigger number. + + Returns + ------- + is_triggered: bool + True if the triggering condition is met + n_triggers: int (Optional) + Count of the total number of triggers in all beamformed traces + """ + + if(trigger_channels is None): + trigger_channels = [channel.get_id() for channel in station.iter_trigger_channels()] + + if(adc_output != 'voltage' and adc_output != 'counts'): + error_msg = 'ADC output type must be "counts" or "voltage". Currently set to:' + str(adc_output) + raise ValueError(error_msg) + + is_triggered = False + trigger_delays = {} + + if(set_not_triggered): + is_triggered = False + trigger_delays = {} + triggered_beams = [] + else: + is_triggered, trigger_delays, trigger_time, trigger_times,\ + maximum_amps, n_triggers, triggered_beams = self.envelope_trigger(station=station, + det=det, + Vrms=Vrms, + threshold=threshold, + trigger_channels=trigger_channels, + phasing_angles=phasing_angles, + ref_index=ref_index, + trigger_adc=trigger_adc, + clock_offset=clock_offset, + adc_output=adc_output, + trigger_filter=trigger_filter, + upsampling_factor=upsampling_factor, + apply_digitization=apply_digitization, + upsampling_method=upsampling_method, + coeff_gain=coeff_gain, + rnog_like=rnog_like + ) + + # Create a trigger object to be returned to the station + trigger = DigitalEnvelopePhasedTrigger( + trigger_name, + threshold, + trigger_channels=trigger_channels, + phasing_angles=phasing_angles, + trigger_delays=trigger_delays, + maximum_amps=maximum_amps + ) + + trigger.set_triggered(is_triggered) + + if is_triggered: + #trigger_time(s)= time(s) from start of trace + start time of trace with respect to moment of first interaction = trigger time from moment of first interaction; time offset to interaction time (channel_trace_start_time) already recognized in self.phased_trigger + trigger.set_trigger_time(trigger_time) + trigger.set_trigger_times(trigger_times) + else: + trigger.set_trigger_time(None) + + station.set_trigger(trigger) + + if return_n_triggers: + return is_triggered, n_triggers + else: + return is_triggered diff --git a/NuRadioReco/modules/trigger/highLowThreshold.py b/NuRadioReco/modules/trigger/highLowThreshold.py index aed46d4d0f..b761c09d97 100644 --- a/NuRadioReco/modules/trigger/highLowThreshold.py +++ b/NuRadioReco/modules/trigger/highLowThreshold.py @@ -253,6 +253,8 @@ def run(self, evt, station, det, self.__t += time.time() - t + return has_triggered + def end(self): from datetime import timedelta dt = timedelta(seconds=self.__t) diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 23ea274cdf..49f1e8bca5 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -188,7 +188,81 @@ def get_stokes(trace_u, trace_v, window_samples=128, squeeze=True): return np.squeeze(stokes) return stokes -def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 ** 7): +def digital_upsampling(trace, adc_sampling_frequency, upsampling_method='fft', upsampling_factor=2, coeff_gain=1, filter_taps=45, adc_output='voltage', rnog_like=False): + """ + Digital upsampling with various methods and settings. + + Parameters + ---------- + trace : 1d array (float or int) + Input trace to upsample + adc_sampling_frequency : float + Original sampling frequency for trace + upsampling_method: str (default 'fft') + Choose between FFT, FIR, or Linear Interpolaion based upsampling methods + upsampling_factor : float (default 2) + The factor which the sampling frequency increases + coeff_gain: int (default 1) + If using the FIR upsampling, this will convert the floating point output of the + scipy filter to a fixed point value by multiplying by this factor and rounding to an int. + If set to 1, this will preserve the float value of the filter coefficients. + filter_taps : int (default 45) + Number of taps in the FIR filter in FIR-based upsampling. + adc_output: string (default "voltage") + - "voltage" : keep upsampled trace as floats + - "counts" : round upsampling trace to ints + rnog_like: bool (default False) + If true, this will apply the RNO-G FLOWER upsampling settings and the fixed point math/rounding done in firmware. + Returns + ------- + upsampled_trace : 1d array (float or int) + Upsampled trace at the new sampling frequency + new_sampling_frequency : float + New sampling frequency + """ + + if (np.abs(int(upsampling_factor) - upsampling_factor) > 1e-3): + warning_msg = "The input upsampling factor does not seem to be close to an integer." + warning_msg += "It has been rounded to {}".format(int(upsampling_factor)) + logger.warning(warning_msg) + + upsampling_factor = int(upsampling_factor) + + if (upsampling_factor <= 1): + error_msg = "Upsampling factor is less or equal to 1. Upsampling will not be performed." + upsampled_trace = trace + new_sampling_freq = adc_sampling_frequency + + else: + new_sampling_freq = adc_sampling_frequency * upsampling_factor + new_len = len(trace) * upsampling_factor + + if upsampling_method=='fft': + upsampled_trace = scipy.signal.resample(trace, new_len) + + elif upsampling_method=='lin': + cur_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/adc_sampling_frequency) + new_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/new_sampling_freq) + upsampled_trace = np.interp(new_t, cur_t, trace) + + elif upsampling_method=='fir': + + if rnog_like: + filter_taps = 45 + coeff_gain = 128 + + upsampled_trace=upsampling_fir(trace, adc_sampling_frequency, int_factor=upsampling_factor, ntaps=filter_taps, coeff_gain=coeff_gain) + + else: + error_msg = 'Interpolation method must be lin, fft, or fir' + raise NotImplementedError(error_msg) + + if adc_output=='counts': + upsampled_trace = np.round(upsampled_trace) + + return upsampled_trace, new_sampling_freq + +def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 ** 7, coeff_gain=128): """ This function performs an upsampling by inserting a number of zeroes between samples and then applying a finite impulse response (FIR) filter. @@ -212,27 +286,15 @@ def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 ** The upsampled trace """ - if (np.abs(int(int_factor) - int_factor) > 1e-3): - warning_msg = "The input upsampling factor does not seem to be close to an integer." - warning_msg += "It has been rounded to {}".format(int(int_factor)) - logger.warning(warning_msg) - - int_factor = int(int_factor) - - if (int_factor <= 1): - error_msg = "Upsampling factor is less or equal to 1. Upsampling will not be performed." - raise ValueError(error_msg) - - zeroed_trace = np.zeros(len(trace) * int_factor) - for i_point, point in enumerate(trace[:-1]): - zeroed_trace[i_point * int_factor] = point - - upsampled_delta_time = 1 / (int_factor * original_sampling_frequency) - upsampled_times = np.arange(0, len(zeroed_trace) * upsampled_delta_time, upsampled_delta_time) + cutoff = 0.5 + up_filt = scipy.signal.firwin(ntaps, original_sampling_frequency*cutoff, pass_zero='lowpass', + fs=original_sampling_frequency*int_factor) + if coeff_gain!=1: + up_filt = np.round(up_filt*coeff_gain)/coeff_gain - cutoff = 1. / int_factor - fir_coeffs = scipy.signal.firwin(ntaps, cutoff, window='boxcar') - upsampled_trace = np.convolve(zeroed_trace, fir_coeffs)[:len(upsampled_times)] * int_factor + zero_padded_sig = np.zeros(len(trace)*int_factor) + zero_padded_sig[::int_factor] = trace[:] + upsampled_trace = np.convolve(zero_padded_sig, up_filt, mode='full')[len(up_filt)//2 : len(zero_padded_sig) + len(up_filt)//2] * int_factor return upsampled_trace From b0f77ce6914a834f9a93253db04351ff919b1661 Mon Sep 17 00:00:00 2001 From: Ryan Krebs Date: Mon, 20 Jan 2025 11:24:44 -0500 Subject: [PATCH 4/7] small things --- NuRadioMC/utilities/Veff.py | 5 ++++- .../modules/trigger/digital_beamformed_envelope_trigger.py | 1 - 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/NuRadioMC/utilities/Veff.py b/NuRadioMC/utilities/Veff.py index b272e487f5..1a70bd7ee2 100644 --- a/NuRadioMC/utilities/Veff.py +++ b/NuRadioMC/utilities/Veff.py @@ -276,7 +276,10 @@ def get_Veff_Aeff_single( else: raise AttributeError(f"attributes do neither contain volume nor area") - Vrms = fin.attrs['Vrms'] + try: + Vrms = fin.attrs['Vrms'] + except: + Vrms = 1 # Solid angle needed for the effective volume calculations out['domega'] = np.abs(phimax - phimin) * np.abs(np.cos(thetamin) - np.cos(thetamax)) diff --git a/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py b/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py index 4154956074..6d90cee1b6 100644 --- a/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py +++ b/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py @@ -76,7 +76,6 @@ def envelope_trigger(self, station, det, upsampling_method='fft', coeff_gain=128, rnog_like=False, - trig_type='power_integration' ): """ simulates phased array trigger for each event From bda57df0394e0e6bf55da0c581b6154332291353 Mon Sep 17 00:00:00 2001 From: Ryan Krebs Date: Mon, 20 Jan 2025 11:26:38 -0500 Subject: [PATCH 5/7] add sim script to compare trigger performances --- .../generate_data_set.py | 349 ++++++++++++++++++ 1 file changed, 349 insertions(+) create mode 100644 NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py diff --git a/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py b/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py new file mode 100644 index 0000000000..e3ef2d2725 --- /dev/null +++ b/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py @@ -0,0 +1,349 @@ +#!/bin/env python3 + +import argparse +import copy +import logging +import numpy as np +import os +import secrets +import datetime as dt +from scipy import constants + +from NuRadioMC.utilities import runner + +from NuRadioMC.EvtGen import generator +from NuRadioMC.simulation import simulation +from NuRadioReco.utilities import units + +from NuRadioReco.detector.RNO_G import rnog_detector + +from NuRadioReco.modules.RNO_G import hardwareResponseIncorporator, triggerBoardResponse +from NuRadioReco.modules.trigger import highLowThreshold +from NuRadioReco.modules.trigger import digital_beamformed_envelope_trigger +from NuRadioReco.modules.trigger import beamformed_power_integration + +import NuRadioReco.modules.channelGenericNoiseAdder + +import matplotlib.pyplot as plt + +root_seed = secrets.randbits(128) +deep_trigger_channels = np.array([0, 1, 2, 3]) + +main_low_angle = np.deg2rad(-60) +main_high_angle = np.deg2rad(60) +pa_channels = [0, 1, 2, 3] + +phasing_angles = np.arcsin(np.linspace(np.sin(main_low_angle), np.sin(main_high_angle), 12)) + +upsampling_factor = 4 +window = 24 #int(16 * units.ns * new_sampling_rate * upsampling_factor) +step = 8 #int(8 * units.ns * new_sampling_rate * upsampling_factor) + + + +def get_vrms_from_temperature_for_trigger_channels(det, station_id, trigger_channels, temperature): + + vrms_per_channel = [] + for channel_id in trigger_channels: + resp = det.get_signal_chain_response(station_id, channel_id, trigger=True) + + freqs = np.linspace(10, 1200, 1000) * units.MHz + filt = resp(freqs) + + # Calculation of Vrms. For details see from elog:1566 and https://en.wikipedia.org/wiki/Johnson%E2%80%93Nyquist_noise + # (last two Eqs. in "noise voltage and power" section) or our wiki https://nu-radio.github.io/NuRadioMC/NuRadioMC/pages/HDF5_structure.html + + # Bandwidth, i.e., \Delta f in equation + integrated_channel_response = np.trapz(np.abs(filt) ** 2, freqs) + + vrms_per_channel.append( + (temperature * 50 * constants.k * integrated_channel_response / units.Hz) ** 0.5 + ) + + return vrms_per_channel + + +def get_fiducial_volume(energy): + # Fiducial volume for a Greenland station. From Martin: https://radio.uchicago.edu/wiki/images/2/26/TriggerSimulation_May2023.pdf + + # key: log10(E), value: radius in km + max_radius_shallow = { + 16.25: 1.5, 16.5: 2.1, 16.75: 2.7, 17.0: 3.1, 17.25: 3.7, 17.5: 3.9, 17.75: 4.4, + 18.00: 4.8, 18.25: 5.1, 18.50: 5.25, 18.75: 5.3, 19.0: 5.6, 100: 6.1, + } + + # key: log10(E), value: depth in km + min_z_shallow = { + 16.25: -0.65, 16.50: -0.8, 16.75: -1.2, 17.00: -1.5, 17.25: -1.7, 17.50: -2.0, + 17.75: -2.1, 18.00: -2.3, 18.25: -2.4, 18.50: -2.55, 100: -2.7, + } + + + def get_limits(dic, E): + # find all energy bins which are higher than E + idx = np.arange(len(dic))[E - 10 ** np.array(list(dic.keys())) * units.eV <= 0] + assert len(idx), f"Energy {E} is too high. Max energy is {10 ** np.amax(dic.keys()):.1e}." + + # take the lowest energy bin which is higher than E + return np.array(list(dic.values()))[np.amin(idx)] * units.km + + r_max = get_limits(max_radius_shallow, energy) + print(f"Maximum radius {r_max}") + z_min = get_limits(min_z_shallow, energy) + print(f"Depth {z_min}") + + volume = { + "fiducial_rmax": r_max, + "fiducial_rmin": 0 * units.km, + "fiducial_zmin": z_min, + "fiducial_zmax": 0 + } + + return volume + + +def RNO_G_HighLow_Thresh(lgRate_per_hz): + # Thresholds calculated using the RNO-G hardware (iglu + flower_lp) + # This applies for the VPol antennas + # parameterization comes from Alan: https://radio.uchicago.edu/wiki/images/e/e6/2023.10.11_Simulating_RNO-G_Trigger.pdf + return (-859 + np.sqrt(39392706 - 3602500 * lgRate_per_hz)) / 1441.0 + + +class mySimulation(simulation.simulation): + + def __init__(self, *args, **kwargs): + # this module is needed in super().__init__ to calculate the vrms + self.rnogHarwareResponse = hardwareResponseIncorporator.hardwareResponseIncorporator() + self.rnogHarwareResponse.begin(trigger_channels=deep_trigger_channels) + + super().__init__(*args, **kwargs) + self.logger = logging.getLogger("NuRadioMC.RNOG_trigger_simulation") + self.deep_trigger_channels = deep_trigger_channels + + + self.highLowThreshold = highLowThreshold.triggerSimulator() + self.rnogADCResponse = triggerBoardResponse.triggerBoardResponse() + self.rnogADCResponse.begin(adc_input_range=2 * units.volt, clock_offset=0.0, adc_output="counts") + + # future TODO: Add noise + self.channel_generic_noise_adder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() + self.channel_generic_noise_adder.begin(seed=root_seed) + self.powerTrigger = beamformed_power_integration.triggerSimulator(log_level=logging.WARNING) + self.envelopeTrigger = digital_beamformed_envelope_trigger.triggerSimulator(log_level=logging.WARNING) + + self.output_mode = {'Channels': self._config['output']['channel_traces'], + 'ElectricFields': self._config['output']['electric_field_traces'], + 'SimChannels': self._config['output']['sim_channel_traces'], + 'SimElectricFields': self._config['output']['sim_electric_field_traces']} + + self.high_low_trigger_thresholds = { + "1Hz": RNO_G_HighLow_Thresh(0) + } + + + def _detector_simulation_filter_amp(self, evt, station, det): + # apply the amplifiers and filters to get to RADIANT-level + self.rnogHarwareResponse.run(evt, station, det, sim_to_data=True) + + def _detector_simulation_trigger(self, evt, station, det): + + vrms_input_to_adc = get_vrms_from_temperature_for_trigger_channels( + det, station.get_id(), self.deep_trigger_channels, 300) + + sampling_rate = det.get_sampling_frequency(station.get_id()) + self.logger.info(f'Radiant sampling rate is {sampling_rate / units.MHz:.1f} MHz') + + # Runs the FLOWER board response + vrms_after_gain, gains = self.rnogADCResponse.run( + evt, station, det, trigger_channels=self.deep_trigger_channels, + vrms=None, digitize_trace=True, apply_adc_gain=True + ) + + # this is only returning the correct value if digitize_trace=True for self.rnogADCResponse.run(..) + flower_sampling_rate = station.get_trigger_channel(self.deep_trigger_channels[0]).get_sampling_rate() + self.logger.info(f'Flower sampling rate is {flower_sampling_rate / units.MHz:.1f} MHz') + + for thresh_key, threshold in self.high_low_trigger_thresholds.items(): + + threshold_high = {channel_id: threshold * vrms for channel_id, vrms in zip(self.deep_trigger_channels, vrms_after_gain)} + threshold_low = {channel_id: -1 * threshold * vrms for channel_id, vrms in zip(self.deep_trigger_channels, vrms_after_gain)} + + self.highLowThreshold.run( + evt, + station, + det, + threshold_high=threshold_high, + threshold_low=threshold_low, + use_digitization=False, #the trace has already been digitized with the rnogADCResponse + high_low_window=6 / flower_sampling_rate, + coinc_window=20 / flower_sampling_rate, + number_concidences=2, + triggered_channels=self.deep_trigger_channels, + trigger_name=f"deep_high_low_{thresh_key}", + ) + + lvrms=0 + for i in pa_channels: + lvrms += vrms_after_gain[i]**2 #coherently summed rms^2 + + self.powerTrigger.run( + evt, + station, + det, + Vrms=None, + threshold=11.84*lvrms, + trigger_channels=pa_channels, + phasing_angles=phasing_angles, + ref_index=1.75, + trigger_name=f"PA_power", + trigger_adc=False, + adc_output=f"counts", + trigger_filter=None, + upsampling_factor=4, + window=window, + step=step, + apply_digitization=False, + upsampling_method='fir', + coeff_gain=256, + rnog_like=True + ) + + lvrms=np.sqrt(lvrms) + + self.envelopeTrigger.run( + evt, + station, + det, + Vrms=None, + threshold=6.61*lvrms, + trigger_channels=pa_channels, + phasing_angles=phasing_angles, + ref_index=1.75, + trigger_name=f"PA_envelope", + trigger_adc=False, + adc_output=f"counts", + trigger_filter=None, + upsampling_factor=1, + apply_digitization=False, + upsampling_method='fir', + coeff_gain=128, + rnog_like=True + ) + + +def task(q, iSim, output_filename, **kwargs): + + defaults = { + "trigger_adc_sampling_frequency": 0.472, + "trigger_adc_nbits": 8, + "trigger_adc_noise_nbits": 2.585, + } + + det = rnog_detector.Detector( + detector_file=args.detectordescription, log_level=logging.INFO, + always_query_entire_description=False, select_stations=args.station_id, + over_write_handset_values=defaults) + + det.update(dt.datetime(2023, 8, 3)) + + volume = get_fiducial_volume(args.energy) + # Simulate fiducial volume around station + pos = det.get_absolute_position(args.station_id) + print(f"Simulating around center x0={pos[0]:.2f}m, y0={pos[1]:.2f}m") + volume.update({"x0": pos[0], "y0": pos[1]}) + + flavor_ids = {"e": [12, -12], "mu": [14, -14], "tau": [16, -16], "all": [12, 14, 16, -12, -14, -16]} + run_proposal = args.proposal and ("cc" in args.interaction_type) and (args.flavor in ["mu", "tau", "all"]) + if run_proposal: + print(f"Using PROPOSAL for simulation of {args.flavor} {args.interaction_type}") + + input_data = generator.generate_eventlist_cylinder( + "on-the-fly", + kwargs["n_events_per_file"], + args.energy, args.energy, + volume, + start_event_id=args.index * args.n_events_per_file + 1, + flavor=flavor_ids[args.flavor], + n_events_per_file=None, + deposited=False, + proposal=run_proposal, + proposal_config="Greenland", + start_file_id=0, + log_level=None, + proposal_kwargs={}, + max_n_events_batch=args.n_events_per_file, + write_events=False, + seed=root_seed + args.index, + interaction_type=args.interaction_type, + ) + + if args.nur_output: + nur_output_filename = output_filename.replace(".hdf5", ".nur") + else: + nur_output_filename = None + + sim = mySimulation( + inputfilename=input_data, + outputfilename=output_filename, + det=det, + evt_time=dt.datetime(2023, 8, 3), + outputfilenameNuRadioReco=None, #nur_output_filename, + config_file=args.config, + trigger_channels=deep_trigger_channels, + # log_level=logging.INFO, + ) + + n_trig=sim.run() + print(f"simulation pass {iSim} with {n_trig} events", flush=True) + q.put(n_trig) + +if __name__ == "__main__": + + ABS_PATH_HERE = str(os.path.dirname(os.path.realpath(__file__))) + def_data_dir = os.path.join(ABS_PATH_HERE, "data") + default_config_path = os.path.join(ABS_PATH_HERE, "../07_RNO_G_simulation/RNO_config.yaml") + + parser = argparse.ArgumentParser(description="Run NuRadioMC simulation") + # Sim steering arguments + parser.add_argument("--config", type=str, default="config.yaml", help="NuRadioMC yaml config file") + parser.add_argument("--detectordescription", '--det', type=str, default=None, help="Path to RNO-G detector description file. If None, query from DB") + parser.add_argument("--station_id", type=int, default=11, help="Set station to be used for simulation") + + # Neutrino arguments + parser.add_argument("--energy", '-e', default=1e18, type=float, help="Neutrino energy [eV]") + parser.add_argument("--flavor", '-f', default="all", type=str, help="the flavor") + parser.add_argument("--interaction_type", '-it', default="ccnc", type=str, help="interaction type cc, nc or ccnc") + + # File meta-variables + parser.add_argument("--index", '-i', default=1, type=int, help="counter to create a unique data-set identifier") + parser.add_argument("--n_events_per_file", '-n', type=int, default=1e4, help="Number of nu-interactions per file") + parser.add_argument("--data_dir", type=str, default=def_data_dir, help="directory name where the library will be created") + parser.add_argument("--proposal", action="store_true", help="Use PROPOSAL for simulation") + parser.add_argument("--nur_output", action="store_true", help="Write nur files.") + parser.add_argument("--n_cpus", type=int, default=1, help="Number of cores for parallel running") + parser.add_argument("--n_triggers", type=int, default=1e2, help="Number of total events to generate") + # Defaults for the trigger simulation which are not yet in the hardware DB + + args = parser.parse_args() + kwargs = args.__dict__ + assert args.station_id is not None, "Please specify a station id with `--station_id`" + + + + output_path = f"{args.data_dir}/station_{args.station_id}/nu_{args.flavor}_{args.interaction_type}" + + if not os.path.exists(output_path): + print("Making dirs", output_path) + os.makedirs(output_path, exist_ok=True) + + output_filename = f"{output_path}/{args.flavor}_{args.interaction_type}_1e{np.log10(args.energy):.2f}eV_{args.index:08d}.hdf5" + + class myrunner(runner.NuRadioMCRunner): + def get_outputfilename(self): + return os.path.join( + self.output_path, f"{np.log10(args.energy):.2f}_{self.kwargs['index']:06d}_{self.i_task:06d}.hdf5" + ) + + # start simulating a library across the chosen number of cpus. Each CPU will only run for 1 day + r = myrunner(args.n_cpus, task, output_path, max_runtime=3600 * 24, n_triggers_max=args.n_triggers, kwargs=kwargs) + r.run() \ No newline at end of file From bf059d43f6dfa10782cff7f25af21aee7fa3d2b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felix=20Schl=C3=BCter?= Date: Tue, 21 Jan 2025 09:48:37 +0100 Subject: [PATCH 6/7] Minor clean up --- NuRadioReco/detector/detector_base.py | 7 +- .../modules/RNO_G/triggerBoardResponse.py | 27 ++++---- NuRadioReco/utilities/trace_utilities.py | 64 +++++++++++-------- 3 files changed, 53 insertions(+), 45 deletions(-) diff --git a/NuRadioReco/detector/detector_base.py b/NuRadioReco/detector/detector_base.py index 4ca809a287..22a6eaaddc 100644 --- a/NuRadioReco/detector/detector_base.py +++ b/NuRadioReco/detector/detector_base.py @@ -878,7 +878,7 @@ def get_amplifier_response(self, station_id, channel_id, frequencies): def get_trigger_amplifier_type(self, station_id, channel_id): """ - returns the type of the amplifier + Returns the type of the amplifier used only for trigger channels. Parameters ---------- @@ -887,7 +887,10 @@ def get_trigger_amplifier_type(self, station_id, channel_id): channel_id: int the channel id - Returns string + Returns + ------- + amp_type: string + The type/name of the amplifier """ res = self.__get_channel(station_id, channel_id) return res['trigger_amp_type'] diff --git a/NuRadioReco/modules/RNO_G/triggerBoardResponse.py b/NuRadioReco/modules/RNO_G/triggerBoardResponse.py index 6aa81b9900..d4920b32a2 100644 --- a/NuRadioReco/modules/RNO_G/triggerBoardResponse.py +++ b/NuRadioReco/modules/RNO_G/triggerBoardResponse.py @@ -1,11 +1,9 @@ -import logging -import numpy as np - from NuRadioReco.modules.base.module import register_run from NuRadioReco.modules.analogToDigitalConverter import analogToDigitalConverter from NuRadioReco.utilities import units, fft -#from scipy.signal import resample, firwin, hilbert +import numpy as np +import logging logger = logging.getLogger("NuRadioReco.triggerBoardResponse") @@ -62,24 +60,23 @@ def get_vrms(self, station, trigger_channels, trace_split=20): Station to use trigger_channels : list Channels that this function should be applied to - trace_split : int (default: 9) + trace_split : int (default: 20) How many chunks each of the waveforms will be split into before calculating the standard deviation Returns ------- - approx_vrms : float - the median RMS voltage of the waveforms - + vrms : list of floats + RMS voltage of the waveforms """ - vrms=[] + vrms = [] for channel_id in trigger_channels: channel = station.get_trigger_channel(channel_id) trace = np.array(channel.get_trace()) trace = trace[: int(trace_split * int(len(trace) / trace_split))].reshape((trace_split, -1)) approx_vrms = np.median(np.std(trace, axis=1)) - logger.debug(f" Ch: {channel_id}\tObser Vrms: {approx_vrms / units.mV:0.3f} mV") + logger.debug(f"\tCh {channel_id}:\tobs. Vrms {approx_vrms / units.mV:0.3f} mV") vrms.append(approx_vrms) self.logger.debug(vrms) @@ -104,17 +101,15 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None, gain_val If set to `None`, this will be estimated using the waveforms gain_values : list (default: None) If set these will be applied to the channel. Otherwise gains are recalculated. + Returns ------- vrms_after_gain : list the RMS voltage of the waveforms after the gain has been applied - ideal_vrms: float the ideal vrms, as measured on the ADC capacitors - ret_gain_values: list gain values applied to each channel - """ if avg_vrms is None: @@ -140,7 +135,7 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None, gain_val vrms_after_gain.append(vrms * gain_values[channel_id]) channel = station.get_trigger_channel(channel_id) channel.set_trace(channel.get_trace() * gain_values[channel_id], channel.get_sampling_rate()) - + ret_gain_values.append(gain_values[channel_id]) else: msg = f"\t Ch: {channel_id}\t Target Vrms: {ideal_vrms / units.mV:0.3f} mV" @@ -168,7 +163,7 @@ def apply_adc_gain(self, station, det, trigger_channels, avg_vrms=None, gain_val self.logger.debug(f"\t Used Vrms: {vrms_after_gain[-1] / units.mV:0.3f} mV" + f"\tADC Gain {gain_to_use}") self.logger.debug(f"\t Eff noise bits: {eff_noise_bits:0.2f}\tRequested: {noise_bits}") - return np.array(vrms_after_gain), ideal_vrms, ret_gain_values + return np.array(vrms_after_gain), ideal_vrms, np.array(ret_gain_values) def digitize_trace(self, station, det, trigger_channels, vrms): for channel_id in trigger_channels: @@ -294,7 +289,7 @@ def end(self): integrated_channel_response = np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) rel_channel_response = np.trapz(np.abs(four_filters_highres[i]) ** 2, fff) bandwidth[i] = integrated_channel_response - Vrms_ratio[i] = np.sqrt(rel_channel_response / (max_freq - min_freq)) + Vrms_ratio[i] = np.sqrt(rel_channel_response / (max_freq - min_freq)) chan_vrms = (noise_temp * 50 * constants.k * integrated_channel_response / units.Hz) ** 0.5 per_channel_vrms.append(chan_vrms) amplitude[i] = chan_vrms / Vrms_ratio[i] diff --git a/NuRadioReco/utilities/trace_utilities.py b/NuRadioReco/utilities/trace_utilities.py index 49f1e8bca5..ed15efdb0f 100644 --- a/NuRadioReco/utilities/trace_utilities.py +++ b/NuRadioReco/utilities/trace_utilities.py @@ -1,10 +1,9 @@ +from NuRadioReco.utilities import units, ice, fft, geometryUtilities as geo_utl + import numpy as np import scipy.constants import scipy.signal -from NuRadioReco.utilities import units -from NuRadioReco.utilities import ice -from NuRadioReco.utilities import geometryUtilities as geo_utl -from NuRadioReco.utilities import fft + import logging logger = logging.getLogger('NuRadioReco.trace_utilities') @@ -20,7 +19,7 @@ def get_efield_antenna_factor(station, frequencies, channels, detector, zenith, Parameters ---------- - + station: Station frequencies: array of complex frequencies of the radio signal for which the antenna response is needed @@ -188,7 +187,10 @@ def get_stokes(trace_u, trace_v, window_samples=128, squeeze=True): return np.squeeze(stokes) return stokes -def digital_upsampling(trace, adc_sampling_frequency, upsampling_method='fft', upsampling_factor=2, coeff_gain=1, filter_taps=45, adc_output='voltage', rnog_like=False): +def digital_upsampling( + trace, adc_sampling_frequency, upsampling_method='fft', + upsampling_factor=2, coeff_gain=1, filter_taps=45, + adc_output='voltage', rnog_like=False): """ Digital upsampling with various methods and settings. @@ -203,16 +205,21 @@ def digital_upsampling(trace, adc_sampling_frequency, upsampling_method='fft', u upsampling_factor : float (default 2) The factor which the sampling frequency increases coeff_gain: int (default 1) - If using the FIR upsampling, this will convert the floating point output of the + If using the FIR upsampling, this will convert the floating point output of the scipy filter to a fixed point value by multiplying by this factor and rounding to an int. If set to 1, this will preserve the float value of the filter coefficients. filter_taps : int (default 45) Number of taps in the FIR filter in FIR-based upsampling. adc_output: string (default "voltage") - - "voltage" : keep upsampled trace as floats - - "counts" : round upsampling trace to ints + Output type of the upsampled trace: + + - "voltage" : keep upsampled trace as floats + - "counts" : round upsampling trace to ints + rnog_like: bool (default False) - If true, this will apply the RNO-G FLOWER upsampling settings and the fixed point math/rounding done in firmware. + If true, this will apply the RNO-G FLOWER upsampling settings and the fixed + point math/rounding done in firmware. + Returns ------- upsampled_trace : 1d array (float or int) @@ -221,37 +228,38 @@ def digital_upsampling(trace, adc_sampling_frequency, upsampling_method='fft', u New sampling frequency """ - if (np.abs(int(upsampling_factor) - upsampling_factor) > 1e-3): + if np.abs(int(upsampling_factor) - upsampling_factor) > 1e-3: warning_msg = "The input upsampling factor does not seem to be close to an integer." warning_msg += "It has been rounded to {}".format(int(upsampling_factor)) logger.warning(warning_msg) upsampling_factor = int(upsampling_factor) - if (upsampling_factor <= 1): + if upsampling_factor <= 1: error_msg = "Upsampling factor is less or equal to 1. Upsampling will not be performed." upsampled_trace = trace new_sampling_freq = adc_sampling_frequency - + else: new_sampling_freq = adc_sampling_frequency * upsampling_factor - new_len = len(trace) * upsampling_factor + new_len = len(trace) * upsampling_factor - if upsampling_method=='fft': + if upsampling_method == 'fft': upsampled_trace = scipy.signal.resample(trace, new_len) - elif upsampling_method=='lin': - cur_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/adc_sampling_frequency) - new_t = np.arange(0, 1/adc_sampling_frequency*len(trace), 1/new_sampling_freq) + elif upsampling_method == 'lin': + cur_t = np.arange(0, 1 / adc_sampling_frequency * len(trace), 1 / adc_sampling_frequency) + new_t = np.arange(0, 1 / adc_sampling_frequency * len(trace), 1 / new_sampling_freq) upsampled_trace = np.interp(new_t, cur_t, trace) - elif upsampling_method=='fir': + elif upsampling_method == 'fir': if rnog_like: filter_taps = 45 coeff_gain = 128 - upsampled_trace=upsampling_fir(trace, adc_sampling_frequency, int_factor=upsampling_factor, ntaps=filter_taps, coeff_gain=coeff_gain) + upsampled_trace = upsampling_fir( + trace, adc_sampling_frequency, int_factor=upsampling_factor, ntaps=filter_taps, coeff_gain=coeff_gain) else: error_msg = 'Interpolation method must be lin, fft, or fir' @@ -262,14 +270,13 @@ def digital_upsampling(trace, adc_sampling_frequency, upsampling_method='fft', u return upsampled_trace, new_sampling_freq -def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 ** 7, coeff_gain=128): +def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2**7, coeff_gain=128): """ This function performs an upsampling by inserting a number of zeroes between samples and then applying a finite impulse response (FIR) filter. Parameters ---------- - trace: array of floats Trace to be upsampled original_sampling_frequency: float @@ -287,14 +294,17 @@ def upsampling_fir(trace, original_sampling_frequency, int_factor=2, ntaps=2 ** """ cutoff = 0.5 - up_filt = scipy.signal.firwin(ntaps, original_sampling_frequency*cutoff, pass_zero='lowpass', - fs=original_sampling_frequency*int_factor) - if coeff_gain!=1: - up_filt = np.round(up_filt*coeff_gain)/coeff_gain + up_filt = scipy.signal.firwin( + ntaps, original_sampling_frequency*cutoff, pass_zero='lowpass', + fs=original_sampling_frequency*int_factor) + + if coeff_gain != 1: + up_filt = np.round(up_filt*coeff_gain) / coeff_gain zero_padded_sig = np.zeros(len(trace)*int_factor) zero_padded_sig[::int_factor] = trace[:] - upsampled_trace = np.convolve(zero_padded_sig, up_filt, mode='full')[len(up_filt)//2 : len(zero_padded_sig) + len(up_filt)//2] * int_factor + upsampled_trace = np.convolve( + zero_padded_sig, up_filt, mode='full')[len(up_filt) // 2 : len(zero_padded_sig) + len(up_filt) // 2] * int_factor return upsampled_trace From 1b147c96ed20ed047da03cb64c7ee7d03727d272 Mon Sep 17 00:00:00 2001 From: Ryan Krebs Date: Tue, 21 Jan 2025 21:33:38 -0500 Subject: [PATCH 7/7] fix module naming --- .../08_RNO_G_trigger_simulation/generate_data_set.py | 8 ++++---- .../analogBeamformedEnvelopeTrigger.py} | 0 .../beamformedPowerIntegrationTrigger.py} | 0 .../digitalBeamformedEnvelopeTrigger.py} | 0 4 files changed, 4 insertions(+), 4 deletions(-) rename NuRadioReco/modules/{trigger/analog_beamformed_envelope_trigger.py => phasedarray/analogBeamformedEnvelopeTrigger.py} (100%) rename NuRadioReco/modules/{trigger/beamformed_power_integration.py => phasedarray/beamformedPowerIntegrationTrigger.py} (100%) rename NuRadioReco/modules/{trigger/digital_beamformed_envelope_trigger.py => phasedarray/digitalBeamformedEnvelopeTrigger.py} (100%) diff --git a/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py b/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py index e3ef2d2725..064f9f78a1 100644 --- a/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py +++ b/NuRadioMC/examples/08_RNO_G_trigger_simulation/generate_data_set.py @@ -19,8 +19,8 @@ from NuRadioReco.modules.RNO_G import hardwareResponseIncorporator, triggerBoardResponse from NuRadioReco.modules.trigger import highLowThreshold -from NuRadioReco.modules.trigger import digital_beamformed_envelope_trigger -from NuRadioReco.modules.trigger import beamformed_power_integration +from NuRadioReco.modules.phasedarray import digitalBeamformedEnvelopeTrigger +from NuRadioReco.modules.phasedarray import beamformedPowerIntegrationTrigger import NuRadioReco.modules.channelGenericNoiseAdder @@ -128,8 +128,8 @@ def __init__(self, *args, **kwargs): # future TODO: Add noise self.channel_generic_noise_adder = NuRadioReco.modules.channelGenericNoiseAdder.channelGenericNoiseAdder() self.channel_generic_noise_adder.begin(seed=root_seed) - self.powerTrigger = beamformed_power_integration.triggerSimulator(log_level=logging.WARNING) - self.envelopeTrigger = digital_beamformed_envelope_trigger.triggerSimulator(log_level=logging.WARNING) + self.powerTrigger = beamformedPowerIntegrationTrigger.triggerSimulator(log_level=logging.WARNING) + self.envelopeTrigger = digitalBeamformedEnvelopeTrigger.triggerSimulator(log_level=logging.WARNING) self.output_mode = {'Channels': self._config['output']['channel_traces'], 'ElectricFields': self._config['output']['electric_field_traces'], diff --git a/NuRadioReco/modules/trigger/analog_beamformed_envelope_trigger.py b/NuRadioReco/modules/phasedarray/analogBeamformedEnvelopeTrigger.py similarity index 100% rename from NuRadioReco/modules/trigger/analog_beamformed_envelope_trigger.py rename to NuRadioReco/modules/phasedarray/analogBeamformedEnvelopeTrigger.py diff --git a/NuRadioReco/modules/trigger/beamformed_power_integration.py b/NuRadioReco/modules/phasedarray/beamformedPowerIntegrationTrigger.py similarity index 100% rename from NuRadioReco/modules/trigger/beamformed_power_integration.py rename to NuRadioReco/modules/phasedarray/beamformedPowerIntegrationTrigger.py diff --git a/NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py b/NuRadioReco/modules/phasedarray/digitalBeamformedEnvelopeTrigger.py similarity index 100% rename from NuRadioReco/modules/trigger/digital_beamformed_envelope_trigger.py rename to NuRadioReco/modules/phasedarray/digitalBeamformedEnvelopeTrigger.py