diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index 24d45def228..d18f7db436f 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -49,6 +49,7 @@ from ..instrument import SubarrayDescription from ..instrument.optics import FocalLengthKind from ..monitoring.interpolation import PointingInterpolator +from ..utils import IndexFinder from .astropy_helpers import read_table from .datalevels import DataLevel from .eventsource import EventSource @@ -230,6 +231,11 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): self.file_.root.configuration.observation.observation_block.col("obs_id") ) pointing_key = "/configuration/telescope/pointing" + # for ctapipe <0.21 + legacy_pointing_key = "/dl1/monitoring/telescope/pointing" + self._legacy_tel_pointing_finders = {} + self._legacy_tel_pointing_tables = {} + self._constant_telescope_pointing = {} if pointing_key in self.file_.root: for h5table in self.file_.root[pointing_key]._f_iter_nodes("Table"): @@ -237,6 +243,17 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): table = QTable(read_table(self.file_, h5table._v_pathname), copy=False) table.add_index("obs_id") self._constant_telescope_pointing[tel_id] = table + elif legacy_pointing_key in self.file_.root: + self.log.info( + "Found file written using ctapipe<0.21, using legacy pointing information" + ) + for node in self.file_.root[legacy_pointing_key]: + tel_id = int(node.name.removeprefix("tel_")) + table = QTable(read_table(self.file_, node._v_pathname), copy=False) + self._legacy_tel_pointing_tables[tel_id] = table + self._legacy_tel_pointing_finders[tel_id] = IndexFinder( + table["time"].mjd + ) self._simulated_shower_distributions = ( self._read_simulated_shower_distributions() @@ -767,12 +784,18 @@ def _fill_telescope_pointing(self, event, tel_pointing_interpolator=None): return for tel_id in event.trigger.tels_with_trigger: - tel_pointing = self._constant_telescope_pointing.get(tel_id) - if tel_pointing is None: - continue - - current = tel_pointing.loc[event.index.obs_id] - event.pointing.tel[tel_id] = TelescopePointingContainer( - altitude=current["telescope_pointing_altitude"], - azimuth=current["telescope_pointing_azimuth"], - ) + if ( + tel_pointing := self._constant_telescope_pointing.get(tel_id) + ) is not None: + current = tel_pointing.loc[event.index.obs_id] + event.pointing.tel[tel_id] = TelescopePointingContainer( + altitude=current["telescope_pointing_altitude"], + azimuth=current["telescope_pointing_azimuth"], + ) + elif (finder := self._legacy_tel_pointing_finders.get(tel_id)) is not None: + index = finder.closest(event.trigger.time.mjd) + row = self._legacy_tel_pointing_tables[tel_id][index] + event.pointing.tel[tel_id] = TelescopePointingContainer( + altitude=row["altitude"], + azimuth=row["azimuth"], + ) diff --git a/src/ctapipe/io/tests/test_hdf5eventsource.py b/src/ctapipe/io/tests/test_hdf5eventsource.py index d1ebf8c3f28..fc8fe6bfede 100644 --- a/src/ctapipe/io/tests/test_hdf5eventsource.py +++ b/src/ctapipe/io/tests/test_hdf5eventsource.py @@ -285,3 +285,17 @@ def test_provenance(dl1_file, provenance): assert inputs[0]["url"] == str(dl1_file) meta = _read_reference_metadata_hdf5(dl1_file) assert inputs[0]["reference_meta"].product.id_ == meta.product.id_ + + +def test_pointing_old_file(): + input_url = "dataset://gamma_diffuse_dl2_train_small.dl2.h5" + + n_read = 0 + with HDF5EventSource(input_url, max_events=5) as source: + for e in source: + assert e.pointing.tel.keys() == set(e.trigger.tels_with_trigger) + for pointing in e.pointing.tel.values(): + assert u.isclose(pointing.altitude, 70 * u.deg) + assert u.isclose(pointing.azimuth, 0 * u.deg) + n_read += 1 + assert n_read == 5 diff --git a/src/ctapipe/tools/tests/test_process.py b/src/ctapipe/tools/tests/test_process.py index d94fb9481cf..49e70144602 100644 --- a/src/ctapipe/tools/tests/test_process.py +++ b/src/ctapipe/tools/tests/test_process.py @@ -499,7 +499,7 @@ def test_only_trigger_and_simulation(tmp_path): [ pytest.param( "dataset://gamma_diffuse_dl2_train_small.dl2.h5", - ["--no-write-images", "--max-events=5"], + ["--no-write-images", "--max-events=20"], id="0.17", ) ], @@ -514,9 +514,24 @@ def test_on_old_file(input_url, args, tmp_path): f"--config={config}", f"--input={input_url}", f"--output={output_path}", + "--write-showers", "--overwrite", *args, ], cwd=tmp_path, raises=True, ) + + with tables.open_file(output_path) as f: + assert "/configuration/telescope/pointing" in f.root + + with TableLoader(output_path) as loader: + events = loader.read_subarray_events() + + # check that we have valid reconstructions and that in case + # we don't, is_valid is False, regression test for #2651 + finite_reco = np.isfinite(events["HillasReconstructor_alt"]) + assert np.any(finite_reco) + np.testing.assert_array_equal( + finite_reco, events["HillasReconstructor_is_valid"] + )