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

NuRadioReco/examples/RNOG/rnog_config.yaml

@@ -0,0 +1,75 @@
+# 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: null
+    # 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
+
+channelAddCableDelay:
+  use: True
+
+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
+
+channelBlockOffsetFitter:
+  use: True
+  kwargs:
+    mode: 'auto'
+
+channelGlitchDetector:
+  use: True
+  begin:
+    kwargs:
+      max_deviation: 2  # 2 V = 2000mV

NuRadioReco/examples/RNOG/rnog_standard_data_processing.py

@@ -0,0 +1,167 @@
+import NuRadioReco.modules.channelAddCableDelay
+import NuRadioReco.modules.channelResampler
+import NuRadioReco.modules.channelBandPassFilter
+import NuRadioReco.modules.channelCWNotchFilter
+
+import NuRadioReco.modules.io.RNO_G.readRNOGDataMattak
+import NuRadioReco.modules.RNO_G.hardwareResponseIncorporator
+import NuRadioReco.modules.RNO_G.channelBlockOffsetFitter
+import NuRadioReco.modules.RNO_G.channelGlitchDetector
+
+import NuRadioReco.modules.io.eventWriter
+
+import NuRadioReco.detector.RNO_G.rnog_detector
+
+from NuRadioReco.utilities import units, logging as nulogging
+
+import argparse
+import logging
+import yaml
+import time
+
+logger = logging.getLogger("NuRadioReco.example.RNOG.rnog_standard_data_processing")
+logger.setLevel(logging.INFO)
+
+
+def use_module(name, config):
+    """
+    Simple helper function returns True if module is supposed to be used
+
+    Parameters
+    ----------
+    name : str
+        Name of the module
+    config : dict
+        (Yaml) configuration
+
+    Returns
+    -------
+    use : bool
+        True if module is supposed to be used
+    """
+    return name in config and config[name]["use"]
+
+
+def process_data(config):
+    """
+    Processing of RNO-G data.
+
+    Parameters
+    ----------
+    config : dict
+        (Yaml) configuration
+    """
+
+    # Initialize modules
+
+    channelAddCableDelay = NuRadioReco.modules.channelAddCableDelay.channelAddCableDelay()
+
+    channelResampler = NuRadioReco.modules.channelResampler.channelResampler()
+    channelResampler.begin()
+
+    channelBandPassFilter = NuRadioReco.modules.channelBandPassFilter.channelBandPassFilter()
+    channelBandPassFilter.begin()
+
+    channelCWNotchFilter = NuRadioReco.modules.channelCWNotchFilter.channelCWNotchFilter()
+    channelCWNotchFilter.begin()
+
+    eventWriter = NuRadioReco.modules.io.eventWriter.eventWriter()
+    eventWriter.begin(filename=config["eventWriter"]["kwargs"]['filename'])
+
+    hardwareResponseIncorporator = NuRadioReco.modules.RNO_G.hardwareResponseIncorporator.hardwareResponseIncorporator()
+    hardwareResponseIncorporator.begin()
+
+    channelBlockOffsetFitter = NuRadioReco.modules.RNO_G.channelBlockOffsetFitter.channelBlockOffsetFitter()
+    channelBlockOffsetFitter.begin()
+
+    channelGlitchDetector = NuRadioReco.modules.RNO_G.channelGlitchDetector.channelGlitchDetector()
+    channelGlitchDetector.begin(config["channelGlitchDetector"]["begin"]["kwargs"]["max_deviation"])
+
+    paths = config["readRNOGDataMattak"]["kwargs"].pop("filenames")
+    readRNOGDataMattak = NuRadioReco.modules.io.RNO_G.readRNOGDataMattak.readRNOGData(log_level=logging.INFO)
+    readRNOGDataMattak.begin(
+        paths, **config["readRNOGDataMattak"]["kwargs"],
+    )
+
+    # Initialize detector class
+    det = NuRadioReco.detector.RNO_G.rnog_detector.Detector(
+        **config["detector"]["kwargs"]
+    )
+
+    # 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(readRNOGDataMattak.run()):
+
+        if (idx + 1) % 50 == 0:
+            print(f'"Processing events: {idx + 1} / {readRNOGDataMattak.get_n_events()}\r', end="")
+
+        t0 = time.time()
+        # Loop over all stations.
+        for station in evt.get_stations():
+
+            det.update(station.get_station_time())
+
+            # Glitch Removal
+            if use_module("channelGlitchDetector", config):
+                channelGlitchDetector.run(evt, station, det)
+
+            # Correcting for block offsets
+            if use_module("channelBlockOffsetFitter", config):
+                channelBlockOffsetFitter.run(
+                    evt, station, det, **config["channelBlockOffsetFitter"]["kwargs"])
+
+            # Add cable delay
+            if use_module("channelAddCableDelay", config):
+                channelAddCableDelay.run(evt, station, det, mode='subtract')
+
+            # Resample
+            if use_module("channelResampler", config):
+                channelResampler.run(evt, station, det, **config["channelResampler"]["kwargs"])
+
+            # Hardware response
+            if use_module("hardwareResponseIncorporator", config):
+                hardwareResponseIncorporator.run(
+                    evt, station, det, sim_to_data=False, **config["hardwareResponseIncorporator"]["kwargs"])
+
+            # CW notch filter
+            if use_module("channelCWNotchFilter", config):
+                channelCWNotchFilter.run(evt, station, det)
+
+            # Bandpass filter
+            if use_module("channelBandPassFilter", config):
+                channelBandPassFilter.run(evt, station, det, **config["channelBandPassFilter"]["kwargs"])
+
+            # Write event - the RNO-G detector class is not stored within the nur files.
+            eventWriter.run(evt, det=None, mode=config["eventWriter"]["kwargs"]['mode'])
+
+        logger.debug("Time for event: %f", time.time() - t0)
+        t_total += time.time() - t0
+
+    readRNOGDataMattak.end()
+
+    logger.info(
+        f"Processed {idx + 1} events:"
+        f"\n\tTotal time: {t_total:.2f}s"
+        f"\n\tTime per event: {t_total / (idx + 1):.2f}s")
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(description='Run standard RNO-G data processing')
+
+    parser.add_argument('config', type=str, help='Yaml config file to steer the data procession and reconstruction')
+    parser.add_argument('--filenames', type=str, nargs="*", help='Specify root data files if not specified in the config file')
+
+    args = parser.parse_args()
+    # nulogging.set_general_log_level(logging.INFO)
+
+    with open(args.config, 'r') as f:
+        config = yaml.safe_load(f)
+
+    if config["readRNOGDataMattak"]["kwargs"]["filenames"] is None:
+        config["readRNOGDataMattak"]["kwargs"]["filenames"] = args.filenames
+
+    process_data(config)

NuRadioReco/examples/RNO_data/read_data_example/run_reconstruction.py → 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))