diff --git a/src/ctapipe_io_zfits/conftest.py b/src/ctapipe_io_zfits/conftest.py index bc39180..e7bd76b 100644 --- a/src/ctapipe_io_zfits/conftest.py +++ b/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 diff --git a/src/ctapipe_io_zfits/dl0.py b/src/ctapipe_io_zfits/dl0.py index efebaca..bc6960b 100644 --- a/src/ctapipe_io_zfits/dl0.py +++ b/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)) } + # /DL0///acada-adh/events///
/ + 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) diff --git a/src/ctapipe_io_zfits/tests/test_dl0.py b/src/ctapipe_io_zfits/tests/test_dl0.py index 1cd26bc..8989c0a 100644 --- a/src/ctapipe_io_zfits/tests/test_dl0.py +++ b/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)) + diff --git a/src/ctapipe_io_zfits/tests/test_multifiles.py b/src/ctapipe_io_zfits/tests/test_multifiles.py index 90a4445..0f4ff35 100644 --- a/src/ctapipe_io_zfits/tests/test_multifiles.py +++ b/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