Merge pull request #52 from cta-observatory/dead-time-calc

Implement the dead time calculation

magicctapipe/scripts/lst1_magic/lst1_magic_dl2_to_dl3.py

@@ -44,16 +44,42 @@
 ORM_LON = -17.89064   # unit: [deg]
 ORM_HEIGHT = 2199.835   # unit: [m]
 
+dead_time_lst = 7.6e-6   # unit: [sec]
+dead_time_magic = 26e-6   # unit: [sec]
+
 MJDREF = Time(0, format='unix', scale='utc').mjd
 
 __all__ = [
+    'calculate_deadc',
     'load_dl2_data_file',
     'create_event_list',
     'create_gti_table',
     'create_pointing_table',
     'dl2_to_dl3',
 ]
 
+def calculate_deadc(time_diffs, dead_time):
+    """
+    Calculates the dead time correction factor.
+
+    Parameters
+    ----------
+    time_diffs: np.ndarray
+        Time differences of event arrival times
+    dead_time: float
+        Dead time due to the read out
+
+    Returns
+    -------
+    deadc: float
+        Dead time correction factor
+    """
+
+    rate = 1 / (time_diffs.mean() - dead_time)
+    deadc = 1 / (1 + rate * dead_time)
+
+    return deadc
+
 
 def load_dl2_data_file(input_file, quality_cuts, irf_type):
     """
@@ -72,6 +98,8 @@ def load_dl2_data_file(input_file, quality_cuts, irf_type):
     -------
     event_table: astropy.table.table.QTable
         Astropy table of DL2 events
+    deadc: float
+        Dead time correction factor for the input data
     """
 
     df_events = pd.read_hdf(input_file, key='events/parameters')
@@ -108,6 +136,35 @@ def load_dl2_data_file(input_file, quality_cuts, irf_type):
     n_events = len(df_events.groupby(['obs_id', 'event_id']).size())
     logger.info(f'--> {n_events} stereo events')
 
+    # Calculate the dead time correction factor.
+    # For MAGIC we select one telescope which has more number of events than the other:
+    logger.info('\nCalculating the dead time correction factor...')
+
+    deadc = 1
+    condition = '(time_diff > 0) & (time_diff < 0.1)'
+
+    time_diffs_lst = df_events.query(f'(tel_id == 1) & {condition}')['time_diff'].to_numpy()
+
+    if len(time_diffs_lst) > 0:
+        deadc_lst = calculate_deadc(time_diffs_lst, dead_time_lst)
+        logger.info(f'LST-1: {deadc_lst}')
+        deadc *= deadc_lst
+
+    time_diffs_m1 = df_events.query(f'(tel_id == 2) & {condition}')['time_diff'].to_numpy()
+    time_diffs_m2 = df_events.query(f'(tel_id == 3) & {condition}')['time_diff'].to_numpy()
+
+    if len(time_diffs_m1) >= len(time_diffs_m2):
+        deadc_magic = calculate_deadc(time_diffs_m1, dead_time_magic)
+        logger.info(f'MAGIC-I: {deadc_magic}')
+    else:
+        deadc_magic = calculate_deadc(time_diffs_m2, dead_time_magic)
+        logger.info(f'MAGIC-II: {deadc_magic}')
+
+    deadc *= deadc_magic
+
+    dead_time_fraction = 100 * (1 - deadc)
+    logger.info(f'--> Total dead time fraction: {dead_time_fraction:.2f}%')
+
     # Compute the mean of the DL2 parameters:
     df_dl2_mean = get_dl2_mean(df_events)
     df_dl2_mean.reset_index(inplace=True)
@@ -125,10 +182,10 @@ def load_dl2_data_file(input_file, quality_cuts, irf_type):
     event_table['reco_dec'] *= u.deg
     event_table['reco_energy'] *= u.TeV
 
-    return event_table
+    return event_table, deadc
 
 
-def create_event_list(event_table, effective_time, elapsed_time,
+def create_event_list(event_table, deadc,
                       source_name=None, source_ra=None, source_dec=None):
     """
     Creates an event list and its header.
@@ -137,10 +194,8 @@ def create_event_list(event_table, effective_time, elapsed_time,
     ----------
     event_table: astropy.table.table.QTable
         Astropy table of the DL2 events surviving gammaness cuts
-    effective_time: float
-        Effective time of the input data
-    elapsed_time: float
-        Elapsed time of the input data
+    deadc: float:
+        Dead time correction factor
     source_name: str
         Name of the observed source
     source_ra:
@@ -158,10 +213,10 @@ def create_event_list(event_table, effective_time, elapsed_time,
 
     time_start = Time(event_table['timestamp'][0], format='unix', scale='utc')
     time_end = Time(event_table['timestamp'][-1], format='unix', scale='utc')
+    time_diffs = np.diff(event_table['timestamp'])
 
-    delta_time = time_end.value - time_start.value
-
-    deadc = effective_time / elapsed_time
+    elapsed_time = np.sum(time_diffs)
+    effective_time = elapsed_time * deadc
 
     event_coords = SkyCoord(ra=event_table['reco_ra'], dec=event_table['reco_dec'], frame='icrs')
 
@@ -203,7 +258,7 @@ def create_event_list(event_table, effective_time, elapsed_time,
     event_header['TIMESYS'] = 'UTC'
     event_header['TIMEREF'] = 'TOPOCENTER'
     event_header['ONTIME'] = elapsed_time
-    event_header['TELAPSE'] = delta_time
+    event_header['TELAPSE'] = elapsed_time
     event_header['DEADC'] = deadc
     event_header['LIVETIME'] = effective_time
     event_header['OBJECT'] = source_name
@@ -340,12 +395,7 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
     logger.info('\nLoading the input DL2 data file:')
     logger.info(input_file_dl2)
 
-    event_table = load_dl2_data_file(input_file_dl2, quality_cuts, irf_type)
-
-    # ToBeUpdated: how to compute the effective time for the software coincidence?
-    # At the moment it does not consider any dead times, which slightly underestimates a source flux.
-    elapsed_time = event_table['timestamp'][-1] - event_table['timestamp'][0]
-    effective_time = elapsed_time
+    event_table, deadc = load_dl2_data_file(input_file_dl2, quality_cuts, irf_type)
 
     # Apply gammaness cuts:
     if 'GH_CUT' in header:
@@ -371,7 +421,7 @@ def dl2_to_dl3(input_file_dl2, input_file_irf, output_dir, config):
     # Create a event list HDU:
     logger.info('\nCreating an event list HDU...')
 
-    event_list, event_header = create_event_list(event_table, effective_time, elapsed_time, **config_dl3)
+    event_list, event_header = create_event_list(event_table, deadc, **config_dl3)
 
     hdu_event = fits.BinTableHDU(event_list, header=event_header, name='EVENTS')
     hdus.append(hdu_event)

magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py

@@ -41,6 +41,7 @@
 from astropy.time import Time
 from ctapipe.containers import EventType
 from ctapipe.instrument import SubarrayDescription
+from lstchain.reco.utils import add_delta_t_key
 from magicctapipe.utils import (
     check_tel_combination,
     save_pandas_to_table,
@@ -105,6 +106,9 @@ def load_lst_data_file(input_file):
     event_data.set_index(['obs_id', 'event_id', 'tel_id'], inplace=True)
     event_data.sort_index(inplace=True)
 
+    # Add the arrival time differences of consecutive events:
+    event_data = add_delta_t_key(event_data)
+
     # Exclude non-reconstructed events:
     event_data.dropna(subset=['intensity', 'time_gradient', 'alt_tel', 'az_tel'], inplace=True)
 
@@ -128,7 +132,8 @@ def load_lst_data_file(input_file):
 
     # Rename the column names:
     event_data.rename(
-        columns={'alt_tel': 'pointing_alt',
+        columns={'delta_t': 'time_diff',
+                 'alt_tel': 'pointing_alt',
                  'az_tel': 'pointing_az',
                  'leakage_pixels_width_1': 'pixels_width_1',
                  'leakage_pixels_width_2': 'pixels_width_2',

magicctapipe/scripts/lst1_magic/magic_calib_to_dl1.py

@@ -76,8 +76,9 @@ class EventInfoContainer(Container):
     tel_id = Field(-1, 'Telescope ID')
     pointing_alt = Field(-1, 'Telescope pointing altitude', u.rad)
     pointing_az = Field(-1, 'Telescope pointing azimuth', u.rad)
-    time_sec = Field(-1, 'Event trigger time second')
-    time_nanosec = Field(-1, 'Event trigger time nanosecond')
+    time_sec = Field(-1, 'Event trigger time second', u.s)
+    time_nanosec = Field(-1, 'Event trigger time nanosecond', u.ns)
+    time_diff = Field(-1, 'Event trigger time difference from the previous event', u.s)
     n_pixels = Field(-1, 'Number of pixels of a cleaned image')
     n_islands = Field(-1, 'Number of islands of a cleaned image')
 
@@ -172,6 +173,8 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False):
             for root_file in event_source.file_list:
                 logger.info(root_file)
 
+        time_diffs = event_source.event_time_diffs
+
     # Configure the MAGIC image cleaning:
     magic_clean = MAGICClean(camera_geom, config_cleaning)
 
@@ -311,8 +314,10 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False):
                 # here we set the integral and fractional parts separately as "time_sec" and "time_nanosec":
                 fractional, integral = np.modf(timestamp)
 
-                time_sec = int(np.round(integral))
-                time_nanosec = int(np.round(fractional * sec2nsec))
+                time_sec = u.Quantity(int(np.round(integral)), u.s)
+                time_nanosec = u.Quantity(int(np.round(fractional * sec2nsec)), u.ns)
+
+                time_diff = time_diffs[event.count]
 
                 # Set the event information to the container:
                 event_info = EventInfoContainer(
@@ -322,6 +327,7 @@ def magic_calib_to_dl1(input_file, output_dir, config, process_run=False):
                     pointing_az=event.pointing.tel[tel_id].azimuth,
                     time_sec=time_sec,
                     time_nanosec=time_nanosec,
+                    time_diff=time_diff,
                     n_pixels=n_pixels,
                     n_islands=n_islands,
                 )
@@ -406,3 +412,4 @@ def main():
 
 if __name__ == '__main__':
     main()
+