From bd7d777cc19206faf029697034628606924116a5 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 31 Oct 2023 17:57:25 +0100 Subject: [PATCH] Implement CEP-2: new event container structure --- ctapipe/calib/camera/calibrator.py | 60 +-- ctapipe/calib/camera/flatfield.py | 28 +- ctapipe/calib/camera/pedestals.py | 28 +- ctapipe/calib/camera/tests/test_calibrator.py | 143 ++++--- ctapipe/calib/camera/tests/test_flatfield.py | 44 +- ctapipe/calib/camera/tests/test_pedestals.py | 22 +- ctapipe/containers.py | 390 ++++++++---------- ctapipe/image/extractor.py | 42 +- ctapipe/image/image_processor.py | 85 ++-- ctapipe/image/muon/processor.py | 33 +- .../image/muon/tests/test_intensity_fit.py | 2 +- ctapipe/image/muon/tests/test_processor.py | 6 +- ctapipe/image/reducer.py | 4 +- ctapipe/image/tests/test_extractor.py | 13 +- ctapipe/image/tests/test_image_processor.py | 50 +-- .../tests/test_sliding_window_correction.py | 4 +- ctapipe/instrument/tests/test_trigger.py | 77 ++-- ctapipe/instrument/trigger.py | 25 +- ctapipe/io/datawriter.py | 373 ++++++++--------- ctapipe/io/eventsource.py | 8 +- ctapipe/io/hdf5eventsource.py | 234 ++++++----- ctapipe/io/hdf5merger.py | 30 +- ctapipe/io/simteleventsource.py | 128 +++--- ctapipe/io/tableloader.py | 34 +- ctapipe/io/tests/test_datawriter.py | 27 +- ctapipe/io/tests/test_event_source.py | 4 +- ctapipe/io/tests/test_hdf5.py | 14 +- ctapipe/io/tests/test_hdf5eventsource.py | 71 ++-- ctapipe/io/tests/test_simteleventsource.py | 81 ++-- ctapipe/io/tests/test_table_loader.py | 5 +- ctapipe/io/tests/test_toysource.py | 3 +- ctapipe/io/toymodel.py | 15 +- ctapipe/reco/hillas_intersection.py | 10 +- ctapipe/reco/hillas_reconstructor.py | 17 +- ctapipe/reco/preprocessing.py | 15 +- ctapipe/reco/reconstructor.py | 29 +- ctapipe/reco/shower_processor.py | 6 +- ctapipe/reco/sklearn.py | 34 +- ctapipe/reco/stereo_combination.py | 37 +- .../reco/tests/test_HillasReconstructor.py | 38 +- .../reco/tests/test_hillas_intersection.py | 13 +- ctapipe/reco/tests/test_preprocessing.py | 10 +- .../reco/tests/test_reconstruction_methods.py | 8 +- ctapipe/reco/tests/test_shower_processor.py | 4 +- ctapipe/reco/tests/test_stereo_combination.py | 28 +- ctapipe/tools/apply_models.py | 9 +- ctapipe/tools/display_dl1.py | 17 +- ctapipe/tools/tests/test_apply_models.py | 12 +- ctapipe/tools/tests/test_process.py | 16 +- ctapipe/tools/tests/test_process_ml.py | 4 +- ctapipe/utils/event_type_filter.py | 2 +- ctapipe/utils/tests/test_event_filter.py | 20 +- ctapipe/visualization/mpl_array.py | 2 +- ctapipe/visualization/tests/test_bokeh.py | 20 +- ctapipe/visualization/tests/test_mpl.py | 33 +- docs/api-reference/containers/index.rst | 2 +- docs/api-reference/io/index.rst | 2 +- docs/user-guide/data_models/dl1.rst | 18 +- test_plugin/ctapipe_test_plugin/__init__.py | 6 +- 59 files changed, 1281 insertions(+), 1214 deletions(-) diff --git a/ctapipe/calib/camera/calibrator.py b/ctapipe/calib/camera/calibrator.py index 0e29f7f228c..445ae8ce9df 100644 --- a/ctapipe/calib/camera/calibrator.py +++ b/ctapipe/calib/camera/calibrator.py @@ -9,7 +9,12 @@ import numpy as np from numba import float32, float64, guvectorize, int64 -from ctapipe.containers import DL0CameraContainer, DL1CameraContainer, PixelStatus +from ctapipe.containers import ( + DL0TelescopeContainer, + DL1TelescopeContainer, + PixelStatus, + TelescopeEventContainer, +) from ctapipe.core import TelescopeComponent from ctapipe.core.traits import ( BoolTelescopeParameter, @@ -153,8 +158,8 @@ def __init__( parent=self, ) - def _check_r1_empty(self, waveforms): - if waveforms is None: + def _check_r1_empty(self, r1): + if r1 is None or r1.waveform is None: if not self._r1_empty_warn: warnings.warn( "Encountered an event with no R1 data. " @@ -165,8 +170,8 @@ def _check_r1_empty(self, waveforms): else: return False - def _check_dl0_empty(self, waveforms): - if waveforms is None: + def _check_dl0_empty(self, dl0): + if dl0 is None or dl0.waveform is None: if not self._dl0_empty_warn: warnings.warn( "Encountered an event with no DL0 data. " @@ -177,12 +182,13 @@ def _check_dl0_empty(self, waveforms): else: return False - def _calibrate_dl0(self, event, tel_id): - r1 = event.r1.tel[tel_id] - - if self._check_r1_empty(r1.waveform): + def r1_to_dl0(self, tel_event: TelescopeEventContainer): + if self._check_r1_empty(tel_event.r1): return + tel_id = tel_event.index.tel_id + r1 = tel_event.r1 + signal_pixels = self.data_volume_reducer( r1.waveform, tel_id=tel_id, @@ -198,7 +204,8 @@ def _calibrate_dl0(self, event, tel_id): # unset dvr bits for removed pixels dl0_pixel_status[~signal_pixels] &= ~np.uint8(PixelStatus.DVR_STATUS) - event.dl0.tel[tel_id] = DL0CameraContainer( + # FIXME: trigger information? + tel_event.dl0 = DL0TelescopeContainer( event_type=r1.event_type, event_time=r1.event_time, waveform=dl0_waveform, @@ -208,22 +215,24 @@ def _calibrate_dl0(self, event, tel_id): calibration_monitoring_id=r1.calibration_monitoring_id, ) - def _calibrate_dl1(self, event, tel_id): - waveforms = event.dl0.tel[tel_id].waveform - if self._check_dl0_empty(waveforms): + def dl0_to_dl1(self, tel_event: TelescopeEventContainer): + if self._check_dl0_empty(tel_event.dl0): return + tel_id = tel_event.index.tel_id + dl0 = tel_event.dl0 + waveforms = dl0.waveform n_pixels, n_samples = waveforms.shape - selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel + selected_gain_channel = dl0.selected_gain_channel broken_pixels = _get_invalid_pixels( n_pixels, - event.mon.tel[tel_id].pixel_status, + tel_event.mon.pixel_status, selected_gain_channel, ) - dl1_calib = event.calibration.tel[tel_id].dl1 - time_shift = event.calibration.tel[tel_id].dl1.time_shift + dl1_calib = tel_event.calibration.dl1 + time_shift = dl1_calib.time_shift readout = self.subarray.tel[tel_id].camera.readout # subtract any remaining pedestal before extraction @@ -239,7 +248,7 @@ def _calibrate_dl1(self, event, tel_id): # - Read into dl1 container directly? # - Don't do anything if dl1 container already filled # - Update on SST review decision - dl1 = DL1CameraContainer( + dl1 = DL1TelescopeContainer( image=waveforms[..., 0].astype(np.float32), peak_time=np.zeros(n_pixels, dtype=np.float32), is_valid=True, @@ -283,7 +292,7 @@ def _calibrate_dl1(self, event, tel_id): ) # store the results in the event structure - event.dl1.tel[tel_id] = dl1 + tel_event.dl1 = dl1 def __call__(self, event): """ @@ -294,13 +303,14 @@ def __call__(self, event): Parameters ---------- event : container - A `~ctapipe.containers.ArrayEventContainer` event container + A `~ctapipe.containers.SubarrayEventContainer` event container """ - # TODO: How to handle different calibrations depending on tel_id? - tel = event.r1.tel or event.dl0.tel or event.dl1.tel - for tel_id in tel.keys(): - self._calibrate_dl0(event, tel_id) - self._calibrate_dl1(event, tel_id) + for tel_event in event.tel.values(): + self.calibrate_tel_event(tel_event) + + def calibrate_tel_event(self, tel_event): + self.r1_to_dl0(tel_event) + self.dl0_to_dl1(tel_event) def shift_waveforms(waveforms, time_shift_samples): diff --git a/ctapipe/calib/camera/flatfield.py b/ctapipe/calib/camera/flatfield.py index 2a6c178f5ea..d461e69a04a 100644 --- a/ctapipe/calib/camera/flatfield.py +++ b/ctapipe/calib/camera/flatfield.py @@ -7,7 +7,7 @@ import numpy as np from astropy import units as u -from ctapipe.containers import DL1CameraContainer +from ctapipe.containers import DL1TelescopeContainer from ctapipe.core import Component from ctapipe.core.traits import Int, List, Unicode from ctapipe.image.extractor import ImageExtractor @@ -108,7 +108,7 @@ def calculate_relative_gain(self, event): Parameters ---------- - event: ctapipe.containers.ArrayEventContainer + event: ctapipe.containers.SubarrayEventContainer Returns: True if the mon.tel[tel_id].flatfield is updated, False otherwise @@ -169,7 +169,7 @@ def __init__(self, **kwargs): self.arrival_times = None # arrival time per event in sample self.sample_masked_pixels = None # masked pixels per event in sample - def _extract_charge(self, event) -> DL1CameraContainer: + def _extract_charge(self, event) -> DL1TelescopeContainer: """ Extract the charge and the time from a calibration event @@ -182,11 +182,11 @@ def _extract_charge(self, event) -> DL1CameraContainer: DL1CameraContainer """ - waveforms = event.r1.tel[self.tel_id].waveform - selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel + waveforms = event.tel[self.tel_id].r1.waveform + selected_gain_channel = event.tel[self.tel_id].r1.selected_gain_channel broken_pixels = _get_invalid_pixels( n_pixels=waveforms.shape[-2], - pixel_status=event.mon.tel[self.tel_id].pixel_status, + pixel_status=event.tel[self.tel_id].mon.pixel_status, selected_gain_channel=selected_gain_channel, ) # Extract charge and time @@ -195,7 +195,7 @@ def _extract_charge(self, event) -> DL1CameraContainer: waveforms, self.tel_id, selected_gain_channel, broken_pixels ) else: - return DL1CameraContainer(image=0, peak_pos=0, is_valid=False) + return DL1TelescopeContainer(image=0, peak_pos=0, is_valid=False) def calculate_relative_gain(self, event): """ @@ -209,23 +209,23 @@ def calculate_relative_gain(self, event): """ # initialize the np array at each cycle - waveform = event.r1.tel[self.tel_id].waveform - container = event.mon.tel[self.tel_id].flatfield + waveform = event.tel[self.tel_id].r1.waveform + container = event.tel[self.tel_id].mon.flatfield # re-initialize counter if self.n_events_seen == self.sample_size: self.n_events_seen = 0 # real data - trigger_time = event.trigger.time + trigger_time = event.dl0.trigger.time if event.meta["origin"] != "hessio": hardware_or_pedestal_mask = np.logical_or( - event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels, - event.mon.tel[self.tel_id].pixel_status.pedestal_failing_pixels, + event.tel[self.tel_id].mon.pixel_status.hardware_failing_pixels, + event.tel[self.tel_id].mon.pixel_status.pedestal_failing_pixels, ) pixel_mask = np.logical_or( hardware_or_pedestal_mask, - event.mon.tel[self.tel_id].pixel_status.flatfield_failing_pixels, + event.tel[self.tel_id].mon.pixel_status.flatfield_failing_pixels, ) else: # patches for MC data @@ -237,7 +237,7 @@ def calculate_relative_gain(self, event): # extract the charge of the event and # the peak position (assumed as time for the moment) - dl1: DL1CameraContainer = self._extract_charge(event) + dl1: DL1TelescopeContainer = self._extract_charge(event) if not dl1.is_valid: return False diff --git a/ctapipe/calib/camera/pedestals.py b/ctapipe/calib/camera/pedestals.py index d8a5576e7b8..4c7ff3385d2 100644 --- a/ctapipe/calib/camera/pedestals.py +++ b/ctapipe/calib/camera/pedestals.py @@ -7,7 +7,7 @@ import numpy as np from astropy import units as u -from ctapipe.containers import DL1CameraContainer +from ctapipe.containers import DL1TelescopeContainer from ctapipe.core import Component from ctapipe.core.traits import Int, List, Unicode from ctapipe.image.extractor import ImageExtractor @@ -134,7 +134,7 @@ def calculate_pedestals(self, event): Parameters ---------- - event: ctapipe.containers.ArrayEventContainer + event: ctapipe.containers.SubarrayEventContainer Returns: True if the mon.tel[tel_id].pedestal is updated, False otherwise @@ -197,24 +197,24 @@ def __init__(self, **kwargs): self.charges = None # charge per event in sample self.sample_masked_pixels = None # pixels tp be masked per event in sample - def _extract_charge(self, event) -> DL1CameraContainer: + def _extract_charge(self, event) -> DL1TelescopeContainer: """ Extract the charge and the time from a pedestal event Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer general event container Returns ------- - DL1CameraContainer + DL1TelescopeContainer """ - waveforms = event.r1.tel[self.tel_id].waveform - selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel + waveforms = event.tel[self.tel_id].r1.waveform + selected_gain_channel = event.tel[self.tel_id].r1.selected_gain_channel broken_pixels = _get_invalid_pixels( n_pixels=waveforms.shape[-2], - pixel_status=event.mon.tel[self.tel_id].pixel_status, + pixel_status=event.tel[self.tel_id].mon.pixel_status, selected_gain_channel=selected_gain_channel, ) @@ -224,7 +224,7 @@ def _extract_charge(self, event) -> DL1CameraContainer: waveforms, self.tel_id, selected_gain_channel, broken_pixels ) else: - return DL1CameraContainer(image=0, peak_pos=0, is_valid=False) + return DL1TelescopeContainer(image=0, peak_pos=0, is_valid=False) def calculate_pedestals(self, event): """ @@ -238,17 +238,17 @@ def calculate_pedestals(self, event): """ # initialize the np array at each cycle - waveform = event.r1.tel[self.tel_id].waveform - container = event.mon.tel[self.tel_id].pedestal + waveform = event.tel[self.tel_id].r1.waveform + container = event.tel[self.tel_id].mon.pedestal # re-initialize counter if self.n_events_seen == self.sample_size: self.n_events_seen = 0 # real data - trigger_time = event.trigger.time + trigger_time = event.dl0.trigger.time if event.meta["origin"] != "hessio": - pixel_mask = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels + pixel_mask = event.tel[self.tel_id].mon.pixel_status.hardware_failing_pixels else: # patches for MC data pixel_mask = np.zeros(waveform.shape[1], dtype=bool) @@ -258,7 +258,7 @@ def calculate_pedestals(self, event): # extract the charge of the event and # the peak position (assumed as time for the moment) - dl1: DL1CameraContainer = self._extract_charge(event) + dl1: DL1TelescopeContainer = self._extract_charge(event) if not dl1.is_valid: return False diff --git a/ctapipe/calib/camera/tests/test_calibrator.py b/ctapipe/calib/camera/tests/test_calibrator.py index b70c6a4aee1..dcaf10e9cd2 100644 --- a/ctapipe/calib/camera/tests/test_calibrator.py +++ b/ctapipe/calib/camera/tests/test_calibrator.py @@ -10,7 +10,14 @@ from traitlets.config import Config from ctapipe.calib.camera.calibrator import CameraCalibrator -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import ( + DL0TelescopeContainer, + DL1TelescopeContainer, + R1TelescopeContainer, + SubarrayEventContainer, + TelescopeEventContainer, + TelescopeEventIndexContainer, +) from ctapipe.image.extractor import ( FullWaveformSum, GlobalPeakWindowSum, @@ -21,11 +28,11 @@ def test_camera_calibrator(example_event, example_subarray): - tel_id = list(example_event.r0.tel)[0] + tel_event = next(iter(example_event.tel.values())) calibrator = CameraCalibrator(subarray=example_subarray) calibrator(example_event) - image = example_event.dl1.tel[tel_id].image - peak_time = example_event.dl1.tel[tel_id].peak_time + image = tel_event.dl1.image + peak_time = tel_event.dl1.peak_time assert image is not None assert peak_time is not None assert image.shape == (1764,) @@ -99,53 +106,68 @@ def test_config(example_subarray): def test_check_r1_empty(example_event, example_subarray): calibrator = CameraCalibrator(subarray=example_subarray) - tel_id = list(example_event.r0.tel)[0] - waveform = example_event.r1.tel[tel_id].waveform.copy() + tel_id, tel_event = next(iter(example_event.tel.items())) + waveform = tel_event.r1.waveform.copy() with pytest.warns(UserWarning): - example_event.r1.tel[tel_id].waveform = None - calibrator._calibrate_dl0(example_event, tel_id) - assert example_event.dl0.tel[tel_id].waveform is None + tel_event.r1.waveform = None + calibrator.r1_to_dl0(tel_event) + assert tel_event.dl0.waveform is None assert calibrator._check_r1_empty(None) is True - assert calibrator._check_r1_empty(waveform) is False + assert calibrator._check_r1_empty(R1TelescopeContainer(waveform=None)) is True + assert calibrator._check_r1_empty(R1TelescopeContainer(waveform=waveform)) is False calibrator = CameraCalibrator( subarray=example_subarray, image_extractor=FullWaveformSum(subarray=example_subarray), ) - event = ArrayEventContainer() - event.dl0.tel[tel_id].waveform = np.full((2048, 128), 2) + event = SubarrayEventContainer() + event.tel[tel_id] = TelescopeEventContainer( + index=TelescopeEventIndexContainer(obs_id=1, event_id=1, tel_id=tel_id), + dl0=DL0TelescopeContainer(waveform=np.full((2048, 128), 2)), + ) with pytest.warns(UserWarning): calibrator(event) - assert (event.dl0.tel[tel_id].waveform == 2).all() - assert (event.dl1.tel[tel_id].image == 2 * 128).all() + assert (event.tel[tel_id].dl0.waveform == 2).all() + assert (event.tel[tel_id].dl1.image == 2 * 128).all() def test_check_dl0_empty(example_event, example_subarray): calibrator = CameraCalibrator(subarray=example_subarray) - tel_id = list(example_event.r0.tel)[0] - calibrator._calibrate_dl0(example_event, tel_id) - waveform = example_event.dl0.tel[tel_id].waveform.copy() + tel_id, tel_event = next(iter(example_event.tel.items())) + + calibrator.r1_to_dl0(tel_event) + waveform = tel_event.dl0.waveform.copy() + with pytest.warns(UserWarning): - example_event.dl0.tel[tel_id].waveform = None - calibrator._calibrate_dl1(example_event, tel_id) - assert example_event.dl1.tel[tel_id].image is None + tel_event.dl0.waveform = None + calibrator.dl0_to_dl1(tel_event) + assert tel_event.dl1.image is None assert calibrator._check_dl0_empty(None) is True - assert calibrator._check_dl0_empty(waveform) is False + assert calibrator._check_dl0_empty(DL0TelescopeContainer(waveform=None)) is True + assert ( + calibrator._check_dl0_empty(DL0TelescopeContainer(waveform=waveform)) is False + ) calibrator = CameraCalibrator(subarray=example_subarray) - event = ArrayEventContainer() - event.dl1.tel[tel_id].image = np.full(2048, 2) + event = SubarrayEventContainer() + tel_event = TelescopeEventContainer( + index=TelescopeEventIndexContainer(obs_id=1, event_id=1, tel_id=tel_id), + dl1=DL1TelescopeContainer(image=np.full(2048, 2)), + ) + event.tel[tel_id] = tel_event with pytest.warns(UserWarning): calibrator(event) - assert (event.dl1.tel[tel_id].image == 2).all() + assert (tel_event.dl1.image == 2).all() def test_dl1_charge_calib(example_subarray): # copy because we mutate the camera, should not affect other tests + rng = np.random.default_rng(1) + tel_id = 1 subarray = deepcopy(example_subarray) - camera = subarray.tel[1].camera + camera = subarray.tel[tel_id].camera # test with a sampling_rate different than 1 to # test if we handle time vs. slices correctly sampling_rate = 2 @@ -155,60 +177,63 @@ def test_dl1_charge_calib(example_subarray): n_samples = 96 mid = n_samples // 2 pulse_sigma = 6 - random = np.random.default_rng(1) x = np.arange(n_samples) # Randomize times and create pulses - time_offset = random.uniform(-10, +10, n_pixels) + time_offset = rng.uniform(-10, +10, n_pixels) y = norm.pdf(x, mid + time_offset[:, np.newaxis], pulse_sigma).astype("float32") camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma)[np.newaxis, :] camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate # Define absolute calibration coefficients - absolute = random.uniform(100, 1000, n_pixels).astype("float32") + absolute = rng.uniform(100, 1000, n_pixels).astype("float32") y *= absolute[:, np.newaxis] # Define relative coefficients - relative = random.normal(1, 0.01, n_pixels) + relative = rng.normal(1, 0.01, n_pixels) y /= relative[:, np.newaxis] # Define pedestal - pedestal = random.uniform(-4, 4, n_pixels) + pedestal = rng.uniform(-4, 4, n_pixels) y += pedestal[:, np.newaxis] - event = ArrayEventContainer() - tel_id = list(subarray.tel.keys())[0] - event.dl0.tel[tel_id].waveform = y - event.dl0.tel[tel_id].selected_gain_channel = np.zeros(len(y), dtype=int) - event.r1.tel[tel_id].selected_gain_channel = np.zeros(len(y), dtype=int) + event = SubarrayEventContainer() + event.tel[tel_id] = TelescopeEventContainer( + index=TelescopeEventIndexContainer(obs_id=1, event_id=1, tel_id=tel_id), + dl0=DL0TelescopeContainer( + waveform=y, + selected_gain_channel=np.zeros(len(y), dtype=int), + ), + ) # Test default calibrator = CameraCalibrator( - subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) + subarray=subarray, + image_extractor=FullWaveformSum(subarray=subarray), ) calibrator(event) - np.testing.assert_allclose(event.dl1.tel[tel_id].image, y.sum(1), rtol=1e-4) + np.testing.assert_allclose(event.tel[tel_id].dl1.image, y.sum(1), rtol=1e-4) - event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal - event.calibration.tel[tel_id].dl1.absolute_factor = absolute - event.calibration.tel[tel_id].dl1.relative_factor = relative + event.tel[tel_id].calibration.dl1.pedestal_offset = pedestal + event.tel[tel_id].calibration.dl1.absolute_factor = absolute + event.tel[tel_id].calibration.dl1.relative_factor = relative # Test without timing corrections calibrator(event) - dl1 = event.dl1.tel[tel_id] + dl1 = event.tel[tel_id].dl1 np.testing.assert_allclose(dl1.image, 1, rtol=1e-5) expected_peak_time = (mid + time_offset) / sampling_rate np.testing.assert_allclose(dl1.peak_time, expected_peak_time, rtol=1e-5) # test with timing corrections - event.calibration.tel[tel_id].dl1.time_shift = time_offset / sampling_rate + event.tel[tel_id].calibration.dl1.time_shift = time_offset / sampling_rate calibrator(event) # more rtol since shifting might lead to reduced integral - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-5) + np.testing.assert_allclose(event.tel[tel_id].dl1.image, 1, rtol=1e-5) np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + event.tel[tel_id].dl1.peak_time, mid / sampling_rate, atol=1 ) # test not applying time shifts @@ -217,9 +242,9 @@ def test_dl1_charge_calib(example_subarray): calibrator.apply_waveform_time_shift = False calibrator(event) - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-4) + np.testing.assert_allclose(event.tel[tel_id].dl1.image, 1, rtol=1e-4) np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, expected_peak_time, atol=1 + event.tel[tel_id].dl1.peak_time, expected_peak_time, atol=1 ) # We now use GlobalPeakWindowSum to see the effect of missing charge @@ -232,9 +257,9 @@ def test_dl1_charge_calib(example_subarray): calibrator(event) # test with timing corrections, should work # higher rtol because we cannot shift perfectly - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=0.01) + np.testing.assert_allclose(event.tel[tel_id].dl1.image, 1, rtol=0.01) np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + event.tel[tel_id].dl1.peak_time, mid / sampling_rate, atol=1 ) # test deactivating timing corrections @@ -243,8 +268,8 @@ def test_dl1_charge_calib(example_subarray): # make sure we chose an example where the time shifts matter # charges should be quite off due to summing around global shift - assert not np.allclose(event.dl1.tel[tel_id].image, 1, rtol=0.1) - assert not np.allclose(event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1) + assert not np.allclose(event.tel[tel_id].dl1.image, 1, rtol=0.1) + assert not np.allclose(event.tel[tel_id].dl1.peak_time, mid / sampling_rate, atol=1) def test_shift_waveforms(): @@ -283,22 +308,22 @@ def test_invalid_pixels(example_event, example_subarray): ) # going to modify this event = deepcopy(example_event) - tel_id = list(event.r0.tel)[0] + tel_id, tel_event = next(iter(event.tel.items())) camera = example_subarray.tel[tel_id].camera sampling_rate = camera.readout.sampling_rate.to_value(u.GHz) - event.mon.tel[tel_id].pixel_status.flatfield_failing_pixels[:, 0] = True - event.r1.tel[tel_id].waveform.fill(0.0) - event.r1.tel[tel_id].waveform[1:, 20] = 1.0 - event.r1.tel[tel_id].waveform[0, 10] = 9999 + tel_event.mon.pixel_status.flatfield_failing_pixels[:, 0] = True + tel_event.r1.waveform.fill(0.0) + tel_event.r1.waveform[1:, 20] = 1.0 + tel_event.r1.waveform[0, 10] = 9999 calibrator = CameraCalibrator( subarray=example_subarray, config=config, ) calibrator(event) - assert np.all(event.dl1.tel[tel_id].image == 1.0) - assert np.all(event.dl1.tel[tel_id].peak_time == 20.0 / sampling_rate) + assert np.all(tel_event.dl1.image == 1.0) + assert np.all(tel_event.dl1.peak_time == 20.0 / sampling_rate) # test we can set the invalid pixel handler to None config.CameraCalibrator.invalid_pixel_handler_type = None @@ -307,5 +332,5 @@ def test_invalid_pixels(example_event, example_subarray): config=config, ) calibrator(event) - assert event.dl1.tel[tel_id].image[0] == 9999 - assert event.dl1.tel[tel_id].peak_time[0] == 10.0 / sampling_rate + assert tel_event.dl1.image[0] == 9999 + assert tel_event.dl1.peak_time[0] == 10.0 / sampling_rate diff --git a/ctapipe/calib/camera/tests/test_flatfield.py b/ctapipe/calib/camera/tests/test_flatfield.py index ec0eeaf028c..a476d653fd3 100644 --- a/ctapipe/calib/camera/tests/test_flatfield.py +++ b/ctapipe/calib/camera/tests/test_flatfield.py @@ -6,7 +6,7 @@ from traitlets.config import Config from ctapipe.calib.camera.flatfield import FlasherFlatFieldCalculator -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import SubarrayEventContainer from ctapipe.instrument import SubarrayDescription @@ -38,49 +38,49 @@ def test_flasherflatfieldcalculator(prod5_sst, reference_location): config=config, ) # create one event - data = ArrayEventContainer() - data.meta["origin"] = "test" - data.trigger.time = Time.now() + event = SubarrayEventContainer() + event.meta["origin"] = "test" + event.dl0.trigger.time = Time.now() # initialize mon and r1 data - data.mon.tel[tel_id].pixel_status.hardware_failing_pixels = np.zeros( + event.tel[tel_id].mon.pixel_status.hardware_failing_pixels = np.zeros( (n_gain, n_pixels), dtype=bool ) - data.mon.tel[tel_id].pixel_status.pedestal_failing_pixels = np.zeros( + event.tel[tel_id].mon.pixel_status.pedestal_failing_pixels = np.zeros( (n_gain, n_pixels), dtype=bool ) - data.mon.tel[tel_id].pixel_status.flatfield_failing_pixels = np.zeros( + event.tel[tel_id].mon.pixel_status.flatfield_failing_pixels = np.zeros( (n_gain, n_pixels), dtype=bool ) - data.r1.tel[tel_id].waveform = np.zeros((n_gain, n_pixels, 40)) - data.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8) + event.tel[tel_id].r1.waveform = np.zeros((n_gain, n_pixels, 40)) + event.tel[tel_id].r1.selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8) # flat-field signal put == delta function of height ff_level at sample 20 - data.r1.tel[tel_id].waveform[:, :, 20] = ff_level - print(data.r1.tel[tel_id].waveform[0, 0, 20]) + event.tel[tel_id].r1.waveform[:, :, 20] = ff_level + print(event.tel[tel_id].r1.waveform[0, 0, 20]) # First test: good event while ff_calculator.n_events_seen < n_events: - if ff_calculator.calculate_relative_gain(data): - assert data.mon.tel[tel_id].flatfield + if ff_calculator.calculate_relative_gain(event): + assert event.tel[tel_id].mon.flatfield - print(data.mon.tel[tel_id].flatfield) - assert np.mean(data.mon.tel[tel_id].flatfield.charge_median) == ff_level - assert np.mean(data.mon.tel[tel_id].flatfield.relative_gain_median) == 1 - assert np.mean(data.mon.tel[tel_id].flatfield.relative_gain_std) == 0 + print(event.tel[tel_id].mon.flatfield) + assert np.mean(event.tel[tel_id].mon.flatfield.charge_median) == ff_level + assert np.mean(event.tel[tel_id].mon.flatfield.relative_gain_median) == 1 + assert np.mean(event.tel[tel_id].mon.flatfield.relative_gain_std) == 0 # Second test: introduce some failing pixels failing_pixels_id = np.array([10, 20, 30, 40]) - data.r1.tel[tel_id].waveform[:, failing_pixels_id, :] = 0 - data.mon.tel[tel_id].pixel_status.pedestal_failing_pixels[ + event.tel[tel_id].r1.waveform[:, failing_pixels_id, :] = 0 + event.tel[tel_id].mon.pixel_status.pedestal_failing_pixels[ :, failing_pixels_id ] = True while ff_calculator.n_events_seen < n_events: - if ff_calculator.calculate_relative_gain(data): + if ff_calculator.calculate_relative_gain(event): # working pixel have good gain - assert data.mon.tel[tel_id].flatfield.relative_gain_median[0, 0] == 1 + assert event.tel[tel_id].mon.flatfield.relative_gain_median[0, 0] == 1 # bad pixels do non influence the gain - assert np.mean(data.mon.tel[tel_id].flatfield.relative_gain_std) == 0 + assert np.mean(event.tel[tel_id].mon.flatfield.relative_gain_std) == 0 diff --git a/ctapipe/calib/camera/tests/test_pedestals.py b/ctapipe/calib/camera/tests/test_pedestals.py index 7f17c9f1703..ea56cff7775 100644 --- a/ctapipe/calib/camera/tests/test_pedestals.py +++ b/ctapipe/calib/camera/tests/test_pedestals.py @@ -11,7 +11,7 @@ PedestalIntegrator, calc_pedestals_from_traces, ) -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import SubarrayEventContainer from ctapipe.instrument import SubarrayDescription @@ -40,24 +40,24 @@ def test_pedestal_integrator(prod5_sst, reference_location): tel_id=tel_id, ) # create one event - data = ArrayEventContainer() - data.meta["origin"] = "test" - data.trigger.time = Time.now() + event = SubarrayEventContainer() + event.meta["origin"] = "test" + event.dl0.trigger.time = Time.now() # fill the values necessary for the pedestal calculation - data.mon.tel[tel_id].pixel_status.hardware_failing_pixels = np.zeros( + event.tel[tel_id].mon.pixel_status.hardware_failing_pixels = np.zeros( (n_gain, n_pixels), dtype=bool ) - data.r1.tel[tel_id].waveform = np.full((2, n_pixels, 40), ped_level) - data.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8) + event.tel[tel_id].r1.waveform = np.full((2, n_pixels, 40), ped_level) + event.tel[tel_id].r1.selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8) while ped_calculator.n_events_seen < n_events: - if ped_calculator.calculate_pedestals(data): - assert data.mon.tel[tel_id].pedestal - assert np.mean(data.mon.tel[tel_id].pedestal.charge_median) == ( + if ped_calculator.calculate_pedestals(event): + assert event.tel[tel_id].mon.pedestal + assert np.mean(event.tel[tel_id].mon.pedestal.charge_median) == ( ped_calculator.extractor.window_width.tel[0] * ped_level ) - assert np.mean(data.mon.tel[tel_id].pedestal.charge_std) == 0 + assert np.mean(event.tel[tel_id].mon.pedestal.charge_std) == 0 def test_calc_pedestals_from_traces(): diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 3fbb88a762c..99266766072 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -12,56 +12,52 @@ from .core import Container, Field, Map __all__ = [ - "ArrayEventContainer", + "BaseHillasParametersContainer", + "BaseTimingParametersContainer", + "CameraHillasParametersContainer", + "CameraTimingParametersContainer", "ConcentrationContainer", - "DL0CameraContainer", - "DL0Container", + "CoreParametersContainer", + "DL0TelescopeContainer", "DL1CameraCalibrationContainer", - "DL1CameraContainer", - "DL1Container", - "DL2Container", - "EventCalibrationContainer", - "EventCameraCalibrationContainer", - "EventIndexContainer", + "DL1TelescopeContainer", + "DL2SubarrayContainer", + "DispContainer", "EventType", "FlatFieldContainer", "HillasParametersContainer", - "CoreParametersContainer", "ImageParametersContainer", + "IntensityStatisticsContainer", "LeakageContainer", - "MonitoringCameraContainer", - "MonitoringContainer", "MorphologyContainer", - "BaseHillasParametersContainer", - "CameraHillasParametersContainer", - "CameraTimingParametersContainer", + "ObservationBlockContainer", + "ObservationBlockState", + "ObservingMode", "ParticleClassificationContainer", + "PeakTimeStatisticsContainer", "PedestalContainer", "PixelStatusContainer", - "R0CameraContainer", - "R0Container", - "R1CameraContainer", - "R1Container", - "ReconstructedContainer", + "R0TelescopeContainer", + "R1TelescopeContainer", + "ReconstructedShowerContainer", "ReconstructedEnergyContainer", "ReconstructedGeometryContainer", - "DispContainer", - "SimulatedCameraContainer", + "SchedulingBlockContainer", "SimulatedShowerContainer", "SimulatedShowerDistribution", "SimulationConfigContainer", - "TelEventIndexContainer", - "BaseTimingParametersContainer", + "SimulationSubarrayContainer", + "SimulationTelescopeContainer", + "StatisticsContainer", + "SubarrayEventContainer", + "SubarrayEventIndexContainer", + "SubarrayTriggerContainer", + "TelescopeCalibrationContainer", + "TelescopeEventContainer", + "TelescopeEventIndexContainer", + "TelescopeMonitoringContainer", "TimingParametersContainer", - "TriggerContainer", "WaveformCalibrationContainer", - "StatisticsContainer", - "IntensityStatisticsContainer", - "PeakTimeStatisticsContainer", - "SchedulingBlockContainer", - "ObservationBlockContainer", - "ObservingMode", - "ObservationBlockState", ] @@ -141,30 +137,6 @@ class CoordinateFrameType(enum.Enum): GALACTIC = 2 -class EventType(enum.Enum): - """Enum of EventTypes as defined in the CTA Data Model [cta_r1event]_""" - - # calibrations are 0-15 - FLATFIELD = 0 - SINGLE_PE = 1 - SKY_PEDESTAL = 2 - DARK_PEDESTAL = 3 - ELECTRONIC_PEDESTAL = 4 - OTHER_CALIBRATION = 15 - - #: For mono-telescope triggers (not used in MC) - MUON = 16 - HARDWARE_STEREO = 17 - - #: ACADA (DAQ) software trigger - DAQ = 24 - - #: Standard Physics stereo trigger - SUBARRAY = 32 - - UNKNOWN = 255 - - class PixelStatus(enum.IntFlag): """ Pixel status information @@ -229,7 +201,31 @@ def is_invalid(pixel_status): return (pixel_status & gain_bits) == 0 -class EventIndexContainer(Container): +class EventType(enum.Enum): + """Enum of EventTypes as defined in the CTA R1 and DL0 Data Model [cta_r1event]_""" + + # calibrations are 0-15 + FLATFIELD = 0 + SINGLE_PE = 1 + SKY_PEDESTAL = 2 + DARK_PEDESTAL = 3 + ELECTRONIC_PEDESTAL = 4 + OTHER_CALIBRATION = 15 + + #: For mono-telescope triggers (not used in MC) + MUON = 16 + HARDWARE_STEREO = 17 + + #: ACADA (DAQ) software trigger + DAQ = 24 + + #: Standard Physics stereo trigger + SUBARRAY = 32 + + UNKNOWN = 255 + + +class SubarrayEventIndexContainer(Container): """index columns to include in event lists, common to all data levels""" default_prefix = "" # don't want to prefix these @@ -237,7 +233,7 @@ class EventIndexContainer(Container): event_id = event_id_field() -class TelEventIndexContainer(Container): +class TelescopeEventIndexContainer(Container): """ index columns to include in telescope-wise event lists, common to all data levels that have telescope-wise information @@ -249,6 +245,25 @@ class TelEventIndexContainer(Container): tel_id = tel_id_field() +class TelescopeTriggerContainer(Container): + default_prefix = "" + time = Field(NAN_TIME, description="Telescope trigger time") + n_trigger_pixels = Field( + -1, description="Number of trigger groups (sectors) listed" + ) + trigger_pixels = Field(None, description="pixels involved in the camera trigger") + event_type = Field(EventType.UNKNOWN, description="Event type") + + +class SubarrayTriggerContainer(Container): + default_prefix = "" + time = Field(NAN_TIME, description="central average time stamp") + tels_with_trigger = Field( + None, description="List of telescope ids that triggered the array event" + ) + event_type = Field(EventType.SUBARRAY, description="Event type") + + class BaseHillasParametersContainer(Container): """ Base container for hillas parameters to @@ -463,7 +478,7 @@ class ImageParametersContainer(Container): ) -class DL1CameraContainer(Container): +class DL1TelescopeContainer(Container): """ Storage of output of camera calibration e.g the final calibrated image in intensity units and the pulse time. @@ -503,15 +518,6 @@ class DL1CameraContainer(Container): ) -class DL1Container(Container): - """DL1 Calibrated Camera Images and associated data""" - - tel = Field( - default_factory=partial(Map, DL1CameraContainer), - description="map of tel_id to DL1CameraContainer", - ) - - class DL1CameraCalibrationContainer(Container): """ Storage of DL1 calibration parameters for the current event @@ -541,7 +547,7 @@ class DL1CameraCalibrationContainer(Container): ) -class R0CameraContainer(Container): +class R0TelescopeContainer(Container): """ Storage of raw data from a single telescope """ @@ -551,18 +557,7 @@ class R0CameraContainer(Container): ) -class R0Container(Container): - """ - Storage of a Merged Raw Data Event - """ - - tel = Field( - default_factory=partial(Map, R0CameraContainer), - description="map of tel_id to R0CameraContainer", - ) - - -class R1CameraContainer(Container): +class R1TelescopeContainer(Container): """ Storage of r1 calibrated data from a single telescope """ @@ -621,20 +616,9 @@ class R1CameraContainer(Container): ) -class R1Container(Container): - """ - Storage of a r1 calibrated Data Event - """ - - tel = Field( - default_factory=partial(Map, R1CameraContainer), - description="map of tel_id to R1CameraContainer", - ) - - -class DL0CameraContainer(Container): +class DL0TelescopeContainer(Container): """ - Storage of data volume reduced dl0 data from a single telescope. + DL0 data from a single telescope. See DL0 Data Model specification: https://redmine.cta-observatory.org/dmsf/files/17552/view @@ -643,6 +627,11 @@ class DL0CameraContainer(Container): event_type = Field(EventType.UNKNOWN, "type of event", type=EventType) event_time = Field(NAN_TIME, "event timestamp") + trigger = Field( + default_factory=TelescopeTriggerContainer, + description="telescope trigger information", + ) + waveform = Field( None, ( @@ -682,14 +671,14 @@ class DL0CameraContainer(Container): ) -class DL0Container(Container): +class DL0SubarrayContainer(Container): """ - Storage of a data volume reduced Event + DL0 data at subarray level """ - tel = Field( - default_factory=partial(Map, DL0CameraContainer), - description="map of tel_id to DL0CameraContainer", + trigger = Field( + default_factory=SubarrayTriggerContainer, + description="subarray trigger information", ) @@ -728,9 +717,9 @@ class SimulatedShowerContainer(Container): ) -class SimulatedCameraContainer(Container): +class SimulationTelescopeContainer(Container): """ - True images and parameters derived from them, analgous to the `DL1CameraContainer` + True images and parameters derived from them, analgous to the `DL1TelescopeContainer` but for simulated data. """ @@ -759,12 +748,11 @@ class SimulatedCameraContainer(Container): ) -class SimulatedEventContainer(Container): +class SimulationSubarrayContainer(Container): shower = Field( default_factory=SimulatedShowerContainer, description="True event information", ) - tel = Field(default_factory=partial(Map, SimulatedCameraContainer)) class SimulationConfigContainer(Container): @@ -851,28 +839,6 @@ class SimulationConfigContainer(Container): ) -class TelescopeTriggerContainer(Container): - default_prefix = "" - time = Field(NAN_TIME, description="Telescope trigger time") - n_trigger_pixels = Field( - -1, description="Number of trigger groups (sectors) listed" - ) - trigger_pixels = Field(None, description="pixels involved in the camera trigger") - - -class TriggerContainer(Container): - default_prefix = "" - time = Field(NAN_TIME, description="central average time stamp") - tels_with_trigger = Field( - None, description="List of telescope ids that triggered the array event" - ) - event_type = Field(EventType.SUBARRAY, description="Event type") - tel = Field( - default_factory=partial(Map, TelescopeTriggerContainer), - description="telescope-wise trigger information", - ) - - class ReconstructedGeometryContainer(Container): """ Standard output of algorithms reconstructing shower geometry @@ -994,16 +960,13 @@ class DispContainer(Container): is_valid = Field(False, "true if the predictions are valid") -class ReconstructedContainer(Container): - """Reconstructed shower info from multiple algorithms""" +class ReconstructedShowerContainer(Container): + """ + Reconstructed Shower Information - # Note: there is a reason why the hiererchy is - # `event.dl2.stereo.geometry[algorithm]` and not - # `event.dl2[algorithm].stereo.geometry` and that is because when writing - # the data, the former makes it easier to only write information that a - # particular reconstructor generates, e.g. only write the geometry in cases - # where energy is not yet computed. Some algorithms will compute all three, - # but most will compute only fill or two of these sub-Contaiers: + May contain reconstructions of multiple algorithms for + each property mapping algorithm name to containers. + """ geometry = Field( default_factory=partial(Map, ReconstructedGeometryContainer), @@ -1019,35 +982,26 @@ class ReconstructedContainer(Container): ) -class TelescopeReconstructedContainer(ReconstructedContainer): - """Telescope-wise reconstructed quantities""" +class DL2SubarrayContainer(ReconstructedShowerContainer): + """ + DL2 subarray-wise information (reconstructed air shower properties) + """ + + +class DL2TelescopeContainer(ReconstructedShowerContainer): + """DL2 telescope-wise quantities""" impact = Field( default_factory=partial(Map, TelescopeImpactParameterContainer), description="map of algorithm to impact parameter info", ) + disp = Field( default_factory=partial(Map, DispContainer), description="map of algorithm to reconstructed disp parameters", ) -class DL2Container(Container): - """Reconstructed Shower information for a given reconstruction algorithm, - including optionally both per-telescope mono reconstruction and per-shower - stereo reconstructions - """ - - tel = Field( - default_factory=partial(Map, TelescopeReconstructedContainer), - description="map of tel_id to single-telescope reconstruction (DL2a)", - ) - stereo = Field( - default_factory=ReconstructedContainer, - description="Stereo Shower reconstruction results", - ) - - class TelescopePointingContainer(Container): """ Container holding pointing information for a single telescope @@ -1060,18 +1014,14 @@ class TelescopePointingContainer(Container): altitude = Field(nan * u.rad, "Altitude", unit=u.rad) -class PointingContainer(Container): - tel = Field( - default_factory=partial(Map, TelescopePointingContainer), - description="Telescope pointing positions", - ) - array_azimuth = Field(nan * u.rad, "Array pointing azimuth", unit=u.rad) - array_altitude = Field(nan * u.rad, "Array pointing altitude", unit=u.rad) - array_ra = Field(nan * u.rad, "Array pointing right ascension", unit=u.rad) - array_dec = Field(nan * u.rad, "Array pointing declination", unit=u.rad) +class SubarrayPointingContainer(Container): + azimuth = Field(nan * u.rad, "Array pointing azimuth", unit=u.rad) + altitude = Field(nan * u.rad, "Array pointing altitude", unit=u.rad) + ra = Field(nan * u.rad, "Array pointing right ascension", unit=u.rad) + dec = Field(nan * u.rad, "Array pointing declination", unit=u.rad) -class EventCameraCalibrationContainer(Container): +class TelescopeCalibrationContainer(Container): """ Container for the calibration coefficients for the current event and camera """ @@ -1082,18 +1032,6 @@ class EventCameraCalibrationContainer(Container): ) -class EventCalibrationContainer(Container): - """ - Container for calibration coefficients for the current event - """ - - # create the camera container - tel = Field( - default_factory=partial(Map, EventCameraCalibrationContainer), - description="map of tel_id to EventCameraCalibrationContainer", - ) - - class MuonRingContainer(Container): """Container for the result of a ring fit in telescope frame""" @@ -1113,6 +1051,8 @@ class MuonRingContainer(Container): class MuonEfficiencyContainer(Container): + """Container for the result of the MuonIntensityFitter""" + width = Field(nan * u.deg, "width of the muon ring in degrees") impact = Field(nan * u.m, "distance of muon impact position from center of mirror") impact_x = Field(nan * u.m, "impact parameter x position") @@ -1126,6 +1066,8 @@ class MuonEfficiencyContainer(Container): class MuonParametersContainer(Container): + """Container for muon-ring-related image parameters""" + containment = Field(nan, "containment of the ring inside the camera") completeness = Field( nan, @@ -1140,7 +1082,7 @@ class MuonParametersContainer(Container): ) -class MuonTelescopeContainer(Container): +class MuonContainer(Container): """ Container for muon analysis """ @@ -1154,15 +1096,6 @@ class MuonTelescopeContainer(Container): ) -class MuonContainer(Container): - """Root container for muon parameters""" - - tel = Field( - default_factory=partial(Map, MuonTelescopeContainer), - description="map of tel_id to MuonTelescopeContainer", - ) - - class FlatFieldContainer(Container): """ Container for flat-field parameters obtained from a set of @@ -1306,7 +1239,7 @@ class WaveformCalibrationContainer(Container): ) -class MonitoringCameraContainer(Container): +class TelescopeMonitoringContainer(Container): """ Container for camera monitoring data """ @@ -1329,18 +1262,6 @@ class MonitoringCameraContainer(Container): ) -class MonitoringContainer(Container): - """ - Root container for monitoring data (MON) - """ - - # create the camera container - tel = Field( - default_factory=partial(Map, MonitoringCameraContainer), - description="map of tel_id to MonitoringCameraContainer", - ) - - class SimulatedShowerDistribution(Container): """ 2D histogram of simulated number of showers simulated as function of energy and @@ -1367,36 +1288,44 @@ class SimulatedShowerDistribution(Container): ) -class ArrayEventContainer(Container): - """Top-level container for all event information""" +class TelescopeEventContainer(Container): + """Container for single-telescope data""" index = Field( - default_factory=EventIndexContainer, description="event indexing information" + default_factory=TelescopeEventIndexContainer, + description="event indexing information", ) - r0 = Field(default_factory=R0Container, description="Raw Data") - r1 = Field(default_factory=R1Container, description="R1 Calibrated Data") + + simulation = Field( + None, + description="Simulated Event Information", + type=SimulationTelescopeContainer, + ) + + r0 = Field(default_factory=R0TelescopeContainer, description="Raw Data") + r1 = Field(default_factory=R1TelescopeContainer, description="R1 Calibrated Data") dl0 = Field( - default_factory=DL0Container, description="DL0 Data Volume Reduced Data" + default_factory=DL0TelescopeContainer, + description="DL0 Data Volume Reduced Data", ) - dl1 = Field(default_factory=DL1Container, description="DL1 Calibrated image") - dl2 = Field(default_factory=DL2Container, description="DL2 reconstruction info") - simulation = Field( - None, description="Simulated Event Information", type=SimulatedEventContainer + dl1 = Field( + default_factory=DL1TelescopeContainer, description="DL1 images and parameters" ) - trigger = Field( - default_factory=TriggerContainer, description="central trigger information" + dl2 = Field( + default_factory=DL2TelescopeContainer, + description="DL2 reconstructed event properties", ) count = Field(0, description="number of events processed") pointing = Field( - default_factory=PointingContainer, + default_factory=TelescopePointingContainer, description="Array and telescope pointing positions", ) calibration = Field( - default_factory=EventCalibrationContainer, + default_factory=TelescopeCalibrationContainer, description="Container for calibration coefficients for the current event", ) mon = Field( - default_factory=MonitoringContainer, + default_factory=TelescopeMonitoringContainer, description="container for event-wise monitoring data (MON)", ) muon = Field( @@ -1404,6 +1333,43 @@ class ArrayEventContainer(Container): ) +class SubarrayEventContainer(Container): + """Top-level container for all event information""" + + count = Field(0, description="number of events processed") + + index = Field( + default_factory=SubarrayEventIndexContainer, + description="event indexing information", + ) + + simulation = Field( + None, + description="Simulated Event Information", + type=SimulationSubarrayContainer, + ) + + dl0 = Field( + default_factory=DL0SubarrayContainer, + description="DL0 subarray wide information", + ) + + dl2 = Field( + default_factory=DL2SubarrayContainer, + description="DL2 reconstructed event properties", + ) + + pointing = Field( + default_factory=SubarrayPointingContainer, + description="Array and telescope pointing positions", + ) + + tel = Field( + default_factory=partial(Map, TelescopeEventContainer), + description="Telescope Events", + ) + + class SchedulingBlockContainer(Container): """Stores information about the scheduling block. This is a simplified version of the SB model, only storing what is necessary for analysis. From diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py index 2cf3eef5519..21564d8a154 100644 --- a/ctapipe/image/extractor.py +++ b/ctapipe/image/extractor.py @@ -33,7 +33,7 @@ from scipy.ndimage import convolve1d from traitlets import Bool, Int -from ctapipe.containers import DL1CameraContainer +from ctapipe.containers import DL1TelescopeContainer from ctapipe.core import TelescopeComponent from ctapipe.core.traits import ( BoolTelescopeParameter, @@ -387,7 +387,7 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): @abstractmethod def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: """ Call the relevant functions to fully extract the charge and time for the particular extractor. @@ -407,7 +407,7 @@ def __call__( Returns ------- - DL1CameraContainer: + DL1TelescopeContainer: extracted images and validity flags """ @@ -419,11 +419,11 @@ class FullWaveformSum(ImageExtractor): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: charge, peak_time = extract_around_peak( waveforms, 0, waveforms.shape[-1], 0, self.sampling_rate_ghz[tel_id] ) - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) class FixedWindowSum(ImageExtractor): @@ -478,7 +478,7 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: charge, peak_time = extract_around_peak( waveforms, self.peak_index.tel[tel_id], @@ -488,7 +488,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) class GlobalPeakWindowSum(ImageExtractor): @@ -557,7 +557,7 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: if self.pixel_fraction.tel[tel_id] == 1.0: # average over pixels then argmax over samples peak_index = waveforms[~broken_pixels].mean(axis=-2).argmax() @@ -579,7 +579,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) class LocalPeakWindowSum(ImageExtractor): @@ -633,7 +633,7 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: peak_index = waveforms.argmax(axis=-1).astype(np.int64) charge, peak_time = extract_around_peak( waveforms, @@ -644,7 +644,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) class SlidingWindowMaxSum(ImageExtractor): @@ -714,13 +714,13 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: charge, peak_time = extract_sliding_window( waveforms, self.window_width.tel[tel_id], self.sampling_rate_ghz[tel_id] ) if self.apply_integration_correction.tel[tel_id]: charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) class NeighborPeakWindowSum(ImageExtractor): @@ -780,7 +780,7 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: neighbors = self.subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse peak_index = neighbor_average_maximum( waveforms, @@ -798,7 +798,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) class BaselineSubtractedNeighborPeakWindowSum(NeighborPeakWindowSum): @@ -814,7 +814,7 @@ class BaselineSubtractedNeighborPeakWindowSum(NeighborPeakWindowSum): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: baseline_corrected = subtract_baseline( waveforms, self.baseline_start, self.baseline_end ) @@ -1263,12 +1263,12 @@ def _apply_second_pass( def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: charge1, pulse_time1, correction1 = self._apply_first_pass(waveforms, tel_id) # FIXME: properly make sure that output is 32Bit instead of downcasting here if self.disable_second_pass: - return DL1CameraContainer( + return DL1TelescopeContainer( image=(charge1 * correction1[selected_gain_channel]).astype("float32"), peak_time=pulse_time1.astype("float32"), is_valid=True, @@ -1284,7 +1284,7 @@ def __call__( broken_pixels, ) # FIXME: properly make sure that output is 32Bit instead of downcasting here - return DL1CameraContainer( + return DL1TelescopeContainer( image=charge2.astype("float32"), peak_time=pulse_time2.astype("float32"), is_valid=is_valid, @@ -1642,7 +1642,7 @@ def clip(x, lo=0.0, hi=np.inf): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: upsampling = self.upsampling.tel[tel_id] integration_window_width = self.window_width.tel[tel_id] integration_window_shift = self.window_shift.tel[tel_id] @@ -1695,4 +1695,4 @@ def __call__( if shift != 0: peak_time -= shift - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1TelescopeContainer(image=charge, peak_time=peak_time, is_valid=True) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index 6fc5e4bd8e5..f84289860d1 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -8,12 +8,13 @@ from ctapipe.coordinates import TelescopeFrame from ..containers import ( - ArrayEventContainer, CameraHillasParametersContainer, CameraTimingParametersContainer, ImageParametersContainer, IntensityStatisticsContainer, PeakTimeStatisticsContainer, + SubarrayEventContainer, + TelescopeEventContainer, TimingParametersContainer, ) from ..core import QualityQuery, TelescopeComponent @@ -116,8 +117,9 @@ def __init__( for tel_id in self.subarray.tel } - def __call__(self, event: ArrayEventContainer): - self._process_telescope_event(event) + def __call__(self, event: SubarrayEventContainer): + for tel_event in event.tel.values(): + self._process_telescope_event(tel_event) def _parameterize_image( self, @@ -209,51 +211,48 @@ def _parameterize_image( # parameterization return default - def _process_telescope_event(self, event): + def _process_telescope_event(self, tel_event: TelescopeEventContainer): """ Loop over telescopes and process the calibrated images into parameters """ - for tel_id, dl1_camera in event.dl1.tel.items(): + tel_id = tel_event.index.tel_id + dl1 = tel_event.dl1 - if self.apply_image_modifier.tel[tel_id]: - dl1_camera.image = self.modify(tel_id=tel_id, image=dl1_camera.image) + if self.apply_image_modifier.tel[tel_id]: + dl1.image = self.modify(tel_id=tel_id, image=dl1.image) - dl1_camera.image_mask = self.clean( - tel_id=tel_id, - image=dl1_camera.image, - arrival_times=dl1_camera.peak_time, - ) + dl1.image_mask = self.clean( + tel_id=tel_id, + image=dl1.image, + arrival_times=dl1.peak_time, + ) + dl1.parameters = self._parameterize_image( + tel_id=tel_id, + image=dl1.image, + signal_pixels=dl1.image_mask, + peak_time=dl1.peak_time, + default=self.default_image_container, + ) - dl1_camera.parameters = self._parameterize_image( - tel_id=tel_id, - image=dl1_camera.image, - signal_pixels=dl1_camera.image_mask, - peak_time=dl1_camera.peak_time, - default=self.default_image_container, + self.log.debug("params: %s", dl1.parameters.as_dict(recursive=True)) + + if ( + tel_event.simulation is not None + and tel_event.simulation.true_image is not None + ): + simulation = tel_event.simulation + simulation.true_parameters = self._parameterize_image( + tel_id, + image=simulation.true_image, + signal_pixels=simulation.true_image > 0, + peak_time=None, # true image from simulation has no peak time + default=DEFAULT_TRUE_IMAGE_PARAMETERS, ) + for container in simulation.true_parameters.values(): + if not container.prefix.startswith("true_"): + container.prefix = f"true_{container.prefix}" - self.log.debug("params: %s", dl1_camera.parameters.as_dict(recursive=True)) - - if ( - event.simulation is not None - and tel_id in event.simulation.tel - and event.simulation.tel[tel_id].true_image is not None - ): - sim_camera = event.simulation.tel[tel_id] - sim_camera.true_parameters = self._parameterize_image( - tel_id, - image=sim_camera.true_image, - signal_pixels=sim_camera.true_image > 0, - peak_time=None, # true image from simulation has no peak time - default=DEFAULT_TRUE_IMAGE_PARAMETERS, - ) - for container in sim_camera.true_parameters.values(): - if not container.prefix.startswith("true_"): - container.prefix = f"true_{container.prefix}" - - self.log.debug( - "sim params: %s", - event.simulation.tel[tel_id].true_parameters.as_dict( - recursive=True - ), - ) + self.log.debug( + "true image parameters: %s", + tel_event.simulation.true_parameters.as_dict(recursive=True), + ) diff --git a/ctapipe/image/muon/processor.py b/ctapipe/image/muon/processor.py index 30b382aef2c..5eb11ac01de 100644 --- a/ctapipe/image/muon/processor.py +++ b/ctapipe/image/muon/processor.py @@ -4,9 +4,9 @@ import numpy as np from ctapipe.containers import ( - ArrayEventContainer, + MuonContainer, MuonParametersContainer, - MuonTelescopeContainer, + SubarrayEventContainer, ) from ctapipe.coordinates import TelescopeFrame from ctapipe.core import QualityQuery, TelescopeComponent @@ -21,7 +21,7 @@ from .intensity_fitter import MuonIntensityFitter from .ring_fitter import MuonRingFitter -INVALID = MuonTelescopeContainer() +INVALID = MuonContainer() INVALID_PARAMETERS = MuonParametersContainer() __all__ = ["MuonProcessor"] @@ -114,23 +114,24 @@ def __init__(self, subarray, **kwargs): self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self) - def __call__(self, event: ArrayEventContainer): - for tel_id in event.dl1.tel: - self._process_telescope_event(event, tel_id) + def __call__(self, event: SubarrayEventContainer): + for tel_event in event.tel.values(): + self._process_telescope_event(tel_event) - def _process_telescope_event(self, event, tel_id): + def _process_telescope_event(self, tel_event): """ Extract and process a ring from a single image. Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer Collection of all event information tel_id: int Telescope ID of the instrument that has measured the image """ - event_index = event.index - event_id = event_index.event_id + index = tel_event.index + tel_id = index.tel_id + event_id = index.event_id if self.subarray.tel[tel_id].optics.n_mirrors != 1: self.log.warning( @@ -139,11 +140,11 @@ def _process_telescope_event(self, event, tel_id): " not supported. Exclude dual mirror telescopes via setting" " 'EventSource.allowed_tels'." ) - event.muon.tel[tel_id] = INVALID + tel_event.muon = INVALID return self.log.debug(f"Processing event {event_id}, telescope {tel_id}") - dl1 = event.dl1.tel[tel_id] + dl1 = tel_event.dl1 image = dl1.image mask = dl1.image_mask if mask is None: @@ -152,7 +153,7 @@ def _process_telescope_event(self, event, tel_id): checks = self.dl1_query(dl1_params=dl1.parameters) if not all(checks): - event.muon.tel[tel_id] = INVALID + tel_event.muon = INVALID return geometry = self.geometries[tel_id] @@ -176,9 +177,7 @@ def _process_telescope_event(self, event, tel_id): checks = self.ring_query(parameters=parameters, ring=ring, mask=mask) if not all(checks): - event.muon.tel[tel_id] = MuonTelescopeContainer( - parameters=parameters, ring=ring - ) + tel_event.muon = MuonContainer(parameters=parameters, ring=ring) return efficiency = self.intensity_fitter( @@ -197,7 +196,7 @@ def _process_telescope_event(self, event, tel_id): f", efficiency={efficiency.optical_efficiency:.2%}" ) - event.muon.tel[tel_id] = MuonTelescopeContainer( + tel_event.muon = MuonContainer( ring=ring, efficiency=efficiency, parameters=parameters ) diff --git a/ctapipe/image/muon/tests/test_intensity_fit.py b/ctapipe/image/muon/tests/test_intensity_fit.py index eeff33131b0..3de6216f5d7 100644 --- a/ctapipe/image/muon/tests/test_intensity_fit.py +++ b/ctapipe/image/muon/tests/test_intensity_fit.py @@ -36,7 +36,7 @@ def test_muon_efficiency_fit(prod5_lst, reference_location): ) center_x = 0.8 * u.deg - center_y = 0.4 * u.deg + center_y = 0.2 * u.deg radius = 1.2 * u.deg ring_width = 0.05 * u.deg impact_parameter = 5 * u.m diff --git a/ctapipe/image/muon/tests/test_processor.py b/ctapipe/image/muon/tests/test_processor.py index e2a7211ed19..513da086fb2 100644 --- a/ctapipe/image/muon/tests/test_processor.py +++ b/ctapipe/image/muon/tests/test_processor.py @@ -20,10 +20,8 @@ def test_processor(dl1_muon_file): for event in source: image_processor(event) muon_processor(event) - for tel_id in event.dl1.tel: - efficiencies.append( - event.muon.tel[tel_id].efficiency.optical_efficiency - ) + for tel_id, tel_event in event.tel.items(): + efficiencies.append(tel_event.muon.efficiency.optical_efficiency) assert len(efficiencies) > 0 # Assert there were events analyzed assert np.any( diff --git a/ctapipe/image/reducer.py b/ctapipe/image/reducer.py index 30ad1eda574..c1dcc3c6c7c 100644 --- a/ctapipe/image/reducer.py +++ b/ctapipe/image/reducer.py @@ -5,7 +5,7 @@ import numpy as np -from ctapipe.containers import DL1CameraContainer +from ctapipe.containers import DL1TelescopeContainer from ctapipe.core import TelescopeComponent from ctapipe.core.traits import ( BoolTelescopeParameter, @@ -183,7 +183,7 @@ def select_pixels(self, waveforms, tel_id=None, selected_gain_channel=None): extractor = self.image_extractors[self.image_extractor_type.tel[tel_id]] # do not treat broken pixels in data volume reduction broken_pixels = np.zeros(camera_geom.n_pixels, dtype=bool) - dl1: DL1CameraContainer = extractor( + dl1: DL1TelescopeContainer = extractor( waveforms, tel_id=tel_id, selected_gain_channel=selected_gain_channel, diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py index 169f9b58413..ec27925e986 100644 --- a/ctapipe/image/tests/test_extractor.py +++ b/ctapipe/image/tests/test_extractor.py @@ -845,18 +845,17 @@ def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): subarray = source.subarray extractor = FlashCamExtractor(subarray) - def is_flashcam(tel_id): + def is_flashcam(kv): + tel_id = kv[0] return subarray.tel[tel_id].camera.name == "FlashCam" for event in source: - for tel_id in filter(is_flashcam, event.trigger.tels_with_trigger): - true_charge = event.simulation.tel[tel_id].true_image + for tel_id, tel_event in filter(is_flashcam, event.tel.items()): + true_charge = tel_event.simulation.true_image - waveforms = event.r1.tel[tel_id].waveform + waveforms = tel_event.r1.waveform selected_gain_channel = np.zeros(waveforms.shape[0], dtype=np.int64) - broken_pixels = event.mon.tel[ - tel_id - ].pixel_status.hardware_failing_pixels[0] + broken_pixels = tel_event.mon.pixel_status.hardware_failing_pixels[0] dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid == True diff --git a/ctapipe/image/tests/test_image_processor.py b/ctapipe/image/tests/test_image_processor.py index 327b17aebf6..5000821aa38 100644 --- a/ctapipe/image/tests/test_image_processor.py +++ b/ctapipe/image/tests/test_image_processor.py @@ -29,16 +29,17 @@ def test_image_processor(example_event, example_subarray): calibrate(example_event) process_images(example_event) - for dl1tel in example_event.dl1.tel.values(): - assert isfinite(dl1tel.image_mask.sum()) - assert isfinite(dl1tel.parameters.hillas.length.value) - dl1tel.parameters.hillas.length.to("deg") - assert isfinite(dl1tel.parameters.timing.slope.value) - assert isfinite(dl1tel.parameters.leakage.pixels_width_1) - assert isfinite(dl1tel.parameters.concentration.cog) - assert isfinite(dl1tel.parameters.morphology.n_pixels) - assert isfinite(dl1tel.parameters.intensity_statistics.max) - assert isfinite(dl1tel.parameters.peak_time_statistics.max) + for tel_event in example_event.tel.values(): + dl1 = tel_event.dl1 + assert isfinite(dl1.image_mask.sum()) + assert isfinite(dl1.parameters.hillas.length.value) + dl1.parameters.hillas.length.to("deg") + assert isfinite(dl1.parameters.timing.slope.value) + assert isfinite(dl1.parameters.leakage.pixels_width_1) + assert isfinite(dl1.parameters.concentration.cog) + assert isfinite(dl1.parameters.morphology.n_pixels) + assert isfinite(dl1.parameters.intensity_statistics.max) + assert isfinite(dl1.parameters.peak_time_statistics.max) process_images.check_image.to_table() @@ -59,16 +60,18 @@ def test_image_processor_camera_frame(example_event, example_subarray): calibrate(event) process_images(event) - for dl1tel in event.dl1.tel.values(): - assert isfinite(dl1tel.image_mask.sum()) - assert isfinite(dl1tel.parameters.hillas.length.value) - dl1tel.parameters.hillas.length.to("meter") - assert isfinite(dl1tel.parameters.timing.slope.value) - assert isfinite(dl1tel.parameters.leakage.pixels_width_1) - assert isfinite(dl1tel.parameters.concentration.cog) - assert isfinite(dl1tel.parameters.morphology.n_pixels) - assert isfinite(dl1tel.parameters.intensity_statistics.max) - assert isfinite(dl1tel.parameters.peak_time_statistics.max) + assert event.tel.keys() == example_event.tel.keys() + for tel_event in event.tel.values(): + dl1 = tel_event.dl1 + assert isfinite(dl1.image_mask.sum()) + assert isfinite(dl1.parameters.hillas.length.value) + dl1.parameters.hillas.length.to("meter") + assert isfinite(dl1.parameters.timing.slope.value) + assert isfinite(dl1.parameters.leakage.pixels_width_1) + assert isfinite(dl1.parameters.concentration.cog) + assert isfinite(dl1.parameters.morphology.n_pixels) + assert isfinite(dl1.parameters.intensity_statistics.max) + assert isfinite(dl1.parameters.peak_time_statistics.max) process_images.check_image.to_table() @@ -76,11 +79,12 @@ def test_image_processor_camera_frame(example_event, example_subarray): # are in the correct frame event = deepcopy(example_event) calibrate(event) - for dl1 in event.dl1.tel.values(): - dl1.image = np.zeros_like(dl1.image) + for tel_event in event.tel.values(): + tel_event.dl1.image[:] = 0 process_images(event) - for dl1 in event.dl1.tel.values(): + for tel_event in event.tel.values(): + dl1 = tel_event.dl1 assert isinstance(dl1.parameters.hillas, CameraHillasParametersContainer) assert isinstance(dl1.parameters.timing, CameraTimingParametersContainer) assert np.isnan(dl1.parameters.hillas.length.value) diff --git a/ctapipe/image/tests/test_sliding_window_correction.py b/ctapipe/image/tests/test_sliding_window_correction.py index 73625511c95..d3a955de6bc 100644 --- a/ctapipe/image/tests/test_sliding_window_correction.py +++ b/ctapipe/image/tests/test_sliding_window_correction.py @@ -8,7 +8,7 @@ from numpy.testing import assert_allclose from traitlets.config.loader import Config -from ctapipe.containers import DL1CameraContainer +from ctapipe.containers import DL1TelescopeContainer from ctapipe.image.extractor import ImageExtractor, SlidingWindowMaxSum from ctapipe.image.toymodel import WaveformModel from ctapipe.instrument import SubarrayDescription @@ -54,7 +54,7 @@ def test_sw_pulse_lst(prod5_lst, reference_location): ) broken_pixels = np.zeros(n_pixels, dtype=bool) - dl1: DL1CameraContainer = extractor( + dl1: DL1TelescopeContainer = extractor( waveform, tel_id, selected_gain_channel, broken_pixels ) print(dl1.image / charge_true) diff --git a/ctapipe/instrument/tests/test_trigger.py b/ctapipe/instrument/tests/test_trigger.py index fec25900e6e..ff8b91fb7fd 100644 --- a/ctapipe/instrument/tests/test_trigger.py +++ b/ctapipe/instrument/tests/test_trigger.py @@ -4,22 +4,25 @@ import pytest from numpy.testing import assert_equal -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import ( + DL0SubarrayContainer, + SubarrayEventContainer, + SubarrayTriggerContainer, + TelescopeEventContainer, +) from ctapipe.io import EventSource -def assert_all_tel_keys(event, expected, ignore=None): - if ignore is None: - ignore = set() - - expected = tuple(expected) - for name, container in event.items(): - if hasattr(container, "tel"): - actual = tuple(container.tel.keys()) - if name not in ignore and actual != expected: - raise AssertionError( - f"Unexpected tel_ids in container {name}:" f"{actual} != {expected}" - ) +def dummy_event(tel_ids): + event = SubarrayEventContainer( + dl0=DL0SubarrayContainer( + trigger=SubarrayTriggerContainer( + tels_with_trigger=tel_ids, + ) + ), + tel={tel_id: TelescopeEventContainer() for tel_id in tel_ids}, + ) + return event @pytest.mark.parametrize("data_type", (list, np.array)) @@ -37,41 +40,35 @@ def test_software_trigger(subarray_prod5_paranal, data_type): ) # only one telescope, no SWAT - event = ArrayEventContainer() - event.trigger.tels_with_trigger = data_type([5]) + event = dummy_event(data_type([5])) assert trigger(event) == False - assert_equal(event.trigger.tels_with_trigger, data_type([])) + assert_equal(event.dl0.trigger.tels_with_trigger, data_type([])) # 1 LST + 1 MST, 1 LST would not have triggered LST hardware trigger # and after LST is removed, we only have 1 telescope, so no SWAT either - event = ArrayEventContainer() - event.trigger.tels_with_trigger = data_type([1, 6]) + event = dummy_event(data_type([1, 6])) assert trigger(event) == False - assert_equal(event.trigger.tels_with_trigger, data_type([])) + assert_equal(event.dl0.trigger.tels_with_trigger, data_type([])) # two MSTs and 1 LST, -> remove single LST - event = ArrayEventContainer() - event.trigger.tels_with_trigger = data_type([1, 5, 6]) + event = dummy_event(data_type([1, 5, 6])) assert trigger(event) == True - assert_equal(event.trigger.tels_with_trigger, data_type([5, 6])) + assert_equal(event.dl0.trigger.tels_with_trigger, data_type([5, 6])) # two MSTs, nothing to change - event = ArrayEventContainer() - event.trigger.tels_with_trigger = data_type([5, 6]) + event = dummy_event(data_type([5, 6])) assert trigger(event) == True - assert_equal(event.trigger.tels_with_trigger, data_type([5, 6])) + assert_equal(event.dl0.trigger.tels_with_trigger, data_type([5, 6])) # three LSTs, nothing to change - event = ArrayEventContainer() - event.trigger.tels_with_trigger = data_type([1, 2, 3]) + event = dummy_event(data_type([1, 2, 3])) assert trigger(event) == True - assert_equal(event.trigger.tels_with_trigger, data_type([1, 2, 3])) + assert_equal(event.dl0.trigger.tels_with_trigger, data_type([1, 2, 3])) # thee LSTs, plus MSTs, nothing to change - event = ArrayEventContainer() - event.trigger.tels_with_trigger = data_type([1, 2, 3, 5, 6, 7]) + event = dummy_event(data_type([1, 2, 3, 5, 6, 7])) assert trigger(event) == True - assert_equal(event.trigger.tels_with_trigger, data_type([1, 2, 3, 5, 6, 7])) + assert_equal(event.dl0.trigger.tels_with_trigger, data_type([1, 2, 3, 5, 6, 7])) @pytest.mark.parametrize("allowed_tels", (None, list(range(1, 20)))) @@ -115,10 +112,10 @@ def test_software_trigger_simtel(allowed_tels): ], ) - for e, expected_tels in zip(source, expected): - trigger(e) - assert_equal(e.trigger.tels_with_trigger, expected_tels) - assert_all_tel_keys(e, expected_tels, ignore={"dl0", "dl1", "dl2", "muon"}) + for event, expected_tels in zip(source, expected): + trigger(event) + assert_equal(event.dl0.trigger.tels_with_trigger, expected_tels) + assert set(event.tel.keys()) == set(expected_tels) def test_software_trigger_simtel_single_lsts(): @@ -162,12 +159,10 @@ def test_software_trigger_simtel_single_lsts(): ], ) - for e, expected_tels in zip(source, expected): - print(e.trigger.tels_with_trigger) - trigger(e) - print(e.trigger.tels_with_trigger, expected_tels) - assert_equal(e.trigger.tels_with_trigger, expected_tels) - assert_all_tel_keys(e, expected_tels, ignore={"dl0", "dl1", "dl2", "muon"}) + for event, expected_tels in zip(source, expected): + trigger(event) + assert_equal(event.dl0.trigger.tels_with_trigger, expected_tels) + assert set(event.tel.keys()) == set(expected_tels) def test_software_trigger_simtel_process(tmp_path): diff --git a/ctapipe/instrument/trigger.py b/ctapipe/instrument/trigger.py index 29e79428f3f..27dcf5a05bb 100644 --- a/ctapipe/instrument/trigger.py +++ b/ctapipe/instrument/trigger.py @@ -1,6 +1,6 @@ import numpy as np -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import SubarrayEventContainer from ctapipe.core import TelescopeComponent from ctapipe.core.traits import Integer, IntTelescopeParameter @@ -64,7 +64,7 @@ def __init__(self, subarray, *args, **kwargs): for type in self.subarray.telescope_types } - def __call__(self, event: ArrayEventContainer) -> bool: + def __call__(self, event: SubarrayEventContainer) -> bool: """ Remove telescope events that have not the required number of telescopes of a given type from the subarray event and decide if the event would @@ -87,7 +87,7 @@ def __call__(self, event: ArrayEventContainer) -> bool: if min_tels == 0: continue - tels_with_trigger = set(event.trigger.tels_with_trigger) + tels_with_trigger = set(event.dl0.trigger.tels_with_trigger) tel_ids = self._ids_by_type[tel_type_str] tels_in_event = tels_with_trigger.intersection(tel_ids) @@ -103,24 +103,17 @@ def __call__(self, event: ArrayEventContainer) -> bool: tels_removed.add(tel_id) # remove any related data - for container in event.values(): - if hasattr(container, "tel"): - tel_map = container.tel - if tel_id in tel_map: - del tel_map[tel_id] + del event.tel[tel_id] if len(tels_removed) > 0: # convert to array with correct dtype to have setdiff1d work correctly tels_removed = np.fromiter(tels_removed, np.uint16, len(tels_removed)) - event.trigger.tels_with_trigger = np.setdiff1d( - event.trigger.tels_with_trigger, tels_removed, assume_unique=True + event.dl0.trigger.tels_with_trigger = np.setdiff1d( + event.dl0.trigger.tels_with_trigger, tels_removed, assume_unique=True ) - if len(event.trigger.tels_with_trigger) < self.min_telescopes: - event.trigger.tels_with_trigger = [] - # remove any related data - for container in event.values(): - if hasattr(container, "tel"): - container.tel.clear() + if len(event.dl0.trigger.tels_with_trigger) < self.min_telescopes: + event.dl0.trigger.tels_with_trigger = [] + event.tel.clear() return False return True diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py index 3a5cc8f00dd..00a9bf3d4bd 100644 --- a/ctapipe/io/datawriter.py +++ b/ctapipe/io/datawriter.py @@ -14,9 +14,9 @@ from traitlets import Dict, Instance from ..containers import ( - ArrayEventContainer, SimulatedShowerDistribution, - TelEventIndexContainer, + SubarrayEventContainer, + TelescopeEventContainer, ) from ..core import Component, Container, Field, Provenance, ToolConfigurationError from ..core.traits import Bool, CaselessStrEnum, Float, Int, Path, Unicode @@ -34,14 +34,6 @@ tables.parameters.NODE_CACHE_SLOTS = 3000 # fixes problem with too many datasets -def _get_tel_index(event, tel_id): - return TelEventIndexContainer( - obs_id=event.index.obs_id, - event_id=event.index.event_id, - tel_id=np.int16(tel_id), - ) - - # define the version of the data model written here. This should be updated # when necessary: # - increase the major number if there is a breaking change to the model @@ -293,46 +285,195 @@ def __enter__(self): def __exit__(self, exc_type, exc_value, exc_traceback): self.finish() - def __call__(self, event: ArrayEventContainer): + def __call__(self, event: SubarrayEventContainer): """ Write a single event to the output file. """ self._at_least_one_event = True self.log.debug("WRITING EVENT %s", event.index) - self._write_subarray_pointing(event) - self._write_trigger(event) + self._write_simulation_subarray(event) + self._write_dl0_subarray(event) + self._write_dl2_subarray(event) + + for tel_id, tel_event in event.tel.items(): + telescope = self._subarray.tel[tel_id] + self.log.debug("WRITING TELESCOPE %s: %s", tel_id, telescope) + self._write_telescope_event(tel_event) + def _write_simulation_subarray(self, event): if event.simulation is not None and event.simulation.shower is not None: self._writer.write( table_name="simulation/event/subarray/shower", containers=[event.index, event.simulation.shower], ) - for tel_id, sim in event.simulation.tel.items(): - table_name = self.table_name(tel_id) - tel_index = _get_tel_index(event, tel_id) - self._writer.write( - f"simulation/event/telescope/impact/{table_name}", - [tel_index, sim.impact], - ) + def _write_dl0_subarray(self, event): + self._writer.write( + table_name="dl0/event/subarray/trigger", + containers=[event.index, event.dl0.trigger], + ) - if self.write_waveforms: - self._write_r1_telescope_events(event) + self._write_subarray_pointing(event) + + def _write_telescope_event(self, tel_event): + table_name = self.table_name(tel_event.index.tel_id) + self._writer.write( + f"simulation/event/telescope/impact/{table_name}", + [tel_event.index, tel_event.simulation.impact], + ) + + if tel_event.simulation is not None: + self._write_simulation_telescope(tel_event) if self.write_raw_waveforms: - self._write_r0_telescope_events(event) + self._write_r0_telescope(tel_event) + + if self.write_waveforms: + self._write_r1_telescope(tel_event) + + # always write dl0, contains needed trigger info etc. + self._write_dl0_telescope(tel_event) + self._write_pointing_telescope(tel_event) - # write telescope event data - self._write_dl1_telescope_events(event) + if self.write_images or self.write_parameters: + self._write_dl1_telescope(tel_event) - # write DL2 info if requested if self.write_showers: - self._write_dl2_telescope_events(event) - self._write_dl2_stereo_event(event) + self._write_dl2_telescope(tel_event) if self.write_muon_parameters: - self._write_muon_telescope_events(event) + self._write_muon_telescope(tel_event) + + def _write_r0_telescope(self, tel_event: TelescopeEventContainer): + table_name = self.table_name(tel_event.index.tel_id) + tel_event.r0.prefix = "" + self._writer.write( + f"r0/event/telescope/{table_name}", + [tel_event.index, tel_event.r0], + ) + + def _write_r1_telescope(self, tel_event: TelescopeEventContainer): + table_name = self.table_name(tel_event.index.tel_id) + tel_event.r1.prefix = "" + self._writer.write( + f"r1/event/telescope/{table_name}", + [tel_event.index, tel_event.r1], + ) + + def _write_dl0_telescope(self, tel_event: TelescopeEventContainer): + self._writer.write( + "dl0/event/telescope/trigger", + [tel_event.index, tel_event.dl0.trigger], + ) + + def _write_pointing_telescope(self, tel_event: TelescopeEventContainer): + tel_index = tel_event.index + tel_id = tel_index.tel_id + current_pointing = (tel_event.pointing.azimuth, tel_event.pointing.altitude) + if current_pointing != self._last_pointing_tel[tel_id]: + table_name = self.table_name(tel_id) + tel_event.pointing.prefix = "" + self._writer.write( + f"dl0/monitoring/telescope/pointing/{table_name}", + [tel_event.dl0.trigger, tel_event.pointing], + ) + self._last_pointing_tel[tel_id] = current_pointing + + def _write_dl1_telescope(self, tel_event: TelescopeEventContainer): + tel_index = tel_event.index + tel_id = tel_index.tel_id + table_name = self.table_name(tel_id) + tel_event.dl1.prefix = "" # don't want a prefix for this container + + if self.write_parameters: + self._writer.write( + table_name=f"dl1/event/telescope/parameters/{table_name}", + containers=[tel_index, *tel_event.dl1.parameters.values()], + ) + + if self.write_images: + if tel_event.dl1.image is None: + raise ValueError( + "DataWriter.write_images is True but event does not contain image" + ) + + self._writer.write( + table_name=f"dl1/event/telescope/images/{table_name}", + containers=[tel_index, tel_event.dl1], + ) + + def _write_simulation_telescope(self, tel_event: TelescopeEventContainer): + tel_index = tel_event.index + tel_id = tel_index.tel_id + table_name = self.table_name(tel_id) + + # always write this, so that at least the sum is included + self._writer.write( + f"simulation/event/telescope/images/{table_name}", + [tel_index, tel_event.simulation], + ) + + has_sim_image = ( + tel_event.simulation is not None + and tel_event.simulation.true_image is not None + ) + if self.write_parameters and has_sim_image: + true_parameters = tel_event.simulation.true_parameters + # only write the available containers, no peak time related + # features for true image available. + self._writer.write( + f"simulation/event/telescope/parameters/{table_name}", + [ + tel_index, + true_parameters.hillas, + true_parameters.leakage, + true_parameters.concentration, + true_parameters.morphology, + true_parameters.intensity_statistics, + ], + ) + + def _write_muon_telescope(self, tel_event: TelescopeEventContainer): + table_name = self.table_name(tel_event.index.tel_id) + muon = tel_event.muon + self._writer.write( + f"dl1/event/telescope/muon/{table_name}", + [tel_event.index, muon.ring, muon.parameters, muon.efficiency], + ) + + def _write_dl2_telescope(self, tel_event: TelescopeEventContainer): + """ + write per-telescope DL2 shower information. + + Currently this writes to a single table per type of shower + reconstruction and per algorithm, with all telescopes combined. + """ + + table_name = self.table_name(tel_event.index.tel_id) + + for container_name, algorithm_map in tel_event.dl2.items(): + for algorithm, container in algorithm_map.items(): + name = f"dl2/event/telescope/{container_name}/{algorithm}/{table_name}" + + self._writer.write( + table_name=name, containers=[tel_event.index, container] + ) + + def _write_dl2_subarray(self, event: SubarrayEventContainer): + """ + write per-telescope DL2 shower information to e.g. + `/dl2/event/stereo/{geometry,energy,classification}/` + """ + for container_name, algorithm_map in event.dl2.items(): + for algorithm, container in algorithm_map.items(): + # note this will only write info if the particular algorithm + # generated it (otherwise the algorithm map is empty, and no + # data will be written) + self._writer.write( + table_name=f"dl2/event/subarray/{container_name}/{algorithm}", + containers=[event.index, container], + ) def finish(self): """called after all events are done""" @@ -431,7 +572,7 @@ def _setup_writer(self): tr_tel_list_to_mask = TelListToMaskTransform(self._subarray) writer.add_column_transform( - table_name="dl1/event/subarray/trigger", + table_name="dl0/event/subarray/trigger", col_name="tels_with_trigger", transform=tr_tel_list_to_mask, ) @@ -439,13 +580,13 @@ def _setup_writer(self): # currently the trigger info is used for the event time, but we dont' # want the other bits of the trigger container in the pointing or other # montitoring containers - writer.exclude("dl1/monitoring/subarray/pointing", "event_type") - writer.exclude("dl1/monitoring/subarray/pointing", "tels_with_trigger") - writer.exclude("dl1/monitoring/subarray/pointing", "n_trigger_pixels") - writer.exclude("/dl1/event/telescope/trigger", "trigger_pixels") - writer.exclude("/dl1/monitoring/telescope/pointing/.*", "n_trigger_pixels") - writer.exclude("/dl1/monitoring/telescope/pointing/.*", "trigger_pixels") - writer.exclude("/dl1/monitoring/event/pointing/.*", "event_type") + writer.exclude("dl0/monitoring/subarray/pointing", "event_type") + writer.exclude("dl0/monitoring/subarray/pointing", "tels_with_trigger") + writer.exclude("dl0/monitoring/subarray/pointing", "n_trigger_pixels") + writer.exclude("/dl0/event/telescope/trigger", "trigger_pixels") + writer.exclude("/dl0/monitoring/telescope/pointing/.*", "n_trigger_pixels") + writer.exclude("/dl0/monitoring/telescope/pointing/.*", "trigger_pixels") + writer.exclude("/dl0/monitoring/event/pointing/.*", "event_type") writer.exclude("/dl1/event/telescope/images/.*", "parameters") writer.exclude("/simulation/event/telescope/images/.*", "true_parameters") @@ -503,13 +644,15 @@ def _setup_writer(self): self._writer = writer self.log.debug("Writer initialized: %s", self._writer) - def _write_subarray_pointing(self, event: ArrayEventContainer): + def _write_subarray_pointing(self, event: SubarrayEventContainer): """store subarray pointing info in a monitoring table""" pnt = event.pointing - current_pointing = (pnt.array_azimuth, pnt.array_altitude) + current_pointing = (pnt.azimuth, pnt.altitude) if current_pointing != self._last_pointing: pnt.prefix = "" - self._writer.write("dl1/monitoring/subarray/pointing", [event.trigger, pnt]) + self._writer.write( + "dl0/monitoring/subarray/pointing", [event.dl0.trigger, pnt] + ) self._last_pointing = current_pointing def _write_scheduling_and_observation_blocks(self): @@ -612,156 +755,6 @@ def table_name(self, tel_id): """construct dataset table names depending on chosen split method""" return f"tel_{tel_id:03d}" - def _write_trigger(self, event: ArrayEventContainer): - """ - Write trigger information - """ - self._writer.write( - table_name="dl1/event/subarray/trigger", - containers=[event.index, event.trigger], - ) - - for tel_id, trigger in event.trigger.tel.items(): - self._writer.write( - "dl1/event/telescope/trigger", (_get_tel_index(event, tel_id), trigger) - ) - - def _write_r1_telescope_events(self, event: ArrayEventContainer): - for tel_id, r1_tel in event.r1.tel.items(): - - tel_index = _get_tel_index(event, tel_id) - table_name = self.table_name(tel_id) - - r1_tel.prefix = "" - self._writer.write(f"r1/event/telescope/{table_name}", [tel_index, r1_tel]) - - def _write_r0_telescope_events(self, event: ArrayEventContainer): - for tel_id, r0_tel in event.r0.tel.items(): - - tel_index = _get_tel_index(event, tel_id) - table_name = self.table_name(tel_id) - - r0_tel.prefix = "" - self._writer.write(f"r0/event/telescope/{table_name}", [tel_index, r0_tel]) - - def _write_dl1_telescope_events(self, event: ArrayEventContainer): - """ - add entries to the event/telescope tables for each telescope in a single - event - """ - - # pointing info - for tel_id, pnt in event.pointing.tel.items(): - current_pointing = (pnt.azimuth, pnt.altitude) - if current_pointing != self._last_pointing_tel[tel_id]: - pnt.prefix = "" - self._writer.write( - f"dl1/monitoring/telescope/pointing/tel_{tel_id:03d}", - [event.trigger.tel[tel_id], pnt], - ) - self._last_pointing_tel[tel_id] = current_pointing - - for tel_id, dl1_camera in event.dl1.tel.items(): - tel_index = _get_tel_index(event, tel_id) - - dl1_camera.prefix = "" # don't want a prefix for this container - telescope = self._subarray.tel[tel_id] - self.log.debug("WRITING TELESCOPE %s: %s", tel_id, telescope) - - table_name = self.table_name(tel_id) - - if self.write_parameters: - self._writer.write( - table_name=f"dl1/event/telescope/parameters/{table_name}", - containers=[tel_index, *dl1_camera.parameters.values()], - ) - - if self.write_images: - if dl1_camera.image is None: - raise ValueError( - "DataWriter.write_images is True but event does not contain image" - ) - - self._writer.write( - table_name=f"dl1/event/telescope/images/{table_name}", - containers=[tel_index, dl1_camera], - ) - - if self._is_simulation: - # always write this, so that at least the sum is included - self._writer.write( - f"simulation/event/telescope/images/{table_name}", - [tel_index, event.simulation.tel[tel_id]], - ) - - has_sim_image = ( - tel_id in event.simulation.tel - and event.simulation.tel[tel_id].true_image is not None - ) - if self.write_parameters and has_sim_image: - true_parameters = event.simulation.tel[tel_id].true_parameters - # only write the available containers, no peak time related - # features for true image available. - self._writer.write( - f"simulation/event/telescope/parameters/{table_name}", - [ - tel_index, - true_parameters.hillas, - true_parameters.leakage, - true_parameters.concentration, - true_parameters.morphology, - true_parameters.intensity_statistics, - ], - ) - - def _write_muon_telescope_events(self, event: ArrayEventContainer): - - for tel_id, muon in event.muon.tel.items(): - table_name = self.table_name(tel_id) - tel_index = _get_tel_index(event, tel_id) - self._writer.write( - f"dl1/event/telescope/muon/{table_name}", - [tel_index, muon.ring, muon.parameters, muon.efficiency], - ) - - def _write_dl2_telescope_events(self, event: ArrayEventContainer): - """ - write per-telescope DL2 shower information. - - Currently this writes to a single table per type of shower - reconstruction and per algorithm, with all telescopes combined. - """ - - for tel_id, dl2_tel in event.dl2.tel.items(): - table_name = self.table_name(tel_id) - - tel_index = _get_tel_index(event, tel_id) - for container_name, algorithm_map in dl2_tel.items(): - for algorithm, container in algorithm_map.items(): - name = ( - f"dl2/event/telescope/{container_name}/{algorithm}/{table_name}" - ) - - self._writer.write( - table_name=name, containers=[tel_index, container] - ) - - def _write_dl2_stereo_event(self, event: ArrayEventContainer): - """ - write per-telescope DL2 shower information to e.g. - `/dl2/event/stereo/{geometry,energy,classification}/` - """ - # pylint: disable=no-self-use - for container_name, algorithm_map in event.dl2.stereo.items(): - for algorithm, container in algorithm_map.items(): - # note this will only write info if the particular algorithm - # generated it (otherwise the algorithm map is empty, and no - # data will be written) - self._writer.write( - table_name=f"dl2/event/subarray/{container_name}/{algorithm}", - containers=[event.index, container], - ) - def _generate_table_indices(self, h5file, start_node): """helper to generate PyTables index tabnles for common columns""" for node in h5file.iter_nodes(start_node): diff --git a/ctapipe/io/eventsource.py b/ctapipe/io/eventsource.py index 8fa293c8bb0..a151709cdec 100644 --- a/ctapipe/io/eventsource.py +++ b/ctapipe/io/eventsource.py @@ -10,10 +10,10 @@ from ctapipe.atmosphere import AtmosphereDensityProfile from ..containers import ( - ArrayEventContainer, ObservationBlockContainer, SchedulingBlockContainer, SimulationConfigContainer, + SubarrayEventContainer, ) from ..core import Provenance, ToolConfigurationError from ..core.component import Component, find_config_in_hierarchy @@ -28,7 +28,7 @@ class EventSource(Component): """ Parent class for EventSources. - EventSources read input files and generate `~ctapipe.containers.ArrayEventContainer` + EventSources read input files and generate `~ctapipe.containers.SubarrayEventContainer` instances when iterated over. A new EventSource should be created for each type of event file read @@ -72,7 +72,7 @@ class EventSource(Component): 0 1 - **NOTE**: EventSource implementations should not reuse the same ArrayEventContainer, + **NOTE**: EventSource implementations should not reuse the same SubarrayEventContainer, as these are mutable and may lead to errors when analyzing multiple events. @@ -300,7 +300,7 @@ def atmosphere_density_profile(self) -> AtmosphereDensityProfile: return None @abstractmethod - def _generator(self) -> Generator[ArrayEventContainer, None, None]: + def _generator(self) -> Generator[SubarrayEventContainer, None, None]: """ Abstract method to be defined in child class. diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py index c7bda1f9b7d..23b9035ed3a 100644 --- a/ctapipe/io/hdf5eventsource.py +++ b/ctapipe/io/hdf5eventsource.py @@ -12,36 +12,42 @@ from ctapipe.instrument.optics import FocalLengthKind from ..containers import ( - ArrayEventContainer, CameraHillasParametersContainer, CameraTimingParametersContainer, ConcentrationContainer, - DL1CameraContainer, - EventIndexContainer, + DL0SubarrayContainer, + DL0TelescopeContainer, + DL1TelescopeContainer, + DL2TelescopeContainer, HillasParametersContainer, ImageParametersContainer, IntensityStatisticsContainer, LeakageContainer, MorphologyContainer, + MuonContainer, MuonEfficiencyContainer, MuonParametersContainer, MuonRingContainer, - MuonTelescopeContainer, ObservationBlockContainer, ParticleClassificationContainer, PeakTimeStatisticsContainer, - R1CameraContainer, + R1TelescopeContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, SchedulingBlockContainer, - SimulatedEventContainer, SimulatedShowerContainer, SimulationConfigContainer, + SimulationSubarrayContainer, + SimulationTelescopeContainer, + SubarrayEventContainer, + SubarrayEventIndexContainer, + SubarrayTriggerContainer, + TelescopeEventContainer, + TelescopeEventIndexContainer, TelescopeImpactParameterContainer, + TelescopePointingContainer, TelescopeTriggerContainer, - TelEventIndexContainer, TimingParametersContainer, - TriggerContainer, ) from ..core import Container, Field from ..core.traits import UseEnum @@ -132,8 +138,8 @@ class HDF5EventSource(EventSource): specifying the file to be read. Looping over the EventSource yields events from the _generate_events - method. An event equals an ArrayEventContainer instance. - See ctapipe.containers.ArrayEventContainer for details. + method. An event equals an SubarrayEventContainer instance. + See ctapipe.containers.SubarrayEventContainer for details. Attributes ---------- @@ -201,13 +207,12 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs): if self.allowed_tels: self._subarray = self._full_subarray.select_subarray(self.allowed_tels) + self._allowed_tels_array = np.array(list(self.allowed_tels)) else: self._subarray = self._full_subarray + self._simulation_configs = self._parse_simulation_configs() - ( - self._scheduling_block, - self._observation_block, - ) = self._parse_sb_and_ob_configs() + self._scheduling_block, self._observation_block = self._parse_sb_and_ob() version = self.file_.root._v_attrs["CTA PRODUCT DATA MODEL VERSION"] self.datamodel_version = tuple(map(int, version.lstrip("v").split("."))) @@ -326,7 +331,7 @@ def simulation_config(self) -> Dict[int, SimulationConfigContainer]: return self._simulation_configs def __len__(self): - n_events = len(self.file_.root.dl1.event.subarray.trigger) + n_events = len(self.file_.root.dl0.event.subarray.trigger) if self.max_events is not None: return min(n_events, self.max_events) return n_events @@ -356,7 +361,7 @@ class ObsIdContainer(Container): else: return {} - def _parse_sb_and_ob_configs(self): + def _parse_sb_and_ob(self): """read Observation and Scheduling block configurations""" sb_reader = HDF5TableReader(self.file_).read( @@ -388,14 +393,14 @@ def _is_hillas_in_camera_frame(self): def _generator(self): """ - Yield ArrayEventContainer to iterate through events. + Yield SubarrayEventContainer to iterate through events. """ self.reader = HDF5TableReader(self.file_) if DataLevel.R1 in self.datalevels: waveform_readers = { table.name: self.reader.read( - f"/r1/event/telescope/{table.name}", R1CameraContainer + f"/r1/event/telescope/{table.name}", R1TelescopeContainer ) for table in self.file_.root.r1.event.telescope } @@ -410,7 +415,7 @@ def _generator(self): image_readers = { table.name: self.reader.read( f"/dl1/event/telescope/images/{table.name}", - DL1CameraContainer, + DL1TelescopeContainer, ignore_columns=ignore_columns, ) for table in self.file_.root.dl1.event.telescope.images @@ -556,48 +561,68 @@ def _generator(self): # Setup iterators for the array events events = HDF5TableReader(self.file_).read( - "/dl1/event/subarray/trigger", - [TriggerContainer, EventIndexContainer], + "/dl0/event/subarray/trigger", + [SubarrayTriggerContainer, SubarrayEventIndexContainer], ignore_columns={"tel"}, ) telescope_trigger_reader = HDF5TableReader(self.file_).read( - "/dl1/event/telescope/trigger", - [TelEventIndexContainer, TelescopeTriggerContainer], + "/dl0/event/telescope/trigger", + [TelescopeEventIndexContainer, TelescopeTriggerContainer], ignore_columns={"trigger_pixels"}, ) - array_pointing_finder = IndexFinder( - self.file_.root.dl1.monitoring.subarray.pointing.col("time") - ) + table_path = "/dl1/monitoring/subarray/pointing" + if table_path not in self.file_.root: + table_path = "/dl0/monitoring/subarray/pointing" + + array_pointing_finder = IndexFinder(self.file_.root[table_path].col("time")) + group_path = "/dl1/monitoring/telescope/pointing" + if group_path not in self.file_.root: + group_path = "/dl0/monitoring/telescope/pointing" tel_pointing_finder = { table.name: IndexFinder(table.col("time")) - for table in self.file_.root.dl1.monitoring.telescope.pointing + for table in self.file_.root[group_path] } counter = 0 for trigger, index in events: - data = ArrayEventContainer( - trigger=trigger, + event = SubarrayEventContainer( + dl0=DL0SubarrayContainer(trigger=trigger), count=counter, index=index, - simulation=SimulatedEventContainer() if self.is_simulation else None, + simulation=SimulationSubarrayContainer() + if self.is_simulation + else None, ) # Maybe take some other metadata, but there are still some 'unknown' # written out by the stage1 tool - data.meta["origin"] = self.file_.root._v_attrs["CTA PROCESS TYPE"] - data.meta["input_url"] = self.input_url - data.meta["max_events"] = self.max_events + event.meta["origin"] = self.file_.root._v_attrs["CTA PROCESS TYPE"] + event.meta["input_url"] = self.input_url + event.meta["max_events"] = self.max_events - data.trigger.tels_with_trigger = self._full_subarray.tel_mask_to_tel_ids( - data.trigger.tels_with_trigger + event.dl0.trigger.tels_with_trigger = ( + self._full_subarray.tel_mask_to_tel_ids( + event.dl0.trigger.tels_with_trigger + ) ) - full_tels_with_trigger = data.trigger.tels_with_trigger.copy() + full_tels_with_trigger = event.dl0.trigger.tels_with_trigger.copy() if self.allowed_tels: - data.trigger.tels_with_trigger = np.intersect1d( - data.trigger.tels_with_trigger, np.array(list(self.allowed_tels)) + event.dl0.trigger.tels_with_trigger = np.intersect1d( + event.dl0.trigger.tels_with_trigger, + self._allowed_tels_array, ) + if self.is_simulation: + event.simulation.shower = next(mc_shower_reader) + + for kind, readers in dl2_readers.items(): + c = getattr(event.dl2, kind) + for algorithm, reader in readers.items(): + c[algorithm] = next(reader) + + self._fill_array_pointing(event, array_pointing_finder) + # the telescope trigger table contains triggers for all telescopes # that participated in the event, so we need to read a row for each # of them, ignoring the ones not in allowed_tels after reading @@ -607,53 +632,43 @@ def _generator(self): if self.allowed_tels and tel_id not in self.allowed_tels: continue - data.trigger.tel[tel_index.tel_id] = tel_trigger - - if self.is_simulation: - data.simulation.shower = next(mc_shower_reader) - - for kind, readers in dl2_readers.items(): - c = getattr(data.dl2.stereo, kind) - for algorithm, reader in readers.items(): - c[algorithm] = next(reader) - - # this needs to stay *after* reading the telescope trigger table - # and after reading all subarray event information, so that we don't - # go out of sync - if len(data.trigger.tels_with_trigger) == 0: - continue - - self._fill_array_pointing(data, array_pointing_finder) - self._fill_telescope_pointing(data, tel_pointing_finder) - - for tel_id in data.trigger.tel.keys(): key = f"tel_{tel_id:03d}" - if self.allowed_tels and tel_id not in self.allowed_tels: - continue + containers = { + "index": tel_index, + "dl0": DL0TelescopeContainer(trigger=tel_trigger), + } - if key in true_impact_readers: - data.simulation.tel[tel_id].impact = next(true_impact_readers[key]) + containers["pointing"] = self._fill_telescope_pointing( + tel_id, tel_trigger.time, tel_pointing_finder + ) if DataLevel.R1 in self.datalevels: - data.r1.tel[tel_id] = next(waveform_readers[key]) + containers["r1"] = next(waveform_readers[key]) + + if self.is_simulation: + containers["simulation"] = SimulationTelescopeContainer() - if self.has_simulated_dl1: - simulated = data.simulation.tel[tel_id] + if key in true_impact_readers: + containers["simulation"].impact = next(true_impact_readers[key]) if DataLevel.DL1_IMAGES in self.datalevels: - data.dl1.tel[tel_id] = next(image_readers[key]) + containers["dl1"] = next(image_readers[key]) if self.has_simulated_dl1: simulated_image_row = next(simulated_image_iterators[key]) - simulated.true_image = simulated_image_row["true_image"] + containers["simulation"].true_image = simulated_image_row[ + "true_image" + ] + else: + containers["dl1"] = DL1TelescopeContainer() if DataLevel.DL1_PARAMETERS in self.datalevels: # Is there a smarter way to unpack this? # Best would probbaly be if we could directly read # into the ImageParametersContainer params = next(param_readers[key]) - data.dl1.tel[tel_id].parameters = ImageParametersContainer( + containers["dl1"].parameters = ImageParametersContainer( hillas=params[0], timing=params[1], leakage=params[2], @@ -664,7 +679,9 @@ def _generator(self): ) if self.has_simulated_dl1: simulated_params = next(simulated_param_readers[key]) - simulated.true_parameters = ImageParametersContainer( + containers[ + "simulation" + ].true_parameters = ImageParametersContainer( hillas=simulated_params[0], leakage=simulated_params[1], concentration=simulated_params[2], @@ -674,14 +691,15 @@ def _generator(self): if self.has_muon_parameters: ring, parameters, efficiency = next(muon_readers[key]) - data.muon.tel[tel_id] = MuonTelescopeContainer( + containers["muon"] = MuonContainer( ring=ring, parameters=parameters, efficiency=efficiency, ) + dl2 = DL2TelescopeContainer() for kind, algorithms in dl2_tel_readers.items(): - c = getattr(data.dl2.tel[tel_id], kind) + c = getattr(dl2, kind) for algorithm, readers in algorithms.items(): c[algorithm] = next(readers[key]) @@ -690,70 +708,74 @@ def _generator(self): prefix = f"{algorithm}_tel_{c[algorithm].default_prefix}" c[algorithm].prefix = prefix - yield data + event.tel[tel_id] = TelescopeEventContainer(**containers, dl2=dl2) + + # this needs to stay *after* reading the telescope trigger table + # and after reading all subarray event information, so that we don't + # go out of sync + if len(event.tel) == 0: + continue + + yield event counter += 1 @lazyproperty def _subarray_pointing_attrs(self): - table = self.file_.root.dl1.monitoring.subarray.pointing + table = self.file_.root.dl0.monitoring.subarray.pointing return get_column_attrs(table) @lru_cache(maxsize=1000) def _telescope_pointing_attrs(self, tel_id): - pointing_group = self.file_.root.dl1.monitoring.telescope.pointing + pointing_group = self.file_.root.dl0.monitoring.telescope.pointing return get_column_attrs(pointing_group[f"tel_{tel_id:03d}"]) - def _fill_array_pointing(self, data, array_pointing_finder): + def _fill_array_pointing(self, event, array_pointing_finder): """ Fill the array pointing information of a given event """ # Only unique pointings are stored, so reader.read() wont work as easily # Thats why we match the pointings based on trigger time - closest_time_index = array_pointing_finder.closest(data.trigger.time.mjd) - table = self.file_.root.dl1.monitoring.subarray.pointing + closest_time_index = array_pointing_finder.closest(event.dl0.trigger.time.mjd) + table = self.file_.root.dl0.monitoring.subarray.pointing array_pointing = table[closest_time_index] - data.pointing.array_azimuth = u.Quantity( - array_pointing["array_azimuth"], - self._subarray_pointing_attrs["array_azimuth"]["UNIT"], + event.pointing.azimuth = u.Quantity( + array_pointing["azimuth"], + self._subarray_pointing_attrs["azimuth"]["UNIT"], ) - data.pointing.array_altitude = u.Quantity( - array_pointing["array_altitude"], - self._subarray_pointing_attrs["array_altitude"]["UNIT"], + event.pointing.altitude = u.Quantity( + array_pointing["altitude"], + self._subarray_pointing_attrs["altitude"]["UNIT"], ) - data.pointing.array_ra = u.Quantity( - array_pointing["array_ra"], - self._subarray_pointing_attrs["array_ra"]["UNIT"], + event.pointing.ra = u.Quantity( + array_pointing["ra"], + self._subarray_pointing_attrs["ra"]["UNIT"], ) - data.pointing.array_dec = u.Quantity( - array_pointing["array_dec"], - self._subarray_pointing_attrs["array_dec"]["UNIT"], + event.pointing.dec = u.Quantity( + array_pointing["dec"], + self._subarray_pointing_attrs["dec"]["UNIT"], ) - def _fill_telescope_pointing(self, data, tel_pointing_finder): + def _fill_telescope_pointing(self, tel_id, time, tel_pointing_finder): """ Fill the telescope pointing information of a given event """ # Same comments as to _fill_array_pointing apply - pointing_group = self.file_.root.dl1.monitoring.telescope.pointing - for tel_id in data.trigger.tel.keys(): - key = f"tel_{tel_id:03d}" - - if self.allowed_tels and tel_id not in self.allowed_tels: - continue + pointing_group = self.file_.root.dl0.monitoring.telescope.pointing + key = f"tel_{tel_id:03d}" - tel_pointing_table = pointing_group[key] - closest_time_index = tel_pointing_finder[key].closest( - data.trigger.tel[tel_id].time.mjd - ) + tel_pointing_table = pointing_group[key] + closest_time_index = tel_pointing_finder[key].closest(time.mjd) - pointing_telescope = tel_pointing_table[closest_time_index] - attrs = self._telescope_pointing_attrs(tel_id) - data.pointing.tel[tel_id].azimuth = u.Quantity( + pointing_telescope = tel_pointing_table[closest_time_index] + attrs = self._telescope_pointing_attrs(tel_id) + return TelescopePointingContainer( + azimuth=u.Quantity( pointing_telescope["azimuth"], attrs["azimuth"]["UNIT"], - ) - data.pointing.tel[tel_id].altitude = u.Quantity( + ), + altitude=u.Quantity( pointing_telescope["altitude"], attrs["altitude"]["UNIT"], - ) + ), + ) diff --git a/ctapipe/io/hdf5merger.py b/ctapipe/io/hdf5merger.py index 8e97d43982b..685ccbacb67 100644 --- a/ctapipe/io/hdf5merger.py +++ b/ctapipe/io/hdf5merger.py @@ -39,6 +39,8 @@ class NodeType(enum.Enum): "/simulation/event/telescope/parameters": NodeType.TEL_GROUP, "/r0/event/telescope": NodeType.TEL_GROUP, "/r1/event/telescope": NodeType.TEL_GROUP, + "/dl0/event/subarray/trigger": NodeType.TABLE, + "/dl0/event/telescope/trigger": NodeType.TABLE, "/dl1/event/subarray/trigger": NodeType.TABLE, "/dl1/event/telescope/trigger": NodeType.TABLE, "/dl1/event/telescope/images": NodeType.TEL_GROUP, @@ -48,6 +50,8 @@ class NodeType(enum.Enum): "/dl2/event/subarray": NodeType.ALGORITHM_GROUP, "/dl1/monitoring/subarray/pointing": NodeType.TABLE, "/dl1/monitoring/telescope/pointing": NodeType.TEL_GROUP, + "/dl0/monitoring/subarray/pointing": NodeType.TABLE, + "/dl0/monitoring/telescope/pointing": NodeType.TEL_GROUP, } @@ -319,13 +323,14 @@ def _append(self, other): self._append_table_group(other, other.root[key]) # DL1 - key = "/dl1/event/subarray/trigger" - if key in other.root: - self._append_table(other, other.root[key]) + for dl in ("dl0", "dl1"): + key = f"/{dl}/event/subarray/trigger" + if key in other.root: + self._append_table(other, other.root[key]) - key = "/dl1/event/telescope/trigger" - if self.telescope_events and key in other.root: - self._append_table(other, other.root[key]) + key = f"/{dl}/event/telescope/trigger" + if self.telescope_events and key in other.root: + self._append_table(other, other.root[key]) key = "/dl1/event/telescope/images" if self.telescope_events and self.dl1_images and key in other.root: @@ -353,13 +358,14 @@ def _append(self, other): self._append_table(other, table) # monitoring - key = "/dl1/monitoring/subarray/pointing" - if self.monitoring and key in other.root: - self._append_table(other, other.root[key]) + for dl in ("dl0", "dl1"): + key = f"/{dl}/monitoring/subarray/pointing" + if self.monitoring and key in other.root: + self._append_table(other, other.root[key]) - key = "/dl1/monitoring/telescope/pointing" - if self.monitoring and self.telescope_events and key in other.root: - self._append_table_group(other, other.root[key]) + key = f"/{dl}/monitoring/telescope/pointing" + if self.monitoring and self.telescope_events and key in other.root: + self._append_table_group(other, other.root[key]) # quality query statistics key = "/dl1/service/image_statistics" diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py index 70e8e93f81a..ee372075354 100644 --- a/ctapipe/io/simteleventsource.py +++ b/ctapipe/io/simteleventsource.py @@ -22,32 +22,38 @@ ) from ..calib.camera.gainselection import GainChannel, GainSelector from ..containers import ( - ArrayEventContainer, CoordinateFrameType, - EventIndexContainer, + DL0SubarrayContainer, + DL0TelescopeContainer, + DL1CameraCalibrationContainer, EventType, ObservationBlockContainer, ObservationBlockState, ObservingMode, PixelStatus, PixelStatusContainer, - PointingContainer, PointingMode, - R0CameraContainer, - R1CameraContainer, + R0TelescopeContainer, + R1TelescopeContainer, SchedulingBlockContainer, SchedulingBlockType, - SimulatedCameraContainer, - SimulatedEventContainer, SimulatedShowerContainer, SimulationConfigContainer, + SimulationSubarrayContainer, + SimulationTelescopeContainer, + SubarrayEventContainer, + SubarrayEventIndexContainer, + SubarrayPointingContainer, + SubarrayTriggerContainer, + TelescopeCalibrationContainer, + TelescopeEventContainer, + TelescopeEventIndexContainer, TelescopeImpactParameterContainer, + TelescopeMonitoringContainer, TelescopePointingContainer, TelescopeTriggerContainer, - TriggerContainer, ) from ..coordinates import CameraFrame, shower_impact_distance -from ..core import Map from ..core.provenance import Provenance from ..core.traits import Bool, ComponentName, Float, Integer, Undefined, UseEnum from ..instrument import ( @@ -750,25 +756,22 @@ def _generate_events(self): obs_id = self.obs_id - trigger = self._fill_trigger_info(array_event) - if trigger.event_type == EventType.SUBARRAY: + array_trigger, tel_trigger = self._fill_trigger_info(array_event) + if array_trigger.event_type == EventType.SUBARRAY: shower = self._fill_simulated_event_information(array_event) else: shower = None - data = ArrayEventContainer( - simulation=SimulatedEventContainer(shower=shower), - pointing=self._fill_array_pointing(), - index=EventIndexContainer(obs_id=obs_id, event_id=event_id), + event = SubarrayEventContainer( count=counter, - trigger=trigger, + simulation=SimulationSubarrayContainer(shower=shower), + dl0=DL0SubarrayContainer(trigger=array_trigger), + pointing=self._fill_array_pointing(), + index=SubarrayEventIndexContainer(obs_id=obs_id, event_id=event_id), ) - data.meta["origin"] = "hessio" - data.meta["input_url"] = self.input_url - data.meta["max_events"] = self.max_events - - telescope_events = array_event["telescope_events"] - tracking_positions = array_event["tracking_positions"] + event.meta["origin"] = "hessio" + event.meta["input_url"] = self.input_url + event.meta["max_events"] = self.max_events photoelectron_sums = array_event.get("photoelectron_sums") if photoelectron_sums is not None: @@ -778,28 +781,30 @@ def _generate_events(self): else: true_image_sums = np.full(self.n_telescopes_original, np.nan) - if data.simulation.shower is not None: + if event.simulation.shower is not None: # compute impact distances of the shower to the telescopes impact_distances = shower_impact_distance( - shower_geom=data.simulation.shower, subarray=self.subarray + shower_geom=event.simulation.shower, subarray=self.subarray ) else: impact_distances = np.full(len(self.subarray), np.nan) * u.m + telescope_events = array_event["telescope_events"] + tracking_positions = array_event["tracking_positions"] for tel_id, telescope_event in telescope_events.items(): adc_samples = telescope_event.get("adc_samples") if adc_samples is None: adc_samples = telescope_event["adc_sums"][:, :, np.newaxis] - n_gains, n_pixels, n_samples = adc_samples.shape + _n_gains, n_pixels, _n_samples = adc_samples.shape true_image = ( array_event.get("photoelectrons", {}) .get(tel_id - 1, {}) .get("photoelectrons", None) ) - if data.simulation is not None: - if data.simulation.shower is not None: + if event.simulation is not None: + if event.simulation.shower is not None: impact_container = TelescopeImpactParameterContainer( distance=impact_distances[ self.subarray.tel_index_array[tel_id] @@ -812,19 +817,17 @@ def _generate_events(self): prefix="true_impact", ) - data.simulation.tel[tel_id] = SimulatedCameraContainer( + simulation = SimulationTelescopeContainer( true_image_sum=true_image_sums[ self.telescope_indices_original[tel_id] ], true_image=true_image, impact=impact_container, ) + else: + simulation = None - data.pointing.tel[tel_id] = self._fill_event_pointing( - tracking_positions[tel_id] - ) - - data.r0.tel[tel_id] = R0CameraContainer(waveform=adc_samples) + r0 = R0TelescopeContainer(waveform=adc_samples) cam_mon = array_event["camera_monitorings"][tel_id] pedestal = cam_mon["pedestal"] / cam_mon["n_ped_slices"] @@ -832,7 +835,7 @@ def _generate_events(self): # fill dc_to_pe and pedestal_per_sample info into monitoring # container - mon = data.mon.tel[tel_id] + mon = TelescopeMonitoringContainer() mon.calibration.dc_to_pe = dc_to_pe mon.calibration.pedestal_per_sample = pedestal mon.pixel_status = self._fill_mon_pixels_status(tel_id) @@ -845,26 +848,43 @@ def _generate_events(self): self.calib_scale, self.calib_shift, ) - pixel_status = self._get_r1_pixel_status( tel_id=tel_id, selected_gain_channel=selected_gain_channel, ) - data.r1.tel[tel_id] = R1CameraContainer( - event_type=trigger.event_type, + r1 = R1TelescopeContainer( waveform=r1_waveform, selected_gain_channel=selected_gain_channel, pixel_status=pixel_status, + event_type=tel_trigger[tel_id].event_type, ) # get time_shift from laser calibration time_calib = array_event["laser_calibrations"][tel_id]["tm_calib"] pix_index = np.arange(n_pixels) - dl1_calib = data.calibration.tel[tel_id].dl1 - dl1_calib.time_shift = time_calib[selected_gain_channel, pix_index] + calibration = TelescopeCalibrationContainer( + dl1=DL1CameraCalibrationContainer( + time_shift=time_calib[selected_gain_channel, pix_index] + ) + ) - yield data + event.tel[tel_id] = TelescopeEventContainer( + index=TelescopeEventIndexContainer( + obs_id=event.index.obs_id, + event_id=event.index.event_id, + tel_id=tel_id, + ), + r0=r0, + r1=r1, + dl0=DL0TelescopeContainer(trigger=tel_trigger[tel_id]), + simulation=simulation, + mon=mon, + pointing=self._fill_event_pointing(tracking_positions[tel_id]), + calibration=calibration, + ) + + yield event def _get_r1_pixel_status(self, tel_id, selected_gain_channel): tel_desc = self.file_.telescope_descriptions[tel_id] @@ -951,7 +971,13 @@ def _fill_trigger_info(self, array_event): central_time = parse_simtel_time(trigger["gps_time"]) - tel = Map(TelescopeTriggerContainer) + array_trigger = SubarrayTriggerContainer( + event_type=event_type, + time=central_time, + tels_with_trigger=tels_with_trigger, + ) + + tel = {} for tel_id, time in zip( trigger["triggered_telescopes"], trigger["trigger_times"] ): @@ -983,25 +1009,21 @@ def _fill_trigger_info(self, array_event): n_trigger_pixels=n_trigger_pixels, trigger_pixels=trigger_pixels, ) - return TriggerContainer( - event_type=event_type, - time=central_time, - tels_with_trigger=tels_with_trigger, - tel=tel, - ) + + return array_trigger, tel def _fill_array_pointing(self): if self.file_.header["tracking_mode"] == 0: az, alt = self.file_.header["direction"] - return PointingContainer( - array_altitude=u.Quantity(alt, u.rad), - array_azimuth=u.Quantity(az, u.rad), + return SubarrayPointingContainer( + altitude=u.Quantity(alt, u.rad), + azimuth=u.Quantity(az, u.rad), ) else: ra, dec = self.file_.header["direction"] - return PointingContainer( - array_ra=u.Quantity(ra, u.rad), - array_dec=u.Quantity(dec, u.rad), + return SubarrayPointingContainer( + ra=u.Quantity(ra, u.rad), + dec=u.Quantity(dec, u.rad), ) def _parse_simulation_header(self): diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py index f7077969e8e..ca3cda7cb6f 100644 --- a/ctapipe/io/tableloader.py +++ b/ctapipe/io/tableloader.py @@ -23,7 +23,10 @@ PARAMETERS_GROUP = "/dl1/event/telescope/parameters" IMAGES_GROUP = "/dl1/event/telescope/images" MUON_GROUP = "/dl1/event/telescope/muon" -TRIGGER_TABLE = "/dl1/event/subarray/trigger" +TRIGGER_TABLE = "/dl0/event/subarray/trigger" +TRIGGER_TABLE_OLD = "/dl1/event/subarray/trigger" +TEL_TRIGGER_TABLE = "/dl0/event/telescope/trigger" +TEL_TRIGGER_TABLE_OLD = "/dl1/event/telescope/trigger" SHOWER_TABLE = "/simulation/event/subarray/shower" TRUE_IMAGES_GROUP = "/simulation/event/telescope/images" TRUE_PARAMETERS_GROUP = "/simulation/event/telescope/parameters" @@ -280,7 +283,19 @@ def __exit__(self, exc_type, exc_value, traceback): def __len__(self): """Number of subarray events in input file""" - return self.h5file.root[TRIGGER_TABLE].shape[0] + return self.h5file.root[self._trigger_table].shape[0] + + @property + def _trigger_table(self): + if TRIGGER_TABLE_OLD in self.h5file.root: + return TRIGGER_TABLE_OLD + return TRIGGER_TABLE + + @property + def _tel_trigger_table(self): + if TEL_TRIGGER_TABLE_OLD in self.h5file.root: + return TEL_TRIGGER_TABLE_OLD + return TEL_TRIGGER_TABLE def _read_telescope_table(self, group, tel_id, start=None, stop=None): key = f"{group}/tel_{tel_id:03d}" @@ -302,7 +317,7 @@ def _get_sort_index(self, start=None, stop=None): """ table = read_table( self.h5file, - TRIGGER_TABLE, + self._trigger_table, start=start, stop=stop, )[["obs_id", "event_id"]] @@ -377,7 +392,14 @@ def read_subarray_events(self, start=None, stop=None, keep_order=True): table: astropy.io.Table Table with primary index columns "obs_id" and "event_id". """ - table = read_table(self.h5file, TRIGGER_TABLE, start=start, stop=stop) + + table = read_table( + self.h5file, + self._trigger_table, + start=start, + stop=stop, + ) + if keep_order: self._add_index_if_needed(table) @@ -462,7 +484,7 @@ def _read_telescope_events_for_id(self, tel_id, start=None, stop=None): table = read_table( self.h5file, - "/dl1/event/telescope/trigger", + self._tel_trigger_table, condition=f"tel_id == {tel_id}", start=trigger_start, stop=trigger_stop, @@ -646,7 +668,7 @@ def _n_telescope_events(self): """ # we need to load the trigger table until "stop" to # know which telescopes participated in which events - table = self.h5file.root[TRIGGER_TABLE] + table = self.h5file.root[self._trigger_table] n_events = table.shape[0] n_telescopes = table.coldescrs["tels_with_trigger"].shape[0] diff --git a/ctapipe/io/tests/test_datawriter.py b/ctapipe/io/tests/test_datawriter.py index abb68cce42c..351a1f42618 100644 --- a/ctapipe/io/tests/test_datawriter.py +++ b/ctapipe/io/tests/test_datawriter.py @@ -28,37 +28,34 @@ def generate_dummy_dl2(event): algos = ["HillasReconstructor", "ImPACTReconstructor"] for algo in algos: - for tel_id in event.dl1.tel: - event.dl2.tel[tel_id].geometry[algo] = ReconstructedGeometryContainer( + for tel_id, tel_event in event.tel.items(): + tel_event.dl2.geometry[algo] = ReconstructedGeometryContainer( alt=70 * u.deg, az=120 * u.deg, prefix=f"{algo}_tel", ) - event.dl2.tel[tel_id].energy[algo] = ReconstructedEnergyContainer( + tel_event.dl2.energy[algo] = ReconstructedEnergyContainer( energy=10 * u.TeV, prefix=f"{algo}_tel", ) - event.dl2.tel[tel_id].classification[ - algo - ] = ParticleClassificationContainer( + tel_event.dl2.classification[algo] = ParticleClassificationContainer( prediction=0.9, prefix=f"{algo}_tel", ) - event.dl2.stereo.geometry[algo] = ReconstructedGeometryContainer( + event.dl2.geometry[algo] = ReconstructedGeometryContainer( alt=72 * u.deg, az=121 * u.deg, telescopes=[1, 2, 4], prefix=algo, ) - - event.dl2.stereo.energy[algo] = ReconstructedEnergyContainer( + event.dl2.energy[algo] = ReconstructedEnergyContainer( energy=10 * u.TeV, telescopes=[1, 2, 4], prefix=algo, ) - event.dl2.stereo.classification[algo] = ParticleClassificationContainer( + event.dl2.classification[algo] = ParticleClassificationContainer( prediction=0.9, telescopes=[1, 2, 4], prefix=algo, @@ -201,13 +198,13 @@ def test_roundtrip(tmpdir: Path): for event in EventSource(output_path): - for tel_id, dl1 in event.dl1.tel.items(): - original_image = events[event.count].dl1.tel[tel_id].image - read_image = dl1.image + for tel_id, tel_event in event.tel.items(): + original_image = events[event.count].tel[tel_id].dl1.image + read_image = tel_event.dl1.image assert np.allclose(original_image, read_image, atol=0.1) - original_peaktime = events[event.count].dl1.tel[tel_id].peak_time - read_peaktime = dl1.peak_time + original_peaktime = events[event.count].tel[tel_id].dl1.peak_time + read_peaktime = tel_event.dl1.peak_time assert np.allclose(original_peaktime, read_peaktime, atol=0.01) diff --git a/ctapipe/io/tests/test_event_source.py b/ctapipe/io/tests/test_event_source.py index c4df13b41e0..f977fa6068c 100644 --- a/ctapipe/io/tests/test_event_source.py +++ b/ctapipe/io/tests/test_event_source.py @@ -4,7 +4,7 @@ from traitlets import TraitError from traitlets.config.loader import Config -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import SubarrayEventContainer from ctapipe.core import Component from ctapipe.io import DataLevel, EventSource, SimTelEventSource from ctapipe.utils import get_dataset_path @@ -25,7 +25,7 @@ class DummyEventSource(EventSource): def _generator(self): for i in range(5): - yield ArrayEventContainer(count=i) + yield SubarrayEventContainer(count=i) @staticmethod def is_compatible(file_path): diff --git a/ctapipe/io/tests/test_hdf5.py b/ctapipe/io/tests/test_hdf5.py index 75f8a19b206..7612cad0e34 100644 --- a/ctapipe/io/tests/test_hdf5.py +++ b/ctapipe/io/tests/test_hdf5.py @@ -11,9 +11,9 @@ from ctapipe.containers import ( HillasParametersContainer, LeakageContainer, - R0CameraContainer, + R0TelescopeContainer, SimulatedShowerContainer, - TelEventIndexContainer, + TelescopeEventIndexContainer, ) from ctapipe.core.container import Container, Field from ctapipe.io import read_table @@ -25,7 +25,7 @@ def test_h5_file(tmp_path_factory): """Test hdf5 file with some tables for the reader tests""" path = tmp_path_factory.mktemp("hdf5") / "test.h5" - r0 = R0CameraContainer() + r0 = R0TelescopeContainer() shower = SimulatedShowerContainer() r0.waveform = np.random.uniform(size=(50, 10)) r0.meta["test_attribute"] = 3.14159 @@ -92,12 +92,12 @@ def test_append_container(tmp_path): with HDF5TableWriter(path, mode="w") as writer: for event_id in range(10): hillas = HillasParametersContainer() - index = TelEventIndexContainer(obs_id=1, event_id=event_id, tel_id=1) + index = TelescopeEventIndexContainer(obs_id=1, event_id=event_id, tel_id=1) writer.write("data", [index, hillas]) with HDF5TableWriter(path, mode="a") as writer: for event_id in range(10): - index = TelEventIndexContainer(obs_id=2, event_id=event_id, tel_id=1) + index = TelescopeEventIndexContainer(obs_id=2, event_id=event_id, tel_id=1) hillas = HillasParametersContainer() writer.write("data", [index, hillas]) @@ -381,8 +381,8 @@ def test_read_container(test_h5_file): # test supplying a single container as well as an # iterable with one entry only simtab = reader.read("/R0/sim_shower", (SimulatedShowerContainer,)) - r0tab1 = reader.read("/R0/tel_001", R0CameraContainer) - r0tab2 = reader.read("/R0/tel_002", R0CameraContainer) + r0tab1 = reader.read("/R0/tel_001", R0TelescopeContainer) + r0tab2 = reader.read("/R0/tel_002", R0TelescopeContainer) # read all 3 tables in sync for _ in range(3): diff --git a/ctapipe/io/tests/test_hdf5eventsource.py b/ctapipe/io/tests/test_hdf5eventsource.py index bae1eb1359b..4d7b6b23a62 100644 --- a/ctapipe/io/tests/test_hdf5eventsource.py +++ b/ctapipe/io/tests/test_hdf5eventsource.py @@ -77,14 +77,14 @@ def test_simulation_info(dl1_file): with HDF5EventSource(input_url=dl1_file) as source: for event in source: assert np.isfinite(event.simulation.shower.energy) - for tel in event.simulation.tel: - assert tel in event.simulation.tel - assert event.simulation.tel[tel].true_image is not None + for tel_event in event.tel.values(): + assert tel_event.simulation is not None + assert tel_event.simulation.true_image is not None reco_lons.append( - event.simulation.tel[tel].true_parameters.hillas.fov_lon.value + tel_event.simulation.true_parameters.hillas.fov_lon.value ) reco_concentrations.append( - event.simulation.tel[tel].true_parameters.concentration.core + tel_event.simulation.true_parameters.concentration.core ) assert not np.isnan(reco_lons).all() assert sum(np.isnan(reco_lons)) == sum(np.isnan(reco_concentrations)) @@ -94,8 +94,8 @@ def test_dl1_a_only_data(dl1_image_file): with HDF5EventSource(input_url=dl1_image_file) as source: assert source.datalevels == (DataLevel.DL1_IMAGES,) for event in source: - for tel in event.dl1.tel: - assert event.dl1.tel[tel].image.any() + for tel_event in event.tel.values(): + assert tel_event.dl1.image.any() def test_dl1_b_only_data(dl1_parameters_file): @@ -104,12 +104,12 @@ def test_dl1_b_only_data(dl1_parameters_file): with HDF5EventSource(input_url=dl1_parameters_file) as source: assert source.datalevels == (DataLevel.DL1_PARAMETERS,) for event in source: - for tel in event.dl1.tel: + for tel_event in event.tel.values(): reco_lons.append( - event.simulation.tel[tel].true_parameters.hillas.fov_lon.value + tel_event.simulation.true_parameters.hillas.fov_lon.value ) reco_concentrations.append( - event.simulation.tel[tel].true_parameters.concentration.core + tel_event.simulation.true_parameters.concentration.core ) assert not np.isnan(reco_lons).all() assert sum(np.isnan(reco_lons)) == sum(np.isnan(reco_concentrations)) @@ -120,14 +120,15 @@ def test_dl1_data(dl1_file): reco_concentrations = [] with HDF5EventSource(input_url=dl1_file) as source: for event in source: - for tel in event.dl1.tel: - assert event.dl1.tel[tel].image.any() + for tel_event in event.tel.values(): + assert tel_event.dl1.image.any() reco_lons.append( - event.simulation.tel[tel].true_parameters.hillas.fov_lon.value + tel_event.simulation.true_parameters.hillas.fov_lon.value ) reco_concentrations.append( - event.simulation.tel[tel].true_parameters.concentration.core + tel_event.simulation.true_parameters.concentration.core ) + assert not np.isnan(reco_lons).all() assert sum(np.isnan(reco_lons)) == sum(np.isnan(reco_concentrations)) @@ -135,23 +136,22 @@ def test_dl1_data(dl1_file): def test_pointing(dl1_file): with HDF5EventSource(input_url=dl1_file) as source: for event in source: - assert np.isclose(event.pointing.array_azimuth.to_value(u.deg), 0) - assert np.isclose(event.pointing.array_altitude.to_value(u.deg), 70) - assert event.pointing.tel - for tel in event.pointing.tel: - assert np.isclose(event.pointing.tel[tel].azimuth.to_value(u.deg), 0) - assert np.isclose(event.pointing.tel[tel].altitude.to_value(u.deg), 70) + assert np.isclose(event.pointing.azimuth.to_value(u.deg), 0) + assert np.isclose(event.pointing.altitude.to_value(u.deg), 70) + for tel_event in event.tel.values(): + assert np.isclose(tel_event.pointing.azimuth.to_value(u.deg), 0) + assert np.isclose(tel_event.pointing.altitude.to_value(u.deg), 70) def test_read_r1(r1_hdf5_file): - print(r1_hdf5_file) with HDF5EventSource(input_url=r1_hdf5_file) as source: e = None assert source.datalevels == (DataLevel.R1,) for e in source: - pass + for tel_event in e.tel.values(): + assert tel_event.r1.waveform is not None assert e is not None assert e.count == 3 @@ -165,8 +165,8 @@ def test_trigger_allowed_tels(dl1_proton_file): i = 0 for i, e in enumerate(s): assert e.count == i - assert set(e.trigger.tels_with_trigger) == e.trigger.tel.keys() - assert len(e.trigger.tels_with_trigger) > 1 + assert set(e.dl0.trigger.tels_with_trigger) == e.tel.keys() + assert len(e.dl0.trigger.tels_with_trigger) > 1 assert i == 1 @@ -182,18 +182,18 @@ def test_read_dl2(dl2_shower_geometry_file): ) e = next(iter(s)) - assert algorithm in e.dl2.stereo.geometry - assert e.dl2.stereo.geometry[algorithm].alt is not None - assert e.dl2.stereo.geometry[algorithm].az is not None - assert e.dl2.stereo.geometry[algorithm].telescopes is not None - assert e.dl2.stereo.geometry[algorithm].prefix == algorithm + assert algorithm in e.dl2.geometry + assert e.dl2.geometry[algorithm].alt is not None + assert e.dl2.geometry[algorithm].az is not None + assert e.dl2.geometry[algorithm].telescopes is not None + assert e.dl2.geometry[algorithm].prefix == algorithm - tel_mask = e.dl2.stereo.geometry[algorithm].telescopes + tel_mask = e.dl2.geometry[algorithm].telescopes tel_ids = s.subarray.tel_mask_to_tel_ids(tel_mask) for tel_id in tel_ids: - assert tel_id in e.dl2.tel - assert algorithm in e.dl2.tel[tel_id].impact - impact = e.dl2.tel[tel_id].impact[algorithm] + assert tel_id in e.tel + assert algorithm in e.tel[tel_id].dl2.impact + impact = e.tel[tel_id].dl2.impact[algorithm] assert impact.prefix == algorithm + "_tel_impact" assert impact.distance is not None @@ -202,7 +202,8 @@ def test_dl1_camera_frame(dl1_camera_frame_file): with HDF5EventSource(dl1_camera_frame_file) as s: tel_id = None for e in s: - for tel_id, dl1 in e.dl1.tel.items(): + for tel_id, tel_event in e.tel.items(): + dl1 = tel_event.dl1 assert isinstance( dl1.parameters.hillas, CameraHillasParametersContainer ) @@ -211,7 +212,7 @@ def test_dl1_camera_frame(dl1_camera_frame_file): ) assert dl1.parameters.hillas.intensity is not None - for tel_id, sim in e.simulation.tel.items(): + sim = tel_event.simulation assert isinstance( sim.true_parameters.hillas, CameraHillasParametersContainer ) diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py index 116a8a91f8a..9a291296ca0 100644 --- a/ctapipe/io/tests/test_simteleventsource.py +++ b/ctapipe/io/tests/test_simteleventsource.py @@ -139,9 +139,9 @@ def test_gamma_file_prod2(): for event in source: if event.count == 0: - assert event.r0.tel.keys() == {38, 47} + assert event.tel.keys() == {38, 47} elif event.count == 1: - assert event.r0.tel.keys() == {11, 21, 24, 26, 61, 63, 118, 119} + assert event.tel.keys() == {11, 21, 24, 26, 61, 63, 118, 119} else: break @@ -175,16 +175,20 @@ def test_pointing(): max_events=3, focal_length_choice="EQUIVALENT", ) as reader: - for e in reader: - assert np.isclose(e.pointing.array_altitude.to_value(u.deg), 70) - assert np.isclose(e.pointing.array_azimuth.to_value(u.deg), 0) - assert np.isnan(e.pointing.array_ra) - assert np.isnan(e.pointing.array_dec) - - # normal run, alle telescopes point to the array direction - for pointing in e.pointing.tel.values(): - assert u.isclose(e.pointing.array_azimuth, pointing.azimuth) - assert u.isclose(e.pointing.array_altitude, pointing.altitude) + for array_event in reader: + assert np.isclose(array_event.pointing.altitude.to_value(u.deg), 70) + assert np.isclose(array_event.pointing.azimuth.to_value(u.deg), 0) + assert np.isnan(array_event.pointing.ra) + assert np.isnan(array_event.pointing.dec) + + # normal run, all telescopes point to the array direction + for tel_event in array_event.tel.values(): + assert u.isclose( + tel_event.pointing.azimuth, array_event.pointing.azimuth + ) + assert u.isclose( + tel_event.pointing.altitude, array_event.pointing.altitude + ) def test_allowed_telescopes(): @@ -198,11 +202,8 @@ def test_allowed_telescopes(): ) as reader: assert not allowed_tels.symmetric_difference(reader.subarray.tel_ids) for event in reader: - assert set(event.r0.tel).issubset(allowed_tels) - assert set(event.r1.tel).issubset(allowed_tels) - assert set(event.dl0.tel).issubset(allowed_tels) - assert set(event.trigger.tels_with_trigger).issubset(allowed_tels) - assert set(event.pointing.tel).issubset(allowed_tels) + assert set(event.tel).issubset(allowed_tels) + assert set(event.dl0.trigger.tels_with_trigger).issubset(allowed_tels) def test_calibration_events(): @@ -231,7 +232,7 @@ def test_calibration_events(): for event, expected_type, expected_id in zip_longest( reader, expected_types, expected_ids ): - assert event.trigger.event_type is expected_type + assert event.dl0.trigger.event_type is expected_type assert event.index.event_id == expected_id @@ -243,11 +244,17 @@ def test_trigger_times(): t0 = Time("2020-05-06T15:30:00") t1 = Time("2020-05-06T15:40:00") - for event in source: - assert t0 <= event.trigger.time <= t1 - for tel_id, trigger in event.trigger.tel.items(): + for array_event in source: + assert t0 <= array_event.dl0.trigger.time <= t1 + for tel_event in array_event.tel.values(): # test single telescope events triggered within 50 ns - assert 0 <= (trigger.time - event.trigger.time).to_value(u.ns) <= 50 + assert ( + 0 + <= (tel_event.dl0.trigger.time - array_event.dl0.trigger.time).to_value( + u.ns + ) + <= 50 + ) def test_true_image(): @@ -256,8 +263,8 @@ def test_true_image(): focal_length_choice="EQUIVALENT", ) as reader: for event in reader: - for tel in event.simulation.tel.values(): - assert np.count_nonzero(tel.true_image) > 0 + for tel in event.tel.values(): + assert np.count_nonzero(tel.simulation.true_image) > 0 def test_instrument(): @@ -389,7 +396,7 @@ def test_calibscale_and_calibshift(prod5_gamma_simtel_path): event = next(iter(source)) # make sure we actually have data - assert len(event.r1.tel) > 0 + assert len(event.tel) > 0 calib_scale = 2.0 @@ -398,10 +405,10 @@ def test_calibscale_and_calibshift(prod5_gamma_simtel_path): ) as source: event_scaled = next(iter(source)) - for tel_id, r1 in event.r1.tel.items(): + for tel_id, tel_event in event.tel.items(): np.testing.assert_allclose( - r1.waveform[0], - event_scaled.r1.tel[tel_id].waveform[0] / calib_scale, + tel_event.r1.waveform[0], + event_scaled.tel[tel_id].r1.waveform[0] / calib_scale, rtol=0.1, ) @@ -412,10 +419,10 @@ def test_calibscale_and_calibshift(prod5_gamma_simtel_path): ) as source: event_shifted = next(iter(source)) - for tel_id, r1 in event.r1.tel.items(): + for tel_id, tel_event in event.tel.items(): np.testing.assert_allclose( - r1.waveform[0], - event_shifted.r1.tel[tel_id].waveform[0] - calib_shift, + tel_event.r1.waveform[0], + event_shifted.tel[tel_id].r1.waveform[0] - calib_shift, rtol=0.1, ) @@ -427,7 +434,7 @@ def test_true_image_sum(): focal_length_choice="EQUIVALENT", ) as s: e = next(iter(s)) - assert np.all(np.isnan(sim.true_image_sum) for sim in e.simulation.tel.values()) + assert np.all(np.isnan(tel.simulation.true_image_sum) for tel in e.tel.values()) with SimTelEventSource( calib_events_path, @@ -436,11 +443,11 @@ def test_true_image_sum(): e = next(iter(s)) true_image_sums = {} - for tel_id, sim_camera in e.simulation.tel.items(): + for tel_id, tel in e.tel.items(): # since the test file contains both sums and individual pixel values # we can compare. - assert sim_camera.true_image_sum == sim_camera.true_image.sum() - true_image_sums[tel_id] = sim_camera.true_image_sum + assert tel.simulation.true_image_sum == tel.simulation.true_image.sum() + true_image_sums[tel_id] = tel.simulation.true_image_sum # check it also works with allowed_tels, since the values # are stored in a flat array in simtel @@ -450,8 +457,8 @@ def test_true_image_sum(): focal_length_choice="EQUIVALENT", ) as s: e = next(iter(s)) - assert e.simulation.tel[2].true_image_sum == true_image_sums[2] - assert e.simulation.tel[3].true_image_sum == true_image_sums[3] + assert e.tel[2].simulation.true_image_sum == true_image_sums[2] + assert e.tel[3].simulation.true_image_sum == true_image_sums[3] def test_extracted_calibevents(): diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py index 864a53ba9f7..b4a5c272568 100644 --- a/ctapipe/io/tests/test_table_loader.py +++ b/ctapipe/io/tests/test_table_loader.py @@ -183,7 +183,7 @@ def test_table_loader_keeps_original_order(dl2_merged_file): from ctapipe.io.tableloader import TableLoader # check that the order is the same as in the file itself - trigger = read_table(dl2_merged_file, "/dl1/event/subarray/trigger") + trigger = read_table(dl2_merged_file, "/dl0/event/subarray/trigger") # check we actually have unsorted input assert not np.all(np.diff(trigger["obs_id"]) >= 0) @@ -278,7 +278,7 @@ def test_chunked(dl2_shower_geometry_file): """Test chunked reading""" from ctapipe.io.tableloader import TableLoader, read_table - trigger = read_table(dl2_shower_geometry_file, "/dl1/event/subarray/trigger") + trigger = read_table(dl2_shower_geometry_file, "/dl0/event/subarray/trigger") n_events = len(trigger) n_read = 0 @@ -415,6 +415,7 @@ def test_order_merged(): path = get_dataset_path("gamma_diffuse_dl2_train_small.dl2.h5") + # FIXME: old file, change to dl0 when updating trigger = read_table(path, "/dl1/event/subarray/trigger") tel_trigger = read_table(path, "/dl1/event/telescope/trigger") with TableLoader( diff --git a/ctapipe/io/tests/test_toysource.py b/ctapipe/io/tests/test_toysource.py index f8e58c4d439..f029938f003 100644 --- a/ctapipe/io/tests/test_toysource.py +++ b/ctapipe/io/tests/test_toysource.py @@ -29,7 +29,8 @@ def test_toyeventsource(subarray): for i, e in enumerate(s): assert e.index.event_id == i - for tel_id, dl1 in e.dl1.tel.items(): + for tel_id, tel_event in e.tel.items(): + dl1 = tel_event.dl1 assert dl1.image.size == subarray.tel[tel_id].camera.geometry.n_pixels assert (i + 1) == s.max_events diff --git a/ctapipe/io/toymodel.py b/ctapipe/io/toymodel.py index 05477f4e41b..bfc150552da 100644 --- a/ctapipe/io/toymodel.py +++ b/ctapipe/io/toymodel.py @@ -9,11 +9,11 @@ import numpy as np from ..containers import ( - ArrayEventContainer, - DL1CameraContainer, - EventIndexContainer, + DL1TelescopeContainer, ObservationBlockContainer, SchedulingBlockContainer, + SubarrayEventContainer, + SubarrayEventIndexContainer, ) from ..core import TelescopeComponent, traits from ..image import toymodel @@ -98,15 +98,12 @@ def _generator(self): def generate_event(self): - event = ArrayEventContainer( - index=EventIndexContainer(obs_id=1, event_id=self.event_id), - trigger=None, - r0=None, + event = SubarrayEventContainer( + index=SubarrayEventIndexContainer(obs_id=1, event_id=self.event_id), dl0=None, dl2=None, simulation=None, count=self.event_id, - calibration=None, ) for tel_id, telescope in self.subarray.tel.items(): @@ -150,6 +147,6 @@ def generate_event(self): ) image, _, _ = model.generate_image(cam, intensity) - event.dl1.tel[tel_id] = DL1CameraContainer(image=image) + event.tel[tel_id].dl1 = DL1TelescopeContainer(image=image) return event diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py index 0906cf0ff4b..3447ab91521 100644 --- a/ctapipe/reco/hillas_intersection.py +++ b/ctapipe/reco/hillas_intersection.py @@ -117,7 +117,7 @@ def __call__(self, event): Parameters ---------- - event : `~ctapipe.containers.ArrayEventContainer` + event : `~ctapipe.containers.SubarrayEventContainer` The event, needs to have dl1 parameters. Will be filled with the corresponding dl2 containers, reconstructed stereo geometry and telescope-wise impact position. @@ -126,20 +126,20 @@ def __call__(self, event): try: hillas_dict = self._create_hillas_dict(event) except (TooFewTelescopesException, InvalidWidthException): - event.dl2.stereo.geometry[self.__class__.__name__] = INVALID + event.dl2.geometry[self.__class__.__name__] = INVALID self._store_impact_parameter(event) return # Due to tracking the pointing of the array will never be a constant array_pointing = SkyCoord( - az=event.pointing.array_azimuth, - alt=event.pointing.array_altitude, + az=event.pointing.azimuth, + alt=event.pointing.altitude, frame=AltAz(), ) telescope_pointings = self._get_telescope_pointings(event) - event.dl2.stereo.geometry[self.__class__.__name__] = self._predict( + event.dl2.geometry[self.__class__.__name__] = self._predict( hillas_dict, array_pointing, telescope_pointings ) diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py index 0f82e2bc51d..68725debc43 100644 --- a/ctapipe/reco/hillas_reconstructor.py +++ b/ctapipe/reco/hillas_reconstructor.py @@ -135,7 +135,7 @@ def __call__(self, event): Parameters ---------- - event : `~ctapipe.containers.ArrayEventContainer` + event : `~ctapipe.containers.SubarrayEventContainer` The event, needs to have dl1 parameters. Will be filled with the corresponding dl2 containers, reconstructed stereo geometry and telescope-wise impact position. @@ -145,7 +145,7 @@ def __call__(self, event): try: hillas_dict = self._create_hillas_dict(event) except (TooFewTelescopesException, InvalidWidthException): - event.dl2.stereo.geometry[self.__class__.__name__] = INVALID + event.dl2.geometry[self.__class__.__name__] = INVALID self._store_impact_parameter(event) return @@ -166,7 +166,7 @@ def __call__(self, event): # store core corrected psi values for tel_id, psi in zip(tel_ids, corrected_psi): - event.dl1.tel[tel_id].parameters.core.psi = u.Quantity( + event.tel[tel_id].dl1.parameters.core.psi = u.Quantity( np.rad2deg(psi), u.deg ) @@ -190,9 +190,7 @@ def __call__(self, event): # az is clockwise, lon counter-clockwise, make sure it stays in [0, 2pi) az = Longitude(-lon) - event.dl2.stereo.geometry[ - self.__class__.__name__ - ] = ReconstructedGeometryContainer( + event.dl2.geometry[self.__class__.__name__] = ReconstructedGeometryContainer( alt=lat, az=az, core_x=core_pos_ground.x, @@ -234,10 +232,9 @@ def initialize_arrays(self, event, hillas_dict): # get one telescope id to check what frame to use altaz = AltAz() - # Due to tracking the pointing of the array will never be a constant array_pointing = SkyCoord( - az=event.pointing.array_azimuth, - alt=event.pointing.array_altitude, + az=event.pointing.azimuth, + alt=event.pointing.altitude, frame=altaz, ) @@ -258,7 +255,7 @@ def initialize_arrays(self, event, hillas_dict): for i, (tel_id, hillas) in enumerate(hillas_dict.items()): tel_ids[i] = tel_id - pointing = event.pointing.tel[tel_id] + pointing = event.tel[tel_id].pointing alt[i] = pointing.altitude.to_value(u.rad) az[i] = pointing.azimuth.to_value(u.rad) diff --git a/ctapipe/reco/preprocessing.py b/ctapipe/reco/preprocessing.py index c815dfb44ba..2822a164783 100644 --- a/ctapipe/reco/preprocessing.py +++ b/ctapipe/reco/preprocessing.py @@ -17,7 +17,7 @@ from ctapipe.coordinates import MissingFrameAttributeWarning, TelescopeFrame -from ..containers import ArrayEventContainer +from ..containers import SubarrayEventContainer LOG = logging.getLogger(__name__) @@ -69,13 +69,13 @@ def table_to_X(table: Table, features: List[str], log=LOG): def collect_features( - event: ArrayEventContainer, tel_id: int, subarray_table=None + event: SubarrayEventContainer, tel_id: int, subarray_table=None ) -> Table: """Loop over all containers with features. Parameters ---------- - event : ArrayEventContainer + event : SubarrayEventContainer The event container from which to collect the features tel_id : int The telscope id for which to collect the features @@ -88,10 +88,11 @@ def collect_features( """ features = {} - features.update(event.trigger.as_dict(recursive=False, flatten=True)) + features.update(event.dl0.trigger.as_dict(recursive=False, flatten=True)) + tel_event = event.tel[tel_id] features.update( - event.dl1.tel[tel_id].parameters.as_dict( + tel_event.dl1.parameters.as_dict( add_prefix=True, recursive=True, flatten=True, @@ -99,7 +100,7 @@ def collect_features( ) features.update( - event.dl2.tel[tel_id].as_dict( + tel_event.dl2.as_dict( add_prefix=True, recursive=True, flatten=True, @@ -108,7 +109,7 @@ def collect_features( ) features.update( - event.dl2.stereo.as_dict( + event.dl2.as_dict( add_prefix=True, recursive=True, flatten=True, diff --git a/ctapipe/reco/reconstructor.py b/ctapipe/reco/reconstructor.py index 28402449ee4..af32ae73020 100644 --- a/ctapipe/reco/reconstructor.py +++ b/ctapipe/reco/reconstructor.py @@ -6,7 +6,7 @@ import numpy as np from astropy.coordinates import AltAz, SkyCoord -from ctapipe.containers import ArrayEventContainer, TelescopeImpactParameterContainer +from ctapipe.containers import SubarrayEventContainer, TelescopeImpactParameterContainer from ctapipe.core import Provenance, QualityQuery, TelescopeComponent from ctapipe.core.traits import List @@ -78,7 +78,7 @@ def __init__(self, subarray, **kwargs): self.quality_query = StereoQualityQuery(parent=self) @abstractmethod - def __call__(self, event: ArrayEventContainer): + def __call__(self, event: SubarrayEventContainer): """ Perform stereo reconstruction on event. @@ -87,7 +87,7 @@ def __call__(self, event: ArrayEventContainer): Parameters ---------- - event : `ctapipe.containers.ArrayEventContainer` + event : `ctapipe.containers.SubarrayEventContainer` The event, needs to have dl1 parameters. Will be filled with the corresponding dl2 containers, reconstructed stereo geometry and telescope-wise impact position. @@ -149,9 +149,9 @@ class HillasGeometryReconstructor(Reconstructor): def _create_hillas_dict(self, event): hillas_dict = { - tel_id: dl1.parameters.hillas - for tel_id, dl1 in event.dl1.tel.items() - if all(self.quality_query(parameters=dl1.parameters)) + tel_id: tel_event.dl1.parameters.hillas + for tel_id, tel_event in event.tel.items() + if all(self.quality_query(parameters=tel_event.dl1.parameters)) } if len(hillas_dict) < 2: @@ -174,16 +174,17 @@ def _create_hillas_dict(self, event): def _get_telescope_pointings(event): return { tel_id: SkyCoord( - alt=event.pointing.tel[tel_id].altitude, - az=event.pointing.tel[tel_id].azimuth, + alt=tel_event.pointing.altitude, + az=tel_event.pointing.azimuth, frame=AltAz(), ) - for tel_id in event.dl1.tel.keys() + for tel_id, tel_event in event.tel.items() } def _store_impact_parameter(self, event): """Compute and store the impact parameter for each reconstruction.""" - geometry = event.dl2.stereo.geometry[self.__class__.__name__] + key = self.__class__.__name__ + geometry = event.dl2.geometry[key] if geometry.is_valid: impact_distances = shower_impact_distance( @@ -195,12 +196,10 @@ def _store_impact_parameter(self, event): impact_distances = u.Quantity(np.full(n_tels, np.nan), u.m) default_prefix = TelescopeImpactParameterContainer.default_prefix - prefix = f"{self.__class__.__name__}_tel_{default_prefix}" - for tel_id in event.trigger.tels_with_trigger: + prefix = f"{key}_tel_{default_prefix}" + for tel_id, tel_event in event.tel.items(): tel_index = self.subarray.tel_indices[tel_id] - event.dl2.tel[tel_id].impact[ - self.__class__.__name__ - ] = TelescopeImpactParameterContainer( + tel_event.dl2.impact[key] = TelescopeImpactParameterContainer( distance=impact_distances[tel_index], prefix=prefix, ) diff --git a/ctapipe/reco/shower_processor.py b/ctapipe/reco/shower_processor.py index af4011b19ad..b9e0b64e416 100644 --- a/ctapipe/reco/shower_processor.py +++ b/ctapipe/reco/shower_processor.py @@ -1,7 +1,7 @@ """ High level processing of showers. """ -from ..containers import ArrayEventContainer +from ..containers import SubarrayEventContainer from ..core import Component, traits from ..instrument import SubarrayDescription from .reconstructor import Reconstructor @@ -69,13 +69,13 @@ def __init__( for reco_type in self.reconstructor_types ] - def __call__(self, event: ArrayEventContainer): + def __call__(self, event: SubarrayEventContainer): """ Apply all configured stereo reconstructors to the given event. Parameters ---------- - event : ctapipe.containers.ArrayEventContainer + event : ctapipe.containers.SubarrayEventContainer Top-level container for all event information. """ for reconstructor in self.reconstructors: diff --git a/ctapipe/reco/sklearn.py b/ctapipe/reco/sklearn.py index 1ae0bc78158..b12fb6f457b 100644 --- a/ctapipe/reco/sklearn.py +++ b/ctapipe/reco/sklearn.py @@ -21,11 +21,11 @@ from ctapipe.exceptions import TooFewEvents from ..containers import ( - ArrayEventContainer, DispContainer, ParticleClassificationContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, + SubarrayEventContainer, ) from ..coordinates import TelescopeFrame from ..core import Component, FeatureGenerator, Provenance, QualityQuery, traits @@ -154,14 +154,14 @@ def __init__(self, subarray=None, models=None, **kwargs): self.prefix = self.model_cls @abstractmethod - def __call__(self, event: ArrayEventContainer) -> None: + def __call__(self, event: SubarrayEventContainer) -> None: """Event-wise prediction for the EventSource-Loop. Fills the event.dl2.[name] container. Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer """ @abstractmethod @@ -364,11 +364,11 @@ class EnergyRegressor(SKLearnRegressionReconstructor): target = "true_energy" property = ReconstructionProperty.ENERGY - def __call__(self, event: ArrayEventContainer) -> None: + def __call__(self, event: SubarrayEventContainer) -> None: """ Apply model for a single event and fill result into the event container """ - for tel_id in event.trigger.tels_with_trigger: + for tel_id, tel_event in event.tel.items(): table = collect_features(event, tel_id, self.instrument_table) table = self.feature_generator(table, subarray=self.subarray) @@ -391,7 +391,7 @@ def __call__(self, event: ArrayEventContainer) -> None: ) container.prefix = f"{self.prefix}_tel" - event.dl2.tel[tel_id].energy[self.prefix] = container + tel_event.dl2.energy[self.prefix] = container self.stereo_combiner(event) @@ -436,8 +436,8 @@ class ParticleClassifier(SKLearnClassificationReconstructor): property = ReconstructionProperty.PARTICLE_TYPE - def __call__(self, event: ArrayEventContainer) -> None: - for tel_id in event.trigger.tels_with_trigger: + def __call__(self, event: SubarrayEventContainer) -> None: + for tel_id, tel_event in event.tel.items(): table = collect_features(event, tel_id, self.instrument_table) table = self.feature_generator(table, subarray=self.subarray) passes_quality_checks = self.quality_query.get_table_mask(table)[0] @@ -458,7 +458,7 @@ def __call__(self, event: ArrayEventContainer) -> None: ) container.prefix = f"{self.prefix}_tel" - event.dl2.tel[tel_id].classification[self.prefix] = container + tel_event.dl2.classification[self.prefix] = container self.stereo_combiner(event) @@ -672,7 +672,7 @@ def _predict(self, key, table): return prediction, valid - def __call__(self, event: ArrayEventContainer) -> None: + def __call__(self, event: SubarrayEventContainer) -> None: """Event-wise prediction for the EventSource-Loop. Fills the event.dl2.tel[tel_id].disp[prefix] container @@ -680,9 +680,9 @@ def __call__(self, event: ArrayEventContainer) -> None: Parameters ---------- - event: ArrayEventContainer + event: SubarrayEventContainer """ - for tel_id in event.trigger.tels_with_trigger: + for tel_id, tel_event in event.tel.items(): table = collect_features(event, tel_id, self.instrument_table) table = self.feature_generator(table, subarray=self.subarray) @@ -696,7 +696,7 @@ def __call__(self, event: ArrayEventContainer) -> None: ) if valid: - hillas = event.dl1.tel[tel_id].parameters.hillas + hillas = tel_event.dl1.parameters.hillas psi = hillas.psi.to_value(u.rad) fov_lon = hillas.fov_lon + disp[0] * np.cos(psi) @@ -705,8 +705,8 @@ def __call__(self, event: ArrayEventContainer) -> None: fov_lon=fov_lon, fov_lat=fov_lat, telescope_pointing=AltAz( - alt=event.pointing.tel[tel_id].altitude, - az=event.pointing.tel[tel_id].azimuth, + alt=tel_event.pointing.altitude, + az=tel_event.pointing.azimuth, ), ).transform_to(AltAz()) @@ -733,8 +733,8 @@ def __call__(self, event: ArrayEventContainer) -> None: disp_container.prefix = f"{self.prefix}_parameter" altaz_container.prefix = f"{self.prefix}_tel" - event.dl2.tel[tel_id].disp[self.prefix] = disp_container - event.dl2.tel[tel_id].geometry[self.prefix] = altaz_container + tel_event.dl2.disp[self.prefix] = disp_container + tel_event.dl2.geometry[self.prefix] = altaz_container self.stereo_combiner(event) diff --git a/ctapipe/reco/stereo_combination.py b/ctapipe/reco/stereo_combination.py index f36e1230891..715e19b4f63 100644 --- a/ctapipe/reco/stereo_combination.py +++ b/ctapipe/reco/stereo_combination.py @@ -11,10 +11,10 @@ from ctapipe.reco.reconstructor import ReconstructionProperty from ..containers import ( - ArrayEventContainer, ParticleClassificationContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, + SubarrayEventContainer, ) from .utils import add_defaults_and_meta @@ -76,7 +76,7 @@ class StereoCombiner(Component): ).tag(config=True) @abstractmethod - def __call__(self, event: ArrayEventContainer) -> None: + def __call__(self, event: SubarrayEventContainer) -> None: """ Fill event container with stereo predictions """ @@ -151,15 +151,14 @@ def _combine_energy(self, event): values = [] weights = [] - for tel_id, dl2 in event.dl2.tel.items(): + for tel_id, tel_event in event.tel.items(): + dl2 = tel_event.dl2 mono = dl2.energy[self.prefix] if mono.is_valid: values.append(mono.energy.to_value(u.TeV)) - if tel_id not in event.dl1.tel: + if tel_event.dl1.parameters is None: raise ValueError("No parameters for weighting available") - weights.append( - self._calculate_weights(event.dl1.tel[tel_id].parameters) - ) + weights.append(self._calculate_weights(tel_event.dl1.parameters)) ids.append(tel_id) if len(values) > 0: @@ -181,7 +180,7 @@ def _combine_energy(self, event): mean = std = np.nan valid = False - event.dl2.stereo.energy[self.prefix] = ReconstructedEnergyContainer( + event.dl2.energy[self.prefix] = ReconstructedEnergyContainer( energy=u.Quantity(mean, u.TeV, copy=False), energy_uncert=u.Quantity(std, u.TeV, copy=False), telescopes=ids, @@ -194,12 +193,14 @@ def _combine_classification(self, event): values = [] weights = [] - for tel_id, dl2 in event.dl2.tel.items(): + for tel_id, tel_event in event.tel.items(): + dl2 = tel_event.dl2 mono = dl2.classification[self.prefix] if mono.is_valid: values.append(mono.prediction) - dl1 = event.dl1.tel[tel_id].parameters - weights.append(self._calculate_weights(dl1) if dl1 else 1) + if tel_event.dl1.parameters is None: + raise ValueError("No parameters for weighting available") + weights.append(self._calculate_weights(tel_event.dl1.parameters)) ids.append(tel_id) if len(values) > 0: @@ -212,7 +213,7 @@ def _combine_classification(self, event): container = ParticleClassificationContainer( prediction=mean, telescopes=ids, is_valid=valid, prefix=self.prefix ) - event.dl2.stereo.classification[self.prefix] = container + event.dl2.classification[self.prefix] = container def _combine_altaz(self, event): ids = [] @@ -220,13 +221,15 @@ def _combine_altaz(self, event): az_values = [] weights = [] - for tel_id, dl2 in event.dl2.tel.items(): + for tel_id, tel_event in event.tel.items(): + dl2 = tel_event.dl2 mono = dl2.geometry[self.prefix] if mono.is_valid: alt_values.append(mono.alt) az_values.append(mono.az) - dl1 = event.dl1.tel[tel_id].parameters - weights.append(self._calculate_weights(dl1) if dl1 else 1) + if tel_event.dl1.parameters is None: + raise ValueError("No parameters for weighting available") + weights.append(self._calculate_weights(tel_event.dl1.parameters)) ids.append(tel_id) if len(alt_values) > 0: # by construction len(alt_values) == len(az_values) @@ -254,7 +257,7 @@ def _combine_altaz(self, event): std = np.nan valid = False - event.dl2.stereo.geometry[self.prefix] = ReconstructedGeometryContainer( + event.dl2.geometry[self.prefix] = ReconstructedGeometryContainer( alt=mean_altaz.alt, az=mean_altaz.az, ang_distance_uncert=u.Quantity(np.rad2deg(std), u.deg, copy=False), @@ -263,7 +266,7 @@ def _combine_altaz(self, event): prefix=self.prefix, ) - def __call__(self, event: ArrayEventContainer) -> None: + def __call__(self, event: SubarrayEventContainer) -> None: """ Calculate the mean prediction for a single array event. """ diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py index fd8ca438495..25814105137 100644 --- a/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/ctapipe/reco/tests/test_HillasReconstructor.py @@ -88,40 +88,40 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev): hillas_reconstructor(event) # test the container is actually there and not only created by Map - assert "HillasReconstructor" in event.dl2.stereo.geometry - result = event.dl2.stereo.geometry["HillasReconstructor"] + assert "HillasReconstructor" in event.dl2.geometry + result = event.dl2.geometry["HillasReconstructor"] assert result.is_valid # copy event container to modify it invalid_event = deepcopy(original_event) # overwrite all image parameters but the last one with dummy ones - for tel_id in list(invalid_event.dl1.tel.keys())[:-1]: - invalid_event.dl1.tel[tel_id].parameters.hillas = HillasParametersContainer() + for tel_id in list(invalid_event.tel)[:-1]: + invalid_event.tel[tel_id].dl1.parameters.hillas = HillasParametersContainer() hillas_reconstructor(invalid_event) # test the container is actually there and not only created by Map - assert "HillasReconstructor" in invalid_event.dl2.stereo.geometry - result = invalid_event.dl2.stereo.geometry["HillasReconstructor"] + assert "HillasReconstructor" in invalid_event.dl2.geometry + result = invalid_event.dl2.geometry["HillasReconstructor"] assert not result.is_valid - tel_id = list(invalid_event.dl1.tel.keys())[-1] + tel_id = next(iter(invalid_event.tel)) # Now use the original event, but overwrite the last width to 0 invalid_event = deepcopy(original_event) - invalid_event.dl1.tel[tel_id].parameters.hillas.width = 0 * u.m + invalid_event.tel[tel_id].dl1.parameters.hillas.width = 0 * u.m hillas_reconstructor(invalid_event) # test the container is actually there and not only created by Map - assert "HillasReconstructor" in invalid_event.dl2.stereo.geometry - result = invalid_event.dl2.stereo.geometry["HillasReconstructor"] + assert "HillasReconstructor" in invalid_event.dl2.geometry + result = invalid_event.dl2.geometry["HillasReconstructor"] assert not result.is_valid # Now use the original event, but overwrite the last width to NaN invalid_event = deepcopy(original_event) - invalid_event.dl1.tel[tel_id].parameters.hillas.width = np.nan * u.m + invalid_event.tel[tel_id].dl1.parameters.hillas.width = np.nan * u.m hillas_reconstructor(invalid_event) # test the container is actually there and not only created by Map - assert "HillasReconstructor" in invalid_event.dl2.stereo.geometry - result = invalid_event.dl2.stereo.geometry["HillasReconstructor"] + assert "HillasReconstructor" in invalid_event.dl2.geometry + result = invalid_event.dl2.geometry["HillasReconstructor"] assert not result.is_valid @@ -146,8 +146,8 @@ def test_reconstruction_against_simulation_telescope_frame( image_processor(event) reconstructor(event) # test the container is actually there and not only created by Map - assert "HillasReconstructor" in event.dl2.stereo.geometry - result = event.dl2.stereo.geometry["HillasReconstructor"] + assert "HillasReconstructor" in event.dl2.geometry + result = event.dl2.geometry["HillasReconstructor"] # get the reconstructed coordinates in the sky reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz()) @@ -183,7 +183,7 @@ def test_reconstruction_against_simulation_camera_frame( calib(event) image_processor(event) reconstructor(event) - result = event.dl2.stereo.geometry[reconstructor.__class__.__name__] + result = event.dl2.geometry[reconstructor.__class__.__name__] # get the reconstructed coordinates in the sky reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz()) @@ -249,11 +249,9 @@ def test_CameraFrame_against_TelescopeFrame(filename): image_processor_camera_frame(event_camera_frame) reconstructor(event_camera_frame) - result_camera_frame = event_camera_frame.dl2.stereo.geometry[ - "HillasReconstructor" - ] + result_camera_frame = event_camera_frame.dl2.geometry["HillasReconstructor"] reconstructor(event_telescope_frame) - result_telescope_frame = event_telescope_frame.dl2.stereo.geometry[ + result_telescope_frame = event_telescope_frame.dl2.geometry[ "HillasReconstructor" ] diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py index 924699b74ed..d2568e7e474 100644 --- a/ctapipe/reco/tests/test_hillas_intersection.py +++ b/ctapipe/reco/tests/test_hillas_intersection.py @@ -273,7 +273,7 @@ def test_reconstruction_works(subarray_and_event_gamma_off_axis_500_gev): ) reconstructor(event) - result = event.dl2.stereo.geometry["HillasIntersection"] + result = event.dl2.geometry["HillasIntersection"] reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz()) assert reco_coord.separation(true_coord) < 0.1 * u.deg @@ -287,14 +287,15 @@ def test_selected_subarray(subarray_and_event_gamma_off_axis_500_gev): allowed_tels = {1, 4} for tel_id in subarray.tel.keys(): if tel_id not in allowed_tels: - event.dl1.tel.pop(tel_id, None) - event.trigger.tels_with_trigger = event.trigger.tels_with_trigger[ - event.trigger.tels_with_trigger != tel_id - ] + event.tel.pop(tel_id, None) + + event.dl0.trigger.tels_with_trigger = list( + set(event.dl0.trigger.tels_with_trigger).intersection(allowed_tels) + ) subarray = subarray.select_subarray(allowed_tels) reconstructor = HillasIntersection(subarray) reconstructor(event) - result = event.dl2.stereo.geometry["HillasIntersection"] + result = event.dl2.geometry["HillasIntersection"] assert result.is_valid diff --git a/ctapipe/reco/tests/test_preprocessing.py b/ctapipe/reco/tests/test_preprocessing.py index 861fb33404e..ef8b4c3651e 100644 --- a/ctapipe/reco/tests/test_preprocessing.py +++ b/ctapipe/reco/tests/test_preprocessing.py @@ -81,20 +81,20 @@ def test_collect_features(example_event, example_subarray): image_processor(event) shower_processor(event) - tel_id = next(iter(event.dl2.tel)) + tel_id, tel_event = next(iter(event.tel.items())) tab = collect_features(event, tel_id=tel_id) k = "HillasReconstructor" - impact = event.dl2.tel[tel_id].impact[k] + impact = tel_event.dl2.impact[k] assert tab[f"{k}_tel_impact_distance"].quantity[0] == impact.distance - geometry = event.dl2.stereo.geometry[k] + geometry = event.dl2.geometry[k] assert tab[f"{k}_az"].quantity[0] == geometry.az - hillas = event.dl1.tel[tel_id].parameters.hillas + hillas = tel_event.dl1.parameters.hillas assert tab["hillas_intensity"].quantity[0] == hillas.intensity - leakage = event.dl1.tel[tel_id].parameters.leakage + leakage = tel_event.dl1.parameters.leakage assert tab["leakage_intensity_width_1"].quantity[0] == leakage.intensity_width_1 tab = collect_features( diff --git a/ctapipe/reco/tests/test_reconstruction_methods.py b/ctapipe/reco/tests/test_reconstruction_methods.py index 73ceccfe635..3a9cf0d5c54 100644 --- a/ctapipe/reco/tests/test_reconstruction_methods.py +++ b/ctapipe/reco/tests/test_reconstruction_methods.py @@ -43,7 +43,7 @@ def test_reconstructors(reconstructors): name = ReconstructorType.__name__ # test the container is actually there and not only created by Map - assert name in event.dl2.stereo.geometry - assert event.dl2.stereo.geometry[name].alt.unit.is_equivalent(u.deg) - assert event.dl2.stereo.geometry[name].az.unit.is_equivalent(u.deg) - assert event.dl2.stereo.geometry[name].core_x.unit.is_equivalent(u.m) + assert name in event.dl2.geometry + assert event.dl2.geometry[name].alt.unit.is_equivalent(u.deg) + assert event.dl2.geometry[name].az.unit.is_equivalent(u.deg) + assert event.dl2.geometry[name].core_x.unit.is_equivalent(u.m) diff --git a/ctapipe/reco/tests/test_shower_processor.py b/ctapipe/reco/tests/test_shower_processor.py index 607568b810d..d786bab82f4 100644 --- a/ctapipe/reco/tests/test_shower_processor.py +++ b/ctapipe/reco/tests/test_shower_processor.py @@ -45,7 +45,7 @@ def test_shower_processor_geometry( process_shower(example_event_copy) for reco_type in reconstructor_types: - DL2a = example_event_copy.dl2.stereo.geometry[reco_type] + DL2a = example_event_copy.dl2.geometry[reco_type] print(DL2a) assert isfinite(DL2a.alt) assert isfinite(DL2a.az) @@ -71,7 +71,7 @@ def test_shower_processor_geometry( process_shower(example_event_copy) for reco_type in reconstructor_types: - DL2a = example_event_copy.dl2.stereo.geometry[reco_type] + DL2a = example_event_copy.dl2.geometry[reco_type] print(DL2a) assert not isfinite(DL2a.alt) assert not isfinite(DL2a.az) diff --git a/ctapipe/reco/tests/test_stereo_combination.py b/ctapipe/reco/tests/test_stereo_combination.py index 7187262af27..3b05cdbc333 100644 --- a/ctapipe/reco/tests/test_stereo_combination.py +++ b/ctapipe/reco/tests/test_stereo_combination.py @@ -5,13 +5,13 @@ from numpy.testing import assert_allclose, assert_array_equal from ctapipe.containers import ( - ArrayEventContainer, + DL2TelescopeContainer, HillasParametersContainer, ImageParametersContainer, ParticleClassificationContainer, - ReconstructedContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, + SubarrayEventContainer, ) from ctapipe.reco.reconstructor import ReconstructionProperty from ctapipe.reco.stereo_combination import StereoMeanCombiner @@ -152,10 +152,10 @@ def test_predict_mean_disp(mono_table): @pytest.mark.parametrize("weights", ["konrad", "intensity", "none"]) def test_mean_prediction_single_event(weights): - event = ArrayEventContainer() + event = SubarrayEventContainer() for tel_id, intensity in zip((25, 125, 130), (100, 200, 400)): - event.dl1.tel[tel_id].parameters = ImageParametersContainer( + event.tel[tel_id].dl1.parameters = ImageParametersContainer( hillas=HillasParametersContainer( intensity=intensity, width=0.1 * u.deg, @@ -163,7 +163,7 @@ def test_mean_prediction_single_event(weights): ) ) - event.dl2.tel[25] = ReconstructedContainer( + event.tel[25].dl2 = DL2TelescopeContainer( energy={ "dummy": ReconstructedEnergyContainer(energy=10 * u.GeV, is_valid=True) }, @@ -176,7 +176,7 @@ def test_mean_prediction_single_event(weights): ) }, ) - event.dl2.tel[125] = ReconstructedContainer( + event.tel[125].dl2 = DL2TelescopeContainer( energy={ "dummy": ReconstructedEnergyContainer(energy=20 * u.GeV, is_valid=True) }, @@ -189,7 +189,7 @@ def test_mean_prediction_single_event(weights): ) }, ) - event.dl2.tel[130] = ReconstructedContainer( + event.tel[130].dl2 = DL2TelescopeContainer( energy={ "dummy": ReconstructedEnergyContainer(energy=0.04 * u.TeV, is_valid=True) }, @@ -222,11 +222,11 @@ def test_mean_prediction_single_event(weights): combine_classification(event) combine_geometry(event) if weights == "none": - assert u.isclose(event.dl2.stereo.energy["dummy"].energy, (70 / 3) * u.GeV) - assert u.isclose(event.dl2.stereo.geometry["dummy"].alt, 63.0738383 * u.deg) - assert u.isclose(event.dl2.stereo.geometry["dummy"].az, 348.0716693 * u.deg) + assert u.isclose(event.dl2.energy["dummy"].energy, (70 / 3) * u.GeV) + assert u.isclose(event.dl2.geometry["dummy"].alt, 63.0738383 * u.deg) + assert u.isclose(event.dl2.geometry["dummy"].az, 348.0716693 * u.deg) elif weights == "intensity": - assert u.isclose(event.dl2.stereo.energy["dummy"].energy, 30 * u.GeV) - assert u.isclose(event.dl2.stereo.geometry["dummy"].alt, 60.9748605 * u.deg) - assert u.isclose(event.dl2.stereo.geometry["dummy"].az, 316.0365515 * u.deg) - assert event.dl2.stereo.classification["dummy"].prediction == 0.6 + assert u.isclose(event.dl2.energy["dummy"].energy, 30 * u.GeV) + assert u.isclose(event.dl2.geometry["dummy"].alt, 60.9748605 * u.deg) + assert u.isclose(event.dl2.geometry["dummy"].az, 316.0365515 * u.deg) + assert event.dl2.classification["dummy"].prediction == 0.6 diff --git a/ctapipe/tools/apply_models.py b/ctapipe/tools/apply_models.py index 2df691d8bf4..32c6ddcc7c6 100644 --- a/ctapipe/tools/apply_models.py +++ b/ctapipe/tools/apply_models.py @@ -227,9 +227,12 @@ def _combine(self, reconstructor, tel_tables, start, stop): # potentially changes the order of the subarray events. # to ensure events are stored in the correct order, # we resort to trigger table order - trigger = read_table( - self.h5file, "/dl1/event/subarray/trigger", start=start, stop=stop - )[["obs_id", "event_id"]] + trigger_table = "/dl1/event/subarray/trigger" + if trigger_table not in self.h5file.root: + trigger_table = "/dl0/event/subarray/trigger" + trigger = read_table(self.h5file, trigger_table, start=start, stop=stop)[ + ["obs_id", "event_id"] + ] trigger["__sort_index__"] = np.arange(len(trigger)) stereo_predictions = _join_subarray_events(trigger, stereo_predictions) diff --git a/ctapipe/tools/display_dl1.py b/ctapipe/tools/display_dl1.py index 7a27305b60b..30baab8dc0a 100644 --- a/ctapipe/tools/display_dl1.py +++ b/ctapipe/tools/display_dl1.py @@ -72,8 +72,8 @@ def _init_figure(self): self.pdf = self._exit_stack.enter_context(PdfPages(self.output_path)) def plot(self, event, tel_id): - image = event.dl1.tel[tel_id].image - peak_time = event.dl1.tel[tel_id].peak_time + image = event.tel[tel_id].dl1.image + peak_time = event.tel[tel_id].dl1.peak_time if self._current_tel != tel_id: self._current_tel = tel_id @@ -192,16 +192,15 @@ def start(self): for event in self.eventsource: self.calibrator(event) - tel_list = event.dl1.tel.keys() - - tel_list = event.dl1.tel.keys() + tels = event.tel.keys() if self.telescope: - if self.telescope not in tel_list: + if self.telescope not in tels: continue - tel_list = [self.telescope] - for tel_id in tel_list: - if all(self.quality_query(dl1=event.dl1.tel[tel_id])): + tels = [self.telescope] + + for tel_id in tels: + if all(self.quality_query(dl1=event.tel[tel_id].dl1)): self.plotter.plot(event, tel_id) def finish(self): diff --git a/ctapipe/tools/tests/test_apply_models.py b/ctapipe/tools/tests/test_apply_models.py index 32b0e4bcc86..df49866af1b 100644 --- a/ctapipe/tools/tests/test_apply_models.py +++ b/ctapipe/tools/tests/test_apply_models.py @@ -2,10 +2,10 @@ import pytest from ctapipe.containers import ( - EventIndexContainer, ParticleClassificationContainer, ReconstructedEnergyContainer, ReconstructedGeometryContainer, + SubarrayEventIndexContainer, ) from ctapipe.core import run_tool from ctapipe.core.tool import ToolConfigurationError @@ -41,7 +41,10 @@ def test_apply_energy_regressor( prefix = "ExtraTreesRegressor" table = read_table(output_path, f"/dl2/event/subarray/energy/{prefix}") for col in "obs_id", "event_id": - assert table[col].description == EventIndexContainer.fields[col].description + assert ( + table[col].description + == SubarrayEventIndexContainer.fields[col].description + ) for name, field in ReconstructedEnergyContainer.fields.items(): colname = f"{prefix}_{name}" @@ -70,7 +73,7 @@ def test_apply_energy_regressor( assert f"{prefix}_tel_is_valid" in tel_events.colnames assert "hillas_intensity" in tel_events.colnames - trigger = read_table(output_path, "/dl1/event/subarray/trigger") + trigger = read_table(output_path, "/dl0/event/subarray/trigger") energy = read_table(output_path, "/dl2/event/subarray/energy/ExtraTreesRegressor") check_equal_array_event_order(trigger, energy) @@ -155,6 +158,7 @@ def test_apply_all( # test file is produced using 0.17, the descriptions don't match # assert table[colname].description == field.description + # FIXME: old file, change to dl0 when updating test file trigger = read_table(output_path, "/dl1/event/subarray/trigger") subarray_tables = ( @@ -167,6 +171,8 @@ def test_apply_all( check_equal_array_event_order(trigger, table) subarray = SubarrayDescription.from_hdf(input_path) + + # FIXME: old file, change to dl0 when updating test file tel_trigger = read_table(output_path, "/dl1/event/telescope/trigger") for tel_id in subarray.tel: tel_keys = ( diff --git a/ctapipe/tools/tests/test_process.py b/ctapipe/tools/tests/test_process.py index 01039649265..2f18d29589c 100644 --- a/ctapipe/tools/tests/test_process.py +++ b/ctapipe/tools/tests/test_process.py @@ -90,20 +90,22 @@ def test_stage_1_dl1(tmp_path, dl1_image_file, dl1_parameters_file): # check tables were written with tables.open_file(dl1b_from_dl1a_file, mode="r") as testfile: + assert testfile.root.dl0 + assert testfile.root.dl0.event.subarray + assert testfile.root.dl0.event.telescope assert testfile.root.dl1 assert testfile.root.dl1.event.telescope - assert testfile.root.dl1.event.subarray assert testfile.root.configuration.instrument.subarray.layout assert testfile.root.configuration.instrument.telescope.optics assert testfile.root.configuration.instrument.telescope.camera.geometry_0 assert testfile.root.configuration.instrument.telescope.camera.readout_0 - assert testfile.root.dl1.monitoring.subarray.pointing.dtype.names == ( + assert testfile.root.dl0.monitoring.subarray.pointing.dtype.names == ( "time", - "array_azimuth", - "array_altitude", - "array_ra", - "array_dec", + "azimuth", + "altitude", + "ra", + "dec", ) # check we can read telescope parameters @@ -457,7 +459,7 @@ def test_muon_reconstruction_simtel(tmp_path): completeness = table["muonparameters_completeness"] for event in source: - muon = event.muon.tel[1] + muon = event.tel[1].muon assert u.isclose(muon.ring.radius, radius[event.count], equal_nan=True) assert np.isclose( muon.parameters.completeness, completeness[event.count], equal_nan=True diff --git a/ctapipe/tools/tests/test_process_ml.py b/ctapipe/tools/tests/test_process_ml.py index a4fb75d428d..f21f3742941 100644 --- a/ctapipe/tools/tests/test_process_ml.py +++ b/ctapipe/tools/tests/test_process_ml.py @@ -106,7 +106,7 @@ def test_process_apply_classification( events = read_table( output, "/dl2/event/subarray/classification/ExtraTreesClassifier" ) - trigger = read_table(output, "/dl1/event/subarray/trigger") + trigger = read_table(output, "/dl0/event/subarray/trigger") assert len(events) == len(trigger) assert "ExtraTreesClassifier_is_valid" in events.colnames assert "ExtraTreesClassifier_prediction" in events.colnames @@ -170,7 +170,7 @@ def test_process_apply_disp( assert "disp_tel_is_valid" in tel_events.colnames events = read_table(output, "/dl2/event/subarray/geometry/disp") - trigger = read_table(output, "/dl1/event/subarray/trigger") + trigger = read_table(output, "/dl0/event/subarray/trigger") assert len(events) == len(trigger) assert "disp_alt" in events.colnames assert "disp_az" in events.colnames diff --git a/ctapipe/utils/event_type_filter.py b/ctapipe/utils/event_type_filter.py index 6369608a4d4..81db188a292 100644 --- a/ctapipe/utils/event_type_filter.py +++ b/ctapipe/utils/event_type_filter.py @@ -31,4 +31,4 @@ def __call__(self, event): if self.allowed_types is None: return True - return event.trigger.event_type in self.allowed_types + return event.dl0.trigger.event_type in self.allowed_types diff --git a/ctapipe/utils/tests/test_event_filter.py b/ctapipe/utils/tests/test_event_filter.py index 6a0fc5990e8..a70e93551d2 100644 --- a/ctapipe/utils/tests/test_event_filter.py +++ b/ctapipe/utils/tests/test_event_filter.py @@ -1,6 +1,6 @@ from traitlets.config import Config -from ctapipe.containers import ArrayEventContainer, EventType +from ctapipe.containers import EventType, SubarrayEventContainer def test_event_filter(): @@ -10,12 +10,12 @@ def test_event_filter(): allowed_types={EventType.SUBARRAY, EventType.FLATFIELD} ) - e = ArrayEventContainer() - e.trigger.event_type = EventType.SUBARRAY + e = SubarrayEventContainer() + e.dl0.trigger.event_type = EventType.SUBARRAY assert event_filter(e) - e.trigger.event_type = EventType.FLATFIELD + e.dl0.trigger.event_type = EventType.FLATFIELD assert event_filter(e) - e.trigger.event_type = EventType.DARK_PEDESTAL + e.dl0.trigger.event_type = EventType.DARK_PEDESTAL assert not event_filter(e) @@ -25,9 +25,9 @@ def test_event_filter_none(): event_filter = EventTypeFilter(allowed_types=None) # all event types should pass - e = ArrayEventContainer() + e = SubarrayEventContainer() for value in EventType: - e.trigger.event_type = value + e.dl0.trigger.event_type = value assert event_filter(e) @@ -53,9 +53,9 @@ def test_event_filter_config(): EventType.SINGLE_PE, } - e = ArrayEventContainer() - e.trigger.event_type = EventType.DARK_PEDESTAL + e = SubarrayEventContainer() + e.dl0.trigger.event_type = EventType.DARK_PEDESTAL assert not event_filter(e) - e.trigger.event_type = EventType.SUBARRAY + e.dl0.trigger.event_type = EventType.SUBARRAY assert event_filter(e) diff --git a/ctapipe/visualization/mpl_array.py b/ctapipe/visualization/mpl_array.py index 78d1ad62adf..5201ff626a8 100644 --- a/ctapipe/visualization/mpl_array.py +++ b/ctapipe/visualization/mpl_array.py @@ -279,7 +279,7 @@ def set_vector_hillas( time_gradient: Dict[int, value of time gradient (no units)] dictionary for value of the time gradient for each telescope angle_offset: Float - This should be the ``event.pointing.array_azimuth`` parameter + This should be the ``event.pointing.azimuth`` parameter """ diff --git a/ctapipe/visualization/tests/test_bokeh.py b/ctapipe/visualization/tests/test_bokeh.py index 2e1dc1c8711..98b9c476fa6 100644 --- a/ctapipe/visualization/tests/test_bokeh.py +++ b/ctapipe/visualization/tests/test_bokeh.py @@ -12,17 +12,17 @@ def test_create_display_without_geometry(example_event, example_subarray): # test we can create it without geometry, and then set all the stuff display = CameraDisplay() - tel_id = next(iter(example_event.r0.tel.keys())) + tel_id, tel_event = next(iter(example_event.tel.items())) display.geometry = example_subarray.tel[tel_id].camera.geometry - display.image = example_event.dl1.tel[tel_id].image + display.image = tel_event.dl1.image def test_camera_display_creation(example_event, example_subarray): """Test we can create a display and check the resulting pixel coordinates""" from ctapipe.visualization.bokeh import CameraDisplay - t = list(example_event.r0.tel.keys())[0] - geom = example_subarray.tel[t].camera.geometry + tel_id = next(iter(example_event.tel)) + geom = example_subarray.tel[tel_id].camera.geometry display = CameraDisplay(geom) assert np.allclose(np.mean(display.datasource.data["xs"], axis=1), geom.pix_x.value) @@ -33,8 +33,8 @@ def test_camera_display_telescope_frame(example_event, example_subarray): """Test we can create a display in telescope frame""" from ctapipe.visualization.bokeh import CameraDisplay - t = list(example_event.r0.tel.keys())[0] - geom = example_subarray.tel[t].camera.geometry.transform_to(TelescopeFrame()) + tel_id = next(iter(example_event.tel)) + geom = example_subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame()) display = CameraDisplay(geom) assert np.allclose(np.mean(display.datasource.data["xs"], axis=1), geom.pix_x.value) @@ -45,8 +45,8 @@ def test_camera_image(example_event, example_subarray, tmp_path): """Test we set an image""" from ctapipe.visualization.bokeh import CameraDisplay - t = list(example_event.r0.tel.keys())[0] - geom = example_subarray.tel[t].camera.geometry + tel_id = next(iter(example_event.tel)) + geom = example_subarray.tel[tel_id].camera.geometry image = np.ones(geom.n_pixels) display = CameraDisplay(geom, image) @@ -64,8 +64,8 @@ def test_camera_enable_pixel_picker(example_event, example_subarray): """Test we can call enable_pixel_picker""" from ctapipe.visualization.bokeh import CameraDisplay - t = list(example_event.r0.tel.keys())[0] - geom = example_subarray.tel[t].camera.geometry + tel_id = next(iter(example_event.tel)) + geom = example_subarray.tel[tel_id].camera.geometry n_pixels = geom.pix_x.value.size image = np.ones(n_pixels) c_display = CameraDisplay(geom, image) diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 4cc0983560e..47a6df0de34 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -14,6 +14,7 @@ from ctapipe.containers import ( CameraHillasParametersContainer, HillasParametersContainer, + TelescopeEventContainer, ) from ctapipe.coordinates.telescope_frame import TelescopeFrame from ctapipe.instrument import PixelShape, SubarrayDescription @@ -157,11 +158,10 @@ def test_camera_display_multiple(prod5_lst_cam, tmp_path): def test_array_display(prod5_mst_nectarcam, reference_location): """check that we can do basic array display functionality""" from ctapipe.containers import ( - ArrayEventContainer, CoreParametersContainer, - DL1CameraContainer, - DL1Container, + DL1TelescopeContainer, ImageParametersContainer, + SubarrayEventContainer, ) from ctapipe.image import timing_parameters from ctapipe.visualization.mpl_array import ArrayDisplay @@ -182,15 +182,15 @@ def test_array_display(prod5_mst_nectarcam, reference_location): # Create a fake event containing telescope-wise information about # the image directions projected on the ground - event = ArrayEventContainer() - event.dl1 = DL1Container() - event.dl1.tel = {1: DL1CameraContainer(), 2: DL1CameraContainer()} - event.dl1.tel[1].parameters = ImageParametersContainer() - event.dl1.tel[2].parameters = ImageParametersContainer() - event.dl1.tel[2].parameters.core = CoreParametersContainer() - event.dl1.tel[1].parameters.core = CoreParametersContainer() - event.dl1.tel[1].parameters.core.psi = u.Quantity(2.0, unit=u.deg) - event.dl1.tel[2].parameters.core.psi = u.Quantity(1.0, unit=u.deg) + event = SubarrayEventContainer() + for tel_id, psi in zip((1, 2), (2 * u.deg, 1 * u.deg)): + event.tel[tel_id] = TelescopeEventContainer( + dl1=DL1TelescopeContainer( + parameters=ImageParametersContainer( + core=CoreParametersContainer(psi=psi), + ) + ) + ) ad = ArrayDisplay(subarray=sub) ad.set_vector_rho_phi(1 * u.m, 90 * u.deg) @@ -230,7 +230,8 @@ def test_array_display(prod5_mst_nectarcam, reference_location): ) gradient_dict = {1: timing_rot20.slope.value, 2: timing_rot20.slope.value} core_dict = { - tel_id: dl1.parameters.core.psi for tel_id, dl1 in event.dl1.tel.items() + tel_id: tel_event.dl1.parameters.core.psi + for tel_id, tel_event in event.tel.items() } ad.set_vector_hillas( hillas_dict=hillas_dict, @@ -346,8 +347,8 @@ def test_overlay_coord(tmp_path, subarray_and_event_gamma_off_axis_500_gev): calib(event) pointing = AltAz( - alt=event.pointing.array_altitude, - az=event.pointing.array_azimuth, + alt=event.pointing.altitude, + az=event.pointing.azimuth, ) # add pointing here, so the transform to CameraFrame / TelescopeFrame works @@ -359,7 +360,7 @@ def test_overlay_coord(tmp_path, subarray_and_event_gamma_off_axis_500_gev): ) geometry = subarray.tel[1].camera.geometry - image = event.dl1.tel[1].image + image = event.tel[1].dl1.image fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5), layout="constrained") diff --git a/docs/api-reference/containers/index.rst b/docs/api-reference/containers/index.rst index 5be56a863e5..15cbdfc15f1 100644 --- a/docs/api-reference/containers/index.rst +++ b/docs/api-reference/containers/index.rst @@ -13,7 +13,7 @@ The `ctapipe.containers` module contains the data model definition of all ctapipe `~ctapipe.core.Container` classes, which provide the container definitions for all ctapipe data. -The base Container for an event is in `ctapipe.containers.ArrayEventContainer`. +The base Container for an event is in `ctapipe.containers.SubarrayEventContainer`. Reference/API diff --git a/docs/api-reference/io/index.rst b/docs/api-reference/io/index.rst index 82f18d7c908..73bdbfe7b9d 100644 --- a/docs/api-reference/io/index.rst +++ b/docs/api-reference/io/index.rst @@ -131,7 +131,7 @@ Writing Output Files: ===================== The `DataWriter` Component allows one to write a series of events (stored in -`ctapipe.containers.ArrayEventContainer`) to a standardized HDF5 format file +`ctapipe.containers.SubarrayEventContainer`) to a standardized HDF5 format file following the data model (see :ref:`datamodels`). This includes all related datasets such as the instrument and simulation configuration information, simulated shower and image information, observed images and parameters and reconstruction diff --git a/docs/user-guide/data_models/dl1.rst b/docs/user-guide/data_models/dl1.rst index f792d495074..d0f8589b2fe 100644 --- a/docs/user-guide/data_models/dl1.rst +++ b/docs/user-guide/data_models/dl1.rst @@ -35,16 +35,16 @@ The following datasets will be written to the group ``/dl1/event/`` in the outp - (group) * - ``subarray/trigger`` - subarray trigger information - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.TriggerContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.SubarrayTriggerContainer` * - ``telescope/`` - Per-telescope Per-event information - (group) * - ``telescope/parameters/tel_{TEL_ID:03d}`` - tables of image parameters (one per telescope) - - :py:class:`~ctapipe.containers.TelEventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` + - :py:class:`~ctapipe.containers.TelescopeEventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` * - ``telescope/images/tel_{TEL_ID:03d}`` - tables of telescope images (one per telescope) - - :py:class:`~ctapipe.containers.TelEventIndexContainer` +, :py:class:`~ctapipe.containers.DL1CameraContainer` + - :py:class:`~ctapipe.containers.TelescopeEventIndexContainer` +, :py:class:`~ctapipe.containers.DL1TelescopeContainer` -------------- DL2 Data Model @@ -65,26 +65,26 @@ output file, where ```` is the identifier of the algorithm (e.g. - Contents * - /geometry - shower geometry reconstruction - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedGeometryContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedGeometryContainer` * - /energy - shower energy reconstruction - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedEnergyContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer`, :py:class:`~ctapipe.containers.ReconstructedEnergyContainer` * - /classification - shower classification parameters - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ParticleClassificationContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer`, :py:class:`~ctapipe.containers.ParticleClassificationContainer` --------------------- Simulation Data Model --------------------- * - ``/simulation/event/subarray/shower`` - true shower parameters from Monte-Carlo simulation - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedShowerContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedShowerContainer` * - ``/simulation/event/telescope/images/tel_{TEL_ID:03d}`` - simulated camera images - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.SimulatedCameraContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.SimulationTelescopeContainer` * - ``/simulation/event/telescope/parameters/tel_{TEL_ID:03d}`` - Parameters derived form the simulated camera images - - :py:class:`~ctapipe.containers.EventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` + - :py:class:`~ctapipe.containers.SubarrayEventIndexContainer` +, :py:class:`~ctapipe.containers.ImageParametersContainer` * - ``/simulation/service/shower_distribution`` - simulated shower distribution histograms - :py:class:`~ctapipe.containers.SimulatedShowerDistribution` diff --git a/test_plugin/ctapipe_test_plugin/__init__.py b/test_plugin/ctapipe_test_plugin/__init__.py index 2c06f79619f..27edc2395b1 100644 --- a/test_plugin/ctapipe_test_plugin/__init__.py +++ b/test_plugin/ctapipe_test_plugin/__init__.py @@ -20,7 +20,7 @@ TelescopeDescription, ) from ctapipe.io import DataLevel, EventSource -from ctapipe.io.datawriter import ArrayEventContainer +from ctapipe.io.datawriter import SubarrayEventContainer from ctapipe.reco import Reconstructor __all__ = [ @@ -90,7 +90,7 @@ def is_compatible(cls, path): def _generator(self): """Foo""" for i in range(10): - yield ArrayEventContainer(count=i) + yield SubarrayEventContainer(count=i) class PluginReconstructor(Reconstructor): @@ -100,6 +100,6 @@ class PluginReconstructor(Reconstructor): True, help="example traitlet to see that it is included in --help" ).tag(config=True) - def __call__(self, event: ArrayEventContainer): + def __call__(self, event: SubarrayEventContainer): """Foo""" event.dl2.geometry["PluginReconstructor"] = ReconstructedGeometryContainer()