script to find forced triggers by the event timing.py

rno_g_data_tools/scripts/find_forced_triggers_by_event_timing.py

+#!/usr/bin/env python
+
+""" Script to identify forced triggers that fired based on event timing
+
+Tries to identify PPS trigger events and SOFT FORCED trigger events
+"""
+
+#TODO: make these variables that can be passed by command line. Also: top level path of data needs to be added
+RUNMIN = 250
+RUNMAX = 300
+STATION = 21
+
+
+
+import uproot
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import astropy
+import astropy.time
+# from astropy.time import Time
+import NuRadioReco.framework.base_trace
+from NuRadioReco.utilities import units
+
+def find_pps(times):
+    """ Find events from the PPS (from GPS) triggers by checking if they arrive at (near) integer second
+  
+    Parameters
+    ----------
+    time: array of astropy.time.Time
+
+    Returns
+    -------
+    is_pps: array of bool
+    """
+    ymdhms = times.ymdhms
+    # get second part of times
+    second = np.array([e[-1] for e in ymdhms])
+    s_int = np.round(second,0)
+    # PPS tick always arrives precisely at the second,
+    # allow < 1e-5 s; sometimes it floats around...example: event 761 in station 22 / run 95
+    is_pps = np.isclose(second, s_int, atol=1e-5, rtol=0)
+    return is_pps
+
+
+def find_forced_triggers(times, mask, dt=astropy.time.TimeDelta(10.0558,format="sec")):
+    """ Find events from the forced software trigger (SOFT) by checking if events occur at the adjacent soft tick times
+  
+    Parameters
+    ----------
+    time: array of astropy.time.Time
+    mask: array of bool
+        Events to take into account, e.g. to ignore identified PPS events.
+    dt: astropy.time.Timedelta (default: 10.0558)
+        delta t between soft events. Typically longer than requested in the run config, at least at the beginning of
+        RNO-G data-taking
+        
+    Returns
+    -------
+    is_pps: array of bool
+    """
+
+    forced_candidates = np.zeros_like(times, dtype=bool)
+    forced_index = np.arange(len(times))
+    mindex = forced_index[mask]
+
+    # define how many soft ticks left and right should be used
+    padding = [5, 5]
+
+    first_time = min(times)
+    last_time = max(times)
+
+    for i_m_idx, m_idx in enumerate(mindex):
+        delta_first_time = times[m_idx] - first_time
+        delta_ticks_before = int(delta_first_time/dt)
+        delta_last_time = last_time - times[m_idx]
+        delta_ticks_after = int(delta_last_time/dt)
+
+        # define how many soft ticks left and right should be used
+        padding = [5, 5]
+        if delta_ticks_before < padding[0]:
+            padding = [delta_ticks_before, padding[0]+padding[1]-delta_ticks_before]
+        if delta_ticks_after < padding[1]:
+            padding = [min(padding[0]+padding[1]-delta_ticks_after, delta_ticks_before), delta_ticks_after]
+        """
+        # shift for first and last events in file
+        # TODO: catch extremely short runs, with less than padding width...
+        if i_m_idx < padding[0]:
+            padding = [i_m_idx, padding[1] + padding[0]-i_m_idx]
+        if len(times)-1 - i_m_idx < padding[1]:
+            padding = [padding[0]+padding[1] - (len(times)-1-i_m_idx), (len(times)-1-i_m_idx)]
+        #print(i_m_idx, padding)
+        """
+        #get the expected tick times for the mask
+        tick_mask = np.arange(-padding[0],padding[1]+1)*dt + times[m_idx]
+        #print(tick_mask)
+        matches = []
+        for tick in tick_mask:
+            delta = times-tick
+            #print(delta.sec)
+
+            # check, if any (!) event is close (<atol seconds) to the tick
+            close_to_tick = np.isclose(delta.sec, 0., atol=0.25, rtol=0.01)
+            n_close_to_tick =  np.count_nonzero(close_to_tick)
+            tick_matched = n_close_to_tick>0
+            matches.append(tick_matched)
+        #print(np.sum(matches))
+        #print("next")
+        #break
+        #if np.sum(matches)>6:
+        ### these are events that are likely a forced trigger
+        #    print(astrotimes[i_m_idx].iso, np.sum(matches))
+        forced_candidates[m_idx] = np.sum(matches)>8#6
+        #else:
+        ### these are the ones that do not have a forced trigger
+        #    print('  ',astrotimes[i_m_idx].iso, np.sum(matches))
+    return forced_candidates
+
+def find_ambiguous_forced_triggers(times, mask):
+    selected_times = times[mask]
+    ambiguous = np.zeros_like(times, dtype=bool)
+    ambiguous_index = np.arange(len(times))
+    for found_t, idx in zip(selected_times, ambiguous_index):
+        delta = selected_times - found_t
+        #print(found_t.iso,np.count_nonzero(np.isclose(delta.sec, 0 , atol=0.25, rtol=0))>1)
+        ambiguous[idx] = np.count_nonzero(np.isclose(delta.sec, 0 , atol=0.25, rtol=0))>1
+    return ambiguous
+
+def find_ambiguous_events(times, mask):
+    ambiguous = np.zeros_like(times, dtype=bool)
+    for i, time in enumerate(times):
+        delta = times - time
+        is_close = np.isclose(delta.sec, 0, atol=0.25, rtol=0)
+        ambiguous[i] = np.count_nonzero(is_close[mask])>1
+    return ambiguous
+
+
+def identify_triggers_from_times(astrotimes, timedelta_soft_sec = 10.0558, check_pps=True, check_soft=True, check_threshold=True):
+    if check_pps:
+        pps_mask = find_pps(astrotimes)
+    else:
+        pps_mask = np.zeros_like(astrotimes.jd, dtype=bool)
+    triggers = pd.DataFrame({"PPS": pps_mask})
+    # typically soft trigger has slightly less than the requested 0.1 Hz
+    dt = astropy.time.TimeDelta(timedelta_soft_sec, format="sec")
+    if check_soft:
+        triggers["soft_candidates"] = find_forced_triggers(astrotimes, ~pps_mask, dt)
+    else:
+        triggers["soft_candidates"] = np.zeros_like(astrotimes.jd, dtype=bool)
+    triggers["maybe_soft"] = triggers["soft_candidates"] & find_ambiguous_events(astrotimes, triggers['soft_candidates'])#(astrotimes) #find_ambiguous_forced_triggers(astrotimes, triggers['soft_candidates'])
+    triggers["is_soft"] = triggers['soft_candidates']&(~triggers["maybe_soft"])
+    triggers["softOrPPS"] = triggers["PPS"]|triggers["is_soft"]
+
+    if not check_threshold:
+        triggers["softOrPPS"] = np.ones_like(astrotimes.jd, dtype=bool)
+        triggers["is_soft"] = ~triggers["PPS"]
+        triggers["maybe_soft"] = np.zeros_like(astrotimes.jd, dtype=bool)
+
+    print("events:  {}, pps: {}, is_soft: {}, maybe_soft: {}".format(len(triggers),
+                np.sum(triggers["PPS"]),
+                np.sum(triggers["is_soft"]),
+                np.sum(triggers["maybe_soft"])
+                ))
+    return triggers
+
+def get_trigger_info(header_file, combined_file, output="trigger_types.npy", check_pps=True, check_soft=True, check_threshold=True):
+    # load the header file and identify triggers based on all events
+    header = uproot.open(header_file)['hdr']
+    times = np.array(header['trigger_time'])
+    if len(times)>10000:
+        print("too many events, probably calibration run...")
+        return 0
+    astrotimes = astropy.time.Time(times, format="unix")
+    triggers = identify_triggers_from_times(astrotimes, check_pps=check_pps, check_soft=check_soft, check_threshold=check_threshold)
+    soft_event_numbers = np.array(header['hdr/event_number'])[np.array(triggers['softOrPPS'])]
+    
+    # load waveform file to see which of the transferred events had forced or threshold triggers
+    waveforms = uproot.open(combined_file)['combined']
+    # select events
+    mask = np.array([evn in soft_event_numbers for evn in np.array(waveforms['waveforms/event_number'])])
+    print("Identified in transferred file {}:  {}".format(combined_file, np.sum(mask)))
+    trigger_dict = triggers.to_dict()
+    trigger_dict["forced_mask"] = mask
+    np.save(output, trigger_dict)
+
+
+
+def plot(station, run, data_file, selection_file, antennas=[12]):
+    trigger_dict = np.load(selection_file,allow_pickle=True)
+    #print(trigger_dict.item()['forced_mask'])
+
+    waveforms = uproot.open(data_file)['combined']
+    
+    traces = np.array(waveforms['waveforms/radiant_data[24][2048]'])
+    trigger_times = np.array(waveforms['header/trigger_time'])
+
+    # select events
+    mask = trigger_dict.item()['forced_mask']
+    traces = traces[mask]
+    trigger_times = trigger_times[mask]
+
+    def convert(data):
+        tr = NuRadioReco.framework.base_trace.BaseTrace()
+        tr.set_trace(data * units.mV, sampling_rate=3.2 * units.GHz)
+        spectrum =  np.abs(tr.get_frequency_spectrum())
+        # set DC to zero
+        spectrum[0] = 0
+        return spectrum
+
+    #plt.plot(spec.T)
+    tr = NuRadioReco.framework.base_trace.BaseTrace()
+    tr.set_trace(np.zeros(2048), sampling_rate=3.2 * units.GHz)
+    for antenna in antennas:
+        spec = np.apply_along_axis(convert, 1, np.array(traces[:,antenna,:]))
+        print(np.shape(spec))
+        fig = plt.figure()
+        plt.pcolormesh(tr.get_frequencies()[:601]*1e3,np.arange(len(spec[:,0])+1), spec[:,:600],norm=matplotlib.colors.LogNorm(vmin=0.1, vmax=50))
+        plt.colorbar()
+        plt.title("station {}, run {}".format(station,run))
+        plt.xlabel("frequencies [MHz]")
+        plt.tight_layout()
+        plt.savefig("noise_station{}_run{}_antenna{}.png".format(station,run,antenna))
+        plt.close(fig)
+        np.save("noise_station{}_run{}_antenna{}_array.npy".format(station,run,antenna), {"spectra":spec[:,:600], "traces":np.array(traces[:,antenna,:]), "times": np.array(trigger_times), "frequencies":tr.get_frequencies()[:600]*1e3, "mask": mask})
+
+
+
+if __name__ == "__main__":
+    run_summary = pd.read_csv("https://www.zeuthen.desy.de/~shallman/rnog_run_summary.csv")
+
+    for item, r_summ in run_summary.iterrows():
+        #print(r_summ)
+        if r_summ["soft_rate [Hz]"] != 0.1:
+            print("only calculating for 0.1 Hz soft rate at the moment")
+            continue
+        has_pps = r_summ["has_pps (PPS signal)"]
+        has_soft = r_summ["has_soft (forced)"]
+        has_threshold = r_summ["has_rf0 (surface)"] | r_summ["has_rf1 (deep)"]
+
+        station = r_summ.station
+        run = r_summ.run
+        
+        if run<RUNMIN or run>RUNMAX or station!=STATION:
+            continue
+        combined_file = r_summ.path
+        header_file = r_summ.path.replace("combined","headers")
+        try:
+            get_trigger_info(header_file, combined_file, output="trigger_flags_station{}_run{}.npy".format(station, run), check_pps=has_pps, check_soft=has_soft, check_threshold=has_threshold)
+            #plot(station, run, data_file=combined_file, selection_file="trigger_flags_station{}_run{}.npy", antennas = [12,13,14])
+        except:
+        #    #plot(station, run, antennas = [12,13])
+            print("station {}, run {} failed".format(station, run))
+            continue