Adding a script for the standard data processing of RNO-G data

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

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)

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)