diff --git a/.pylintrc b/.pylintrc index d09e823fb17..5dbceccf948 100644 --- a/.pylintrc +++ b/.pylintrc @@ -425,6 +425,7 @@ good-names=i, x, y, n, + f, ex, _ diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py index c03cc8f1be4..8ce9ce9c867 100644 --- a/ctapipe/conftest.py +++ b/ctapipe/conftest.py @@ -137,9 +137,69 @@ def prod5_proton_simtel_path(): @pytest.fixture(scope="session") def dl1_tmp_path(tmp_path_factory): + """Temporary directory for global dl1 test data""" return tmp_path_factory.mktemp("dl1") +@pytest.fixture(scope="session") +def dl2_tmp_path(tmp_path_factory): + """Temporary directory for global dl2 test data""" + return tmp_path_factory.mktemp("dl2") + + +@pytest.fixture(scope="session") +def dl2_shower_geometry_file(dl2_tmp_path, prod5_gamma_simtel_path): + """ + File containing both parameters and shower geometry from a gamma simulation set. + """ + from ctapipe.core import run_tool + from ctapipe.tools.process import ProcessorTool + + output = dl2_tmp_path / "gamma.training.h5" + + # prevent running process multiple times in case of parallel tests + with FileLock(output.with_suffix(output.suffix + ".lock")): + if output.is_file(): + return output + + argv = [ + f"--input={prod5_gamma_simtel_path}", + f"--output={output}", + "--write-images", + "--write-stereo-shower", + "--max-events=20", + ] + assert run_tool(ProcessorTool(), argv=argv, cwd=dl2_tmp_path) == 0 + return output + + +@pytest.fixture(scope="session") +def dl2_shower_geometry_file_type(dl2_tmp_path, prod5_gamma_simtel_path): + """ + File containing both parameters and shower geometry from a gamma simulation set. + """ + from ctapipe.core import run_tool + from ctapipe.tools.process import ProcessorTool + + output = dl2_tmp_path / "gamma_by_type.training.h5" + + # prevent running process multiple times in case of parallel tests + with FileLock(output.with_suffix(output.suffix + ".lock")): + if output.is_file(): + return output + + argv = [ + f"--input={prod5_gamma_simtel_path}", + f"--output={output}", + "--write-images", + "--write-stereo-shower", + "--max-events=20", + "--DataWriter.split_datasets_by=tel_type", + ] + assert run_tool(ProcessorTool(), argv=argv, cwd=dl2_tmp_path) == 0 + return output + + @pytest.fixture(scope="session") def dl1_file(dl1_tmp_path, prod5_gamma_simtel_path): """ @@ -166,6 +226,32 @@ def dl1_file(dl1_tmp_path, prod5_gamma_simtel_path): return output +@pytest.fixture(scope="session") +def dl1_by_type_file(dl1_tmp_path, prod5_gamma_simtel_path): + """ + DL1 file containing both images and parameters from a gamma simulation set. + """ + from ctapipe.tools.process import ProcessorTool + from ctapipe.core import run_tool + + output = dl1_tmp_path / "gamma_by_type.dl1.h5" + + # prevent running stage1 multiple times in case of parallel tests + with FileLock(output.with_suffix(output.suffix + ".lock")): + if output.is_file(): + return output + + argv = [ + f"--input={prod5_gamma_simtel_path}", + f"--output={output}", + "--write-images", + "--max-events=20", + "--DataWriter.split_datasets_by=tel_type", + ] + assert run_tool(ProcessorTool(), argv=argv, cwd=dl1_tmp_path) == 0 + return output + + @pytest.fixture(scope="session") def dl1_image_file(dl1_tmp_path, prod5_gamma_simtel_path): """ diff --git a/ctapipe/core/tool.py b/ctapipe/core/tool.py index 9e8b2447c8c..03985b9e818 100644 --- a/ctapipe/core/tool.py +++ b/ctapipe/core/tool.py @@ -52,8 +52,8 @@ class Tool(Application): ``name``, ``description`` and ``examples`` class attributes as strings. The ``aliases`` attribute can be set to cause a lower-level `~ctapipe.core.Component` parameter to become a high-level command-line - parameter (See example below). The `setup()`, `start()`, and - `finish()` methods should be defined in the sub-class. + parameter (See example below). The `setup`, `start`, and + `finish` methods should be defined in the sub-class. Additionally, any `ctapipe.core.Component` used within the `Tool` should have their class in a list in the ``classes`` attribute, @@ -61,7 +61,7 @@ class Tool(Application): tool. Once a tool is constructed and the virtual methods defined, the - user can call the `run()` method to setup and start it. + user can call the `run` method to setup and start it. .. code:: python @@ -254,14 +254,14 @@ def setup(self): @abstractmethod def start(self): """main body of tool (override in subclass). This is automatically - called after `initialize()` when the `run()` is called. + called after `Tool.initialize` when the `Tool.run` is called. """ pass @abstractmethod def finish(self): """finish up (override in subclass). This is called automatically - after `start()` when `run()` is called.""" + after `Tool.start` when `Tool.run` is called.""" self.log.info("Goodbye") def run(self, argv=None): diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py index 6100d0fb393..98ebba40335 100644 --- a/ctapipe/instrument/subarray.py +++ b/ctapipe/instrument/subarray.py @@ -1,6 +1,7 @@ """ Description of Arrays or Subarrays of telescopes """ +from typing import Dict, List, Union from pathlib import Path import warnings @@ -74,7 +75,7 @@ class SubarrayDescription: def __init__(self, name, tel_positions=None, tel_descriptions=None): self.name = name self.positions = tel_positions or dict() - self.tels = tel_descriptions or dict() + self.tels: Dict[int, TelescopeDescription] = tel_descriptions or dict() if self.positions.keys() != self.tels.keys(): raise ValueError("Telescope ids in positions and descriptions do not match") @@ -360,17 +361,17 @@ def peek(self): plt.tight_layout() @lazyproperty - def telescope_types(self): + def telescope_types(self) -> List[TelescopeDescription]: """ list of telescope types in the array""" return list({tel for tel in self.tel.values()}) @lazyproperty - def camera_types(self): + def camera_types(self) -> List[CameraDescription]: """ list of camera types in the array """ return list({tel.camera for tel in self.tel.values()}) @lazyproperty - def optics_types(self): + def optics_types(self) -> List[OpticsDescription]: """ list of optics types in the array """ return list({tel.optics for tel in self.tel.values()}) @@ -391,6 +392,41 @@ def get_tel_ids_for_type(self, tel_type): return [id for id, descr in self.tels.items() if str(descr) == tel_str] + def get_tel_ids( + self, telescopes: List[Union[int, str, TelescopeDescription]] + ) -> List[int]: + """ + Convert a list of telescope ids and telescope descriptions to + a list of unique telescope ids. + + Parameters + ---------- + telescopes: List[Union[int, str, TelescopeDescription]] + List of Telescope IDs and descriptions. + Supported inputs for telescope descriptions are instances of + `~ctapipe.instrument.TelescopeDescription` as well as their + string representation. + + Returns + ------- + tel_ids: List[int] + List of unique telescope ids matching ``telescopes`` + """ + ids = set() + + valid_tel_types = {str(tel_type) for tel_type in self.telescope_types} + + for telescope in telescopes: + if isinstance(telescope, int): + ids.add(telescope) + + if isinstance(telescope, str) and telescope not in valid_tel_types: + raise ValueError("Invalid telescope type input.") + + ids.update(self.get_tel_ids_for_type(telescope)) + + return sorted(ids) + def __eq__(self, other): if not isinstance(other, self.__class__): return False diff --git a/ctapipe/instrument/telescope.py b/ctapipe/instrument/telescope.py index 7adb622ca48..69c48ef84ea 100644 --- a/ctapipe/instrument/telescope.py +++ b/ctapipe/instrument/telescope.py @@ -27,7 +27,7 @@ class TelescopeDescription: Describes a Cherenkov Telescope and its associated `~ctapipe.instrument.OpticsDescription` and `~ctapipe.instrument.CameraDescription` - Parameters + Attributes ---------- name: str Telescope name @@ -40,9 +40,25 @@ class TelescopeDescription: """ def __init__( - self, name, tel_type, optics: OpticsDescription, camera: CameraDescription + self, + name: str, + tel_type: str, + optics: OpticsDescription, + camera: CameraDescription, ): + if not isinstance(name, str): + raise TypeError("`name` must be a str") + + if not isinstance(tel_type, str): + raise TypeError("`tel_type` must be a str") + + if not isinstance(optics, OpticsDescription): + raise TypeError("`optics` must be an instance of `OpticsDescription`") + + if not isinstance(camera, CameraDescription): + raise TypeError("`camera` must be an instance of `CameraDescription`") + self.name = name self.type = tel_type self.optics = optics diff --git a/ctapipe/instrument/tests/test_subarray.py b/ctapipe/instrument/tests/test_subarray.py index 2404969d4b9..62ca7afa8c5 100644 --- a/ctapipe/instrument/tests/test_subarray.py +++ b/ctapipe/instrument/tests/test_subarray.py @@ -1,10 +1,12 @@ """ Tests for SubarrayDescriptions """ +from copy import deepcopy + import numpy as np from astropy import units as u from astropy.coordinates import SkyCoord -from ctapipe.coordinates import TelescopeFrame -from copy import deepcopy +import pytest +from ctapipe.coordinates import TelescopeFrame from ctapipe.instrument import ( CameraDescription, OpticsDescription, @@ -188,3 +190,24 @@ def test_hdf_duplicate_string_repr(tmp_path): assert ( read.tel[1].optics.num_mirror_tiles == read.tel[2].optics.num_mirror_tiles + 1 ) + + +def test_get_tel_ids(example_subarray): + """Test for SubarrayDescription.get_tel_ids""" + subarray = example_subarray + sst = TelescopeDescription.from_name("SST-ASTRI", "CHEC") + + telescopes = [1, 2, "MST_MST_FlashCam", sst] + tel_ids = subarray.get_tel_ids(telescopes) + + true_tel_ids = ( + subarray.get_tel_ids_for_type("MST_MST_FlashCam") + + subarray.get_tel_ids_for_type(sst) + + [1, 2] + ) + + assert sorted(tel_ids) == sorted(true_tel_ids) + + # test invalid telescope type + with pytest.raises(Exception): + tel_ids = subarray.get_tel_ids(["It's a-me, Mario!"]) diff --git a/ctapipe/io/__init__.py b/ctapipe/io/__init__.py index 17e57976e1a..c0e8b10e6c8 100644 --- a/ctapipe/io/__init__.py +++ b/ctapipe/io/__init__.py @@ -2,6 +2,7 @@ from .eventsource import EventSource from .hdf5tableio import HDF5TableReader, HDF5TableWriter from .tableio import TableWriter, TableReader +from .tableloader import TableLoader from .datalevels import DataLevel from .astropy_helpers import read_table from .datawriter import DataWriter, DATA_MODEL_VERSION @@ -21,6 +22,7 @@ "HDF5TableReader", "TableWriter", "TableReader", + "TableLoader", "EventSeeker", "EventSource", "SimTelEventSource", diff --git a/ctapipe/io/astropy_helpers.py b/ctapipe/io/astropy_helpers.py index 647643d9fee..caab438861f 100644 --- a/ctapipe/io/astropy_helpers.py +++ b/ctapipe/io/astropy_helpers.py @@ -6,7 +6,7 @@ from pathlib import Path import tables -from astropy.table import Table +from astropy.table import Table, join import numpy as np from .tableio import ( @@ -18,7 +18,7 @@ from contextlib import ExitStack -__all__ = ["read_table"] +__all__ = ["read_table", "join_allow_empty"] def read_table(h5file, path, start=None, stop=None, step=None, condition=None) -> Table: @@ -119,3 +119,41 @@ def _parse_hdf5_attrs(table): other_attrs[attr] = str(value) if isinstance(value, np.str_) else value return column_transforms, column_descriptions, other_attrs + + +def join_allow_empty(left, right, keys, join_type="left", **kwargs): + """ + Join two astropy tables, allowing both sides to be empty tables. + + See https://github.com/astropy/astropy/issues/12012 for why + this is necessary. + + This behaves as `~astropy.table.join`, with the only difference of + allowing empty tables to be joined. + """ + + left_empty = len(left) == 0 + right_empty = len(right) == 0 + + if join_type == "inner": + if left_empty: + return left.copy() + if right_empty: + return right.copy() + + elif join_type == "left": + if left_empty or right_empty: + return left.copy() + + elif join_type == "right": + if left_empty or right_empty: + return right.copy() + + elif join_type == "outer": + if left_empty: + return right.copy() + + if right_empty: + return left.copy() + + return join(left, right, keys, join_type=join_type, **kwargs) diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py new file mode 100644 index 00000000000..101e87a13fe --- /dev/null +++ b/ctapipe/io/tableloader.py @@ -0,0 +1,343 @@ +""" +Class and related functions to read DL1 (a,b) and/or DL2 (a) data +from an HDF5 file produced with ctapipe-process. +""" + +import re +from collections import defaultdict +from typing import Dict + +import numpy as np +from astropy.table import join, vstack, Table +import tables + +from ..core import Component, traits, Provenance +from ..instrument import SubarrayDescription +from .astropy_helpers import read_table, join_allow_empty + +__all__ = ["TableLoader"] + +PARAMETERS_GROUP = "/dl1/event/telescope/parameters" +IMAGES_GROUP = "/dl1/event/telescope/images" +GEOMETRY_GROUP = "/dl2/event/subarray/geometry" +TRIGGER_TABLE = "/dl1/event/subarray/trigger" +SHOWER_TABLE = "/simulation/event/subarray/shower" +TRUE_IMAGES_GROUP = "/simulation/event/telescope/images" +TRUE_PARAMETERS_GROUP = "/simulation/event/telescope/parameters" + +by_id_RE = re.compile(r"tel_\d+") + + +SUBARRAY_EVENT_KEYS = ["obs_id", "event_id"] +TELESCOPE_EVENT_KEYS = ["obs_id", "event_id", "tel_id"] + + +def _empty_telescope_events_table(): + """ + Create a new astropy table with correct column names and dtypes + for telescope event based data. + """ + return Table(names=TELESCOPE_EVENT_KEYS, dtype=[np.int64, np.int64, np.int16]) + + +def _empty_subarray_events_table(): + """ + Create a new astropy table with correct column names and dtypes + for subarray event based data. + """ + return Table(names=SUBARRAY_EVENT_KEYS, dtype=[np.int64, np.int64]) + + +def _get_structure(h5file): + """ + Check if data is stored by telescope id or by telescope type in h5file. + """ + + if PARAMETERS_GROUP in h5file: + group = h5file.root[PARAMETERS_GROUP] + elif IMAGES_GROUP in h5file: + group = h5file.root[IMAGES_GROUP] + else: + raise ValueError(f"No DL1 parameters data in h5file: {h5file}") + + key = next(iter(group._v_children)) # pylint: disable=protected-access + + if by_id_RE.fullmatch(key): + return "by_id" + + return "by_type" + + +class TableLoader(Component): + """ + Load telescope-event or subarray-event data from ctapipe HDF5 files + + This class provides high-level access to data stored in ctapipe HDF5 files, + such as created by the ctapipe-process tool (`~ctapipe.tools.process.ProcessorTool`). + + The following `TableLoader` methods load data from all relevant tables, + depending on the options, and joins them into single tables: + * `TableLoader.read_subarray_events` + * `TableLoader.read_telescope_events` + + `TableLoader.read_telescope_events_by_type` retuns a dict with a table per + telescope type, which is needed for e.g. DL1 image data that might have + different shapes for each of the telescope types as tables do not support + variable length columns. + """ + + input_url = traits.Path(directory_ok=False, exists=True).tag(config=True) + + load_dl1_images = traits.Bool(False, help="load extracted images").tag(config=True) + load_dl1_parameters = traits.Bool( + True, help="load reconstructed image parameters" + ).tag(config=True) + load_dl2_geometry = traits.Bool( + False, help="load reconstructed shower geometry information" + ).tag(config=True) + load_simulated = traits.Bool(False, help="load simulated shower information").tag( + config=True + ) + load_true_images = traits.Bool(False, help="load simulated shower images").tag( + config=True + ) + load_true_parameters = traits.Bool( + False, help="load image parameters obtained from true images" + ).tag(config=True) + load_trigger = traits.Bool(True, help="load subarray trigger information").tag( + config=True + ) + load_instrument = traits.Bool( + False, help="join subarray instrument information to each event" + ).tag(config=True) + + def __init__(self, input_url=None, **kwargs): + # enable using input_url as posarg + if input_url not in {None, traits.Undefined}: + kwargs["input_url"] = input_url + super().__init__(**kwargs) + + self.subarray = SubarrayDescription.from_hdf(self.input_url) + self.h5file = tables.open_file(self.input_url, mode="r") + + Provenance().add_input_file(self.input_url, role="Event data") + + try: + self.structure = _get_structure(self.h5file) + except ValueError: + self.structure = None + + self.instrument_table = None + if self.load_instrument: + table = self.subarray.to_table() + optics = self.subarray.to_table(kind="optics") + optics["optics_index"] = np.arange(len(optics)) + optics.remove_columns(["name", "description", "type"]) + table = join( + table, + optics, + keys="optics_index", + # conflicts for TAB_VER, TAB_TYPE, not needed here, ignore + metadata_conflicts="silent", + ) + + table.remove_columns(["optics_index", "camera_index"]) + self.instrument_table = table + + def close(self): + """Close the underlying hdf5 file""" + self.h5file.close() + + def __del__(self): + self.close() + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_value, traceback): + self.close() + + def _read_telescope_table(self, group, tel_id): + if self.structure == "by_id": + key = f"{group}/tel_{tel_id:03d}" + condition = None + else: + key = f"{group}/{self.subarray.tel[tel_id]!s}" + condition = f"tel_id == {tel_id}" + + if key in self.h5file: + table = read_table(self.h5file, key, condition=condition) + else: + table = _empty_telescope_events_table() + + return table + + def read_subarray_events(self): + """Read subarray-based event information. + + Returns + ------- + table: astropy.io.Table + Table with primary index columns "obs_id" and "event_id". + """ + table = _empty_subarray_events_table() + + if self.load_trigger: + trigger = read_table(self.h5file, TRIGGER_TABLE) + table = join_allow_empty(table, trigger, SUBARRAY_EVENT_KEYS, "outer") + + if self.load_simulated and SHOWER_TABLE in self.h5file: + showers = read_table(self.h5file, SHOWER_TABLE) + table = join_allow_empty(table, showers, SUBARRAY_EVENT_KEYS, "outer") + + if self.load_dl2_geometry: + shower_geometry_group = self.h5file.root[GEOMETRY_GROUP] + + for reconstructor in shower_geometry_group._v_children: + geometry = read_table(self.h5file, f"{GEOMETRY_GROUP}/{reconstructor}") + + # rename DL2 columns to explicit reconstructor + # TBD: we could skip this if only 1 reconstructor is present + # or simply find another way to deal with multiple reconstructions + for col in set(geometry.colnames) - set(SUBARRAY_EVENT_KEYS): + geometry.rename_column(col, f"{reconstructor}_{col}") + + table = join_allow_empty(table, geometry, SUBARRAY_EVENT_KEYS, "outer") + + return table + + def _read_telescope_events_for_id(self, tel_id): + """Read telescope-based event information for a single telescope. + + This is the most low-level function doing the actual reading. + + Parameters + ---------- + tel_id: int + Telescope identification number. + + Returns + ------- + table: astropy.io.Table + Table with primary index columns "obs_id", "event_id" and "tel_id". + """ + + if tel_id is None: + raise ValueError("Please, specify a telescope ID.") + + table = _empty_telescope_events_table() + + if self.load_dl1_parameters: + parameters = self._read_telescope_table(PARAMETERS_GROUP, tel_id) + table = join_allow_empty( + table, parameters, join_type="outer", keys=TELESCOPE_EVENT_KEYS + ) + + if self.load_dl1_images: + images = self._read_telescope_table(IMAGES_GROUP, tel_id) + table = join_allow_empty( + table, images, join_type="outer", keys=TELESCOPE_EVENT_KEYS + ) + + if self.load_true_images: + true_images = self._read_telescope_table(TRUE_IMAGES_GROUP, tel_id) + table = join_allow_empty( + table, true_images, join_type="outer", keys=TELESCOPE_EVENT_KEYS + ) + + if self.load_true_parameters: + true_parameters = self._read_telescope_table(TRUE_PARAMETERS_GROUP, tel_id) + + for col in set(true_parameters.colnames) - set(TELESCOPE_EVENT_KEYS): + true_parameters.rename_column(col, f"true_{col}") + + table = join_allow_empty( + table, true_parameters, join_type="outer", keys=TELESCOPE_EVENT_KEYS + ) + + if self.load_instrument: + table = join_allow_empty( + table, self.instrument_table, keys=["tel_id"], join_type="left" + ) + + return table + + def _read_telescope_events_for_ids(self, tel_ids): + + table = vstack( + [self._read_telescope_events_for_id(tel_id) for tel_id in tel_ids] + ) + + return table + + def _join_subarray_info(self, table): + subarray_events = self.read_subarray_events() + table = join_allow_empty( + table, subarray_events, keys=SUBARRAY_EVENT_KEYS, join_type="left" + ) + return table + + def read_telescope_events(self, telescopes=None): + """ + Read telescope-based event information. + + If the corresponding traitlets are True, also subarray event information + is joined onto the table. + + Parameters + ---------- + telescopes: Optional[List[Union[int, str, TelescopeDescription]]] + A list containing any combination of telescope IDs and/or + telescope descriptions. If None, all available telescopes are read. + + Returns + ------- + events: astropy.io.Table + Table with primary index columns "obs_id", "event_id" and "tel_id". + """ + + if telescopes is None: + tel_ids = self.subarray.tel.keys() + else: + tel_ids = self.subarray.get_tel_ids(telescopes) + + table = self._read_telescope_events_for_ids(tel_ids) + + if any([self.load_trigger, self.load_simulated, self.load_dl2_geometry]): + table = self._join_subarray_info(table) + + return table + + def read_telescope_events_by_type(self, telescopes=None) -> Dict[str, Table]: + """Read telescope-based event information. + + Parameters + ---------- + telescopes: List[Union[int, str, TelescopeDescription]] + Any list containing a combination of telescope IDs or telescope_descriptions. + + Returns + ------- + tables: dict(astropy.io.Table) + Dictionary of tables organized by telescope types + Table with primary index columns "obs_id", "event_id" and "tel_id". + """ + + if telescopes is None: + tel_ids = self.subarray.tel.keys() + else: + tel_ids = self.subarray.get_tel_ids(telescopes) + + by_type = defaultdict(list) + + for tel_id in tel_ids: + key = str(self.subarray.tel[tel_id]) + by_type[key].append(self._read_telescope_events_for_id(tel_id)) + + by_type = {k: vstack(ts) for k, ts in by_type.items()} + + if any([self.load_trigger, self.load_simulated, self.load_dl2_geometry]): + for key, table in by_type.items(): + by_type[key] = self._join_subarray_info(table) + + return by_type diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py new file mode 100644 index 00000000000..74ea39d9bd8 --- /dev/null +++ b/ctapipe/io/tests/test_table_loader.py @@ -0,0 +1,185 @@ +import pytest +import tables +import numpy as np + + +@pytest.fixture(params=["by_type", "by_id"]) +def test_file(request, dl1_file, dl1_by_type_file): + """Test dl1 files in both structures""" + if request.param == "by_type": + f = dl1_by_type_file + else: + f = dl1_file + + return request.param, f + + +@pytest.fixture(params=["by_type", "by_id"]) +def test_file_dl2(request, dl2_shower_geometry_file, dl2_shower_geometry_file_type): + """Test dl2 files in both structures""" + if request.param == "by_type": + f = dl2_shower_geometry_file + else: + f = dl2_shower_geometry_file_type + + return request.param, f + + +def test_get_structure(test_file): + """Test _get_structure""" + from ctapipe.io.tableloader import _get_structure + + expected, path = test_file + + with tables.open_file(path, "r") as f: + assert _get_structure(f) == expected + + +def test_telescope_events_for_tel_id(test_file): + """Test loading data for a single telescope""" + from ctapipe.io.tableloader import TableLoader + + _, dl1_file = test_file + + loader = TableLoader(dl1_file, load_dl1_parameters=True, load_trigger=True) + + with loader as table_loader: + table = table_loader.read_telescope_events([25]) + assert "hillas_length" in table.colnames + assert "time" in table.colnames + assert "event_type" in table.colnames + assert np.all(table["tel_id"] == 25) + + with TableLoader(dl1_file, load_dl1_images=True) as table_loader: + table = table_loader.read_telescope_events([25]) + assert "image" in table.colnames + assert np.all(table["tel_id"] == 25) + + assert not table_loader.h5file.isopen + + +def test_load_instrument(test_file): + """Test joining instrument data onto telescope events""" + from ctapipe.io.tableloader import TableLoader + + _, dl1_file = test_file + + with TableLoader(dl1_file, load_instrument=True) as table_loader: + expected = table_loader.subarray.tel[25].optics.equivalent_focal_length + table = table_loader.read_telescope_events([25]) + assert "equivalent_focal_length" in table.colnames + assert np.all(table["equivalent_focal_length"] == expected) + + +def test_load_simulated(test_file): + """Test joining simulation info onto telescope events""" + from ctapipe.io.tableloader import TableLoader + + _, dl1_file = test_file + + with TableLoader(dl1_file, load_simulated=True) as table_loader: + table = table_loader.read_subarray_events() + assert "true_energy" in table.colnames + + table = table_loader.read_telescope_events([25]) + assert "true_energy" in table.colnames + + +def test_true_images(test_file): + """Test joining true images onto telescope events""" + from ctapipe.io.tableloader import TableLoader + + _, dl1_file = test_file + + with TableLoader( + dl1_file, load_dl1_parameters=False, load_true_images=True + ) as table_loader: + table = table_loader.read_telescope_events(["MST_MST_NectarCam"]) + assert "true_image" in table.colnames + + +def test_true_parameters(test_file): + """Test joining true parameters onto telescope events""" + from ctapipe.io.tableloader import TableLoader + + _, dl1_file = test_file + + with TableLoader( + dl1_file, load_dl1_parameters=False, load_true_parameters=True + ) as table_loader: + table = table_loader.read_telescope_events() + assert "true_hillas_intensity" in table.colnames + + +def test_read_subarray_events(test_file_dl2): + """Test reading subarray events""" + + from ctapipe.io.tableloader import TableLoader + + _, dl2_file = test_file_dl2 + + with TableLoader( + dl2_file, load_dl2_geometry=True, load_simulated=True, load_trigger=True + ) as table_loader: + table = table_loader.read_subarray_events() + assert "HillasReconstructor_alt" in table.colnames + assert "true_energy" in table.colnames + assert "time" in table.colnames + + +def test_read_telescope_events_type(test_file_dl2): + """Test reading telescope events for a given telescope type""" + + from ctapipe.io.tableloader import TableLoader + + _, dl2_file = test_file_dl2 + + with TableLoader( + dl2_file, + load_dl1_images=False, + load_dl1_parameters=False, + load_dl2_geometry=True, + load_simulated=True, + load_true_images=True, + load_trigger=False, + load_instrument=True, + ) as table_loader: + + table = table_loader.read_telescope_events(["MST_MST_FlashCam"]) + + assert "HillasReconstructor_alt" in table.colnames + assert "true_energy" in table.colnames + assert "true_image" in table.colnames + assert set(table["tel_id"].data).issubset([25, 125, 130]) + assert "equivalent_focal_length" in table.colnames + + +def test_read_telescope_events_by_type(test_file_dl2): + """Test reading telescope events for by types""" + + from ctapipe.io.tableloader import TableLoader + + _, dl2_file = test_file_dl2 + + with TableLoader( + dl2_file, + load_dl1_images=False, + load_dl1_parameters=False, + load_dl2_geometry=True, + load_simulated=True, + load_true_images=True, + load_trigger=False, + load_instrument=True, + ) as table_loader: + + tables = table_loader.read_telescope_events_by_type([25, 130]) + + for tel_type in ["MST_MST_NectarCam", "MST_MST_FlashCam"]: + + table = tables[tel_type] + + assert "HillasReconstructor_alt" in table.colnames + assert "true_energy" in table.colnames + assert "true_image" in table.colnames + assert set(table["tel_id"].data).issubset([25, 125, 130]) + assert "equivalent_focal_length" in table.colnames diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py index a99ac8cdac5..e9c9610d2b0 100644 --- a/ctapipe/tools/process.py +++ b/ctapipe/tools/process.py @@ -21,6 +21,8 @@ DataLevel.DL1_PARAMETERS, ] +__all__ = ["ProcessorTool"] + class ProcessorTool(Tool): """ @@ -209,6 +211,9 @@ def _write_processing_statistics(self): ) def start(self): + """ + Process events + """ self.log.info("(re)compute DL1: %s", self.should_compute_dl1) self.log.info("(re)compute DL2: %s", self.should_compute_dl2) self.event_source.subarray.info(printer=self.log.info) @@ -233,6 +238,9 @@ def start(self): self.write(event) def finish(self): + """ + Last steps after processing events. + """ self.write.write_simulation_histograms(self.event_source) self.write.finish() self._write_processing_statistics() diff --git a/docs/conf.py b/docs/conf.py index 7c1afe236f6..1894a3d84a0 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -106,6 +106,7 @@ def setup(app): ("py:class", "traitlets.config.application.Application"), ("py:obj", "traitlets.config.boolean_flag"), ("py:obj", "traitlets.TraitError"), + ("py:obj", "-v"), ("py:meth", "MetaHasDescriptors.__init__"), ("py:meth", "HasTraits.__new__"), ("py:meth", "BaseDescriptor.instance_init"), diff --git a/docs/ctapipe_api/io/index.rst b/docs/ctapipe_api/io/index.rst index 715f3c2c05d..d2987d5dbf0 100644 --- a/docs/ctapipe_api/io/index.rst +++ b/docs/ctapipe_api/io/index.rst @@ -165,6 +165,9 @@ Reference/API .. automodapi:: ctapipe.io.tableio :no-inheritance-diagram: +.. automodapi:: ctapipe.io.tableloader + :no-inheritance-diagram: + .. automodapi:: ctapipe.io.hdf5tableio :no-inheritance-diagram: diff --git a/docs/ctapipe_api/tools/index.rst b/docs/ctapipe_api/tools/index.rst index 6dcec106cac..a399dbbd18d 100644 --- a/docs/ctapipe_api/tools/index.rst +++ b/docs/ctapipe_api/tools/index.rst @@ -83,3 +83,6 @@ Reference/API .. automodapi:: ctapipe.tools :no-inheritance-diagram: + +.. automodapi:: ctapipe.tools.process + :no-inheritance-diagram: