diff --git a/NuRadioReco/examples/RNO_data/calculate_cable_delays/coax_delays/plots_coax_delays.py b/NuRadioReco/examples/RNOG/calculate_cable_delays/coax_delays/plots_coax_delays.py similarity index 100% rename from NuRadioReco/examples/RNO_data/calculate_cable_delays/coax_delays/plots_coax_delays.py rename to NuRadioReco/examples/RNOG/calculate_cable_delays/coax_delays/plots_coax_delays.py diff --git a/NuRadioReco/examples/RNO_data/calculate_cable_delays/fiber_delays/plots_fiber_delays.py b/NuRadioReco/examples/RNOG/calculate_cable_delays/fiber_delays/plots_fiber_delays.py similarity index 100% rename from NuRadioReco/examples/RNO_data/calculate_cable_delays/fiber_delays/plots_fiber_delays.py rename to NuRadioReco/examples/RNOG/calculate_cable_delays/fiber_delays/plots_fiber_delays.py diff --git a/NuRadioReco/examples/RNO_data/calculate_cable_delays/gps_positions/gps_positions.py b/NuRadioReco/examples/RNOG/calculate_cable_delays/gps_positions/gps_positions.py similarity index 100% rename from NuRadioReco/examples/RNO_data/calculate_cable_delays/gps_positions/gps_positions.py rename to NuRadioReco/examples/RNOG/calculate_cable_delays/gps_positions/gps_positions.py diff --git a/NuRadioReco/examples/RNOG/data_analysis_example.py b/NuRadioReco/examples/RNOG/data_analysis_example.py new file mode 100644 index 0000000000..74cf12973c --- /dev/null +++ b/NuRadioReco/examples/RNOG/data_analysis_example.py @@ -0,0 +1,145 @@ +import argparse +import logging +import time +import os +import numpy as np +from matplotlib import pyplot as plt +import NuRadioReco.modules.channelResampler +import NuRadioReco.modules.RNO_G.dataProviderRNOG +import NuRadioReco.modules.io.eventWriter +import NuRadioReco.modules.channelResampler +import NuRadioReco.detector.RNO_G.rnog_detector +from NuRadioReco.utilities import units, logging as nulogging + +from NuRadioReco.examples.RNOG.processing import process_event + +logger = logging.getLogger("NuRadioReco.example.RNOG.rnog_standard_data_processing") +logger.setLevel(nulogging.LOGGING_STATUS) + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description='Run standard RNO-G data processing') + + parser.add_argument('--filenames', type=str, nargs="*", + help='Specify root data files if not specified in the config file') + parser.add_argument('--outputfile', type=str, nargs=1, default=None) + + args = parser.parse_args() + nulogging.set_general_log_level(logging.WARNING) + + if args.outputfile is None: + if len(args.filenames) > 1: + raise ValueError("Please specify an output file") + + path = args.filenames[0] + + if path.endswith(".root"): + args.outputfile = path.replace(".root", ".nur") + elif os.path.isdir(path): + args.outputfile = os.path.join(path, "output.nur") + else: + args.outputfile = args.outputfile[0] + + logger.status(f"writing output to {args.outputfile}") + + # Initialize detector class + det = NuRadioReco.detector.RNO_G.rnog_detector.Detector() + + # Initialize io modules + dataProviderRNOG = NuRadioReco.modules.RNO_G.dataProviderRNOG.dataProvideRNOG() + dataProviderRNOG.begin(files=args.filenames, det=det) + + eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter() + eventWriter.begin(filename=args.outputfile) + + # initialize additional modules + channelResampler = NuRadioReco.modules.channelResampler.channelResampler() + channelResampler.begin() + + # For time logging + t_total = 0 + + # Loop over all events (the reader module has options to select events - + # see class documentation or module arguements in config file) + for idx, evt in enumerate(dataProviderRNOG.run()): + + if (idx + 1) % 50 == 0: + logger.info(f'"Processing events: {idx + 1}\r') + + t0 = time.time() + # perform standard RNO-G data processing + process_event(evt, det) + + ######################## + ######################## + # + # add more usercode here + # + ######################## + ######################## + + # the following code is just an example of how to access station and channel information: + + station = evt.get_station() # this assumes that there is only one station per event. If multiple stations are present, use evt.get_stations() and iterate over them + + # the following code is just an example of how to access reconstructed parameters + # (the ones that were calculated in the processing steps by the channelSignalReconstructor) + from NuRadioReco.framework.parameters import channelParameters as chp + max_SNR = 0 + for channel in station.iter_channels(): + SNR = channel[chp.SNR]['peak_2_peak_amplitude'] # alternatively, channel.get_parameter(chp.SNR) + max_SNR = max(max_SNR, SNR) + signal_amplitude = channel[chp.maximum_amplitude_envelope] + signal_time = channel[chp.signal_time] + + print(f"Channel {channel.get_id()}: SNR={SNR:.1f}, signal amplitude={signal_amplitude/units.mV:.2f}mV, signal time={signal_time/units.ns:.2f}ns") + + + # the following code is just an example of how to access the channel waveforms and plot them + # we do it only for high-SNR events + if max_SNR > 5: + # iterate over some channels in station and plot them + fig, ax = plt.subplots(4, 2) + for channel_id in range(4): # iterate over the first 4 channels + channel = station.get_channel(channel_id) + ax[channel_id, 0].plot(channel.get_times()/units.ns, channel.get_trace()/units.mV) # plot the timetrace in nanoseconds vs. microvolts + ax[channel_id, 1].plot(channel.get_frequencies()/units.MHz, np.abs(channel.get_frequency_spectrum())/units.mV) # plot the frequency spectrum in MHz vs. microvolts + ax[channel_id, 0].set_xlabel('Time [ns]') + ax[channel_id, 0].set_ylabel('Voltage [mV]') + ax[channel_id, 1].set_xlabel('Frequency [MHz]') + ax[channel_id, 1].set_ylabel('Voltage [mV]') + plt.tight_layout() + fig.savefig(f'channel_traces_{idx}.png') + plt.close(fig) + + # before saving events to disk, it is advisable to downsample back to the two-times the maximum frequency available in the data, i.e., back to the Nquist frequency + # this will save disk space and make the data processing faster. The preprocessing applied a 600MHz low-pass filter, so we can downsample to 2GHz without losing information + channelResampler.run(evt, station, det, sampling_rate=2 * units.GHz) + + # it is advisable to only save the full waveform information for events that pass certain analysis cuts + # this will save disk space and make the data processing faster + # Here, we implement a simple SNR cut as an example + interesting_event = False + if(max_SNR > 5): + interesting_event = True #determined by some analysis cuts + # Write event - the RNO-G detector class is not stored within the nur files. + if interesting_event: + # save full waveform information + print("saving full waveform information") + eventWriter.run(evt, det=None, mode={'Channels':True, "ElectricFields":True}) + else: + # only save meta information but no traces to save disk space + print("saving only meta information") + eventWriter.run(evt, det=None, mode={'Channels':False, "ElectricFields":False}) + + logger.debug("Time for event: %f", time.time() - t0) + t_total += time.time() - t0 + + dataProviderRNOG.end() + eventWriter.end() + + logger.status( + f"Processed {idx + 1} events:" + f"\n\tTotal time: {t_total:.2f}s" + f"\n\tTime per event: {t_total / (idx + 1):.2f}s") diff --git a/NuRadioReco/examples/RNOG/make_glitch_summary_plot.py b/NuRadioReco/examples/RNOG/make_glitch_summary_plot.py new file mode 100644 index 0000000000..3a5ffe838e --- /dev/null +++ b/NuRadioReco/examples/RNOG/make_glitch_summary_plot.py @@ -0,0 +1,105 @@ +import argparse, numpy as np +from NuRadioReco.modules.io import eventReader +import NuRadioReco.framework.parameters as parameters + +def plot(plot_data, outpath, fs=13, zoomfact=1.1): + + channel_colors = { + 0: "tab:blue", + 1: "tab:orange", + 2: "tab:green", + 3: "tab:red" + } + + num_bins = 10 + + channels = [key for key in plot_data.keys() if isinstance(key, int)] + num_channels = len(channels) + + import matplotlib as mpl + mpl.use('Agg') + import matplotlib.pyplot as plt + from matplotlib.gridspec import GridSpec + from matplotlib.ticker import AutoMinorLocator + + fig = plt.figure(figsize = (6, 5), layout = "constrained") + gs = GridSpec(num_channels, 2, figure = fig, width_ratios = [0.8, 0.2]) + + for ind, channel in enumerate(channels): + color = channel_colors.get(channel, "black") + + evt_nums = plot_data["evt_num"] + channel_ts = plot_data[channel] + yscale = zoomfact * np.max(np.abs(channel_ts)) + + ax = fig.add_subplot(gs[ind, 0]) + ax_hist = fig.add_subplot(gs[ind, 1], sharey = ax) + ax.set_ylim(-yscale, yscale) + ax.set_xlim(min(evt_nums), max(evt_nums)) + ax_hist.set_xlim(0.0, 5 * len(channel_ts) / num_bins) + + ax.text(0.02, 0.95, f"CH {channel}", fontsize = fs, transform = ax.transAxes, color = color, ha = "left", va = "top") + ax.set_ylabel("Test stat.", fontsize = fs) + + ax.scatter(evt_nums, channel_ts, s = 2, color = color) + ax_hist.hist(channel_ts, histtype = "step", bins = num_bins, orientation = "horizontal", color = color) + + ax.tick_params(axis = "x", direction = "in", which = "both", bottom = True, top = True, labelbottom = ind == len(channels) - 1, + labelsize = fs) + ax.tick_params(axis = "y", direction = "in", which = "both", left = True, right = True, labelsize = fs) + + ax.axhline(0.0, color = "gray", ls = "dashed", lw = 1.0) + ax_hist.axhline(0.0, color = "gray", ls = "dashed", lw = 1.0) + + ax_hist.tick_params(axis = "x", direction = "in", which = "both", bottom = True, top = True, labelbottom = ind == len(channels) - 1, + labelsize = fs) + ax_hist.tick_params(axis = "y", direction = "in", which = "both", left = True, right = True, labelsize = fs, + labelleft = False) + + ax.fill_between(evt_nums, 0.0, yscale, facecolor = "tab:red", alpha = 0.25, zorder = 0) + ax.fill_between(evt_nums, 0.0, -yscale, facecolor = "tab:green", alpha = 0.25, zorder = 0) + + ax_hist.fill_between(evt_nums, 0.0, yscale, facecolor = "tab:red", alpha = 0.25, zorder = 0) + ax_hist.fill_between(evt_nums, 0.0, -yscale, facecolor = "tab:green", alpha = 0.25, zorder = 0) + + ax_hist.text(0.95, 0.95, "Glitching", ha = "right", va = "top", color = "tab:red", fontsize = 7, transform = ax_hist.transAxes) + ax_hist.text(0.95, 0.05, "No Glitching", ha = "right", va = "bottom", color = "tab:green", fontsize = 7, transform = ax_hist.transAxes) + + ax.set_xlabel("Event Number", fontsize = fs) + ax_hist.set_xlabel("Evts. / bin", fontsize = fs) + + fig.savefig(outpath) + plt.close() + +def make_glitch_summary_plot(nur_path, plot_path, channels=[0, 1, 2, 3]): + + reader = eventReader.eventReader() + reader.begin(nur_path) + + plot_data = {"evt_num": []} + + for evt in reader.run(): + station = evt.get_station() + plot_data["evt_num"].append(evt.get_id()) + + for channel in channels: + ch_data = station.get_channel(channel) + ch_data.add_parameter_type(parameters.channelParametersRNOG) + ch_number = ch_data.get_id() + + ts = ch_data.get_parameter(parameters.channelParametersRNOG.glitch_test_statistic) + + if ch_number not in plot_data: + plot_data[ch_number] = [] + plot_data[ch_number].append(ts) + + plot(plot_data, plot_path) + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + parser.add_argument("--nur", type = str, action = "store", dest = "nur_path") + parser.add_argument("--plot", type = str, action = "store", dest = "plot_path") + args = vars(parser.parse_args()) + + make_glitch_summary_plot(**args) diff --git a/NuRadioReco/examples/RNOG/processing.py b/NuRadioReco/examples/RNOG/processing.py new file mode 100644 index 0000000000..4236a895e7 --- /dev/null +++ b/NuRadioReco/examples/RNOG/processing.py @@ -0,0 +1,89 @@ +import logging +import NuRadioReco.modules.channelResampler +import NuRadioReco.modules.channelBandPassFilter +import NuRadioReco.modules.channelCWNotchFilter +import NuRadioReco.modules.channelSignalReconstructor + +import NuRadioReco.modules.RNO_G.dataProviderRNOG +import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator +import NuRadioReco.modules.io.eventWriter + +import NuRadioReco.detector.RNO_G.rnog_detector +from NuRadioReco.utilities import units + + +logger = logging.getLogger("NuRadioReco.example.RNOG.rnog_standard_data_processing") +logger.setLevel(logging.INFO) + +channelResampler = NuRadioReco.modules.channelResampler.channelResampler() +channelResampler.begin() + +channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter() +channelBandPassFilter.begin() + +channelCWNotchFilter = NuRadioReco.modules.channelCWNotchFilter.channelCWNotchFilter() +channelCWNotchFilter.begin() + +hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator() +hardwareResponseIncorporator.begin() + +channelSignalReconstructor = NuRadioReco.modules.channelSignalReconstructor.channelSignalReconstructor(log_level=logging.WARNING) +channelSignalReconstructor.begin() + + +def process_event(evt, det): + """ + Recommended preprocessing for RNO-G events + + Parameters + ---------- + evt : NuRadioReco.event.Event + Event to process + det : NuRadioReco.detector.detector.Detector + Detector object + """ + + # loop over all stations in the event. Typically we only have a single station per event. + for station in evt.get_stations(): + # The RNO-G detector changed over time (e.g. because certain hardware components were replaced). + # The time-dependent detector description has this information but needs to be updated with the + # current time of the event. + det.update(station.get_station_time()) + + # The first step is to upsample the data to a higher sampling rate. This will e.g. allow to + # determine the maximum amplitude and time of the signal more accurately. Studies showed that + # upsampling to 5 GHz is good enough for most task. + # In general, we need to find a good compromise between the data size (and thus processing time) + # and the required accuracy. + # Also remember to always downsample the data again before saving it to disk to avoid unnecessary + # large files. + channelResampler.run(evt, station, det, sampling_rate=5 * units.GHz) + + + # Our antennas are only sensitive in a certain frequency range. Hence, we should apply a bandpass filter + # in the range where the antennas are sensitive. This will reduce the noise in the data and make the + # signal more visible. The optimal frequency range and filter type depends on the antenna type and + # the expected signal. For the RNO-G antennas, a bandpass filter between 100 MHz and 600 MHz is a good + # general choice. + channelBandPassFilter.run( + evt, station, det, + passband=[0.1 * units.GHz, 0.6 * units.GHz], + filter_type='butter', order=10) + + # The signal chain amplifies and disperses the signal. This module will correct for the effects of the + # analog signal chain, i.e., everything between the antenna and the ADC. This will typically increase the + # signal-to-noise ratio and make the signal more visible. + hardwareResponseIncorporator.run(evt, station, det, sim_to_data=False, mode='phase_only') + + # The antennas often pick up continuous wave (CW) signals from noise various sources. These signals can be + # very strong and can make it difficult to see other signals. This module will remove the CW signals from + # the data by dynamically identifying and removing the contaminated frequency bins. + # An alternative module is the channelSineWaveFilter, which only removes the noise contribution from the CW + # but not the thermal noise of that frequency. However, this is more computationally expensive. + channelCWNotchFilter.run(evt, station, det) + + # The data is now preprocessed and ready for further analysis. The next steps depend on the analysis + # you want to perform. For example, you can now search for signals in the data, determine the arrival + # direction of the signal, or reconstruct the energy of the signal. + # The channelSignalReconstructor module is a good starting point for the signal reconstruction. + channelSignalReconstructor.run(evt, station, det) \ No newline at end of file diff --git a/NuRadioReco/examples/RNO_data/read_data_example/pulser_data_21.root b/NuRadioReco/examples/RNOG/pulser_data_21.root similarity index 100% rename from NuRadioReco/examples/RNO_data/read_data_example/pulser_data_21.root rename to NuRadioReco/examples/RNOG/pulser_data_21.root diff --git a/NuRadioReco/examples/RNO_data/read_data_example/read_rnog.py b/NuRadioReco/examples/RNOG/read_rnog.py similarity index 100% rename from NuRadioReco/examples/RNO_data/read_data_example/read_rnog.py rename to NuRadioReco/examples/RNOG/read_rnog.py diff --git a/NuRadioReco/examples/RNOG/rnog_config.yaml b/NuRadioReco/examples/RNOG/rnog_config.yaml new file mode 100644 index 0000000000..3c984d372f --- /dev/null +++ b/NuRadioReco/examples/RNOG/rnog_config.yaml @@ -0,0 +1,61 @@ +# Config file for the standard RNO-G data processing + +# Yaml instructions: +# - Use always NuRadioReco default units as specifed in `NuRadioReco.utilities.units` +# - If you want to use floats with scientific notation, use the following format: 1.0e-9 (not 1e-9!) +# - `null` is the equivalent of `None` in python + +readRNOGDataMattak: + kwargs: + filenames: null + # Use mattak to calibrate data on the fly. Requires calibration files to be available + read_calibrated_data: False + # Do not apply a baseline correction in the read module, do it afterwards + apply_baseline_correction: none + # Only used when read_calibrated_data==False, performs a linear voltage calibration with hardcoded values + convert_to_voltage: True + # Can be used instead of defining a selector (only for triggers) + select_triggers: null + # If true, and if the RunTable database is available select runs based on the following criteria + select_runs: False + # Only use runs of a certain run type + run_types: ["physics"] + # Only use runs with a maximum trigger rate of 1 Hz + max_trigger_rate: 2.0e-9 + # This might be necessary to read older files which have an invalid value stored as sampling rate + overwrite_sampling_rate: 3.2 + +detector: + kwargs: + # Specify a path to the detector description file. If null (None in python), + #the detector description will be queried from the hardware database. + detector_file: null + # If you want to query information only from a specific station, + # you can specify the station ID here. + select_stations: null + +channelResampler: + use: True + kwargs: + sampling_rate: 5 # GHz + +hardwareResponseIncorporator: + use: True + kwargs: + mode: 'phase_only' + +channelCWNotchFilter: + use: False + kwargs: + +channelBandPassFilter: + use: True + kwargs: + passband: [0.1, 0.6] # GHz + filter_type: 'butter' + order: 10 + +eventWriter: + kwargs: + mode: null + filename: rnog_output_test.nur diff --git a/NuRadioReco/examples/RNO_data/read_data_example/run_reconstruction.py b/NuRadioReco/examples/RNOG/run_reconstruction.py similarity index 71% rename from NuRadioReco/examples/RNO_data/read_data_example/run_reconstruction.py rename to NuRadioReco/examples/RNOG/run_reconstruction.py index f676bd9a93..3432312301 100644 --- a/NuRadioReco/examples/RNO_data/read_data_example/run_reconstruction.py +++ b/NuRadioReco/examples/RNOG/run_reconstruction.py @@ -1,16 +1,23 @@ -import NuRadioReco +import NuRadioReco.modules.io.RNO_G.readRNOGDataMattak +import NuRadioReco.modules.channelBandPassFilter +import NuRadioReco.modules.sphericalWaveFitter +import NuRadioReco.modules.channelAddCableDelay + +from NuRadioReco.detector import detector +from NuRadioReco.utilities import units + import matplotlib.pyplot as plt -from NuRadioReco.modules.io.RNO_G.readRNOGDataMattak import readRNOGData import pandas as pd import numpy as np -from NuRadioReco.utilities import units -from NuRadioReco.modules import channelBandPassFilter -from NuRadioReco.detector import detector import datetime -from NuRadioReco.modules import sphericalWaveFitter -from NuRadioReco.modules import channelAddCableDelay +import os -""" An example to show how to read RNO-G data and how to perform simple reconstructions. The data used is pulser data and a simple (brute force, not optimized) spherical wave reconstruction is performed to obtain the pulser position """"" + +""" +An example to show how to read RNO-G data and how to perform simple reconstructions. +The data used is pulser data and a simple (brute force, not optimized) spherical +wave reconstruction is performed to obtain the pulser position. +""" """ Initiazize modules needed""" @@ -23,7 +30,9 @@ """ Specify the detector. """ -det = detector.Detector(json_filename = "../../../detector/RNO_G/RNO_season_2021.json") +json_file_path = os.path.dirname(__file__) + "/../../detector/RNO_G/RNO_season_2021.json" + +det = detector.Detector(json_filename=json_file_path) det.update(datetime.datetime(2022, 10, 1)) """ Get positions for the pulsers from the detector file as starting positions for the fit """ @@ -49,15 +58,12 @@ sphericalWaveFitter.run(event, station, det, start_pulser_position = rel_pulser_position, n_index = 1.78, debug =True) if plots: - fig = plt.figure() - i = 1 - for channel in station.iter_channels(): + fig, axs = plt.subplots(3, 2) + for ax, channel in zip(axs.flatten(), station.iter_channels()): if channel.get_id() in use_channels: - ax = fig.add_subplot(3, 2, i) ax.plot(channel.get_trace()) ax.title.set_text("channel id {}".format(channel.get_id())) ax.grid() - i+= 1 fig.tight_layout() fig.savefig("trace_{}.pdf".format(i_event)) diff --git a/NuRadioReco/framework/event.py b/NuRadioReco/framework/event.py index bb601e0c3f..2757a45260 100644 --- a/NuRadioReco/framework/event.py +++ b/NuRadioReco/framework/event.py @@ -5,9 +5,12 @@ import NuRadioReco.framework.sim_emitter import NuRadioReco.framework.hybrid_information import NuRadioReco.framework.particle -import NuRadioReco.framework.parameters as parameters import NuRadioReco.framework.parameter_storage +from NuRadioReco.framework.parameters import ( + eventParameters as evp, channelParameters as chp, showerParameters as shp, + particleParameters as pap, generatorAttributes as gta) + from NuRadioReco.utilities import io_utilities, version import astropy.time @@ -24,7 +27,7 @@ class Event(NuRadioReco.framework.parameter_storage.ParameterStorage): def __init__(self, run_number, event_id): - super().__init__([parameters.eventParameters, parameters.generatorAttributes]) + super().__init__([evp, gta]) self.__run_number = run_number self._id = event_id @@ -311,9 +314,9 @@ def get_parent(self, particle_or_shower): returns the parent of a particle or a shower """ if isinstance(particle_or_shower, NuRadioReco.framework.base_shower.BaseShower): - par_id = particle_or_shower[parameters.showerParameters.parent_id] + par_id = particle_or_shower[shp.parent_id] elif isinstance(particle_or_shower, NuRadioReco.framework.particle.Particle): - par_id = particle_or_shower[parameters.particleParameters.parent_id] + par_id = particle_or_shower[pap.parent_id] else: raise ValueError("particle_or_shower needs to be an instance of NuRadioReco.framework.base_shower.BaseShower or NuRadioReco.framework.particle.Particle") if par_id is None: @@ -348,12 +351,12 @@ def get_interaction_products(self, parent_particle, showers=True, particles=True # iterate over sim_showers to look for parent id if showers is True: for shower in self.get_showers(): - if shower[parameters.showerParameters.parent_id] == parent_id: + if shower[shp.parent_id] == parent_id: yield shower # iterate over secondary particles to look for parent id if particles is True: for particle in self.get_particles(): - if particle[parameters.particleParameters.parent_id] == parent_id: + if particle[pap.parent_id] == parent_id: yield particle def add_shower(self, shower): @@ -562,14 +565,53 @@ def get_hybrid_information(self): """ return self.__hybrid_information + def has_glitch(self): + """ + Returns true if any channel in any station has a glitch + """ + for station in self.get_stations(): + for channel in station.iter_channels(): + if channel.has_parameter(chp.glitch) and channel.get_parameter(chp.glitch): + return True + + return False + + def avg_SNR_all_stations(self): + """ + Returns the average SNR in every station + """ + if (len(list(self.get_stations())) == 1): + for station in self.get_stations(): + return station.get_average_signal_to_noise_ratio() + else: + avg_snrs = {} + for station in self.get_stations(): + avg_snrs[station.get_id()] = station.get_average_signal_to_noise_ratio() + + return avg_snrs + + def avg_RPR_all_stations(self): + """ + Returns the average RPR in every station + """ + if (len(list(self.get_stations())) == 1): + for station in self.get_stations(): + return station.get_average_root_power_ratio() + else: + avg_rprs = {} + for station in self.get_stations(): + avg_rprs[station.get_id()] = station.get_average_root_power_ratio() + + return avg_rprs + def serialize(self, mode): stations_pkl = [] try: commit_hash = version.get_NuRadioMC_commit_hash() - self.set_parameter(parameters.eventParameters.hash_NuRadioMC, commit_hash) + self.set_parameter(evp.hash_NuRadioMC, commit_hash) except: logger.warning("Event is serialized without commit hash!") - self.set_parameter(parameters.eventParameters.hash_NuRadioMC, None) + self.set_parameter(evp.hash_NuRadioMC, None) for station in self.get_stations(): stations_pkl.append(station.serialize(mode)) diff --git a/NuRadioReco/framework/parameters.py b/NuRadioReco/framework/parameters.py index ee5f5d781c..fa201df8fe 100644 --- a/NuRadioReco/framework/parameters.py +++ b/NuRadioReco/framework/parameters.py @@ -63,6 +63,26 @@ class channelParameters(Enum): Vrms_NuRadioMC_simulation = 19 #: the noise rms used in the MC simulation bandwidth_NuRadioMC_simulation = 20 #: the integrated channel response (=bandwidth for signal chains without amplification) used in the MC simulation Vrms_trigger_NuRadioMC_simulation = 21 #: the noise rms of the trigger channels (optional) used in the MC simulation + root_power_ratio = 22 #: the root power ratio (float) + +class channelParametersRNOG(Enum): + # RNO-G specific channel parameters + + # FS: I did not start with a negative parameter on the 1, hence I chose 100 + glitch = 100 #: True if channel is likely to have a glitch. See 'NuRadioReco.modules.RNO_G.channelGlitchDetector' + glitch_test_statistic = 101 #: Numerical value delivered by the glitch detector. Positive values indicate a likely glitch. + +class eventParametersRNOG(Enum): + max_corr_coords = 4 #: azimuth, zenith corresponding to max correlation + max_corr = 5 #: maximum correlation value obtained through reconstruction + csw_snr = 6 #: signal to noise ratio of coherently summed waveform + csw_rpr = 7 #: root power ratio of coherently summed waveform + csw_hilbert_snr = 8 #: signal to noise ratio of hilberted coherently summed waveform + csw_impulsivity = 9 #: coherently summed waveform's impulsivity and other statistical measures (dict) + surf_corr_ratio = 10 #: ratio comparing max surface correlation to max overall correlation + max_surf_corr = 11 #: maximum surface correlation + avg_snr = 12 #: average signal to noise ratio across all channels ("peak to peak") + avg_rpr = 13 #: average root power ratio across all channels class electricFieldParameters(Enum): ray_path_type = 1 #: the type of the ray tracing solution ('direct', 'refracted' or 'reflected') diff --git a/NuRadioReco/framework/station.py b/NuRadioReco/framework/station.py index 262c1043ca..6fd5c2dfb4 100644 --- a/NuRadioReco/framework/station.py +++ b/NuRadioReco/framework/station.py @@ -261,6 +261,40 @@ def get_magnetic_field_vector(self, time=None): logger.error(f"Reference reconstruction not set / unknown: {self.__reference_reconstruction}") raise ValueError(f"Reference reconstruction not set / unknown: {self.__reference_reconstruction}") + def get_average_signal_to_noise_ratio(self, channels_to_include = [0, 1, 2, 3, 5, 6, 7, 22, 23]): + """ + Returns the average SNR across all channels + """ + snr_all = [] + for channel in self.iter_channels(channels_to_include): + if (channel.has_parameter(chp.SNR)): + snr_ch = channel.get_parameter(chp.SNR) + if (snr_ch == np.inf): + return np.inf + else: + snr_all.append(snr_ch['peak_amplitude']) + else: + raise LookupError(f"Channel {channel.get_id()} has no parameter SNR") + + return np.mean(snr_all) + + def get_average_root_power_ratio(self, channels_to_include = [0, 1, 2, 3, 5, 6, 7, 22, 23]): + """ + Returns the average RPR across all channels + """ + rpr_all = [] + for channel in self.iter_channels(channels_to_include): + if (channel.has_parameter(chp.root_power_ratio)): + rpr_ch = channel.get_parameter(chp.root_power_ratio) + if (rpr_ch == np.inf): + return np.inf + else: + rpr_all.append(rpr_ch) + else: + raise LookupError(f"Channel {channel.get_id()} has no parameter RPR") + + return np.mean(rpr_all) + def serialize(self, mode): save_efield_traces = 'ElectricFields' in mode and mode['ElectricFields'] is True base_station_pkl = NuRadioReco.framework.base_station.BaseStation.serialize( diff --git a/NuRadioReco/modules/RNO_G/channelGlitchDetector.py b/NuRadioReco/modules/RNO_G/channelGlitchDetector.py new file mode 100644 index 0000000000..1e3e7cdc80 --- /dev/null +++ b/NuRadioReco/modules/RNO_G/channelGlitchDetector.py @@ -0,0 +1,142 @@ +from NuRadioReco.modules.base.module import register_run +from NuRadioReco.utilities import units +from NuRadioReco.framework.parameters import channelParametersRNOG as chp + +import numpy as np +import logging +import collections + +class channelGlitchDetector: + """ + This module detects scrambled data (digitizer "glitches") in the channels. + + RNO-G data (in particular the RNO-G data recorded with the LAB4D digitizers on the + first-generation radiants) is known to have "glitches" in the channels. + When a glitch is present in a channel, the 64-sample readout blocks were scrambled + by the readout electronics which results in sudden unphyiscal jumps in the signal. + These jumps are detected with this module (but not corrected). + """ + + def __init__(self, cut_value=0.0, glitch_fraction_warn_level=0.1, log_level=logging.NOTSET): + """ + Parameters + ---------- + cut_value : float + This module calculates a test statistic that is sensitive to the presence + of digitizer glitches. This parameter marks the critical value at which the test is + said to have detected a glitch. This is a free parameter; increasing its value results + in a lower false-positive rate (i.e. healthy events are incorrectly marked as glitching), + but also a lower true-positive rate (i.e. glitching events are correctly marked as such). + The default value of 0.0 is a good starting point. + + glitch_fraction_warn_level : float + Print warning messages at the end of a run if a channel shows glitching in more than a fraction of + `glitch_fraction_warn_level` of all events. + + log_level: enum + Set verbosity level of logger. If logging.DEBUG, set mattak to verbose (unless specified in mattak_kwargs). + (Default: logging.NOTSET, ie adhere to general log level) + + """ + + self.logger = logging.getLogger('NuRadioReco.RNO_G.channelGlitchDetector') + + self.ts_cut_value = cut_value + self.glitch_fraction_warn_level = glitch_fraction_warn_level + + # Total sampling buffer of the LAB4D: 2048 samples + self._lab4d_readout_size = 2048 + + # Length of a sampling block in the LAB4D: 64 samples + self._lab4d_sampling_blocksize = 64 + + def begin(self): + # Per-run glitching statistics + self.events_checked = 0 + self.events_glitching_per_channel = collections.defaultdict(int) + + def end(self): + # Print glitch statistic summary + for ch_id, events_glitching in self.events_glitching_per_channel.items(): + glitching_fraction = events_glitching / self.events_checked + if glitching_fraction > self.glitch_fraction_warn_level: + self.logger.warning(f"Channel {ch_id} shows glitches in {events_glitching} / {self.events_checked} = {100*glitching_fraction:.2f}% of events!") + else: + self.logger.info(f"Channel {ch_id} shows glitches in {events_glitching} / {self.events_checked} = {100*glitching_fraction:.2f}% of events!") + + @register_run() + def run(self, event, station, det=None): + """ Run over channel traces and sets `channelParameter.glitch`. + + Parameters + ---------- + event : `NuRadioReco.framework.event.Event` + The event object + station : `NuRadioReco.framework.station.Station` + The station object + det : `NuRadioReco.detector.detector.Detector` (default: None) + Detector object, not used! + """ + + def diff_sq(eventdata): + """ + Returns sum of squared differences of samples across seams of 128-sample chunks. + + `eventdata`: channel waveform + """ + + block_size = self._lab4d_sampling_blocksize + twice_block_size = 2 * block_size + + runsum = 0.0 + for chunk in range(len(eventdata) // twice_block_size - 1): + runsum += (eventdata[chunk * twice_block_size + block_size - 1] - eventdata[chunk * twice_block_size + block_size]) ** 2 + return np.sum(runsum) + + def unscramble(trace): + """ + Applies an unscrambling operation to the passed `trace`. + Note: the first and last sampling block are unusable and hence replaced by zeros in the returned waveform. + + Parameters + ---------- + `trace`: channel waveform + """ + + readout_size = self._lab4d_readout_size + block_size = self._lab4d_sampling_blocksize + twice_block_size = 2 * block_size + + new_trace = np.zeros_like(trace) + + for i_section in range(len(trace) // block_size): + section_start = i_section * block_size + section_end = i_section * block_size + block_size + if i_section % 2 == 0: + new_trace[(section_start + twice_block_size) % readout_size :\ + (section_end + twice_block_size) % readout_size] = trace[section_start:section_end] + elif i_section > 1: + new_trace[(section_start - twice_block_size) % readout_size :\ + (section_end - twice_block_size) % readout_size] = trace[section_start:section_end] + new_trace[0:block_size] = 0 + + return new_trace + + # update event counter + self.events_checked += 1 + + for ch in station.iter_channels(): + ch_id = ch.get_id() + + trace = ch.get_trace() + trace_us = unscramble(trace) + + # glitching test statistic and boolean discriminate + glitch_ts = (diff_sq(trace) - diff_sq(trace_us)) / np.var(trace) + glitch_disc = glitch_ts > self.ts_cut_value + + ch.set_parameter(chp.glitch, glitch_disc) + ch.set_parameter(chp.glitch_test_statistic, glitch_ts) + + # update glitching statistics + self.events_glitching_per_channel[ch_id] += glitch_disc diff --git a/NuRadioReco/modules/RNO_G/dataProviderRNOG.py b/NuRadioReco/modules/RNO_G/dataProviderRNOG.py new file mode 100644 index 0000000000..e43d04b3ba --- /dev/null +++ b/NuRadioReco/modules/RNO_G/dataProviderRNOG.py @@ -0,0 +1,50 @@ +import NuRadioReco.modules.io.RNO_G.readRNOGDataMattak +import NuRadioReco.modules.RNO_G.channelGlitchDetector +import NuRadioReco.modules.RNO_G.channelBlockOffsetFitter + +import NuRadioReco.modules.channelAddCableDelay + +import logging +logger = logging.getLogger('NuRadioReco.RNO_G.dataProviderRNOG') + +class dataProvideRNOG: + + def __init__(self): + + self.channelGlitchDetector = NuRadioReco.modules.RNO_G.channelGlitchDetector.channelGlitchDetector() + self.channelBlockOffsetFitter = NuRadioReco.modules.RNO_G.channelBlockOffsetFitter.channelBlockOffsets() + self.reader = NuRadioReco.modules.io.RNO_G.readRNOGDataMattak.readRNOGData() + + self.channelCableDelayAdder = NuRadioReco.modules.channelAddCableDelay.channelAddCableDelay() + + + def begin(self, files, reader_kwargs={}, det=None): + self.files = files + + self.channelGlitchDetector.begin() + self.channelBlockOffsetFitter.begin() + self.reader.begin(self.files, **reader_kwargs) + self.channelCableDelayAdder.begin() + + assert det is not None, "Detector object is None, please provide a detector object." + self.detector = det + + def end(self): + self.reader.end() + self.channelGlitchDetector.end() + + def run(self): + + for event in self.reader.run(): + + # This will throw an error if the event has more than one station + station = event.get_station() + self.detector.update(station.get_station_time()) + + self.channelGlitchDetector.run(event, station, self.detector) + + self.channelBlockOffsetFitter.run(event, station, self.detector) + + self.channelCableDelayAdder.run(event, station, self.detector, mode='subtract') + + yield event diff --git a/NuRadioReco/modules/channelBeamFormingDirectionFitter.py b/NuRadioReco/modules/channelBeamFormingDirectionFitter.py new file mode 100644 index 0000000000..e3ec69bc17 --- /dev/null +++ b/NuRadioReco/modules/channelBeamFormingDirectionFitter.py @@ -0,0 +1,607 @@ +import matplotlib.pyplot as plt +from NuRadioReco.modules.base.module import register_run +from NuRadioReco.utilities import units +import NuRadioReco.detector.detector +import uproot +import numpy as np +import itertools +import argparse +from plotting import plotter +import datetime +import sys, os +from scipy import signal +import yaml +from scipy.interpolate import RegularGridInterpolator +from NuRadioReco.modules.RNO_G.channelBlockOffsetFitter import fit_block_offsets +import pandas as pd +import copy +import json +from NuRadioReco.modules.base import module +from NuRadioReco.framework.parameters import stationParameters as stnp +from NuRadioReco.framework.parameters import channelParameters as chp + +class channelBeamFormingDirectionFitter(): + def __init__(self): + self.begin() + + def begin(self): + # here you could define some variables that are the same for all events + pass + + @register_run() + def run(self, evt, station, det, config): + + if isinstance(config, str): + config = self.load_config(config) + + map_info = MapInfo( + config["limits"], + config["stepsize"], + config["channels"], + config["calibrate"], + config["events"], + ) + + # I removed the StationInfo class as all this information is already available in the NuRadio station object + + + channel_ids = map_info.used_channels + positions = Positions( + station.get_id(), + map_info, + config["coords"], + config["distance"], + config["rec_type"], + det=det + ) + + ant_locs = positions.get_ant_locs()[1] + + delay_matrices = positions.get_t_delay_matrices( + ant_locs, config["rec_type"] + ) + + corr_matrix_sum = None + max_corr_vals = [] + rec_coord0s = [] + rec_coord1s = [] + + + # upsample = 1 # upsampling is best done by the channelResampler module that is executed before this module + + # the current sampling rate can be requested from the station object + sampling_rate = station.get_channel(0).get_sampling_rate() #assuming that there is a channel 0 + + # the padding is being done in the channelStopFilter module, therefor we don't need to pad the trace here + # pad_amount = int( + # max(delays.values()) + # * sampling_rate + # ) + + times = station.get_channel(0).get_times() + + # this doesn't seem to do anyting. The `times` array from above is already the correct time array + # new_times = np.linspace( + # times[0], times[-1], len(times) + # ) + # dt = times[1] - times[0] + + # as the traces are already padded in a previous module, we don't need to pad them here: + # shifted time array only calculated once since time axis is the same + # for all waveforms in a run + # shifted_time_array = np.array( + # [ + # float(n) * dt + # for n in range(1, 2 * pad_amount + len(times) + 1) + # ] + # ) + + shifted_volt_array = [] + for ch in map_info.used_channels: + # waveform = station.get_channel(ch).get_trace() + # it is better to use the channelResampler module to resample the waveform + # and the channelStopFilter module to pad the trace with zeros. + # then, the waveform can be shifted in an optimized way (using the Fourier shift theorem) by + channel_copy = copy.copy(station.get_channel(ch)) + channel_copy.apply_time_shift(det.get_cable_delay(station.get_id(), channel_copy.get_id())) + shifted_waveform = channel_copy.get_trace() + # waveform = signal.resample_poly(waveform, upsample, 1) + # padded_wf = np.pad(waveform, pad_amount) + # element_shift = int( + # self.station_info.delays[ch] * self.sampling_rate * upsample + # ) + # shifted_waveform = np.roll(padded_wf, -element_shift) + + scaled_waveform = shifted_waveform / np.max(shifted_waveform) + scaled_waveform = scaled_waveform / np.std(scaled_waveform) + scaled_waveform = scaled_waveform - np.mean(scaled_waveform) + + shifted_volt_array.append(scaled_waveform) + + + v_array_pairs = list(itertools.combinations(shifted_volt_array, 2)) + + corr_matrix, max_corr = self.correlator( + times, v_array_pairs, delay_matrices + ) + + + rec_coord0, rec_coord1 = positions.get_rec_locs_from_corr_map( + corr_matrix + ) + + rec_coord0s.append(rec_coord0) + rec_coord1s.append(rec_coord1) + max_corr_vals.append(max_corr) + surface_corr = [] + + if config["coords"] == "cylindrical": + num_rows_to_10m = int(config["limits"][-1] / 10) + 1 + surface_corr.append( + self.get_surf_corr_ratio(corr_matrix, num_rows_to_10m) + ) + + # to safe output parameters, add new parameters to station or channel level in `NuRadioReco/framework/parameters.py` + # then set the parameter like this: + # e.g. station.set_parameter(stnp.beamforming_direction, rec_coord0) + # the eventWriter module will automatically write them to disk. + # you can also use the `channelParameters` class to set parameters on the channel level. + # if you want to access the parameters in a later module, you can use the `get_parameter` method of the station or channel object. + # e.g. rec_coord0 = station.get_parameter(stnp.beamforming_direction) + + + def end(self): + pass + + def load_config(self, config_file): + with open(config_file, "r") as f: + return yaml.safe_load(f) + + def get_surf_corr_ratio(self, corr_map, num_rows_for_10m): + surf_corr_ratio = np.max(corr_map[:num_rows_for_10m]) / np.max(corr_map) + return surf_corr_ratio + + def load_interpolator(self, table_filename): + time_differences_fine_3d = np.load(table_filename) + + x_range = np.arange(-200, 201, 1) + y_range = np.arange(-200, 201, 1) + z_range = np.arange(-100, 1, 1) + + interpolator = RegularGridInterpolator( + (x_range, y_range, z_range), time_differences_fine_3d, method="linear" + ) + + return interpolator + + def corr_index_from_t(self, time_delay, times): + center = len(times) + dt = times[1] - times[0] + + time_delay_clean = np.nan_to_num(time_delay, nan=0.0) + + index = np.rint(time_delay_clean / dt + (center - 1)).astype(int) + + return index + + def get_max_val_indices(self, matrix): + max_value = matrix.argmax() + row_by_col = matrix.shape + best_row_index, best_col_index = np.unravel_index(max_value, row_by_col) + return best_col_index, best_row_index + + def correlator(self, times, v_array_pairs, delay_matrices): + amplitude_corrector = 1.0 / (float(len(v_array_pairs[0][0]))) + + volt_corrs = [ + amplitude_corrector + * np.asarray(signal.correlate(v_array_pair[0], v_array_pair[1])) + for v_array_pair in v_array_pairs + ] + + pair_corr_matrices = [] + + for volt_corr, time_delay in zip(volt_corrs, delay_matrices): + indices = self.corr_index_from_t(time_delay, times) + mask = np.isnan(time_delay) + corr_matrix = np.zeros_like(time_delay) + + for i in range(indices.shape[0]): + for j in range(indices.shape[1]): + if not mask[i, j]: + corr_matrix[i, j] = volt_corr[indices[i, j]] + + pair_corr_matrices.append(corr_matrix) + + mean_corr_matrix = np.mean(np.array(pair_corr_matrices), axis=0) + max_corr = np.max(np.array(mean_corr_matrix)) + + return mean_corr_matrix, max_corr + + +def get_surf_corr_ratio(corr_map, num_rows_for_10m): + # print(f"slice of corr: {corr_map[:num_rows_for_10m]}") + surf_corr_ratio = np.max(corr_map[:num_rows_for_10m]) / np.max(corr_map) + + return surf_corr_ratio + + +def load_interpolator(table_filename): + time_differences_fine_3d = np.load(table_filename) + + # Define the grid used for generating lookup tables + x_range = np.arange(-200, 201, 1) + y_range = np.arange(-200, 201, 1) + z_range = np.arange(-100, 1, 1) + + # Create the interpolator + interpolator = RegularGridInterpolator( + (x_range, y_range, z_range), time_differences_fine_3d, method="linear" + ) + + return interpolator + + +class MapInfo: + + def __init__( + self, + limits, + stepsize, + used_channels, + cal_locs_file, + selected_events, + ): + self.limits = limits + self.stepsize = stepsize + self.used_channels = used_channels + self.cal_locs_file = cal_locs_file + self.selected_events = selected_events + + +class StationInfo: + + def __init__(self, station, det): + self.det = det + self.station = station + det.update(datetime.datetime(2023, 10, 1)) + self.channels = range(0, 24) + self.delays = dict( + zip( + self.channels, + [ + self.det.get_cable_delay(int(self.station), int(channel)) + for channel in self.channels + ], + ) + ) + + +class Positions: + + def __init__(self, station_id, map_info, coord_system, dist, rec_type, det): + self.rec_type = rec_type + self.phis = None + self.map_info = map_info + self.station_id = station_id + self.coord_system = coord_system + self.cal_locs_file = self.map_info.cal_locs_file + self.pulser_id = None + self.distance = dist + self.det = det + + def get_t_delay_matrices(self, ant_locs, rec_type): + + ch_pairs = list(itertools.combinations(self.map_info.used_channels, 2)) + print(f"ch pairs used: {ch_pairs}\n") + posn_pairs = list(itertools.combinations(ant_locs, 2)) + + z_offset = -self.get_xyz_origin()[2] + + src_enu_matrix_and_coord_vecs = self.get_source_enu_matrix(rec_type) + src_posn_enu_matrix = src_enu_matrix_and_coord_vecs[0] + + # for src_posn_enu_matrix in src_posn_enu_matrices: + time_delay_matrices = [] + for ch_pair, posn_pair in zip(ch_pairs, np.array(posn_pairs)): + # my_interp = load_interpolator( + # f"/storage/group/szw5718/default/rno-g/data/interp_tables/station{self.station_info.station}/station{self.station_info.station}_2023json_nocal_time_differences_3d_{ch_pair[0]}_{ch_pair[1]}_200_1deg_grid_noraytracingasdfasdfdstest.npy" + # ) + # my_interp = load_interpolator( + # f"/storage/group/szw5718/default/rno-g/data/interp_tables/station{self.station_info.station}/station{self.station_info.station}_2023json_nocal_time_differences_3d_{ch_pair[0]}_{ch_pair[1]}_200_1deg_grid_noraytracingasdfasdfdstest.npy" + # ) + if self.station_id == 23: + my_interp = load_interpolator( + f"/storage/group/szw5718/default/rno-g/data/interp_tables/station{self.station_info.station}/station{self.station_info.station}_2023json_nocal_time_differences_3d_{ch_pair[0]}_{ch_pair[1]}_200_1deg_grid_noraytracing.npy" + ) + elif self.station_id == 13: + my_interp = load_interpolator( + f"/storage/group/szw5718/default/rno-g/data/interp_tables/station{self.station_info.station}/station{self.station_info.station}_2023json_nocal_time_differences_3d_{ch_pair[0]}_{ch_pair[1]}_200_1deg_grid_noraytracingasdfasdfdstest.npy" + ) + else: + print("Station not found, returning None") + return None + + time_delay_matrix = my_interp(src_posn_enu_matrix) + time_delay_matrices.append(time_delay_matrix) + + return time_delay_matrices + + def get_xyz_origin(self): + + if self.cal_locs_file: + xyz_ch_1_loc = np.array(self.get_ant_locs()[0]["1"]) + xyz_ch_2_loc = np.array(self.get_ant_locs()[0]["2"]) + else: + xyz_ch_1_loc = self.det.get_relative_position(self.station_id, 1) + xyz_ch_2_loc = self.det.get_relative_position(self.station_id, 2) + + xyz_origin_loc = (xyz_ch_1_loc + xyz_ch_2_loc) / 2.0 + + return xyz_origin_loc + + def get_ant_locs(self): + + channel_ids = [0, 1, 2, 3, 5, 6, 7, 9, 10, 22, 23] + if self.cal_locs_file: + data = np.load(self.cal_locs_file) + + all_xyz_ant_locs = {} + for channel_idx, channel in enumerate(data): + all_xyz_ant_locs[str(channel_ids[channel_idx])] = [ + channel[0][231], + channel[1][231], + channel[2][231], + ] + else: + all_xyz_ant_locs = {} + for ch in channel_ids: + all_xyz_ant_locs[str(ch)] = ( + self.det.get_relative_position(self.station_id, int(ch) + ).tolist() + ) + + used_xyz_ant_locs = [ + all_xyz_ant_locs[str(ch)] for ch in self.map_info.used_channels + ] + + return all_xyz_ant_locs, used_xyz_ant_locs + + def generate_coord_arrays(self): + """ + From specified limits and step size, create arrays for each pulser + coordinate in its proper units. + """ + left, right, bottom, top = self.map_info.limits + + coord0_vec = np.arange( + left, right + self.map_info.stepsize, self.map_info.stepsize + ) + coord1_vec = np.arange( + bottom, top - self.map_info.stepsize, -self.map_info.stepsize + ) + + # else: + # coord1_vec = np.arange(bottom, top + self.map_info.stepsize, self.map_info.stepsize) + + # print(f"coord0 vec: {coord0_vec}") + # print(f"coord1 vec: {coord1_vec}") + + if coord0_vec[-1] < right: + coord0_vec = np.append(coord0_vec, right) + if coord1_vec[-1] < top: + coord1_vec = np.append(coord1_vec, top) + + if self.rec_type == "surface": + coord0_vec = [coord0 * units.m for coord0 in coord0_vec] + coord1_vec = [coord1 * units.deg for coord1 in coord1_vec] + else: + if self.coord_system == "cylindrical": + if self.rec_type == "phiz": + coord0_vec = [coord0 * units.deg for coord0 in coord0_vec] + coord1_vec = [coord1 * units.m for coord1 in coord1_vec] + elif self.rec_type == "rhoz": + coord0_vec = [coord0 * units.m for coord0 in coord0_vec] + coord1_vec = [coord1 * units.m for coord1 in coord1_vec] + elif self.coord_system == "spherical": + coord0_vec = [coord0 * units.deg for coord0 in coord0_vec] + coord1_vec = [coord1 * units.deg for coord1 in coord1_vec] + + # print(f"coord1_vec: {[angle/units.deg for angle in coord1_vec]}") + + else: + sys.exit( + "Coordinate system not found. Please specify either " + "'cylindrical' or 'spherical' coordinates." + ) + + return coord0_vec, coord1_vec + + def get_surface_coords(self, coord0_vec, coord1_vec): + + surface_coord_dict = defaultdict(list) + x = [] + y = [] + z = [-self.get_xyz_origin()[2]] * (len(coord0_vec) * len(coord1_vec)) + for R in coord0_vec: + for phi in coord1_vec: + x.append(R * np.cos(phi / units.rad)) + y.append(R * np.sin(phi / units.rad)) + + surface_coord_dict[(R, phi)] = [x, y, z] + + return surface_coord_dict + + def get_source_enu_matrix(self, rec_type): + """ + Returns a matrix of potential source locations in enu coords with + respect to the station origin. It must be in this + coord system for the time difference calculator to work properly. + """ + coord0_vec, coord1_vec = self.generate_coord_arrays() + + if self.coord_system == "cylindrical": + if rec_type == "phiz": + rho_to_pulser = self.distance + + phi_grid, z_grid = np.meshgrid(coord0_vec, coord1_vec) + rho_grid = np.full_like(phi_grid, rho_to_pulser) + + x_grid, y_grid, z_grid = self.get_enu_wrt_new_origin( + [rho_grid, phi_grid, z_grid] + ) + + elif rec_type == "rhoz": + phi = 0 * units.deg + + rho_grid, z_grid = np.meshgrid(coord0_vec, coord1_vec) + + # rho_grid = rho_grid[:len(coord1_vec)-1, :len(coord0_vec)-1] + # z_grid = z_grid[:len(coord1_vec)-1, :len(coord0_vec)-1] + + phi_grid = np.full_like(rho_grid, phi) + + x_grid, y_grid, z_grid = self.get_enu_wrt_new_origin( + [rho_grid, phi_grid, z_grid] + ) + + else: + rval = self.distance + phi_grid, theta_grid = np.meshgrid(coord0_vec, coord1_vec) + r_grid = np.full_like(phi_grid, rval) + + x_grid, y_grid, z_grid = self.get_enu_wrt_new_origin( + [r_grid, phi_grid, theta_grid] + ) + + src_xyz_loc_matrix = np.stack((x_grid, y_grid, z_grid), axis=-1) + + x_values = src_xyz_loc_matrix[:, :, 0] + y_values = src_xyz_loc_matrix[:, :, 1] + z_values = src_xyz_loc_matrix[:, :, 2] + + # Compute min and max for x, y, z + min_x, max_x = np.min(x_values), np.max(x_values) + min_y, max_y = np.min(y_values), np.max(y_values) + min_z, max_z = np.min(z_values), np.max(z_values) + + # Print results rounded to 2 decimal places + print(f"Min X: {min_x:.2f}, Max X: {max_x:.2f}") + print(f"Min Y: {min_y:.2f}, Max Y: {max_y:.2f}") + print(f"Min Z: {min_z:.2f}, Max Z: {max_z:.2f}") + + # print("\n") + # print(src_xyz_loc_matrix) + # print(np.shape(src_xyz_loc_matrix)) + # print("\n") + + # print(src_xyz_loc_matrix[0][0]) + + # print(f"\nsrc loc matrix: {src_xyz_loc_matrix}\n") + + # print(f"\nsize of src matrix: {np.shape(src_xyz_loc_matrix)}") + + return src_xyz_loc_matrix, [coord0_vec, coord1_vec] + + def get_enu_wrt_new_origin(self, coords): + """ + Converts enu to enu of xyz_origin coord system + """ + enu_origin = self.get_xyz_origin() + + if self.coord_system == "cylindrical": + rhos = coords[0] + phis = coords[1] + zs = coords[2] + + eastings = rhos * np.cos(phis) + northings = rhos * np.sin(phis) + elevations = zs + + else: + rs = coords[0] + phis = coords[1] + thetas = coords[2] + + eastings = rs * np.sin(thetas) * np.cos(phis) + northings = rs * np.sin(thetas) * np.sin(phis) + elevations = rs * np.cos(thetas) + + # prob should re-enable this? + + if self.coord_system == "spherical": + # prob should re-enable for the rest too? + eastings = eastings + enu_origin[0] + northings = northings + enu_origin[1] + elevations = elevations + enu_origin[2] + + # print(f"\n\nenu origin: {enu_origin}\n\n") + # print(f"\n\nenu origin z: {enu_origin[2]}\n\n") + + return eastings, northings, elevations + + def get_rec_locs_from_corr_map(self, corr_matrix): + + coord0_vec, coord1_vec = self.generate_coord_arrays() + + rec_pulser_loc0_idx, rec_pulser_loc1_idx = get_max_val_indices( + corr_matrix + ) + coord0_best = coord0_vec[rec_pulser_loc0_idx] + coord1_best = coord1_vec[rec_pulser_loc1_idx] + + return coord0_best, coord1_best + + +def corr_index_from_t(time_delay, times): + """ + Finds the index in the correlation corresponding to the given time delay. + """ + center = len(times) + dt = times[1] - times[0] + + time_delay_clean = np.nan_to_num(time_delay, nan=0.0) + + index = np.rint(time_delay_clean / dt + (center - 1)).astype(int) + + return index + + +def get_max_val_indices(matrix): + + max_value = matrix.argmax() + row_by_col = matrix.shape + best_row_index, best_col_index = np.unravel_index(max_value, row_by_col) + + return best_col_index, best_row_index + + +def correlator(times, v_array_pairs, delay_matrices): + + amplitude_corrector = 1.0 / (float(len(v_array_pairs[0][0]))) + + volt_corrs = [ + amplitude_corrector + * np.asarray(signal.correlate(v_array_pair[0], v_array_pair[1])) + for v_array_pair in v_array_pairs + ] + + pair_corr_matrices = [] + + for volt_corr, time_delay in zip(volt_corrs, delay_matrices): + indices = corr_index_from_t(time_delay, times) + mask = np.isnan(time_delay) + corr_matrix = np.zeros_like(time_delay) + + for i in range(indices.shape[0]): + for j in range(indices.shape[1]): + if not mask[i, j]: + corr_matrix[i, j] = volt_corr[indices[i, j]] + + pair_corr_matrices.append(corr_matrix) + + mean_corr_matrix = np.mean(np.array(pair_corr_matrices), axis=0) + max_corr = np.max(np.array(mean_corr_matrix)) + + return mean_corr_matrix, max_corr + diff --git a/NuRadioReco/modules/channelSignalReconstructor.py b/NuRadioReco/modules/channelSignalReconstructor.py index a5fbc3f1d1..c6890fa8af 100644 --- a/NuRadioReco/modules/channelSignalReconstructor.py +++ b/NuRadioReco/modules/channelSignalReconstructor.py @@ -1,6 +1,7 @@ from NuRadioReco.modules.base.module import register_run import numpy as np from scipy import signal +from scipy.ndimage import uniform_filter1d import time from NuRadioReco.utilities import units @@ -62,7 +63,7 @@ def begin( self.__noise_window_length = noise_window_length self.__debug = debug - def get_SNR(self, station_id, channel, det, stored_noise=False, rms_stage=None): + def get_SNR_and_RPR(self, station_id, channel, det, stored_noise=False, rms_stage=None): """ Parameters ---------- @@ -79,6 +80,8 @@ def get_SNR(self, station_id, channel, det, stored_noise=False, rms_stage=None): ------- SNR: dict dictionary of various SNR parameters + RPR: float + root power ratio of a channel """ trace = channel.get_trace() @@ -119,32 +122,53 @@ def get_SNR(self, station_id, channel, det, stored_noise=False, rms_stage=None): plt.ylabel("Power") plt.legend() - # Calculating SNR - SNR = {} - if (noise_rms == 0) or (noise_int == 0): - logger.info("RMS of noise is zero, calculating an SNR is not useful. All SNRs are set to infinity.") - SNR['peak_2_peak_amplitude'] = np.inf - SNR['peak_amplitude'] = np.inf - SNR['integrated_power'] = np.inf + # Calculating SNR - signal to noise ratio + snr = {} + if noise_rms == 0 or noise_int == 0: + logger.info("RMS of noise is zero, calculating an snr is not useful. All snrs are set to infinity.") + snr['peak_2_peak_amplitude'] = np.inf + snr['peak_amplitude'] = np.inf + snr['integrated_power'] = np.inf else: - SNR['integrated_power'] = np.sum(np.square(trace[signal_window_mask])) - noise_int - if SNR['integrated_power'] <= 0: - logger.debug("Integrated signal {0} smaller than noise {1}, power SNR 0".format(SNR['integrated_power'], noise_int)) - SNR['integrated_power'] = 0. + snr['integrated_power'] = np.sum(np.square(trace[signal_window_mask])) - noise_int + if snr['integrated_power'] <= 0: + logger.debug("Integrated signal {0} smaller than noise {1}, power snr 0".format(snr['integrated_power'], noise_int)) + snr['integrated_power'] = 0. else: - SNR['integrated_power'] /= (noise_int) - SNR['integrated_power'] = np.sqrt(SNR['integrated_power']) + snr['integrated_power'] /= (noise_int) + snr['integrated_power'] = np.sqrt(snr['integrated_power']) # calculate amplitude values - SNR['peak_2_peak_amplitude'] = np.max(trace[signal_window_mask]) - np.min(trace[signal_window_mask]) - SNR['peak_2_peak_amplitude'] /= noise_rms - SNR['peak_2_peak_amplitude'] /= 2 + snr['peak_2_peak_amplitude'] = np.max(trace[signal_window_mask]) - np.min(trace[signal_window_mask]) + snr['peak_2_peak_amplitude'] /= noise_rms + snr['peak_2_peak_amplitude'] /= 2 - SNR['peak_amplitude'] = np.max(np.abs(trace[signal_window_mask])) / noise_rms + snr['peak_amplitude'] = np.max(np.abs(trace[signal_window_mask])) / noise_rms # SCNR - SNR['Seckel_2_noise'] = 5 + snr['Seckel_2_noise'] = 5 + + # Calculating RPR (Root Power Ratio) + if noise_rms == 0: + root_power_ratio = np.inf + else: + wf_len = len(trace) + channel_wf = trace ** 2 + + # Calculate the smoothing window size based on sampling rate + dt = times[1] - times[0] + sum_win = 25 # Smoothing window in ns + sum_win_idx = int(np.round(sum_win / dt)) # Convert window size to sample points + + channel_wf = np.sqrt(uniform_filter1d(channel_wf, size=sum_win_idx, mode='constant')) + + # Find the maximum value of the smoothed waveform + max_bin = np.argmax(channel_wf) + max_val = channel_wf[max_bin] + + root_power_ratio = max_val / noise_rms + if self.__debug: plt.figure() @@ -154,7 +178,7 @@ def get_SNR(self, station_id, channel, det, stored_noise=False, rms_stage=None): plt.legend() plt.show() - return SNR, noise_rms + return snr, noise_rms, root_power_ratio @register_run() def run(self, evt, station, det, stored_noise=False, rms_stage='amp'): @@ -176,12 +200,14 @@ def run(self, evt, station, det, stored_noise=False, rms_stage='amp'): trace = channel.get_trace() h = np.abs(signal.hilbert(trace)) max_amplitude = np.max(np.abs(trace)) - logger.info(f"event {evt.get_run_number()}.{evt.get_id()} station {station.get_id()} channel {channel.get_id()} max amp = {max_amplitude:.6g} max amp env {h.max():.6g}") - if(logger.level >= logging.DEBUG): - tmp = "" - for amp in trace: - tmp += f"{amp:.6g}, " - logger.debug(tmp) + + logger.info(f"Event {evt.get_run_number()}.{evt.get_id()}, station.channel " + f"{station.get_id()}. {channel.get_id()}: max amp = " + f"{max_amplitude:.6g} max amp env {h.max():.6g}") + + if logger.level >= logging.DEBUG: + logger.debug(", ".join([f"{x:.6g}" for x in trace])) + channel[chp.signal_time] = times[np.argmax(h)] max_amplitude_station = max(max_amplitude_station, max_amplitude) channel[chp.maximum_amplitude] = max_amplitude @@ -189,9 +215,11 @@ def run(self, evt, station, det, stored_noise=False, rms_stage='amp'): channel[chp.P2P_amplitude] = np.max(trace) - np.min(trace) # Use noise precalculated from forced triggers - signal_to_noise, noise_rms = self.get_SNR(station.get_id(), channel, det, stored_noise=stored_noise, rms_stage=rms_stage) + signal_to_noise, noise_rms, root_power_ratio = self.get_SNR_and_RPR( + station.get_id(), channel, det, stored_noise=stored_noise, rms_stage=rms_stage) channel[chp.SNR] = signal_to_noise channel[chp.noise_rms] = noise_rms + channel[chp.root_power_ratio] = root_power_ratio station[stnp.channels_max_amplitude] = max_amplitude_station diff --git a/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py b/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py index 9f8d56f284..0126d0dc12 100644 --- a/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py +++ b/NuRadioReco/modules/io/RNO_G/readRNOGDataMattak.py @@ -14,7 +14,7 @@ import NuRadioReco.framework.station import NuRadioReco.framework.channel import NuRadioReco.framework.trigger -import NuRadioReco.framework.parameters +from NuRadioReco.framework.parameters import channelParameters, channelParametersRNOG from NuRadioReco.utilities import units import mattak.Dataset @@ -214,8 +214,8 @@ def __init__(self, run_table_path=None, load_run_table=True, log_level=logging.N "Runs can not be filtered.") except ImportError: self.logger.warn( - "import run_table failed. You can still use readRNOGData, but runs can not be filtered. " - "To install the run table, run\n\n" + "\nImport `from rnog_data.runtable import RunTable` failed. You can still use readRNOGData, " + "but runs can not be filtered. To install the run table, run:\n\n" "\tpip install git+ssh://git@github.com/RNO-G/rnog-runtable.git\n" ) else: @@ -391,7 +391,7 @@ def begin(self, self.__n_runs = 0 # Set verbose for mattak - verbose = mattak_kwargs.pop("verbose", self.logger.level >= logging.DEBUG) + verbose = mattak_kwargs.pop("verbose", self.logger.level <= logging.DEBUG) for dir_file in dirs_files: @@ -775,7 +775,6 @@ def _get_event(self, event_info, waveforms): evt: NuRadioReco.framework.event """ - trigger_time = event_info.triggerTime # only overwrite sampling rate if the stored value is invalid if self._overwrite_sampling_rate is not None and event_info.sampleRate in [0, None]: @@ -799,6 +798,9 @@ def _get_event(self, event_info, waveforms): readout_delays = event_info.readoutDelay for channel_id, wf in enumerate(waveforms): channel = NuRadioReco.framework.channel.Channel(channel_id) + # this allows to store `channelParametersRNOG` via `set_parameter` + channel.add_parameter_type(channelParametersRNOG) + if self._read_calibrated_data: channel.set_trace(wf * units.V, sampling_rate * units.GHz) else: @@ -812,7 +814,7 @@ def _get_event(self, event_info, waveforms): time_offset = get_time_offset(event_info.triggerType) + readout_delays[channel_id] channel.set_trace_start_time(-time_offset) # relative to event/trigger time if block_offsets is not None: - channel.set_parameter(NuRadioReco.framework.parameters.channelParameters.block_offsets, block_offsets.T[channel_id]) + channel.set_parameter(channelParameters.block_offsets, block_offsets.T[channel_id]) station.add_channel(channel) @@ -838,12 +840,13 @@ def run(self): dataset.setEntries((0, dataset.N())) # read all event infos of the entire dataset (= run) - for evtinfo, wf in dataset.iterate(calibrated=self._read_calibrated_data, - selectors=self._select_events, - max_entries_in_mem=self._max_in_mem): + for evtinfo, wf in dataset.iterate( + calibrated=self._read_calibrated_data, selectors=self._select_events, + max_entries_in_mem=self._max_in_mem): + t0 = time.time() evt = self._get_event(evtinfo, wf) - + self._time_run += time.time() - t0 yield evt @@ -936,14 +939,13 @@ def get_event(self, run_nr, event_id): return evt - def end(self): if self.__counter: self.logger.info( f"\n\tRead {self.__counter} events ({self.__skipped} events are skipped (filtered), {self.__invalid} invalid events)" f"\n\tTime to initialize data sets : {self._time_begin:.2f}s" - f"\n\tTime to read all events : {self._time_run:.2f}s" - f"\n\tTime to per event : {self._time_run / self.__counter:.2f}s" + f"\n\tTime to read all events : {self._time_run:.4f}s" + f"\n\tTime to per event : {self._time_run / self.__counter:.4f}s" f"\n\tRead {self.__n_runs} runs, skipped {self.__skipped_runs} runs.") else: self.logger.warning( @@ -951,6 +953,8 @@ def end(self): f"\n\tTime to initialize data sets : {self._time_begin:.2f}s" f"\n\tTime to read all events : {self._time_run:.2f}s") + def get_n_events(self): + return self._n_events_total ### we create a wrapper for readRNOGData to mirror the interface of the .nur reader class _readRNOGData_eventbrowser(readRNOGData): @@ -982,9 +986,6 @@ def get_event_i(self, i): def get_event(self, event_id): return super().get_event(*event_id) - def get_n_events(self): - return self._n_events_total - def get_detector(self): """Not implemented in mattak reader""" return None diff --git a/NuRadioReco/test/RNO_G/test_read_rnog_data.sh b/NuRadioReco/test/RNO_G/test_read_rnog_data.sh index 767ebab98c..4ec78a8b59 100755 --- a/NuRadioReco/test/RNO_G/test_read_rnog_data.sh +++ b/NuRadioReco/test/RNO_G/test_read_rnog_data.sh @@ -6,4 +6,4 @@ wget https://rnog-data.zeuthen.desy.de/rnog_share/forced_triggers/station23_run3 mkdir -p tests/data/station23/run325 mv station23_run325.root tests/data/station23/run325/combined.root -python3 NuRadioReco/examples/RNO_data/read_data_example/read_rnog.py tests/data/station23/run325/combined.root test.nur \ No newline at end of file +python3 NuRadioReco/examples/RNOG/read_rnog.py tests/data/station23/run325/combined.root test.nur \ No newline at end of file