Merge pull request #191 from cta-observatory/cta_r1_main

Add support for new CTA r1 data format (EVBv6)

environment.yml

@@ -6,9 +6,9 @@ dependencies:
   - astropy=5.2
   - python=3.11 # nail the python version, so conda does not try upgrading / dowgrading
   - ctapipe=0.19
+  - protozfits=2.2
   - eventio
   - corsikaio
-  - protozfits=2
   - zeromq
   - ipython
   - numba

setup.cfg

@@ -32,7 +32,7 @@ zip_safe = False
 install_requires=
     astropy~=5.2
     ctapipe~=0.19.0
-    protozfits~=2.0
+    protozfits~=2.2
     numpy>=1.20
 
 [options.package_data]

src/ctapipe_io_lst/anyarray_dtypes.py

@@ -50,16 +50,18 @@
     ('unknown', np.uint8),  # called arbitraryInformation in C-Struct
 ]).newbyteorder('<')
 
+
 SWAT_DTYPE = np.dtype([
-    ('timestamp', np.uint64),
-    ('counter1', np.uint32),
-    ('counter2', np.uint32),
-    ('event_type', np.uint8),
-    ('camera_flag', np.uint8),
-    ('camera_event_num', np.uint32),
-    ('array_flag', np.uint8),
-    ('event_num', np.uint32),
-]).newbyteorder('<')
+    ("assigned_event_id", np.uint32),
+    ("trigger_id", np.uint64),
+    ("trigger_type", np.uint8),
+    ("trigger_time_s", np.uint32),
+    ("trigger_time_qns", np.uint32),
+    ("readout_requested", np.bool_),
+    ("data_available", np.bool_),
+    ("hardware_stereo_trigger_mask", np.uint16),
+    ("negative_flag", np.uint8),
+], align=True).newbyteorder('<')
 
 
 def parse_tib_10MHz_counter(counter):

src/ctapipe_io_lst/containers.py

@@ -1,6 +1,7 @@
 """
 Container structures for data that should be read or written to disk
 """
+import numpy as np
 from ctapipe.core import Container, Field, Map
 from ctapipe.containers import ArrayEventContainer
 from functools import partial
@@ -21,41 +22,58 @@ class LSTServiceContainer(Container):
     """
 
     # Data from the CameraConfig table
+    # present in CTA R1
     telescope_id = Field(-1, "telescope id")
-    cs_serial = Field(None, "serial number of the camera server")
-    configuration_id = Field(None, "id of the CameraConfiguration")
-    date = Field(None, "NTP start of run date")
+    local_run_id = Field(-1, "local run id")
+    date = Field(None, "config_time_s in data model")
+    configuration_id = Field(None, "camera_config_id in the date model")
+    pixel_ids = Field([], "pixel_id_map in the data model")
+    module_ids = Field([], "module_id_map in the data model")
+
+    num_modules = Field(-1, "number of modules")
     num_pixels = Field(-1, "number of pixels")
+    num_channels = Field(-1, "number of channels")
     num_samples = Field(-1, "num samples")
-    pixel_ids = Field([], "id of the pixels in the waveform array")
+
     data_model_version = Field(None, "data model version")
+    calibration_service_id = Field(-1, "calibration service id")
+    calibration_algorithm_id = Field(-1, "calibration service id")
 
-    idaq_version = Field(0, "idaq version")
+    # in debug in CTA R1 debug
+    cs_serial = Field(None, "serial number of the camera server")
+    idaq_version = Field(0, "idaq/evb version")
     cdhs_version = Field(0, "cdhs version")
+    tdp_type = Field(None, "tdp type")
+    tdp_action = Field(None, "tdp action")
+    ttype_pattern = Field(None, "ttype_pattern")
+
+    # only in old R1
     algorithms = Field(None, "algorithms")
     pre_proc_algorithms = Field(None, "pre processing algorithms")
-    module_ids = Field([], "module ids")
-    num_modules = Field(-1, "number of modules")
 
 
 class LSTEventContainer(Container):
     """
     Container for Fields that are specific to each LST event
     """
+    # present in CTA R1 but not in ctapipe R1CameraEvent
+    pixel_status = Field(None, "status of the pixels (n_pixels)", dtype=np.uint8)
+    first_capacitor_id = Field(None, "first capacitor id")
+    calibration_monitoring_id = Field(None, "calibration id of applied pre-calibration")
+    local_clock_counter = Field(None, "Dragon local 133 MHz counter (n_modules)")
 
-    # Data from the CameraEvent table
-    configuration_id = Field(None, "id of the CameraConfiguration")
-    event_id = Field(None, "local id of the event")
-    tel_event_id = Field(None, "global id of the event")
-    pixel_status = Field([], "status of the pixels (n_pixels)")
-    ped_id = Field(None, "tel_event_id of the event used for pedestal substraction")
-    module_status = Field([], "status of the modules (n_modules)")
+    # in debug event
+    module_status = Field(None, "status of the modules (n_modules)")
     extdevices_presence = Field(None, "presence of data for external devices")
-
-    tib_event_counter = Field(-1, "TIB event counter")
-    tib_pps_counter = Field(-1, "TIB pps counter")
-    tib_tenMHz_counter = Field(-1, "TIB 10 MHz counter")
-    tib_stereo_pattern = Field(0, "TIB stereo pattern")
+    chips_flags = Field(None, "chips flags")
+    charges_hg = Field(None, "charges of high gain channel")
+    charges_lg = Field(None, "charges of low gain channel")
+    tdp_action = Field(None, "tdp action")
+
+    tib_event_counter = Field(np.uint32(0), "TIB event counter", dtype=np.uint32)
+    tib_pps_counter = Field(np.uint16(0), "TIB pps counter", dtype=np.uint16)
+    tib_tenMHz_counter = Field(np.uint32(0), "TIB 10 MHz counter", dtype=np.uint32)
+    tib_stereo_pattern = Field(np.uint16(0), "TIB stereo pattern", dtype=np.uint16)
     tib_masked_trigger = Field(0, "TIB trigger mask")
 
     ucts_event_counter =  Field(-1, "UCTS event counter")
@@ -71,26 +89,31 @@ class LSTEventContainer(Container):
     ucts_num_in_bunch = Field(-1, "UCTS num in bunch (for debugging)")
     ucts_cdts_version = Field(-1, "UCTS CDTS version")
 
-    swat_timestamp = Field(-1, "SWAT timestamp")
-    swat_counter1 = Field(-1, "SWAT event counter 1")
-    swat_counter2 = Field(-1, "SWAT event counter 2")
-    swat_event_type = Field(0, "SWAT event type")
-    swat_camera_flag = Field(-1, "SWAT camera flag ")
-    swat_camera_event_num = Field(-1, "SWAT camera event number")
-    swat_array_flag = Field(-1, "SWAT array negative flag")
-    swat_array_event_num = Field(-1, "SWAT array event number")
+    swat_assigned_event_id = Field(np.uint32(0), "SWAT assigned event id")
+    swat_trigger_id = Field(np.uint64(0), "SWAT trigger id")
+    swat_trigger_type = Field(np.uint8(0), "SWAT trigger type")
+    swat_trigger_time_s = Field(np.uint32(0), "SWAT trigger_time_s")
+    swat_trigger_time_qns = Field(np.uint32(0), "SWAT trigger_time_qns")
+    swat_readout_requested = Field(np.bool_(False), "SWAT readout requested")
+    swat_data_available = Field(np.bool_(False), "SWAT data available")
+    swat_hardware_stereo_trigger_mask = Field(np.uint16(0), "SWAT hardware stereo trigger mask")
+    swat_negative_flag = Field(np.uint8(0), "SWAT negative flag")
 
     pps_counter= Field(None, "Dragon pulse per second counter (n_modules)")
     tenMHz_counter = Field(None, "Dragon 10 MHz counter (n_modules)")
     event_counter = Field(None, "Dragon event counter (n_modules)")
     trigger_counter = Field(None, "Dragon trigger counter (n_modules)")
-    local_clock_counter = Field(None, "Dragon local 133 MHz counter (n_modules)")
 
-    chips_flags = Field(None, "chips flags")
-    first_capacitor_id = Field(None, "first capacitor id")
+    # Only in old R1
+    configuration_id = Field(None, "id of the CameraConfiguration")
+    event_id = Field(None, "global id of the event")
+    tel_event_id = Field(None, "local id of the event")
+    ped_id = Field(None, "tel_event_id of the event used for pedestal substraction")
+
     drs_tag_status = Field(None, "DRS tag status")
     drs_tag = Field(None, "DRS tag")
 
+    # custom here
     ucts_jump = Field(False, "A ucts jump happened in the current event")
 
 

src/ctapipe_io_lst/event_time.py

@@ -19,7 +19,57 @@
     TimeUnixTai.epoch_scale = 'tai'
 
 
+EPOCH = Time(0, format='unix_tai', scale='tai')
+DAY_TO_S = 86400
 CENTRAL_MODULE = 132
+S_TO_QNS = 4e9
+QNS_TO_S = 0.25e-9
+
+
+def cta_high_res_to_time(seconds, quarter_nanoseconds):
+    '''Convert cta high resolution timestamp to astropy Time'''
+    # unix_tai accepts two floats for maximum precision
+    # we can just pass integral and fractional part
+    fractional_seconds = quarter_nanoseconds * QNS_TO_S
+    return Time(
+        seconds,
+        fractional_seconds, 
+        format='unix_tai',
+        # this is only for displaying iso timestamp, not any actual precision
+        precision=9,
+    )
+
+
+def to_seconds(days):
+    '''Returns whole and fractional seconds from a number in days'''
+    seconds = days * DAY_TO_S
+    return np.divmod(seconds, 1)
+
+
+def time_to_cta_high(time):
+    '''Convert astropy Time to cta high precision timestamp'''
+    # make sure we are in TAI
+    time = time.tai
+
+    # internally, astropy always uses jd values
+    # jd1 is integral and jd2 is in [-0.5, 0.5]
+    # we get the integral and fractional seconds from both jd values
+    # relative to the epoch
+    seconds_jd1, fractional_seconds_jd1 = to_seconds(time.jd1 - EPOCH.jd1)
+    seconds_jd2, fractional_seconds_jd2 = to_seconds(time.jd2 - EPOCH.jd2)
+
+    # add up the integral number of seconds
+    seconds = seconds_jd1 + seconds_jd2
+
+    # convert fractional seconds to quarter nanoseconds
+    fractional_seconds = fractional_seconds_jd1 + fractional_seconds_jd2
+    quarter_nanoseconds = np.round(fractional_seconds * S_TO_QNS)
+
+    # convert to uint32
+    seconds = np.asanyarray(seconds, dtype=np.uint32)
+    quarter_nanoseconds = np.asanyarray(quarter_nanoseconds, dtype=np.uint32)
+    return seconds, quarter_nanoseconds
+
 
 
 # fix for https://github.com/ipython/traitlets/issues/637

src/ctapipe_io_lst/multifiles.py

@@ -1,5 +1,4 @@
 import re
-import warnings
 from pathlib import Path
 from collections import namedtuple, defaultdict
 from itertools import count
@@ -160,9 +159,10 @@ def _load_next_subrun(self, stream):
             self._files.pop(stream).close()
 
         Provenance().add_input_file(str(path), "R0")
-        self._files[stream] = File(str(path))
+        file_ = File(str(path))
+        self._files[stream] = file_
         self.log.info("Opened file %s", path)
-        self._events_tables[stream] = self._files[stream].Events
+        self._events_tables[stream] = file_.Events
         self._headers[stream] = self._events_tables[stream].header
         dvr_applied = self._headers[stream].get("LSTDVR", False)
         if self.dvr_applied is None:
@@ -175,7 +175,13 @@ def _load_next_subrun(self, stream):
         self._events.put_nowait(NextEvent(event.event_id, event, stream))
 
         # make sure we have a camera config
-        config = next(self._files[stream].CameraConfig)
+        if hasattr(file_, "CameraConfig"):
+            config = next(file_.CameraConfig)
+            self.cta_r1 = False
+        else:
+            # new files use CameraConfiguration
+            self.cta_r1 = True
+            config = next(file_.CameraConfiguration)
 
         if self.camera_config is None:
             self.camera_config = config