Merge pull request #1252 from cta-observatory/teff_calculation

Update eff. time computation

lstchain/reco/tests/test_utils.py

@@ -3,6 +3,7 @@
 from astropy.time import Time
 from astropy.coordinates import SkyCoord
 from astropy.table import QTable
+from ctapipe.containers import EventType
 import numpy as np
 import pandas as pd
 
@@ -144,13 +145,20 @@ def test_get_obstime_real():
     n_cosmics = np.random.poisson(cosmics_rate * t_obs)
 
     timestamps = np.random.uniform(0, t_obs, n_cosmics)
+    event_types = EventType.SUBARRAY.value * np.ones(n_cosmics)
     timestamps = np.append(timestamps, np.arange(t0_pedestal, t_obs, 1 / pedestal_rate))
+    n_pedestals = len(timestamps) - n_cosmics
+    event_types = np.append(event_types,  EventType.SKY_PEDESTAL.value * np.ones(n_pedestals))
     timestamps = np.append(
         timestamps, np.arange(t0_flatfield, t_obs, 1 / flatfield_rate)
     )
+    n_flatfield = len(timestamps) - n_cosmics - n_pedestals
+    event_types = np.append(event_types, EventType.FLATFIELD.value * np.ones(n_flatfield))
+
     # sort events by timestamp:
+    sorted_event_types = event_types[np.argsort(timestamps)]
     timestamps.sort()
-
+    
     # time to previous event:
     delta_t = np.insert(np.diff(timestamps), 0, 0)
 
@@ -169,6 +177,7 @@ def test_get_obstime_real():
         {
             "delta_t": delta_t[recorded_events][cut],
             "dragon_time": timestamps[recorded_events][cut],
+            "event_type": sorted_event_types[recorded_events][cut],
         }
     )
     t_eff, t_elapsed = utils.get_effective_time(events)
@@ -179,7 +188,8 @@ def test_get_obstime_real():
     # now test with a QTable:
     a = delta_t[recorded_events][cut] * u.s
     b = timestamps[recorded_events][cut] * u.s
-    events = QTable([a, b], names=("delta_t", "dragon_time"))
+    c = sorted_event_types[recorded_events][cut]
+    events = QTable([a, b, c], names=("delta_t", "dragon_time", "event_type"))
     t_eff, t_elapsed = utils.get_effective_time(events)
     print(t_obs, t_elapsed, true_t_eff, t_eff)
     # test accuracy to 0.05%:

lstchain/reco/utils.py

@@ -17,6 +17,7 @@
 from astropy.coordinates import AltAz, SkyCoord
 from astropy.time import Time
 from ctapipe.coordinates import CameraFrame
+from ctapipe.containers import EventType
 from ctapipe_io_lst import OPTICS
 from ctapipe_io_lst.constants import LST1_LOCATION
 
@@ -654,8 +655,19 @@ def get_effective_time(events):
     t_eff: astropy Quantity (in seconds, if input has no units)
     t_elapsed: astropy Quantity (ditto)
     """
-    timestamp = np.array(events["dragon_time"])
-    delta_t = np.array(events["delta_t"])
+
+    # For consistency with the event selection applied in the DL2 to DL3 stage
+    # we require the events to be tagged as "physics triggers". In this way we
+    # count as livetime only periods in which the telescope is recording showers
+    # and properly tagging them. Without this filter the effective time was in 
+    # *rare* occasions (example: run 7199) a few % larger than it should be - it was 
+    # including periods in which showers were tagged "UNKNOWN", which were not 
+    # present in the DL3 file. For most runs the effect is zero.
+
+    typemask = events["event_type"] == EventType.SUBARRAY.value
+
+    timestamp = np.array(events["dragon_time"][typemask])
+    delta_t = np.array(events["delta_t"][typemask])
 
     if not isinstance(timestamp, u.Quantity):
         timestamp *= u.s