Implement CEP-2: new event container structure

docs/api-reference/containers/index.rst

@@ -14,7 +14,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

docs/api-reference/io/index.rst

@@ -133,7 +133,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

docs/changes/2204.api.rst

@@ -0,0 +1,21 @@
+CEP-002 has been implemented. This is a major API breaking change,
+as the whole event container structure of ctapipe has been updated.
+
+The main change is the inversion of the hierarchy of telescopes
+and data levels.
+
+Where e.g. before the data level came first and then below that
+were mappings for each telescope:
+
+.. code:: python
+
+   event.dl1.tel[tel_id].parameters.hillas
+
+will now become accessing first the telescope event,
+and then the dl1 container inside:
+
+.. code:: python
+
+   event.tel[tel_id].dl1.parameters.hillas
+
+Fore more details, see :ref:`cep-002`.

docs/user-guide/data_models/dl1.rst

@@ -34,16 +34,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
@@ -64,27 +64,27 @@ output file, where ``<algorithm>`` 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`

examples/tutorials/calibrated_data_exploration.py

@@ -1,7 +1,6 @@
 """
 Explore Calibrated Data
 =======================
-
 """
 
 import numpy as np
@@ -12,7 +11,6 @@
 from ctapipe.calib import CameraCalibrator
 from ctapipe.image import hillas_parameters, tailcuts_clean
 from ctapipe.io import EventSource
-from ctapipe.utils.datasets import get_dataset_path
 from ctapipe.visualization import ArrayDisplay, CameraDisplay
 
 # %matplotlib inline
@@ -21,71 +19,53 @@
 print(ctapipe.__version__)
 print(ctapipe.__file__)
 
-
 ######################################################################
 # Let’s first open a raw event file and get an event out of it:
 #
 
-filename = get_dataset_path("gamma_prod5.simtel.zst")
-source = EventSource(filename, max_events=2)
+# this is using datasets from ctapipe's test data server
+source = EventSource("dataset://gamma_prod5.simtel.zst", max_events=2)
 
 for event in source:
     print(event.index.event_id)
 
-######################################################################
-filename
-
-######################################################################
+# %%
 source
 
-######################################################################
+# %%
 event
 
-######################################################################
-print(event.r1)
-
-
-######################################################################
-# Perform basic calibration:
-# --------------------------
-#
-# Here we will use a ``CameraCalibrator`` which is just a simple wrapper
-# that runs the three calibraraton and trace-integration phases of the
-# pipeline, taking the data from levels:
-#
-# **R0** → **R1** → **DL0** → **DL1**
-#
-# You could of course do these each separately, by using the classes
-# ``R1Calibrator``, ``DL0Reducer``, and ``DL1Calibrator``. Note that we
-# have not specified any configuration to the ``CameraCalibrator``, so it
-# will be using the default algorithms and thresholds, other than
-# specifying that the product is a “HESSIOR1Calibrator” (hopefully in the
-# near future that will be automatic).
-#
+# %%
+# to see for which telescopes we have data
+print(event.tel.keys())
+# get the first one:
+tel_id = next(iter(event.tel))
 
+# %%
+event.tel[tel_id]
 
+# %%
 calib = CameraCalibrator(subarray=source.subarray)
 calib(event)
 
-
 ######################################################################
-# Now the *r1*, *dl0* and *dl1* containers are filled in the event
+# Now the *r1*, *dl0* and *dl1* containers are filled in the telescope events
 #
-# -  **r1.tel[x]**: contains the “r1-calibrated” waveforms, after
+# -  **r1**: contains the “r1-calibrated” waveforms, after
 #    gain-selection, pedestal subtraciton, and gain-correction
-# -  **dl0.tel[x]**: is the same but with optional data volume reduction
+# -  **dl0**: is the same but with optional data volume reduction
 #    (some pixels not filled), in this case this is not performed by
 #    default, so it is the same as r1
-# -  **dl1.tel[x]**: contains the (possibly re-calibrated) waveforms as
+# -  **dl1**: contains the (possibly re-calibrated) waveforms as
 #    dl0, but also the time-integrated *image* that has been calculated
 #    using a ``ImageExtractor`` (a ``NeighborPeakWindowSum`` by default)
 #
 
-for tel_id in event.dl1.tel:
+for tel_id, tel_event in event.tel.items():
     print("TEL{:03}: {}".format(tel_id, source.subarray.tel[tel_id]))
-    print("  - r0  wave shape  : {}".format(event.r0.tel[tel_id].waveform.shape))
-    print("  - r1  wave shape  : {}".format(event.r1.tel[tel_id].waveform.shape))
-    print("  - dl1 image shape : {}".format(event.dl1.tel[tel_id].image.shape))
+    print("  - r0  wave shape  : {}".format(tel_event.r0.waveform.shape))
+    print("  - r1  wave shape  : {}".format(tel_event.r1.waveform.shape))
+    print("  - dl1 image shape : {}".format(tel_event.dl1.image.shape))
 
 
 ######################################################################
@@ -95,56 +75,46 @@
 # Let’s look at the image
 #
 
-
-tel_id = sorted(event.r1.tel.keys())[1]
 sub = source.subarray
 geometry = sub.tel[tel_id].camera.geometry
-image = event.dl1.tel[tel_id].image
+image = event.tel[tel_id].dl1.image
 
-######################################################################
+# %%
 disp = CameraDisplay(geometry, image=image)
 
-######################################################################
+# %%
 mask = tailcuts_clean(
     geometry,
     image,
     picture_thresh=10,
     boundary_thresh=5,
     min_number_picture_neighbors=2,
 )
-cleaned = image.copy()
-cleaned[~mask] = 0
-disp = CameraDisplay(geometry, image=cleaned)
+disp = CameraDisplay(geometry, image=image)
+disp.highlight_pixels(mask)
 
-######################################################################
-params = hillas_parameters(geometry, cleaned)
+# %%
+params = hillas_parameters(geometry[mask], image[mask])
 print(params)
 params
 
-######################################################################
-params = hillas_parameters(geometry, cleaned)
-
+# %%
 plt.figure(figsize=(10, 10))
 disp = CameraDisplay(geometry, image=image)
 disp.add_colorbar()
 disp.overlay_moments(params, color="red", lw=3)
-disp.highlight_pixels(mask, color="white", alpha=0.3, linewidth=2)
-
-plt.xlim(params.x.to_value(u.m) - 0.5, params.x.to_value(u.m) + 0.5)
-plt.ylim(params.y.to_value(u.m) - 0.5, params.y.to_value(u.m) + 0.5)
-
-######################################################################
-source.metadata
+disp.highlight_pixels(mask, color="xkcd:light blue", linewidth=5)
 
+plt.xlim(params.x.to_value(u.m) - 0.05, params.x.to_value(u.m) + 0.05)
+plt.ylim(params.y.to_value(u.m) - 0.05, params.y.to_value(u.m) + 0.05)
 
 ######################################################################
 # More complex image processing:
 # ------------------------------
 #
 # Let’s now explore how stereo reconstruction works.
 #
-# first, look at a summed image from multiple telescopes
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# First, look at a summed image from multiple telescopes
 #
 # For this, we want to use a ``CameraDisplay`` again, but since we can’t
 # sum and display images with different cameras, we’ll just sub-select
@@ -153,53 +123,42 @@
 # These are the telescopes that are in this event:
 #
 
-tels_in_event = set(
-    event.dl1.tel.keys()
-)  # use a set here, so we can intersect it later
+tels_in_event = set(event.tel.keys())
 tels_in_event
 
-######################################################################
+# %%
 cam_ids = set(sub.get_tel_ids_for_type("MST_MST_NectarCam"))
 cam_ids
 
-######################################################################
+# %%
 cams_in_event = tels_in_event.intersection(cam_ids)
 first_tel_id = list(cams_in_event)[0]
 tel = sub.tel[first_tel_id]
 print("{}s in event: {}".format(tel, cams_in_event))
 
-
 ######################################################################
 # Now, let’s sum those images:
 #
 
-image_sum = np.zeros_like(
-    tel.camera.geometry.pix_x.value
-)  # just make an array of 0's in the same shape as the camera
+image_sum = np.zeros(tel.camera.geometry.n_pixels)
 
 for tel_id in cams_in_event:
-    image_sum += event.dl1.tel[tel_id].image
-
+    image_sum += event.tel[tel_id].dl1.image
 
 ######################################################################
-# And finally display the sum of those images
-#
+# And finally display the sum of those images:
 
 plt.figure(figsize=(8, 8))
 
 disp = CameraDisplay(tel.camera.geometry, image=image_sum)
 disp.overlay_moments(params, with_label=False)
 plt.title("Sum of {}x {}".format(len(cams_in_event), tel))
 
-
 ######################################################################
 # let’s also show which telescopes those were. Note that currently
-# ArrayDisplay’s value field is a vector by ``tel_index``, not ``tel_id``,
-# so we have to convert to a tel_index. (this may change in a future
-# version to be more user-friendly)
+# ArrayDisplay’s value field is a vector by ``tel_index``.
 #
 
-
 nectarcam_subarray = sub.select_subarray(cam_ids, name="NectarCam")
 
 hit_pattern = np.zeros(shape=nectarcam_subarray.n_tels)