Merge branch 'issue1-trigger-time-shifted-by-42sec-in-root-files' into 'main'

Issue1 trigger time shifted by 42sec in root files

See merge request !1

rno_g_data_tools/offline_fixes/offlinefix_trigger_times.py

+import ROOT
+import sys
+import os
+from array import array
+import numpy as np
+import astropy
+from astropy.time import Time
+from pathlib import Path
+
+import logging
+logger = logging.getLogger("RNOG_ROOT_OFFLINE_FIX")
+logger.setLevel(logging.DEBUG)
+# create console handler and set level to debug
+ch = logging.StreamHandler()
+ch.setLevel(logging.DEBUG)
+# create formatter
+formatter = logging.Formatter('%(name)s - %(levelname)s - %(message)s')
+
+# add formatter to ch
+ch.setFormatter(formatter)
+# add ch to logger
+logger.addHandler(ch)
+
+# load libmattak.so at install location
+ROOT.gSystem.Load("/".join([os.getenv("MATTAK"),"libmattak.so"]))
+
+def sysclk_time(sysclk, sysclk_last_pps, sysclk_last_last_pps, ignore_overflow=False):
+    """Calculates the subsecond part of the trigger time from the 100MHz sysclk cycle count
+    Parameters
+    ----------
+    sysclk : int
+        The sysclk at trigger time
+    sysclk_last_pps : int
+        The sysclk of the last GPS PPS pulse
+    sysclk_last_last_pps : int
+        The sysclk of the 2nd last GPS PPS pulse
+    ignore_overflow : bool, optional
+        sysclk are 32 bit integers increasing with ~100 MHz, i.e. they will overflow after ~42s. If set to True the overflow is ignored, the calculated sysclk_time will be wrong! (default False)
+
+    Returns
+    -------
+    float
+        subsecond part of the trigger time
+    """
+    sysclk_diff = float(sysclk) - float(sysclk_last_pps)
+    sysclk_diff_last = float(sysclk_last_pps) - float(sysclk_last_last_pps)
+
+    if not ignore_overflow:
+        if sysclk_diff < 0:
+            sysclk_diff += 2**32
+        if sysclk_diff_last < 0:
+            sysclk_diff_last += 2**32
+
+    if sysclk_diff_last == 0:
+        logger.error("Same sysclk last: (%s) / last last: (%s), assuming 100MHz", sysclk_last_pps, sysclk_last_last_pps)
+        sysclk_diff_last = 1e8
+    return sysclk_diff/sysclk_diff_last
+
+
+def fix_trigger_time(input_filename, output_filename, tree_name="hdr", write_mode="RECREATE"):
+    """Fixes wrongly calculated trigger times
+
+    These were  induced by wrong calculation at the beginning of data taking and did not take into account the sysclk overflow correctly
+
+    Parameters
+    ----------
+    input_filename : str
+        The .root file to be fixed
+    output_filename : str
+        The output file with fixed trigger times.
+        The output is only written if any event had a wrong trigger time
+    tree_name : str, optional
+        Name of the input root tree. (default "hdr")
+        For 'hdr', assume the 'hdr' is also the branch name
+        For 'combined', assume 'header' is the branch name
+    write_mode : str, optional
+        root write mode (default "RECREATE")
+
+    Returns
+    -------
+    int
+        number of trigger times that were fixed
+    """
+
+    if tree_name == "hdr":
+        branch_name = "hdr"
+    if tree_name == "combined":
+        branch_name = "header"
+
+    outfile = ROOT.TFile(output_filename, write_mode)
+    logger.info("Fixing trigger times for file {}".format(input_filename))
+    logger.debug("Output file: {}".format(output_filename))
+
+    hdr = ROOT.mattak.Header()
+
+    infile = ROOT.TFile(input_filename, "READONLY")
+    try:
+        oldtree = infile.Get(tree_name)
+    except:
+        logger.error("tree {} not found in input file, skipping file".format(tree_name))
+        return np.nan
+
+    try:
+        oldtree.SetBranchAddress(branch_name, hdr)
+    except:
+        logger.error("tree {} in input file is empty/not found, skipping file".format(branch_name))
+        return np.nan
+
+    outfile.cd()
+    newtree = oldtree.CloneTree(0)
+    n_fixed = 0
+    for i in range(oldtree.GetEntries()):
+        oldtree.GetEntry(i)
+
+        time_before = hdr.trigger_time
+        
+        #subsecond part
+        hdr.trigger_time = sysclk_time(hdr.sysclk, hdr.sysclk_last_pps, hdr.sysclk_last_last_pps)
+        
+        # get second part from readout time
+        readout_time_subsecs = (hdr.readout_time - int(hdr.readout_time))
+        readout_time_fullsecs = int(hdr.readout_time)
+        if readout_time_subsecs <  hdr.trigger_time:
+            hdr.trigger_time += readout_time_fullsecs - 1
+        else:
+            hdr.trigger_time += readout_time_fullsecs
+
+        if not np.isclose(time_before, hdr.trigger_time, atol=1e-3, rtol=0):
+            try:
+                logger.info(" ".join(["Fixed time:", astropy.time.Time(time_before, format="unix").iso, "->", astropy.time.Time(hdr.trigger_time, format="unix").iso]))
+            except:
+                logger.error("Problematic time: {} -> {}".format(time_before, hdr.trigger_time))
+            n_fixed += 1
+        newtree.Fill()
+    if n_fixed > 0:
+        newtree.Write()
+        outfile.Close()
+        logger.warning("Fixed {} event times".format(n_fixed))
+    else:
+        outfile.Close()
+        os.remove(output_filename)
+        logger.warning("No wrong timestamps found. Not saving file.")
+    return n_fixed
+
+if __name__ == '__main__':
+    INPUT_DIR="/pnfs/ifh.de/acs/radio/diskonly/data/inbox/"
+    FIXED_DIR="/lustre/fs22/group/radio/shallman/RNO_G_DATA/inbox/"
+
+    logger.warning("INPUT DIRECTORY: {}".format(INPUT_DIR))
+    logger.warning("OUTPUT DIRECTORY: {}".format(FIXED_DIR))
+
+    # set to small number for testing   
+    nmax = 10#0000
+    
+    # fix header files
+    for i, filename_headers in enumerate(Path(INPUT_DIR).rglob("headers.root")):
+        if i>nmax:
+            break
+        # get directory
+        filename_repaired = Path(str(filename_headers).replace(INPUT_DIR, FIXED_DIR))
+        if not os.path.exists(filename_repaired.parent):
+            logger.info("Generating directories")
+            Path(filename_repaired.parent).mkdir(parents=True, exist_ok=True)
+        n_fixed = fix_trigger_time(str(filename_headers), str(filename_repaired))
+
+    # fix combined root files
+    for i, filename_headers in enumerate(Path(INPUT_DIR).rglob("combined.root")):
+        if i>nmax:
+            break
+        filename_repaired = Path(str(filename_headers).replace(INPUT_DIR, FIXED_DIR))
+        if not os.path.exists(filename_repaired.parent):
+            logger.info("Generating directories")
+            Path(filename_repaired.parent).mkdir(parents=True, exist_ok=True)
+        n_fixed = fix_trigger_time(str(filename_headers), str(filename_repaired), tree_name="combined")

rno_g_data_tools/offline_fixes/setenv.sh

+export MATTAK=/afs/ifh.de/group/radio/scratch/shallman/software/mattak/