Add telescope event source

src/ctapipe_io_zfits/conftest.py

@@ -62,7 +62,7 @@ def dl0_base(acada_base):
     dl0 = acada_base / "DL0"
     dl0.mkdir(exist_ok=True)
 
-    lst_base = dl0 / "LSTN-01" / acada_user / "acada-adh"
+    lst_base = dl0 / "TEL001" / acada_user / "acada-adh"
     lst_events = lst_base / "events"
     lst_monitoring = lst_base / "monitoring"
     array_triggers = dl0 / "array" / acada_user / "acada-adh" / "triggers"
@@ -80,7 +80,7 @@ def dl0_base(acada_base):
 @pytest.fixture(scope="session")
 def dummy_dl0(dl0_base):
     trigger_dir = dl0_base / "array" / acada_user / "acada-adh/triggers/2023/08/01/"
-    lst_event_dir = dl0_base / "LSTN-01" / acada_user / "acada-adh/events/2023/08/01/"
+    lst_event_dir = dl0_base / "TEL001" / acada_user / "acada-adh/events/2023/08/01/"
     subarray_id = 1
     sb_id = 123
     obs_id = 456
@@ -207,3 +207,9 @@ def open_next_event_file(sdh_id):
             time = time + 0.001 * u.s
 
     return trigger_path
+
+
+@pytest.fixture(scope="session")
+def dummy_tel_file(dummy_dl0, dl0_base):
+    name = "TEL001_SDH001_20230802T021531_SBID0000000000000000123_OBSID0000000000000000456_TEL_SHOWER_CHUNK000.fits.fz"
+    return dl0_base / "TEL001/ctao-acada-n/acada-adh/events/2023/08/01/" / name

src/ctapipe_io_zfits/dl0.py

@@ -15,6 +15,7 @@
     ObservationBlockContainer,
     SchedulingBlockContainer,
     TriggerContainer,
+    TelescopeTriggerContainer,
 )
 from ctapipe.core.traits import Integer
 from ctapipe.instrument import SubarrayDescription
@@ -34,6 +35,45 @@
 ARRAY_ELEMENTS = get_array_elements_by_id()
 
 
+def _is_compatible(input_url, extname, allowed_protos):
+    from astropy.io import fits
+
+    # this allows us to just use try/except around the opening of the fits file,
+    stack = ExitStack()
+
+    with stack:
+        try:
+            hdul = stack.enter_context(fits.open(input_url))
+        except Exception as e:
+            log.debug(f"Error trying to open input file as fits: {e}")
+            return False
+
+        if extname not in hdul:
+            log.debug("FITS file does not contain '%s' HDU", extname)
+            return False
+
+        header = hdul[extname].header
+
+    if header["XTENSION"] != "BINTABLE":
+        log.debug("%s HDU is not a bintable", extname)
+        return False
+
+    if not header.get("ZTABLE", False):
+        log.debug("ZTABLE is not in header or False")
+        return False
+
+    proto_class = header.get("PBFHEAD")
+    if proto_class is None:
+        log.debug("Missing PBFHEAD key")
+        return False
+
+    if proto_class not in allowed_protos:
+        log.debug(f"Unsupported PBFHEAD: {proto_class} not in {allowed_protos}")
+        return False
+
+    return True
+
+
 class ProtozfitsDL0EventSource(EventSource):
     """
     DL0 Protozfits EventSource.
@@ -76,19 +116,34 @@ def __init__(self, input_url=None, **kwargs):
             self.sb_id: SchedulingBlockContainer(sb_id=np.uint64(self.sb_id))
         }
 
+        # <prefix>/DL0/<ae-id>/<acada-user>/acada-adh/events/<YYYY>/<MM>/<DD>/ 
+        self._dl0_base = self.input_url.parents[7]
+        self._acada_user = self.input_url.parents[5].name
+        self._date_dirs = self.input_url.parent.relative_to(self.input_url.parents[3])
+
         self._open_telescope_files()
 
+    def _get_tel_events_directory(self, tel_id):
+        return (
+            self._dl0_base / f"TEL{tel_id:03d}" / self._acada_user
+            / "acada-adh/events" / self._date_dirs
+        )
+
+    @classmethod
+    def is_compatible(cls, input_url):
+        return _is_compatible(
+            input_url,
+            extname="SubarrayEvents",
+            allowed_protos={"DL0v1.Subarray.Event"},
+        )
+
     def _open_telescope_files(self):
         self._telescope_files = {}
         for tel_id in self.subarray.tel:
-            name = ARRAY_ELEMENTS[tel_id]["name"]
 
             # get the directory, where we should look for files
-            tel_dir = Path(
-                str(self.input_url.parent)
-                .replace("triggers", "events")
-                .replace("array", name)
-            )
+            tel_dir = self._get_tel_events_directory(tel_id)
+
             try:
                 first_file = sorted(
                     tel_dir.glob(f"*_SBID*{self.sb_id}_OBSID*{self.obs_id}*.fits.fz")
@@ -176,45 +231,124 @@ def _generator(self):
 
             yield array_event
 
+
+class ProtozfitsDL0TelescopeEventSource(EventSource):
+    """
+    DL0 Protozfits EventSource.
+
+    The ``input_url`` must be the subarray trigger file, the source
+    will then look for the other data files according to the filename and
+    directory schema layed out in the draft of the ACADA - DPPS ICD.
+    """
+    subarray_id = Integer(default_value=1).tag(config=True)
+
     @classmethod
     def is_compatible(cls, input_url):
-        from astropy.io import fits
+        return _is_compatible(
+            input_url,
+            extname="Events",
+            allowed_protos={"DL0v1.Telescope.Event"},
+        )
+
+    def __init__(self, input_url=None, **kwargs):
+        # this enables passing input_url as posarg, kwarg and via the config/parent
+        if input_url is not None:
+            kwargs["input_url"] = input_url
 
-        # this allows us to just use try/except around the opening of the fits file,
-        stack = ExitStack()
+        super().__init__(**kwargs)
 
-        with stack:
-            try:
-                hdul = stack.enter_context(fits.open(input_url))
-            except Exception as e:
-                log.debug(f"Error trying to open input file as fits: {e}")
-                return False
+        # we will open a lot of files, this helps keeping it clean
+        self._exit_stack = ExitStack()
+        self._subarray = build_subarray_description(self.subarray_id)
 
-            if "SubarrayEvents" not in hdul:
-                log.debug("FITS file does not contain an Events HDU, returning False")
-                return False
+        self._multi_file = self._exit_stack.enter_context(
+            MultiFiles(self.input_url)
+        )
+        self.sb_id = self._multi_file.data_stream.sb_id
+        self.obs_id = self._multi_file.data_stream.obs_id
+        self.tel_id = self._multi_file.data_stream.tel_id
 
-            header = hdul["SubarrayEvents"].header
+        self._observation_blocks = {
+            self.obs_id: ObservationBlockContainer(
+                obs_id=np.uint64(self.obs_id), sb_id=np.uint64(self.sb_id)
+            )
+        }
+        self._scheduling_blocks = {
+            self.sb_id: SchedulingBlockContainer(sb_id=np.uint64(self.sb_id))
+        }
 
-        if header["XTENSION"] != "BINTABLE":
-            log.debug("Events HDU is not a bintable")
-            return False
 
-        if not header.get("ZTABLE", False):
-            log.debug("ZTABLE is not in header or False")
-            return False
+    def close(self):
+        self._exit_stack.__exit__(None, None, None)
 
-        if header.get("ORIGIN", "") != "CTA":
-            log.debug("ORIGIN != CTA")
-            return False
+    def __exit__(self, exc_type, exc_value, traceback):
+        self._exit_stack.__exit__(exc_type, exc_value, traceback)
 
-        proto_class = header.get("PBFHEAD")
-        if proto_class is None:
-            log.debug("Missing PBFHEAD key")
-            return False
+    @property
+    def is_simulation(self) -> bool:
+        return False
 
-        if proto_class != "DL0v1.Subarray.Event":
-            log.debug(f"Unsupported PBFHEAD: {proto_class}")
-            return False
+    @property
+    def datalevels(self) -> Tuple[DataLevel]:
+        return (DataLevel.DL0, )
+
+    @property
+    def subarray(self) -> SubarrayDescription:
+        return self._subarray
+
+    @property
+    def observation_blocks(self) -> Dict[int, ObservationBlockContainer]:
+        return self._observation_blocks
 
-        return True
+    @property
+    def scheduling_blocks(self) -> Dict[int, SchedulingBlockContainer]:
+        return self._scheduling_blocks
+
+    def _fill_event(self, zfits_event) -> ArrayEventContainer:
+        tel_id = self.tel_id
+        # until ctapipe allows telescope event sources
+        # we have to fill an arrayevent with just one telescope here
+        time = cta_high_res_to_time(
+            zfits_event.event_time_s,
+            zfits_event.event_time_qns
+        )
+        array_event = ArrayEventContainer(
+            index=EventIndexContainer(
+                obs_id=self.obs_id,
+                event_id=zfits_event.event_id,
+            ),
+            trigger=TriggerContainer(
+                tels_with_trigger=[self.tel_id],
+                time=time,
+            ),
+        )
+        array_event.trigger.tel[tel_id] = TelescopeTriggerContainer(time=time)
+
+        n_channels = zfits_event.num_channels
+        n_pixels = zfits_event.num_pixels_survived
+        n_samples = zfits_event.num_samples
+        shape = (n_channels, n_pixels, n_samples)
+        # FIXME: ctapipe  can only handle a single gain
+        waveform = zfits_event.waveform.reshape(shape)[0]
+        offset = self._multi_file.data_stream.waveform_offset
+        scale = self._multi_file.data_stream.waveform_scale
+        waveform = (waveform.astype(np.float32) - offset) / scale
+
+        array_event.dl0.tel[tel_id] = DL0CameraContainer(
+            pixel_status=zfits_event.pixel_status,
+            event_type=EventType(zfits_event.event_type),
+            selected_gain_channel=np.zeros(n_pixels, dtype=np.int8),
+            event_time=cta_high_res_to_time(
+                zfits_event.event_time_s,
+                zfits_event.event_time_qns,
+            ),
+            waveform=waveform,
+            first_cell_id=zfits_event.first_cell_id,
+            # module_hires_local_clock_counter=zfits_event.module_hires_local_clock_counter,
+        )
+
+        return array_event
+
+    def _generator(self):
+        for event in self._multi_file:
+            yield self._fill_event(event)

src/ctapipe_io_zfits/tests/test_dl0.py

@@ -4,7 +4,7 @@
 from astropy.time import Time
 from ctapipe.core.tool import run_tool
 from ctapipe.instrument import SubarrayDescription
-from ctapipe.io import EventSource
+from ctapipe.io import EventSource, TableLoader
 from ctapipe.tools.process import ProcessorTool
 
 
@@ -60,3 +60,40 @@ def test_process(dummy_dl0, tmp_path):
             ],
             raises=True,
         )
+
+
+
+def test_telescope_event_source(dummy_tel_file):
+    from ctapipe_io_zfits.dl0 import ProtozfitsDL0TelescopeEventSource
+
+    assert ProtozfitsDL0TelescopeEventSource.is_compatible(dummy_tel_file)
+
+    with EventSource(dummy_tel_file) as source:
+        assert isinstance(source, ProtozfitsDL0TelescopeEventSource)
+
+        for event in source:
+            assert event.dl0.tel.keys() == {1}
+
+
+def test_process_tel_events(dummy_tel_file, tmp_path):
+    path = tmp_path / "dummy.dl1.h5"
+
+    with pytest.warns(UserWarning, match="Encountered an event with no R1 data"):
+        run_tool(
+            ProcessorTool(),
+            [
+                f"--input={dummy_tel_file}",
+                f"--output={path}",
+                "--write-images",
+                "--write-parameters",
+            ],
+            raises=True,
+        )
+
+    with TableLoader(path) as loader:
+        events = loader.read_telescope_events()
+        assert len(events) == 100
+        np.testing.assert_array_equal(events["obs_id"], 456)
+        np.testing.assert_array_equal(events["tel_id"], 1)
+        np.testing.assert_array_equal(events["event_id"], np.arange(1, 101))
+

src/ctapipe_io_zfits/tests/test_multifiles.py

@@ -4,7 +4,7 @@
 def test_multifiles(dummy_dl0, dl0_base):
     from ctapipe_io_zfits.multifile import MultiFiles
 
-    directory = dl0_base / "LSTN-01/ctao-acada-n/acada-adh/events/2023/08/01"
+    directory = dl0_base / "TEL001/ctao-acada-n/acada-adh/events/2023/08/01"
     filename = "TEL001_SDH001_20230802T021531_SBID0000000000000000123_OBSID0000000000000000456_TEL_SHOWER_CHUNK000.fits.fz"  # noqa
     path = directory / filename