From 5474dcabdc0e7fe78049de0c30edc5836c0b5e57 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 30 Apr 2024 12:34:49 +0200 Subject: [PATCH] Implement CEP-2: new event container structure --- docs/api-reference/containers/index.rst | 2 +- docs/api-reference/io/index.rst | 2 +- docs/changes/2204.api.rst | 23 + docs/user-guide/data_models/dl1.rst | 18 +- .../tutorials/calibrated_data_exploration.py | 99 ++-- examples/tutorials/coordinates_example.py | 69 +-- examples/tutorials/ctapipe_handson.py | 171 +++---- examples/tutorials/ctapipe_overview.py | 482 ++++++++---------- examples/tutorials/raw_data_exploration.py | 88 ++-- .../tutorials/tableloader_and_eventsources.py | 44 +- examples/tutorials/theta_square.py | 5 +- examples/visualization/array_display.py | 43 +- examples/visualization/camera_display.py | 4 +- src/ctapipe/calib/camera/calibrator.py | 66 ++- src/ctapipe/calib/camera/flatfield.py | 33 +- src/ctapipe/calib/camera/pedestals.py | 30 +- .../calib/camera/tests/test_calibrator.py | 136 +++-- .../calib/camera/tests/test_flatfield.py | 42 +- .../calib/camera/tests/test_pedestals.py | 20 +- src/ctapipe/conftest.py | 2 +- src/ctapipe/containers.py | 359 ++++++------- src/ctapipe/image/cleaning.py | 37 +- src/ctapipe/image/extractor.py | 43 +- src/ctapipe/image/image_processor.py | 94 ++-- src/ctapipe/image/muon/processor.py | 46 +- .../image/muon/tests/test_processor.py | 6 +- src/ctapipe/image/reducer.py | 10 +- src/ctapipe/image/tests/test_extractor.py | 17 +- .../image/tests/test_image_processor.py | 50 +- src/ctapipe/instrument/tests/test_trigger.py | 107 ++-- src/ctapipe/instrument/trigger.py | 25 +- src/ctapipe/io/datawriter.py | 356 +++++++------ src/ctapipe/io/eventsource.py | 8 +- src/ctapipe/io/hdf5eventsource.py | 209 ++++---- src/ctapipe/io/hdf5merger.py | 30 +- src/ctapipe/io/simteleventsource.py | 131 +++-- src/ctapipe/io/tableloader.py | 33 +- src/ctapipe/io/tests/test_datawriter.py | 27 +- src/ctapipe/io/tests/test_event_source.py | 4 +- src/ctapipe/io/tests/test_hdf5.py | 14 +- src/ctapipe/io/tests/test_hdf5eventsource.py | 98 ++-- .../io/tests/test_simteleventsource.py | 77 +-- src/ctapipe/io/tests/test_table_loader.py | 5 +- src/ctapipe/io/tests/test_toysource.py | 3 +- src/ctapipe/io/toymodel.py | 15 +- src/ctapipe/reco/hillas_intersection.py | 10 +- src/ctapipe/reco/hillas_reconstructor.py | 17 +- src/ctapipe/reco/impact.py | 31 +- src/ctapipe/reco/preprocessing.py | 15 +- src/ctapipe/reco/reconstructor.py | 29 +- src/ctapipe/reco/shower_processor.py | 6 +- src/ctapipe/reco/sklearn.py | 41 +- src/ctapipe/reco/stereo_combination.py | 37 +- .../reco/tests/test_HillasReconstructor.py | 38 +- src/ctapipe/reco/tests/test_ImPACT.py | 8 +- .../reco/tests/test_hillas_intersection.py | 13 +- src/ctapipe/reco/tests/test_preprocessing.py | 10 +- .../reco/tests/test_reconstruction_methods.py | 8 +- .../reco/tests/test_shower_processor.py | 7 +- .../reco/tests/test_stereo_combination.py | 28 +- src/ctapipe/tools/apply_models.py | 9 +- src/ctapipe/tools/display_dl1.py | 17 +- src/ctapipe/tools/tests/test_apply_models.py | 12 +- src/ctapipe/tools/tests/test_process.py | 6 +- src/ctapipe/tools/tests/test_process_ml.py | 4 +- src/ctapipe/utils/event_type_filter.py | 2 +- src/ctapipe/utils/tests/test_event_filter.py | 20 +- src/ctapipe/visualization/mpl_array.py | 2 +- src/ctapipe/visualization/tests/test_bokeh.py | 20 +- src/ctapipe/visualization/tests/test_mpl.py | 33 +- test_plugin/ctapipe_test_plugin/__init__.py | 6 +- 71 files changed, 1787 insertions(+), 1825 deletions(-) create mode 100644 docs/changes/2204.api.rst diff --git a/docs/api-reference/containers/index.rst b/docs/api-reference/containers/index.rst index 2c2c4a2e3af..1b575da4f30 100644 --- a/docs/api-reference/containers/index.rst +++ b/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 diff --git a/docs/api-reference/io/index.rst b/docs/api-reference/io/index.rst index 9d4bd7f3dbe..9edaed81023 100644 --- a/docs/api-reference/io/index.rst +++ b/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 diff --git a/docs/changes/2204.api.rst b/docs/changes/2204.api.rst new file mode 100644 index 00000000000..28b0073896a --- /dev/null +++ b/docs/changes/2204.api.rst @@ -0,0 +1,23 @@ +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. + +Before, the data level came first and then below that +were mappings for each telescope: + +.. code:: python + + event.dl1.tel[tel_id].parameters.hillas + +Whereas now, the subarray event has a collection of telescope events, +which then contain the different data levels: + +.. code:: python + + event.tel[tel_id].dl1.parameters.hillas + +Fore more details, see :ref:`cep-002`. + +With this, also the container structure itself has been largely refactored. diff --git a/docs/user-guide/data_models/dl1.rst b/docs/user-guide/data_models/dl1.rst index 0e555449545..17af1902e5e 100644 --- a/docs/user-guide/data_models/dl1.rst +++ b/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,13 +64,13 @@ 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 @@ -78,13 +78,13 @@ 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/examples/tutorials/calibrated_data_exploration.py b/examples/tutorials/calibrated_data_exploration.py index 40d9d4bb202..c050c93528e 100644 --- a/examples/tutorials/calibrated_data_exploration.py +++ b/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,20 +19,16 @@ 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 @@ -42,50 +36,36 @@ event ###################################################################### -print(event.r1) - +# to see for which telescopes we have data +print(event.tel.keys()) +# get the first one: +tel_id = next(iter(event.tel)) ###################################################################### -# 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). -# - +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,11 +75,9 @@ # 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) @@ -112,30 +90,23 @@ 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: @@ -143,8 +114,7 @@ # # 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,9 +123,7 @@ # 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 ###################################################################### @@ -168,22 +136,17 @@ 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)) @@ -191,15 +154,11 @@ 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) diff --git a/examples/tutorials/coordinates_example.py b/examples/tutorials/coordinates_example.py index e2f1ca57a79..2319ffb57ba 100644 --- a/examples/tutorials/coordinates_example.py +++ b/examples/tutorials/coordinates_example.py @@ -4,17 +4,13 @@ """ -import copy - import astropy.units as u import matplotlib.pyplot as plt from astropy.coordinates import AltAz, EarthLocation, SkyCoord from astropy.time import Time from ctapipe.coordinates import CameraFrame, NominalFrame, TelescopeFrame -from ctapipe.instrument import SubarrayDescription from ctapipe.io import EventSource -from ctapipe.utils import get_dataset_path from ctapipe.visualization import CameraDisplay # %matplotlib inline @@ -24,36 +20,28 @@ plt.rcParams["figure.figsize"] = (12, 8) plt.rcParams["font.size"] = 16 - ###################################################################### # Open test dataset # ----------------- # -filename = get_dataset_path("gamma_prod5.simtel.zst") -source = EventSource(filename) - -events = [copy.deepcopy(event) for event in source] -event = events[4] - -layout = set(source.subarray.tel_ids) - +with EventSource("dataset://gamma_prod5.simtel.zst") as source: + events = [event for event in source] + event = events[4] + subarray = source.subarray + layout = set(subarray.tel_ids) ###################################################################### # Choose event with LST -# ~~~~~~~~~~~~~~~~~~~~~ -# - ###################################################################### # This ensures that the telescope is not “parked” (as it would be in an # event where it is not triggered) but is actually pointing to a source. # -print(f"Telescope with data: {event.r1.tel.keys()}") +print(f"Telescope with data: {event.tel.keys()}") tel_id = 3 - ###################################################################### # AltAz # ----- @@ -69,21 +57,19 @@ # oriented East of North (i.e., N=0°, E=90°). # - obstime = Time("2013-11-01T03:00") location = EarthLocation.of_site("Roque de los Muchachos") altaz = AltAz(location=location, obstime=obstime) array_pointing = SkyCoord( - alt=event.pointing.array_azimuth, - az=event.pointing.array_altitude, + alt=event.pointing.azimuth, + az=event.pointing.altitude, frame=altaz, ) print(array_pointing) - ###################################################################### # CameraFrame # ----------- @@ -93,28 +79,28 @@ # The camera frame is a 2d cartesian frame, describing position of objects # in the focal plane of the telescope. # -# The frame is defined as in H.E.S.S., starting at the horizon, the -# telescope is pointed to magnetic north in azimuth and then up to zenith. +# The frame is defined as in ``sim_telarray``, starting at the horizon, the +# telescope is pointed North in azimuth and then up to zenith. # -# Now, x points north and y points west, so in this orientation, the -# camera coordinates line up with the CORSIKA ground coordinate system. +# Now, x points North and y points West, so in this orientation, the +# camera coordinates line up with the ground coordinate system. # # MAGIC and FACT use a different camera coordinate system: Standing at the # dish, looking at the camera, x points right, y points up. To transform -# MAGIC/FACT to ctapipe, do x’ = -y, y’ = -x. +# MAGIC/FACT to ctapipe, do x’ = -y, y’ = -x. This frame is implemented in +# ctapipe as `~ctapipe.coordinates.EngineeringCameraFrame`. # # **Typical usage**: Position of pixels in the focal plane. -# -geometry = source.subarray.tel[tel_id].camera.geometry +geometry = subarray.tel[tel_id].camera.geometry pix_x = geometry.pix_x pix_y = geometry.pix_y -focal_length = source.subarray.tel[tel_id].optics.equivalent_focal_length +focal_length = subarray.tel[tel_id].optics.equivalent_focal_length ###################################################################### telescope_pointing = SkyCoord( - alt=event.pointing.tel[tel_id].altitude, - az=event.pointing.tel[tel_id].azimuth, + alt=event.tel[tel_id].pointing.altitude, + az=event.tel[tel_id].pointing.azimuth, frame=altaz, ) @@ -135,7 +121,6 @@ plt.ylabel(f"y / {cam_coords.y.unit}") plt.axis("square") - ###################################################################### # The implementation of the coordinate system with astropy makes it easier # to use time of the observation and location of the observing site, to @@ -143,7 +128,6 @@ # and how they might be visible in the camera. # - location = EarthLocation.of_site("Roque de los Muchachos") obstime = Time("2018-11-01T04:00") @@ -161,7 +145,6 @@ ) -subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst") cam = subarray.tel[1].camera.geometry fig, ax = plt.subplots() display = CameraDisplay(cam, ax=ax) @@ -188,7 +171,6 @@ plt.show() - ###################################################################### # TelescopeFrame # -------------- @@ -215,8 +197,6 @@ telescope_coords = cam_coords.transform_to(telescope_frame) ###################################################################### -wrap_angle = telescope_pointing.az + 180 * u.deg - plt.axis("equal") plt.scatter( telescope_coords.fov_lon.deg, telescope_coords.fov_lat.deg, alpha=0.2, color="gray" @@ -239,7 +219,6 @@ plt.xlabel("fov_lon / {}".format(telescope_coords.altaz.az.unit)) plt.ylabel("fov_lat / {}".format(telescope_coords.altaz.alt.unit)) - ###################################################################### # NominalFrame # ------------ @@ -254,7 +233,6 @@ # this frame # - ###################################################################### # Let’s play a bit with 3 LSTs with divergent pointing # @@ -344,16 +322,14 @@ # **Typical usage**: positions of telescopes on the ground (x, y, z) # -source.subarray.peek() - +subarray.peek() ###################################################################### # In case a layout is selected, the following line will produce a # different output from the picture above. # -source.subarray.select_subarray(layout, name="Prod3b layout").peek() - +subarray.select_subarray(layout, name="Prod5 layout").peek() ###################################################################### # .. figure:: ground_frame.png @@ -361,19 +337,16 @@ # # Ground Frame - ###################################################################### # In this image all the telescope from the ``gamma_test.simtel.gz`` file # are plotted as spheres in the GroundFrame. # - ###################################################################### # TiltedGroundFrame # ----------------- # - ###################################################################### # Tilted ground coordinate frame. # @@ -386,14 +359,12 @@ # This frame is used for the reconstruction of the shower core position. # - ###################################################################### # .. figure:: tilted_ground_frame.png # :alt: Tilted Ground Frame # # Tilted Ground Frame - ###################################################################### # This image picture both the telescopes in the GroundFrame (red) and in # the TiltedGroundFrame (green) are displayed: in this case since the diff --git a/examples/tutorials/ctapipe_handson.py b/examples/tutorials/ctapipe_handson.py index 48219067be9..faddc5cad09 100644 --- a/examples/tutorials/ctapipe_handson.py +++ b/examples/tutorials/ctapipe_handson.py @@ -2,28 +2,25 @@ Getting Started with ctapipe ============================ -This hands-on was presented at the Paris CTA Consoritum meeting (K. -Kosack) - +This hands-on was initially presented at the Paris CTA Consoritum meeting by K. +Kosack and has been updated since to reflect changes in ctapipe itself. """ - ###################################################################### -# Part 1: load and loop over data -# ------------------------------- -# +# Part 1: load and loop over events +# --------------------------------- import glob +import astropy.units as u import numpy as np -import pandas as pd from ipywidgets import interact from matplotlib import pyplot as plt from ctapipe import utils from ctapipe.calib import CameraCalibrator from ctapipe.image import hillas_parameters, tailcuts_clean -from ctapipe.io import EventSource, HDF5TableWriter +from ctapipe.io import EventSource, HDF5TableWriter, read_table from ctapipe.visualization import CameraDisplay # %matplotlib inline @@ -32,26 +29,32 @@ path = utils.get_dataset_path("gamma_prod5.simtel.zst") ###################################################################### -source = EventSource(path, max_events=5) - -for event in source: - print(event.count, event.index.event_id, event.simulation.shower.energy) +with EventSource(path, max_events=5) as source: + for event in source: + print( + f"{event.count} {event.index.event_id} {event.simulation.shower.energy:.3f} {len(event.tel)}" + ) ###################################################################### event ###################################################################### -event.r1 +event.tel.keys() ###################################################################### -for event in EventSource(path, max_events=5): - print(event.count, event.r1.tel.keys()) +tel_id = 118 +event.tel[tel_id] ###################################################################### -event.r0.tel[3] +with EventSource(path, max_events=5) as source: + for event in source: + print(event.count, event.tel.keys()) + +###################################################################### +event.tel[tel_id].r0 ###################################################################### -r0tel = event.r0.tel[3] +r0tel = event.tel[tel_id].r0 ###################################################################### r0tel.waveform @@ -59,7 +62,6 @@ ###################################################################### r0tel.waveform.shape - ###################################################################### # note that this is (:math:`N_{channels}`, :math:`N_{pixels}`, # :math:`N_{samples}`) @@ -68,17 +70,20 @@ plt.pcolormesh(r0tel.waveform[0]) ###################################################################### -brightest_pixel = np.argmax(r0tel.waveform[0].sum(axis=1)) -print(f"pixel {brightest_pixel} has sum {r0tel.waveform[0,1535].sum()}") +waveform_sums = r0tel.waveform[0].sum(axis=1) +brightest_pixel = np.argmax(waveform_sums) +print(f"pixel {brightest_pixel} has sum {waveform_sums[brightest_pixel]}") ###################################################################### plt.plot(r0tel.waveform[0, brightest_pixel], label="channel 0 (high-gain)") plt.plot(r0tel.waveform[1, brightest_pixel], label="channel 1 (low-gain)") plt.legend() - ###################################################################### -@interact +n_channels, n_pixels, n_samples = r0tel.waveform.shape + + +@interact(chan=(0, n_channels - 1), pix_id=(0, n_pixels - 1)) def view_waveform(chan=0, pix_id=brightest_pixel): plt.plot(r0tel.waveform[chan, pix_id]) @@ -87,7 +92,6 @@ def view_waveform(chan=0, pix_id=brightest_pixel): # try making this compare 2 waveforms # - ###################################################################### # Part 2: Explore the instrument description # ------------------------------------------ @@ -95,9 +99,7 @@ def view_waveform(chan=0, pix_id=brightest_pixel): # This is all well and good, but we don’t really know what camera or # telescope this is… how do we get instrumental description info? # -# Currently this is returned *inside* the event (it will soon change to be -# separate in next version or so) -# +# This information is provided by the ``SubarrayDescription`` of the event source: subarray = source.subarray @@ -111,16 +113,16 @@ def view_waveform(chan=0, pix_id=brightest_pixel): subarray.to_table() ###################################################################### -subarray.tel[2] +subarray.tel[tel_id] ###################################################################### -subarray.tel[2].camera +subarray.tel[tel_id].camera ###################################################################### -subarray.tel[2].optics +subarray.tel[tel_id].optics ###################################################################### -tel = subarray.tel[2] +tel = subarray.tel[tel_id] ###################################################################### tel.camera @@ -143,138 +145,113 @@ def view_waveform(chan=0, pix_id=brightest_pixel): ###################################################################### disp = CameraDisplay(tel.camera.geometry) -disp.image = r0tel.waveform[ - 0, :, 10 -] # display channel 0, sample 0 (try others like 10) +# display channel 0, sample 20 (try others like 15) +disp.image = r0tel.waveform[0, :, 20] ###################################################################### # \*\* aside: \*\* show demo using a CameraDisplay in interactive mode in # ipython rather than notebook # - ###################################################################### # Part 3: Apply some calibration and trace integration # ---------------------------------------------------- # - calib = CameraCalibrator(subarray=subarray) ###################################################################### -for event in EventSource(path, max_events=5): - calib(event) # fills in r1, dl0, and dl1 - print(event.dl1.tel.keys()) - -###################################################################### -event.dl1.tel[3] +with EventSource(path, max_events=5) as source: + for event in source: + calib(event) # fills in dl0, and dl1 ###################################################################### -dl1tel = event.dl1.tel[3] +dl1tel = event.tel[tel_id].dl1 ###################################################################### -dl1tel.image.shape # note this will be gain-selected in next version, so will be just 1D array of 1855 +dl1tel.image.shape ###################################################################### dl1tel.peak_time ###################################################################### -CameraDisplay(tel.camera.geometry, image=dl1tel.image) +d = CameraDisplay(tel.camera.geometry, image=dl1tel.image) +d.add_colorbar() ###################################################################### CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time) - ###################################################################### # Now for Hillas Parameters # - image = dl1tel.image -mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5) +mask = tailcuts_clean( + tel.camera.geometry, + image, + picture_thresh=10, + boundary_thresh=5, + keep_isolated_pixels=False, + min_number_picture_neighbors=2, +) mask ###################################################################### -CameraDisplay(tel.camera.geometry, image=mask) - -###################################################################### -cleaned = image.copy() -cleaned[~mask] = 0 +d = CameraDisplay(tel.camera.geometry, image=image) +d.highlight_pixels(mask) ###################################################################### -disp = CameraDisplay(tel.camera.geometry, image=cleaned) -disp.cmap = plt.cm.coolwarm -disp.add_colorbar() -plt.xlim(0.5, 1.0) -plt.ylim(-1.0, 0.0) - -###################################################################### -params = hillas_parameters(tel.camera.geometry, cleaned) +params = hillas_parameters(tel.camera.geometry[mask], image[mask]) print(params) ###################################################################### -disp = CameraDisplay(tel.camera.geometry, image=cleaned) -disp.cmap = plt.cm.coolwarm +disp = CameraDisplay(tel.camera.geometry, image=image) +disp.highlight_pixels(mask, linewidth=2) disp.add_colorbar() -plt.xlim(0.5, 1.0) -plt.ylim(-1.0, 0.0) -disp.overlay_moments(params, color="white", lw=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) +disp.overlay_moments(params, color="white", lw=2, n_sigma=2) ###################################################################### # Part 4: Let’s put it all together: # ---------------------------------- # # - loop over events, selecting only telescopes of the same type -# (e.g. LST:LSTCam) +# (e.g. LST_LST_LSTCam) # - for each event, apply calibration/trace integration # - calculate Hillas parameters # - write out all hillas parameters to a file that can be loaded with # Pandas # - ###################################################################### -# first let’s select only those telescopes with LST:LSTCam -# +# first let’s select only those telescopes with LST_LST_LSTCam subarray.telescope_types ###################################################################### subarray.get_tel_ids_for_type("LST_LST_LSTCam") - ###################################################################### # Now let’s write out program # data = utils.get_dataset_path("gamma_prod5.simtel.zst") -source = EventSource(data) # remove the max_events limit to get more stats - -###################################################################### -for event in source: - calib(event) - - for tel_id, tel_data in event.dl1.tel.items(): - tel = source.subarray.tel[tel_id] - mask = tailcuts_clean(tel.camera.geometry, tel_data.image) - if np.count_nonzero(mask) > 0: - params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask]) - ###################################################################### with HDF5TableWriter(filename="hillas.h5", group_name="dl1", overwrite=True) as writer: - source = EventSource(data, allowed_tels=[1, 2, 3, 4], max_events=10) - for event in source: - calib(event) - - for tel_id, tel_data in event.dl1.tel.items(): - tel = source.subarray.tel[tel_id] - mask = tailcuts_clean(tel.camera.geometry, tel_data.image) - params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask]) - writer.write("hillas", params) + with EventSource(data) as source: + for event in source: + calib(event) + for tel_id, tel_event in event.tel.items(): + tel = source.subarray.tel[tel_id] + image = tel_event.dl1.image + mask = tailcuts_clean(tel.camera.geometry, image) + if np.count_nonzero(mask) > 0: + params = hillas_parameters(tel.camera.geometry[mask], image[mask]) + writer.write("hillas", params) ###################################################################### # We can now load in the file we created and plot it @@ -284,12 +261,16 @@ def view_waveform(chan=0, pix_id=brightest_pixel): ###################################################################### -hillas = pd.read_hdf("hillas.h5", key="/dl1/hillas") +hillas = read_table("hillas.h5", "/dl1/hillas") hillas ###################################################################### -_ = hillas.hist(figsize=(8, 8)) +fig, axs = plt.subplots(4, 3, figsize=(8, 8), layout="constrained") +axs = axs.ravel() +for col, ax in zip(hillas.colnames, axs): + ax.hist(hillas[col]) + ax.set_title(col) ###################################################################### # If you do this yourself, chose a larger file to loop over more events to diff --git a/examples/tutorials/ctapipe_overview.py b/examples/tutorials/ctapipe_overview.py index 38ed0ac3873..323bb0a25a5 100644 --- a/examples/tutorials/ctapipe_overview.py +++ b/examples/tutorials/ctapipe_overview.py @@ -2,18 +2,13 @@ Analyzing Events Using ctapipe ============================== -""" - - -###################################################################### -# Initially presented @ LST Analysis Bootcamp -# in Padova, 26.11.2018 -# by Maximilian Nöthe (@maxnoe) & Kai A. Brügge (@mackaiver) +Initially presented @ LST Analysis Bootcamp in Padova, 26.11.2018 +by Maximilian Linhoff (@maxnoe) & Kai A. Brügge (@mackaiver). - -import tempfile +Updated since to stay compatible with current ctapipe. +""" import timeit -from copy import deepcopy +import tempfile import astropy.units as u import matplotlib.pyplot as plt @@ -34,218 +29,197 @@ timing_parameters, toymodel, ) -from ctapipe.image.cleaning import tailcuts_clean +from ctapipe.image.cleaning import TailcutsImageCleaner from ctapipe.io import DataWriter, EventSource, TableLoader from ctapipe.reco import ShowerProcessor from ctapipe.utils.datasets import get_dataset_path from ctapipe.visualization import ArrayDisplay, CameraDisplay # %matplotlib inline - -###################################################################### +############################################################################### plt.rcParams["figure.figsize"] = (12, 8) plt.rcParams["font.size"] = 14 plt.rcParams["figure.figsize"] - -###################################################################### +############################################################################### # General Information # ------------------- # - -###################################################################### +############################################################################### # Design -# ~~~~~~ +# ^^^^^^ # -# - DL0 → DL3 analysis +# - DL0 → DL3 analysis # -# - Currently some R0 → DL2 code to be able to analyze simtel files +# - Currently some R0 → DL2 code to be able to analyze simtel files # -# - ctapipe is built upon the Scientific Python Stack, core dependencies -# are +# - ctapipe is built upon the Scientific Python Stack, core dependencies +# are # -# - numpy -# - scipy -# - astropy -# - numba +# - numpy +# - scipy +# - astropy +# - numba # - -###################################################################### +############################################################################### # Development -# ~~~~~~~~~~~~ +# ^^^^^^^^^^^ # -# - ctapipe is developed as Open Source Software (BSD 3-Clause License) -# at https://github.com/cta-observatory/ctapipe +# - ctapipe is developed as Open Source Software (BSD 3-Clause License) +# at https://github.com/cta-observatory/ctapipe # -# - We use the “Github-Workflow”: +# - We use the “Github-Workflow”: # -# - Few people (e.g. @kosack, @maxnoe) have write access to the main -# repository -# - Contributors fork the main repository and work on branches -# - Pull Requests are merged after Code Review and automatic execution -# of the test suite +# - Few people (e.g. @kosack, @maxnoe) have write access to the main +# repository +# - Contributors fork the main repository and work on branches +# - Pull Requests are merged after Code Review and automatic execution +# of the test suite # -# - Early development stage ⇒ backwards-incompatible API changes might -# and will happen +# - Early development stage ⇒ backwards-incompatible API changes might +# and will happen # - -###################################################################### +############################################################################### # What’s there? -# ~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^ +# +# - Reading simtel simulation files +# - Simple calibration, cleaning and feature extraction functions +# - Camera and Array plotting +# - Coordinate frames and transformations +# - Stereo-reconstruction using line intersections # -# - Reading simtel simulation files -# - Simple calibration, cleaning and feature extraction functions -# - Camera and Array plotting -# - Coordinate frames and transformations -# - Stereo-reconstruction using line intersections # - -###################################################################### +############################################################################### # What’s still missing? -# ~~~~~~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^^^^^^ # -# - Good integration with machine learning techniques -# - IRF calculation -# - Documentation, e.g. formal definitions of coordinate frames +# - IRF calculation +# - Documentation, e.g. formal definitions of coordinate frames # - -###################################################################### +############################################################################### # What can you do? -# ~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^ # -# - Report issues +# - Report issues # -# - Hard to get started? Tell us where you are stuck -# - Tell user stories -# - Missing features +# - Hard to get started? Tell us where you are stuck +# - Tell user stories +# - Missing features # # - Start contributing # -# - ctapipe needs more workpower -# - Implement new reconstruction features +# - ctapipe needs more workpower +# - Implement new reconstruction features # - -###################################################################### -# A simple hillas analysis +############################################################################### +# A simple Hillas analysis # ------------------------ # - -###################################################################### # Reading in simtel files -# ~~~~~~~~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^^^^^^^^ # - - input_url = get_dataset_path("gamma_prod5.simtel.zst") - +############################################################################### # EventSource() automatically detects what kind of file we are giving it, -# if already supported by ctapipe -source = EventSource(input_url, max_events=5) - -print(type(source)) +# if already supported by ctapipe or an installed plugin +with EventSource(input_url, max_events=5) as source: + print(type(source)) -###################################################################### -for event in source: - print( - "Id: {}, E = {:1.3f}, Telescopes: {}".format( - event.count, event.simulation.shower.energy, len(event.r0.tel) + for event in source: + print( + "Id: {}, E = {:1.3f}, Telescopes: {}".format( + event.count, event.simulation.shower.energy, len(event.tel) + ) ) - ) - -###################################################################### -# Each event is a ``DataContainer`` holding several ``Field``\ s of data, +############################################################################### +# Each event is a ``SubarrayEventContainer`` holding several ``Field`` s of data, # which can be containers or just numbers. Let’s look a one event: -# - event -###################################################################### +############################################################################### source.subarray.camera_types -###################################################################### -len(event.r0.tel), len(event.r1.tel) +############################################################################### +# Telescope-wise data is stored in ``TelescopeEventContainer``s +# under ``.tel``, which is a mapping by telescope id: +len(event.tel), event.tel.keys() +############################################################################### +event.tel[3] -###################################################################### +############################################################################### # Data calibration -# ~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^ # # The ``CameraCalibrator`` calibrates the event (obtaining the ``dl1`` # images). -# - - calibrator = CameraCalibrator(subarray=source.subarray) -###################################################################### +############################################################################### calibrator(event) - -###################################################################### +############################################################################### # Event displays -# ~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^ # # Let’s use ctapipe’s plotting facilities to plot the telescope images -# - -event.dl1.tel.keys() - -###################################################################### tel_id = 130 -###################################################################### +############################################################################### geometry = source.subarray.tel[tel_id].camera.geometry -dl1 = event.dl1.tel[tel_id] +dl1 = event.tel[tel_id].dl1 -geometry, dl1 +print(geometry) +print(dl1) -###################################################################### +############################################################################### dl1.image - -###################################################################### +############################################################################### display = CameraDisplay(geometry) - -# right now, there might be one image per gain channel. -# This will change as soon as display.image = dl1.image display.add_colorbar() - -###################################################################### +############################################################################### # Image Cleaning -# ~~~~~~~~~~~~~~ -# - - -# unoptimized cleaning levels -cleaning_level = { - "CHEC": (2, 4, 2), - "LSTCam": (3.5, 7, 2), - "FlashCam": (3.5, 7, 2), - "NectarCam": (4, 8, 2), -} - -###################################################################### -boundary, picture, min_neighbors = cleaning_level[geometry.name] - -clean = tailcuts_clean( - geometry, +# ^^^^^^^^^^^^^^ +# +# ctapipe allows most configuration options to be configured by telescope type. +# The example below configures some of the parameters for the +# `~ctapipe.image.TailcutsImageCleaner`: + +############################################################################### +print(source.subarray.telescope_types) +############################################################################### +cleaning = TailcutsImageCleaner( + source.subarray, + picture_threshold_pe=[ + ("type", "*", 7), # global default + ("type", "MST_MST_NectarCam", 8), + ("type", "SST_ASTRI_CHEC", 4), + ], + boundary_threshold_pe=[ + ("type", "*", 3.5), + ("type", "MST_MST_NectarCam", 4), + ("type", "SST_ASTRI_CHEC", 2), + ], +) +############################################################################### +image_mask = cleaning( + tel_id, dl1.image, - boundary_thresh=boundary, - picture_thresh=picture, - min_number_picture_neighbors=min_neighbors, ) -###################################################################### +############################################################################### fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5)) d1 = CameraDisplay(geometry, ax=ax1) @@ -253,6 +227,7 @@ ax1.set_title("Image") d1.image = dl1.image +d1.highlight_pixels(image_mask, color="w") d1.add_colorbar(ax=ax1) ax2.set_title("Pulse Time") @@ -261,66 +236,57 @@ d2.add_colorbar(ax=ax2) d2.set_limits_minmax(-20, 20) -d1.highlight_pixels(clean, color="red", linewidth=1) - -###################################################################### +############################################################################### # Image Parameters -# ~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^ # - - -hillas = hillas_parameters(geometry[clean], dl1.image[clean]) +hillas = hillas_parameters(geometry[image_mask], dl1.image[image_mask]) print(hillas) -###################################################################### -display = CameraDisplay(geometry) - -# set "unclean" pixels to 0 -cleaned = dl1.image.copy() -cleaned[~clean] = 0.0 - -display.image = cleaned +############################################################################### +plt.figure() +display = CameraDisplay(geometry, image=dl1.image) +display.highlight_pixels(image_mask, color="w") display.add_colorbar() +display.overlay_moments(hillas, color="xkcd:red", n_sigma=2) -display.overlay_moments(hillas, color="xkcd:red") - -###################################################################### -timing = timing_parameters(geometry, dl1.image, dl1.peak_time, hillas, clean) - +############################################################################### +timing = timing_parameters(geometry, dl1.image, dl1.peak_time, hillas, image_mask) print(timing) -###################################################################### +############################################################################### long, trans = camera_to_shower_coordinates( geometry.pix_x, geometry.pix_y, hillas.x, hillas.y, hillas.psi ) -plt.plot(long[clean], dl1.peak_time[clean], "o") -plt.plot(long[clean], timing.slope * long[clean] + timing.intercept) +plt.figure() +plt.plot(long[image_mask], dl1.peak_time[image_mask], "o") +plt.plot(long[image_mask], timing.slope * long[image_mask] + timing.intercept) -###################################################################### -leakage = leakage_parameters(geometry, dl1.image, clean) +############################################################################### +leakage = leakage_parameters(geometry, dl1.image, image_mask) print(leakage) -###################################################################### +############################################################################### disp = CameraDisplay(geometry) disp.image = dl1.image disp.highlight_pixels(geometry.get_border_pixel_mask(1), linewidth=2, color="xkcd:red") -###################################################################### -n_islands, island_id = number_of_islands(geometry, clean) +############################################################################### +n_islands, island_id = number_of_islands(geometry, image_mask) print(n_islands) -###################################################################### +############################################################################### conc = concentration_parameters(geometry, dl1.image, hillas) print(conc) -###################################################################### +############################################################################### # Putting it all together / Stereo reconstruction -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ # # All these steps are now unified in several components configurable # through the config system, mainly: @@ -358,50 +324,51 @@ ) input_url = get_dataset_path("gamma_prod5.simtel.zst") -source = EventSource(input_url) - -calibrator = CameraCalibrator(subarray=source.subarray) -image_processor = ImageProcessor( - subarray=source.subarray, config=image_processor_config -) -shower_processor = ShowerProcessor(subarray=source.subarray) -horizon_frame = AltAz() - f = tempfile.NamedTemporaryFile(suffix=".hdf5") -with DataWriter(source, output_path=f.name, overwrite=True, write_dl2=True) as writer: - for event in source: - energy = event.simulation.shower.energy - n_telescopes_r1 = len(event.r1.tel) - event_id = event.index.event_id - print(f"Id: {event_id}, E = {energy:1.3f}, Telescopes (R1): {n_telescopes_r1}") - - calibrator(event) - image_processor(event) - shower_processor(event) - - stereo = event.dl2.stereo.geometry["HillasReconstructor"] - if stereo.is_valid: - print(" Alt: {:.2f}°".format(stereo.alt.deg)) - print(" Az: {:.2f}°".format(stereo.az.deg)) - print(" Hmax: {:.0f}".format(stereo.h_max)) - print(" CoreX: {:.1f}".format(stereo.core_x)) - print(" CoreY: {:.1f}".format(stereo.core_y)) - print(" Multiplicity: {:d}".format(len(stereo.telescopes))) - - # save a nice event for plotting later - if event.count == 3: - plotting_event = deepcopy(event) - - writer(event) - - -###################################################################### +with EventSource(input_url) as source: + calibrator = CameraCalibrator(subarray=source.subarray) + image_processor = ImageProcessor( + subarray=source.subarray, config=image_processor_config + ) + shower_processor = ShowerProcessor(subarray=source.subarray) + horizon_frame = AltAz() + + with DataWriter( + source, output_path=f.name, overwrite=True, write_dl2=True + ) as writer: + for event in source: + energy = event.simulation.shower.energy + n_telescopes = len(event.tel) + event_id = event.index.event_id + print(f"Id: {event_id}, E = {energy:1.3f}, Telescopes: {n_telescopes}") + + calibrator(event) + image_processor(event) + shower_processor(event) + + stereo = event.dl2.geometry["HillasReconstructor"] + if stereo.is_valid: + print(" Alt: {:.2f}°".format(stereo.alt.deg)) + print(" Az: {:.2f}°".format(stereo.az.deg)) + print(" Hmax: {:.0f}".format(stereo.h_max)) + print(" CoreX: {:.1f}".format(stereo.core_x)) + print(" CoreY: {:.1f}".format(stereo.core_y)) + print(" Multiplicity: {:d}".format(len(stereo.telescopes))) + + # save a nice event for plotting later + if event.count == 3: + plotting_event = event + + writer(event) + + +############################################################################### loader = TableLoader(f.name) events = loader.read_subarray_events() -###################################################################### +############################################################################### theta = angular_separation( events["HillasReconstructor_az"].quantity, events["HillasReconstructor_alt"].quantity, @@ -414,21 +381,18 @@ None -###################################################################### +############################################################################### # ArrayDisplay # ------------ # +angle_offset = plotting_event.pointing.azimuth -angle_offset = plotting_event.pointing.array_azimuth - -plotting_hillas = { - tel_id: dl1.parameters.hillas for tel_id, dl1 in plotting_event.dl1.tel.items() -} - -plotting_core = { - tel_id: dl1.parameters.core.psi for tel_id, dl1 in plotting_event.dl1.tel.items() -} +plotting_hillas = {} +plotting_core = {} +for tel_id, tel_event in plotting_event.tel.items(): + plotting_hillas[tel_id] = tel_event.dl1.parameters.hillas + plotting_core[tel_id] = tel_event.dl1.parameters.core.psi disp = ArrayDisplay(source.subarray) @@ -444,8 +408,8 @@ label="True Impact", ) plt.scatter( - plotting_event.dl2.stereo.geometry["HillasReconstructor"].core_x, - plotting_event.dl2.stereo.geometry["HillasReconstructor"].core_y, + plotting_event.dl2.geometry["HillasReconstructor"].core_x, + plotting_event.dl2.geometry["HillasReconstructor"].core_y, s=200, c="r", marker="x", @@ -453,13 +417,10 @@ ) plt.legend() -# plt.xlim(-400, 400) -# plt.ylim(-400, 400) - -###################################################################### +############################################################################### # Reading the LST dl1 data -# ~~~~~~~~~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^^^^^^^^^ # loader = TableLoader(f.name) @@ -470,7 +431,7 @@ true_parameters=False, ) -###################################################################### +############################################################################### plt.scatter( np.log10(dl1_table["true_energy"].quantity / u.TeV), np.log10(dl1_table["hillas_intensity"]), @@ -480,7 +441,7 @@ None -###################################################################### +############################################################################### # Isn’t python slow? # ------------------ # @@ -497,41 +458,42 @@ # **But: “Premature Optimization is the root of all evil” — Donald Knuth** # # So profile to find exactly what is slow. -# + +############################################################################### # Why use python then? -# ~~~~~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^^^^^ # # - Python works very well as *glue* for libraries of all kinds of # languages # - Python has a rich ecosystem for data science, physics, algorithms, # astronomy -# + +############################################################################### # Example: Number of Islands -# ~~~~~~~~~~~~~~~~~~~~~~~~~~ +# ^^^^^^^^^^^^^^^^^^^^^^^^^^ # # Find all groups of pixels, that survived the cleaning -# geometry = loader.subarray.tel[1].camera.geometry -###################################################################### +############################################################################### # Let’s create a toy images with several islands; # -np.random.seed(42) +rng = np.random.default_rng(42) image = np.zeros(geometry.n_pixels) for i in range(9): model = toymodel.Gaussian( - x=np.random.uniform(-0.8, 0.8) * u.m, - y=np.random.uniform(-0.8, 0.8) * u.m, - width=np.random.uniform(0.05, 0.075) * u.m, - length=np.random.uniform(0.1, 0.15) * u.m, - psi=np.random.uniform(0, 2 * np.pi) * u.rad, + x=rng.uniform(-0.8, 0.8) * u.m, + y=rng.uniform(-0.8, 0.8) * u.m, + width=rng.uniform(0.05, 0.075) * u.m, + length=rng.uniform(0.1, 0.15) * u.m, + psi=rng.uniform(0, 2 * np.pi) * u.rad, ) new_image, sig, bg = model.generate_image( @@ -539,23 +501,20 @@ ) image += new_image -###################################################################### -clean = tailcuts_clean( - geometry, - image, - picture_thresh=10, - boundary_thresh=5, - min_number_picture_neighbors=2, +############################################################################### +image_mask = cleaning( + tel_id=1, + image=image, ) -###################################################################### +############################################################################### disp = CameraDisplay(geometry) disp.image = image -disp.highlight_pixels(clean, color="xkcd:red", linewidth=1.5) +disp.highlight_pixels(image_mask, color="xkcd:red", linewidth=1.5) disp.add_colorbar() -###################################################################### +############################################################################### def num_islands_python(camera, clean): """A breadth first search to find connected islands of neighboring pixels in the cleaning set""" @@ -594,11 +553,10 @@ def num_islands_python(camera, clean): return n_islands, island_ids -###################################################################### -n_islands, island_ids = num_islands_python(geometry, clean) - +############################################################################### +n_islands, island_ids = num_islands_python(geometry, image_mask) -###################################################################### +############################################################################### cmap = plt.get_cmap("Paired") cmap = ListedColormap(cmap.colors[:n_islands]) cmap.set_under("k") @@ -609,11 +567,12 @@ def num_islands_python(camera, clean): disp.set_limits_minmax(0.5, n_islands + 0.5) disp.add_colorbar() -###################################################################### -timeit.timeit(lambda: num_islands_python(geometry, clean), number=1000) / 1000 +############################################################################### +timeit.timeit(lambda: num_islands_python(geometry, image_mask), number=1000) / 1000 + +############################################################################### -###################################################################### def num_islands_scipy(geometry, clean): neighbors = geometry.neighbor_matrix_sparse @@ -626,27 +585,22 @@ def num_islands_scipy(geometry, clean): return num_islands, island_ids -###################################################################### -n_islands_s, island_ids_s = num_islands_scipy(geometry, clean) +############################################################################### +n_islands_s, island_ids_s = num_islands_scipy(geometry, image_mask) -###################################################################### +############################################################################### disp = CameraDisplay(geometry) disp.image = island_ids_s disp.cmap = cmap disp.set_limits_minmax(0.5, n_islands_s + 0.5) disp.add_colorbar() -###################################################################### -timeit.timeit(lambda: num_islands_scipy(geometry, clean), number=10000) / 10000 +############################################################################### +timeit.timeit(lambda: num_islands_scipy(geometry, image_mask), number=10000) / 10000 -###################################################################### +############################################################################### # **A lot less code, and a factor 3 speed improvement** # - - -###################################################################### # Finally, current ctapipe implementation is using numba: # - -###################################################################### -timeit.timeit(lambda: number_of_islands(geometry, clean), number=100000) / 100000 +timeit.timeit(lambda: number_of_islands(geometry, image_mask), number=100000) / 100000 diff --git a/examples/tutorials/raw_data_exploration.py b/examples/tutorials/raw_data_exploration.py index 4610996478f..3cdc219f742 100644 --- a/examples/tutorials/raw_data_exploration.py +++ b/examples/tutorials/raw_data_exploration.py @@ -10,11 +10,11 @@ """ - ###################################################################### # Setup: # +import astropy.units as u import numpy as np from matplotlib import pyplot as plt from scipy import signal @@ -25,7 +25,6 @@ # %matplotlib inline - ###################################################################### # To read SimTelArray format data, ctapipe uses the ``pyeventio`` library # (which is installed automatically along with ctapipe). The following @@ -42,7 +41,6 @@ source = EventSource(get_dataset_path("gamma_prod5.simtel.zst"), max_events=5) - ###################################################################### # Explore the contents of an event # -------------------------------- @@ -53,17 +51,14 @@ # so we can advance through events one-by-one event_iterator = iter(source) - event = next(event_iterator) - ###################################################################### # the event is just a class with a bunch of data items in it. You can see # a more compact representation via: # -event.r0 - +event ###################################################################### # printing the event structure, will currently print the value all items @@ -74,15 +69,14 @@ print(event.simulation.shower) ###################################################################### -print(event.r0.tel.keys()) - +print(event.tel.keys()) ###################################################################### # note that the event has 3 telescopes in it: Let’s try the next one: # event = next(event_iterator) -print(event.r0.tel.keys()) +print(event.tel.keys()) ###################################################################### @@ -90,10 +84,11 @@ # them: # -teldata = event.r0.tel[26] -print(teldata) -teldata +tel_event = event.tel[26] +tel_event +###################################################################### +print(tel_event.r0) ###################################################################### # Note that some values are unit quantities (``astropy.units.Quantity``) @@ -104,10 +99,10 @@ event.simulation.shower.energy ###################################################################### -event.simulation.shower.energy.to("GeV") +event.simulation.shower.energy.to(u.GeV) ###################################################################### -event.simulation.shower.energy.to("J") +event.simulation.shower.energy.to(u.J) ###################################################################### event.simulation.shower.alt @@ -115,7 +110,6 @@ ###################################################################### print("Altitude in degrees:", event.simulation.shower.alt.deg) - ###################################################################### # Look for signal pixels in a camera # ---------------------------------- @@ -127,24 +121,25 @@ # if we see which pixels contain Cherenkov light signals: # -plt.pcolormesh(teldata.waveform[0]) # note the [0] is for channel 0 +plt.pcolormesh(tel_event.r0.waveform[0]) # note the [0] is for channel 0 plt.colorbar() plt.xlabel("sample number") plt.ylabel("Pixel_id") - ###################################################################### # Let’s zoom in to see if we can identify the pixels that have the # Cherenkov signal in them # -plt.pcolormesh(teldata.waveform[0]) +plt.pcolormesh(tel_event.r0.waveform[0]) plt.colorbar() plt.ylim(700, 750) plt.xlabel("sample number") plt.ylabel("pixel_id") -print("waveform[0] is an array of shape (N_pix,N_slice) =", teldata.waveform[0].shape) - +print( + "waveform is an array of shape (n_channels, n_pixels, n_samples) =", + tel_event.r0.waveform.shape, +) ###################################################################### # Now we can really see that some pixels have a signal in them! @@ -152,10 +147,9 @@ # Lets look at a 1D plot of pixel 270 in channel 0 and see the signal: # -trace = teldata.waveform[0][719] +trace = tel_event.r0.waveform[0][719] plt.plot(trace, drawstyle="steps") - ###################################################################### # Great! It looks like a *standard Cherenkov signal*! # @@ -164,7 +158,9 @@ for pix_id in range(718, 723): plt.plot( - teldata.waveform[0][pix_id], label="pix {}".format(pix_id), drawstyle="steps" + tel_event.r0.waveform[0][pix_id], + label="pix {}".format(pix_id), + drawstyle="steps", ) plt.legend() @@ -183,7 +179,7 @@ # for pix_id in range(718, 723): - plt.plot(teldata.waveform[0][pix_id], "+-") + plt.plot(tel_event.r0.waveform[0][pix_id], "+-") plt.fill_betweenx([0, 1600], 19, 24, color="red", alpha=0.3, label="Ped window") plt.fill_betweenx([0, 1600], 5, 9, color="green", alpha=0.3, label="Signal window") plt.legend() @@ -199,34 +195,37 @@ # which is the high-gain channel): # -data = teldata.waveform[0] +data = tel_event.r0.waveform[0] peds = data[:, 19:24].mean(axis=1) # mean of samples 20 to 29 for all pixels -sums = data[:, 5:9].sum(axis=1) / (13 - 8) # simple sum integration +sums = data[:, 5:9].sum(axis=1) / (9 - 5) # simple sum integration ###################################################################### -phist = plt.hist(peds, bins=50, range=[0, 150]) +phist = plt.hist(peds, bins=50) plt.title("Pedestal Distribution of all pixels for a single event") - ###################################################################### # let’s now take a look at the pedestal-subtracted sums and a # pedestal-subtracted signal: # +peds + +###################################################################### plt.plot(sums - peds) plt.xlabel("pixel id") plt.ylabel("Pedestal-subtracted Signal") - ###################################################################### # Now, we can clearly see that the signal is centered at 0 where there is # no Cherenkov light, and we can also clearly see the shower around pixel -# 250. +# 700. # # we can also subtract the pedestals from the traces themselves, which would be needed to compare peaks properly -for ii in range(270, 280): - plt.plot(data[ii] - peds[ii], drawstyle="steps", label="pix{}".format(ii)) +for pixel in range(700, 710): + plt.plot( + data[pixel] - peds[pixel], drawstyle="steps", label="pixel {}".format(pixel) + ) plt.legend() @@ -246,13 +245,14 @@ camgeom = source.subarray.tel[24].camera.geometry ###################################################################### -title = "CT24, run {} event {} ped-sub".format(event.index.obs_id, event.index.event_id) +title = "CT24, obs_id {} event_id {} ped-sub".format( + event.index.obs_id, event.index.event_id +) + disp = CameraDisplay(camgeom, title=title) + disp.image = sums - peds -disp.cmap = plt.cm.RdBu_r disp.add_colorbar() -disp.set_limits_percent(95) # autoscale - ###################################################################### # It looks like a nice signal! We have plotted our pedestal-subtracted @@ -264,19 +264,19 @@ # pedestals for each) # -for tel in event.r0.tel.keys(): +for tel_id, tel_event in event.tel.items(): plt.figure() - camgeom = source.subarray.tel[tel].camera.geometry - title = "CT{}, run {} event {}".format( - tel, event.index.obs_id, event.index.event_id + camgeom = source.subarray.tel[tel_id].camera.geometry + title = "obs_id {}, event_id {}, tel_id {}".format( + event.index.obs_id, + event.index.event_id, + tel_id, ) disp = CameraDisplay(camgeom, title=title) - disp.image = event.r0.tel[tel].waveform[0].sum(axis=1) - disp.cmap = plt.cm.RdBu_r + disp.image = tel_event.r0.waveform[0].sum(axis=1) disp.add_colorbar() disp.set_limits_percent(95) - ###################################################################### # some signal processing… # ----------------------- diff --git a/examples/tutorials/tableloader_and_eventsources.py b/examples/tutorials/tableloader_and_eventsources.py index 815654a40ca..428f0830900 100644 --- a/examples/tutorials/tableloader_and_eventsources.py +++ b/examples/tutorials/tableloader_and_eventsources.py @@ -6,7 +6,7 @@ at CEA Paris-Saclay (J. Hackfeld). """ -# %% +###################################################################### # Introduction # ------------ # ``ctapipe`` provides basically two different ways of accessing its data products: @@ -57,7 +57,7 @@ # **It is recommended to use the** ``TableLoader`` **when loading data into Tables, # because its much faster than EventSources!** -# %% +###################################################################### # Code Examples # ------------- # First import some classes/modules and get the example dataset path: @@ -70,7 +70,7 @@ simtel_path = utils.get_dataset_path("gamma_prod5.simtel.zst") -# %% +###################################################################### # EventSource(s) # -------------- # The already implemented EventSources are: @@ -80,25 +80,25 @@ for source in sources.values(): print(f"{source.__name__:>{maxlen}s} -- {source.__module__}.{source.__qualname__}") -# %% +###################################################################### # ``EventSource`` will create an appropriate subclass if a compatible source is found for # the given ``input_url`` : source = EventSource(input_url=simtel_path, max_events=5) print(source) -# %% +###################################################################### # You can now loop over the ``ctapipe.containers.ArrayEventContainer`` generated by the # source. for event in source: print(event.count) -# %% +###################################################################### print(repr(event)) -# %% +###################################################################### # Every time a new loop is started through the source, # it tries to restart from the first event, which might not be supported # by the event source. It is encouraged to use ``EventSource`` in a **context manager** @@ -108,10 +108,10 @@ for event in source: print( f"Event Count: {event.count}," - f"Tels with trigger: {event.trigger.tels_with_trigger}" + f"Tels with trigger: {event.dl0.trigger.tels_with_trigger}" ) -# %% +###################################################################### # You can hand in the ID's of the telescopes to be included in the data # with the ``allowed_tels`` attribute. If given, only this subset of telescopes # will be present in the generated events. If None, all available telescopes are used. @@ -124,10 +124,10 @@ for event in source: print( f"Event Count: {event.count}," - f"Tels with trigger: {event.trigger.tels_with_trigger}" + f"Tels with trigger: {event.dl0.trigger.tels_with_trigger}" ) -# %% +###################################################################### # If you want to calibrate your data in the event loop and write it to an .h5 file with # the ``DataWriter`` : # @@ -146,7 +146,7 @@ calibrate(event) write_data(event) -# %% +###################################################################### # Alternatively doing it with ``ctapipe-process`` would look like this: # # .. code-block:: bash @@ -154,7 +154,7 @@ # ! ctapipe-process -i {simtel_path} -o events.dl1.h5 --overwrite --progress # -# %% +###################################################################### # TableLoader # ----------- # @@ -162,7 +162,7 @@ loader = TableLoader(input_url="events.dl1.h5") -# %% +###################################################################### # Alternatively using a config file: tableloader_config = Config( @@ -175,7 +175,7 @@ loader = TableLoader(config=tableloader_config) -# %% +###################################################################### # Reading subarray-based event information: subarray_events = loader.read_subarray_events( @@ -188,7 +188,7 @@ subarray_events -# %% +###################################################################### # Reading subarray-based event information in chunks: subarray_events_chunked = loader.read_subarray_events_chunked( @@ -201,7 +201,7 @@ for chunk in subarray_events_chunked: print(" \n", chunk) -# %% +###################################################################### # Reading just LST events: lst_events = loader.read_telescope_events( @@ -221,7 +221,7 @@ lst_events -# %% +###################################################################### # Loading telescope events by type returns a dict with the different telescope types: telescope_events_by_type = loader.read_telescope_events_by_type( @@ -242,7 +242,7 @@ for tel_type, table in telescope_events_by_type.items(): print(f"Telescope Type: {tel_type} \n", table, "\n") -# %% +###################################################################### # Loading telescope events by ID returns a dict with the different telescope IDs: telescope_events_by_id = loader.read_telescope_events_by_id( @@ -262,14 +262,14 @@ for tel_id, table in telescope_events_by_id.items(): print(f"Telescope ID: {tel_id} \n", table, "\n") -# %% +###################################################################### # - ``read_telescope_events_chunked`` # - ``read_telescope_events_by_type_chunked`` # - ``read_telescope_events_by_id_chunked`` # # are also available. -# %% +###################################################################### # Reading e.g. simulation- or observation-information: simulation_configuration = loader.read_simulation_configuration() @@ -277,7 +277,7 @@ simulation_configuration -# %% +###################################################################### # Now you have ``astropy.table`` s including all the relevant data you need # for your analyses. # diff --git a/examples/tutorials/theta_square.py b/examples/tutorials/theta_square.py index 4387ac1beeb..8adcd0d5bcf 100644 --- a/examples/tutorials/theta_square.py +++ b/examples/tutorials/theta_square.py @@ -4,7 +4,6 @@ This is a basic example to analyze some events and make a :math:`\Theta^2` plot - """ # %matplotlib inline @@ -43,7 +42,7 @@ image_processor(event) shower_processor(event) - reco_result = event.dl2.stereo.geometry["HillasReconstructor"] + reco_result = event.dl2.geometry["HillasReconstructor"] # get angular offset between reconstructed shower direction and MC # generated shower direction @@ -55,7 +54,6 @@ # Appending all estimated off angles off_angles.append(off_angle.to(u.deg).value) - ###################################################################### # calculate theta square for angles which are not nan # @@ -63,7 +61,6 @@ off_angles = np.array(off_angles) thetasquare = off_angles[np.isfinite(off_angles)] ** 2 - ###################################################################### # Plot the results # ---------------- diff --git a/examples/visualization/array_display.py b/examples/visualization/array_display.py index f0c6fb6a85a..0145cf94324 100644 --- a/examples/visualization/array_display.py +++ b/examples/visualization/array_display.py @@ -11,6 +11,7 @@ import numpy as np from astropy import units as u from astropy.coordinates import SkyCoord +from matplotlib.colors import ListedColormap from ctapipe.calib import CameraCalibrator from ctapipe.coordinates import EastingNorthingFrame @@ -49,7 +50,6 @@ disp = ArrayDisplay(subarray) - ###################################################################### # You can specify the Frame you want as long as it is compatible with # ``GroundFrame``. ``EastingNorthingFrame`` is probably the most useful. @@ -59,7 +59,6 @@ disp = ArrayDisplay(subarray, frame=EastingNorthingFrame()) disp.add_labels() - ###################################################################### # Using color to show information # ------------------------------- @@ -69,24 +68,19 @@ # the ``values`` attribute, like a trigger pattern # -plt.set_cmap("rainbow") # the array display will use the current colormap for values + +cmap = ListedColormap(["xkcd:red", "xkcd:blue"]) ad = ArrayDisplay(subarray) +ad.telescopes.set_cmap(cmap) ad.telescopes.set_linewidth(0) # to turn off the telescope borders trigger_pattern = np.zeros(subarray.n_tels) -trigger_pattern[ - [ - 1, - 4, - 5, - 6, - ] -] = 1 +trigger_pattern[[1, 4, 5, 6]] = 1 + ad.values = trigger_pattern # display certain telescopes in a color ad.add_labels() - ###################################################################### # or for example, you could use color to represent the telescope distance # to the impact point @@ -94,18 +88,22 @@ shower_impact = SkyCoord(200 * u.m, -200 * u.m, 0 * u.m, frame=EastingNorthingFrame()) -plt.set_cmap("rainbow") # the array display will use the current colormap for values + ad = ArrayDisplay(subarray) + +ad.telescopes.set_cmap("inferno") ad.telescopes.set_linewidth(0) # to turn off the telescope borders -plt.scatter(shower_impact.easting, shower_impact.northing, marker="+", s=200) + +ad.axes.scatter(shower_impact.easting, shower_impact.northing, marker="+", s=200) distances = np.hypot( subarray.tel_coords.cartesian.x - shower_impact.cartesian.x, subarray.tel_coords.cartesian.y - shower_impact.cartesian.y, ) + ad.values = distances -plt.colorbar(ad.telescopes, label="Distance (m)") +plt.colorbar(ad.telescopes, label="Distance / m") ###################################################################### # Overlaying vectors @@ -120,15 +118,14 @@ # set of parameters. # -np.random.seed(0) -phis = np.random.uniform(0, 180.0, size=subarray.n_tels) * u.deg -rhos = np.ones(subarray.n_tels) * 200 * u.m +rng = np.random.default_rng(0) +phis = rng.uniform(0, 180.0, size=subarray.n_tels) * u.deg +rhos = np.full(subarray.n_tels, 200) * u.m ad = ArrayDisplay(subarray, frame=EastingNorthingFrame(), tel_scale=2) ad.set_vector_rho_phi(rho=rhos, phi=phis) - ###################################################################### # Overlaying Image Axes # --------------------- @@ -154,11 +151,11 @@ def plot_event(event, subarray, ax): true and reconstructed impact position overlaid """ - event.pointing.array_azimuth + event.pointing.azimuth disp = ArrayDisplay(subarray, axes=ax) - hillas_dict = {tid: tel.parameters.hillas for tid, tel in event.dl1.tel.items()} - core_dict = {tid: tel.parameters.core.psi for tid, tel in event.dl1.tel.items()} + hillas_dict = {tid: tel.dl1.parameters.hillas for tid, tel in event.tel.items()} + core_dict = {tid: tel.dl1.parameters.core.psi for tid, tel in event.tel.items()} disp.set_line_hillas( hillas_dict, @@ -166,7 +163,7 @@ def plot_event(event, subarray, ax): 500, ) - reco_shower = event.dl2.stereo.geometry["HillasReconstructor"] + reco_shower = event.dl2.geometry["HillasReconstructor"] ax.scatter( event.simulation.shower.core_x, diff --git a/examples/visualization/camera_display.py b/examples/visualization/camera_display.py index 56eece592f6..f631cdaec4d 100644 --- a/examples/visualization/camera_display.py +++ b/examples/visualization/camera_display.py @@ -343,9 +343,9 @@ def update(frame): ) as source: event = next(iter(source)) -tel_id = list(event.r0.tel.keys())[0] +tel_id = next(iter(event.tel)) geom = source.subarray.tel[tel_id].camera.geometry -waveform = event.r0.tel[tel_id].waveform +waveform = event.tel[tel_id].r0.waveform n_chan, n_pix, n_samp = waveform.shape diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index baf3d2f1057..e70932ac8b1 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -8,7 +8,13 @@ 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, + TelescopeTriggerContainer, +) from ctapipe.core import TelescopeComponent from ctapipe.core.traits import ( BoolTelescopeParameter, @@ -162,8 +168,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: self.log.debug( "Encountered an event with no R1 data. " @@ -174,8 +180,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: self.log.warning( "Encountered an event with no DL0 data. " @@ -186,12 +192,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, @@ -207,7 +214,14 @@ 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: copy trigger from existing dl0 if there. + # FIXME: how to deal with trigger information at R1 for simulations? + if tel_event.dl0 is not None and tel_event.dl0.trigger is not None: + trigger = tel_event.dl0.trigger + else: + trigger = TelescopeTriggerContainer() + + tel_event.dl0 = DL0TelescopeContainer( event_type=r1.event_type, event_time=r1.event_time, waveform=dl0_waveform, @@ -215,25 +229,28 @@ def _calibrate_dl0(self, event, tel_id): pixel_status=dl0_pixel_status, first_cell_id=r1.first_cell_id, calibration_monitoring_id=r1.calibration_monitoring_id, + trigger=trigger, ) - 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_channels, 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_channels, n_pixels, - event.mon.tel[tel_id].pixel_status, + tel_event.mon.pixel_status, selected_gain_channel, ) pixel_index = _get_pixel_index(n_pixels) - dl1_calib = event.calibration.tel[tel_id].dl1 + dl1_calib = tel_event.calibration.dl1 readout = self.subarray.tel[tel_id].camera.readout # subtract any remaining pedestal before extraction @@ -250,7 +267,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=np.squeeze(waveforms).astype(np.float32), peak_time=np.zeros(n_pixels, dtype=np.float32), is_valid=True, @@ -309,7 +326,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): """ @@ -320,13 +337,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/src/ctapipe/calib/camera/flatfield.py b/src/ctapipe/calib/camera/flatfield.py index c3f74678aae..e6aefab5974 100644 --- a/src/ctapipe/calib/camera/flatfield.py +++ b/src/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 @@ -100,14 +100,14 @@ def __init__(self, subarray, **kwargs): self.log.info(f"extractor {self.extractor}") @abstractmethod - def calculate_relative_gain(self, event): + def calculate_relative_gain(self, event) -> bool: """ Calculate the flat-field statistics and fill the mon.tel[tel_id].flatfield container Parameters ---------- - event: ctapipe.containers.ArrayEventContainer + event: ctapipe.containers.SubarrayEventContainer Returns: True if the mon.tel[tel_id].flatfield is updated, False otherwise @@ -168,7 +168,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 @@ -181,13 +181,14 @@ def _extract_charge(self, event) -> DL1CameraContainer: DL1CameraContainer """ - waveforms = event.r1.tel[self.tel_id].waveform + tel_event = event.tel[self.tel_id] + waveforms = tel_event.r1.waveform n_channels, n_pixels, _ = waveforms.shape - selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel + selected_gain_channel = tel_event.r1.selected_gain_channel broken_pixels = _get_invalid_pixels( n_channels=n_channels, n_pixels=n_pixels, - pixel_status=event.mon.tel[self.tel_id].pixel_status, + pixel_status=tel_event.mon.pixel_status, selected_gain_channel=selected_gain_channel, ) # Extract charge and time @@ -196,7 +197,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): """ @@ -210,21 +211,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 + tel_event = event.tel[self.tel_id] + waveform = tel_event.r1.waveform + container = tel_event.mon.flatfield # re-initialize counter if self.n_events_seen == self.sample_size: self.n_events_seen = 0 - trigger_time = event.trigger.time + # real data + trigger_time = event.dl0.trigger.time 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, + tel_event.mon.pixel_status.hardware_failing_pixels, + tel_event.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, + tel_event.mon.pixel_status.flatfield_failing_pixels, ) if self.n_events_seen == 0: @@ -233,7 +236,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/src/ctapipe/calib/camera/pedestals.py b/src/ctapipe/calib/camera/pedestals.py index ba8ee336838..2ffd834ee71 100644 --- a/src/ctapipe/calib/camera/pedestals.py +++ b/src/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,26 +197,27 @@ 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 + tel_event = event.tel[self.tel_id] + waveforms = tel_event.r1.waveform n_channels, n_pixels, _ = waveforms.shape - selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel + selected_gain_channel = tel_event.r1.selected_gain_channel broken_pixels = _get_invalid_pixels( n_channels=n_channels, n_pixels=n_pixels, - pixel_status=event.mon.tel[self.tel_id].pixel_status, + pixel_status=tel_event.mon.pixel_status, selected_gain_channel=selected_gain_channel, ) # Extract charge and time @@ -225,7 +226,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 +239,18 @@ def calculate_pedestals(self, event): event : general event container """ + tel_event = event.tel[self.tel_id] # initialize the np array at each cycle - waveform = event.r1.tel[self.tel_id].waveform - container = event.mon.tel[self.tel_id].pedestal + waveform = tel_event.r1.waveform + container = tel_event.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 - pixel_mask = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels + trigger_time = tel_event.dl0.trigger.time + pixel_mask = tel_event.mon.pixel_status.hardware_failing_pixels if self.n_events_seen == 0: self.time_start = trigger_time @@ -256,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/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index b5a3875f62e..b5e77caf7f1 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/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,41 +106,55 @@ 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() 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((1, 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((1, 2048, 128), 2)), + ) 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() + 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 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): + rng = np.random.default_rng(1) + tel_id = 1 # copy because we mutate the camera, should not affect other tests 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 @@ -143,7 +164,6 @@ 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) camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma) @@ -156,31 +176,35 @@ def test_dl1_charge_calib(example_subarray): gain_channel = [1, 2] for n_channels in gain_channel: # Randomize times and create pulses - time_offset = random.uniform(-10, +10, (n_channels, n_pixels)) + time_offset = rng.uniform(-10, +10, (n_channels, n_pixels)) y = norm.pdf(x, mid + time_offset[..., np.newaxis], pulse_sigma).astype( "float32" ) # Define absolute calibration coefficients - absolute = random.uniform(100, 1000, (n_channels, n_pixels)).astype("float32") + absolute = rng.uniform(100, 1000, (n_channels, n_pixels)).astype("float32") y *= absolute[..., np.newaxis] # Define relative coefficients - relative = random.normal(1, 0.01, (n_channels, n_pixels)) + relative = rng.normal(1, 0.01, (n_channels, n_pixels)) y /= relative[..., np.newaxis] # Define pedestal - pedestal = random.uniform(-4, 4, (n_channels, n_pixels)) + pedestal = rng.uniform(-4, 4, (n_channels, n_pixels)) y += pedestal[..., np.newaxis] - event = ArrayEventContainer() - tel_id = list(subarray.tel.keys())[0] - event.dl0.tel[tel_id].waveform = y selected_gain_channel = None if n_channels == 1: selected_gain_channel = np.zeros(n_pixels, dtype=int) - event.dl0.tel[tel_id].selected_gain_channel = selected_gain_channel - event.r1.tel[tel_id].selected_gain_channel = selected_gain_channel + + 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=selected_gain_channel, + ), + ) # Test default calibrator = CameraCalibrator( @@ -188,16 +212,16 @@ def test_dl1_charge_calib(example_subarray): ) calibrator(event) np.testing.assert_allclose( - event.dl1.tel[tel_id].image, y.sum(-1).squeeze(), rtol=1e-4 + event.tel[tel_id].dl1.image, y.sum(-1).squeeze(), 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 @@ -206,13 +230,13 @@ def test_dl1_charge_calib(example_subarray): ) # 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 @@ -221,9 +245,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.squeeze(), atol=1 + event.tel[tel_id].dl1.peak_time, expected_peak_time.squeeze(), atol=1 ) # We now use GlobalPeakWindowSum to see the effect of missing charge @@ -236,9 +260,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 @@ -247,9 +271,9 @@ 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.tel[tel_id].dl1.image, 1, rtol=0.1) assert not np.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 ) @@ -305,22 +329,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 @@ -329,8 +353,8 @@ 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 event.tel[tel_id].dl1.image[0] == 9999 + assert event.tel[tel_id].dl1.peak_time[0] == 10.0 / sampling_rate def test_no_gain_selection(prod5_gamma_simtel_path): @@ -341,15 +365,15 @@ def test_no_gain_selection(prod5_gamma_simtel_path): tested_n_channels = set() - for tel_id in event.r1.tel: + for tel_id, tel_event in event.tel.items(): readout = source.subarray.tel[tel_id].camera.readout tested_n_channels.add(readout.n_channels) calibrator = CameraCalibrator(subarray=source.subarray) calibrator(event) - image = event.dl1.tel[tel_id].image - peak_time = event.dl1.tel[tel_id].peak_time + image = tel_event.dl1.image + peak_time = tel_event.dl1.peak_time assert image.ndim == 2 assert peak_time.ndim == 2 assert image.shape == (readout.n_channels, readout.n_pixels) diff --git a/src/ctapipe/calib/camera/tests/test_flatfield.py b/src/ctapipe/calib/camera/tests/test_flatfield.py index 7c68f613af9..2c61a71fd39 100644 --- a/src/ctapipe/calib/camera/tests/test_flatfield.py +++ b/src/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,47 +38,47 @@ 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)) + event.tel[tel_id].r1.waveform = np.zeros((n_gain, n_pixels, 40)) # 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/src/ctapipe/calib/camera/tests/test_pedestals.py b/src/ctapipe/calib/camera/tests/test_pedestals.py index af37677d79f..7a9d85abb40 100644 --- a/src/ctapipe/calib/camera/tests/test_pedestals.py +++ b/src/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,23 +40,23 @@ 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) + event.tel[tel_id].r1.waveform = np.full((2, n_pixels, 40), ped_level) 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/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index 4d743e5d217..2bb8f388d61 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -659,7 +659,7 @@ def dl1_mon_pointing_file(dl1_file, dl1_tmp_path): path = dl1_tmp_path / "dl1_mon_ponting.dl1.h5" shutil.copy(dl1_file, path) - events = read_table(path, "/dl1/event/subarray/trigger") + events = read_table(path, "/dl0/event/subarray/trigger") subarray = SubarrayDescription.from_hdf(path) # create some dummy monitoring data diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 9a5b39f4da8..f4f032cc446 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -13,56 +13,66 @@ from .core import Container, Field, Map __all__ = [ - "ArrayEventContainer", + "BaseHillasParametersContainer", + "BaseTimingParametersContainer", + "CameraHillasParametersContainer", + "CameraTimingParametersContainer", "ConcentrationContainer", - "DL0CameraContainer", - "DL0Container", + "CoordinateFrameType", + "CoreParametersContainer", + "DispContainer", + "DL0SubarrayContainer", + "DL0TelescopeContainer", "DL1CameraCalibrationContainer", - "DL1CameraContainer", - "DL1Container", - "DL2Container", - "EventCalibrationContainer", - "EventCameraCalibrationContainer", - "EventIndexContainer", + "DL1TelescopeContainer", + "DL2SubarrayContainer", + "DL2TelescopeContainer", "EventType", "FlatFieldContainer", "HillasParametersContainer", - "CoreParametersContainer", "ImageParametersContainer", + "IntensityStatisticsContainer", "LeakageContainer", - "MonitoringCameraContainer", - "MonitoringContainer", "MorphologyContainer", - "BaseHillasParametersContainer", - "CameraHillasParametersContainer", - "CameraTimingParametersContainer", + "MuonContainer", + "MuonEfficiencyContainer", + "MuonParametersContainer", + "MuonRingContainer", + "ObservationBlockContainer", + "ObservationBlockState", + "ObservingMode", "ParticleClassificationContainer", + "PeakTimeStatisticsContainer", "PedestalContainer", + "PixelStatus", "PixelStatusContainer", - "R0CameraContainer", - "R0Container", - "R1CameraContainer", - "R1Container", - "ReconstructedContainer", + "PointingMode", + "R0TelescopeContainer", + "R1TelescopeContainer", "ReconstructedEnergyContainer", "ReconstructedGeometryContainer", - "DispContainer", - "SimulatedCameraContainer", + "ReconstructedShowerContainer", + "SchedulingBlockContainer", + "SchedulingBlockType", "SimulatedShowerContainer", "SimulatedShowerDistribution", "SimulationConfigContainer", - "TelEventIndexContainer", - "BaseTimingParametersContainer", + "SimulationSubarrayContainer", + "SimulationTelescopeContainer", + "StatisticsContainer", + "SubarrayEventContainer", + "SubarrayEventIndexContainer", + "SubarrayPointingContainer", + "SubarrayTriggerContainer", + "TelescopeCalibrationContainer", + "TelescopeEventContainer", + "TelescopeEventIndexContainer", + "TelescopeImpactParameterContainer", + "TelescopeMonitoringContainer", + "TelescopePointingContainer", + "TelescopeTriggerContainer", "TimingParametersContainer", - "TriggerContainer", "WaveformCalibrationContainer", - "StatisticsContainer", - "IntensityStatisticsContainer", - "PeakTimeStatisticsContainer", - "SchedulingBlockContainer", - "ObservationBlockContainer", - "ObservingMode", - "ObservationBlockState", ] @@ -238,10 +248,10 @@ class TelescopeConfigurationIndexContainer(Container): tel_id = tel_id_field() -class EventIndexContainer(Container): +class SubarrayEventIndexContainer(Container): """Index columns to include in event lists. - Common to all data levels + Common to all data levels. """ default_prefix = "" @@ -249,7 +259,7 @@ class EventIndexContainer(Container): event_id = event_id_field() -class TelEventIndexContainer(Container): +class TelescopeEventIndexContainer(Container): """ index columns to include in telescope-wise event lists @@ -262,6 +272,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 @@ -476,7 +505,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. @@ -516,15 +545,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 @@ -555,7 +575,7 @@ class DL1CameraCalibrationContainer(Container): ) -class R0CameraContainer(Container): +class R0TelescopeContainer(Container): """ Storage of raw data from a single telescope """ @@ -565,18 +585,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 """ @@ -635,20 +644,9 @@ class R1CameraContainer(Container): ) -class R1Container(Container): +class DL0TelescopeContainer(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): - """ - 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 @@ -657,6 +655,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, ( @@ -696,14 +699,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", ) @@ -742,10 +745,9 @@ class SimulatedShowerContainer(Container): ) -class SimulatedCameraContainer(Container): +class SimulationTelescopeContainer(Container): """ - True images and parameters derived from them, analogous to the `DL1CameraContainer` - but for simulated data. + Per-telescope simulations truths and derived parameters. """ default_prefix = "" @@ -773,12 +775,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): @@ -865,28 +866,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 @@ -1014,16 +993,13 @@ class DispContainer(Container): ) -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 container. + """ geometry = Field( default_factory=partial(Map, ReconstructedGeometryContainer), @@ -1039,35 +1015,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 @@ -1081,18 +1048,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 """ @@ -1103,18 +1066,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""" @@ -1134,6 +1085,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") @@ -1147,6 +1100,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, @@ -1161,7 +1116,7 @@ class MuonParametersContainer(Container): ) -class MuonTelescopeContainer(Container): +class MuonContainer(Container): """ Container for muon analysis """ @@ -1175,15 +1130,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 @@ -1327,7 +1273,7 @@ class WaveformCalibrationContainer(Container): ) -class MonitoringCameraContainer(Container): +class TelescopeMonitoringContainer(Container): """ Container for camera monitoring data """ @@ -1350,18 +1296,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 @@ -1388,36 +1322,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", + ) + + simulation = Field( + None, + description="Simulated Event Information", + type=SimulationTelescopeContainer, ) - r0 = Field(default_factory=R0Container, description="Raw Data") - r1 = Field(default_factory=R1Container, description="R1 Calibrated Data") + + 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( @@ -1425,6 +1367,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. diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py index e96b2bef48f..1bba48e7128 100644 --- a/src/ctapipe/image/cleaning.py +++ b/src/ctapipe/image/cleaning.py @@ -28,9 +28,7 @@ import numpy as np -from ctapipe.image.statistics import n_largest - -from ..containers import MonitoringCameraContainer +from ..containers import TelescopeMonitoringContainer from ..core import TelescopeComponent from ..core.traits import ( BoolTelescopeParameter, @@ -38,6 +36,7 @@ IntTelescopeParameter, ) from .morphology import brightest_island, largest_island, number_of_islands +from .statistics import n_largest def tailcuts_clean( @@ -647,9 +646,9 @@ def __call__( self, tel_id: int, image: np.ndarray, - arrival_times: np.ndarray = None, + arrival_times: np.ndarray | None = None, *, - monitoring: MonitoringCameraContainer = None, + monitoring: TelescopeMonitoringContainer | None = None, ) -> np.ndarray: """ Identify pixels with signal, and reject those with pure noise. @@ -663,9 +662,9 @@ def __call__( image pixel data corresponding to the camera geometry arrival_times : np.ndarray image of arrival time (not used in this method) - monitoring : ctapipe.containers.MonitoringCameraContainer - MonitoringCameraContainer to make use of additional parameters - from monitoring data e.g. pedestal std. + monitoring: `ctapipe.containers.TelescopeMonitoringContainer` + Monitoring information to make use of additional parameters + e.g. pedestal std. Returns ------- @@ -703,9 +702,9 @@ def __call__( self, tel_id: int, image: np.ndarray, - arrival_times: np.ndarray = None, + arrival_times: np.ndarray | None = None, *, - monitoring: MonitoringCameraContainer = None, + monitoring: TelescopeMonitoringContainer | None = None, ) -> np.ndarray: """ Apply standard picture-boundary cleaning. See `ImageCleaner.__call__()` @@ -790,16 +789,16 @@ def __call__( self, tel_id: int, image: np.ndarray, - arrival_times: np.ndarray = None, + arrival_times: np.ndarray | None = None, *, - monitoring: MonitoringCameraContainer = None, + monitoring: TelescopeMonitoringContainer | None = None, ) -> np.ndarray: """ Apply NSB image cleaning used by lstchain. See `ImageCleaner.__call__()` """ pedestal_std = None if monitoring is not None: - pedestal_std = monitoring.tel[tel_id].pedestal.charge_std + pedestal_std = monitoring.pedestal.charge_std return nsb_image_cleaning( self.subarray.tel[tel_id].camera.geometry, @@ -829,9 +828,9 @@ def __call__( self, tel_id: int, image: np.ndarray, - arrival_times: np.ndarray = None, + arrival_times: np.ndarray | None = None, *, - monitoring: MonitoringCameraContainer = None, + monitoring: TelescopeMonitoringContainer | None = None, ) -> np.ndarray: """ Apply MARS-style image cleaning. See `ImageCleaner.__call__()` @@ -861,9 +860,9 @@ def __call__( self, tel_id: int, image: np.ndarray, - arrival_times: np.ndarray = None, + arrival_times: np.ndarray | None = None, *, - monitoring: MonitoringCameraContainer = None, + monitoring: TelescopeMonitoringContainer | None = None, ) -> np.ndarray: """Apply FACT-style image cleaning. see ImageCleaner.__call__()""" @@ -896,9 +895,9 @@ def __call__( self, tel_id: int, image: np.ndarray, - arrival_times: np.ndarray = None, + arrival_times: np.ndarray | None = None, *, - monitoring: MonitoringCameraContainer = None, + monitoring: TelescopeMonitoringContainer | None = None, ) -> np.ndarray: """ Apply MAGIC-like image cleaning with timing information. See `ImageCleaner.__call__()` diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index ce772744c58..f08e0d4b449 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/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, @@ -428,7 +428,7 @@ def _apply_correction(charge, correction, selected_gain_channel): @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. @@ -451,7 +451,7 @@ def __call__( Returns ------- - DL1CameraContainer: + DL1TelescopeContainer: extracted images and validity flags """ @@ -463,7 +463,7 @@ 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] ) @@ -473,7 +473,7 @@ def __call__( charge = charge[0] peak_time = peak_time[0] - 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): @@ -510,7 +510,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], @@ -527,7 +527,7 @@ def __call__( charge = charge[0] peak_time = peak_time[0] - 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): @@ -578,7 +578,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.mean( @@ -614,7 +614,7 @@ def __call__( charge = charge[0] peak_time = peak_time[0] - 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): @@ -650,7 +650,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, @@ -668,7 +668,7 @@ def __call__( charge = charge[0] peak_time = peak_time[0] - 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): @@ -716,7 +716,7 @@ 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] ) @@ -730,7 +730,7 @@ def __call__( charge = charge[0] peak_time = peak_time[0] - 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): @@ -772,7 +772,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__( charge = charge[0] peak_time = peak_time[0] - 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,18 +1263,17 @@ def _apply_second_pass( def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels - ) -> DL1CameraContainer: + ) -> DL1TelescopeContainer: if waveforms.shape[-3] != 1: raise AttributeError( "The data needs to be gain selected to use the TwoPassWindowSum." ) waveforms = waveforms[0, :, :] - 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, @@ -1290,7 +1289,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, @@ -1647,7 +1646,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] @@ -1705,4 +1704,4 @@ def __call__( charge = charge[0] peak_time = peak_time[0] - 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/src/ctapipe/image/image_processor.py b/src/ctapipe/image/image_processor.py index bb810283a71..b3db8eb60af 100644 --- a/src/ctapipe/image/image_processor.py +++ b/src/ctapipe/image/image_processor.py @@ -9,12 +9,13 @@ from ctapipe.coordinates import TelescopeFrame from ..containers import ( - ArrayEventContainer, CameraHillasParametersContainer, CameraTimingParametersContainer, ImageParametersContainer, IntensityStatisticsContainer, PeakTimeStatisticsContainer, + SubarrayEventContainer, + TelescopeEventContainer, TimingParametersContainer, ) from ..core import QualityQuery, TelescopeComponent @@ -62,6 +63,7 @@ class ImageQualityQuery(QualityQuery): class ImageProcessor(TelescopeComponent): """ Takes DL1/Image data and cleans and parametrizes the images into DL1/parameters. + Should be run after CameraCalibrator to produce all DL1 information. """ @@ -117,8 +119,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, @@ -210,51 +213,50 @@ 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 + Process the calibrated image into parameters for a single telescope event. """ - for tel_id, dl1_camera in event.dl1.tel.items(): - if self.apply_image_modifier.tel[tel_id]: - dl1_camera.image = self.modify(tel_id=tel_id, image=dl1_camera.image) - - dl1_camera.image_mask = self.clean( - tel_id=tel_id, - image=dl1_camera.image, - arrival_times=dl1_camera.peak_time, - monitoring=event.mon.tel[tel_id], - ) + tel_id = tel_event.index.tel_id + dl1 = tel_event.dl1 + + if self.apply_image_modifier.tel[tel_id]: + dl1.image = self.modify(tel_id=tel_id, image=dl1.image) + + dl1.image_mask = self.clean( + tel_id=tel_id, + image=dl1.image, + arrival_times=dl1.peak_time, + monitoring=tel_event.mon, + ) - 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, + 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, + ) + + 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/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index 30b382aef2c..8e87e19756a 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/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"] @@ -79,7 +79,10 @@ class RingQuery(QualityQuery): class MuonProcessor(TelescopeComponent): """ - Takes cleaned images and extracts muon rings. Should be run after ImageProcessor. + Takes cleaned images and extracts muon rings. + + Relies on image parameters for quality query, thus should be run after + the ImageProcessor. """ completeness_threshold = FloatTelescopeParameter( @@ -114,23 +117,20 @@ 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): + """ + Process all telescope events of the given subarray event. + """ + 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 - 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 +139,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 +152,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 +176,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 +195,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/src/ctapipe/image/muon/tests/test_processor.py b/src/ctapipe/image/muon/tests/test_processor.py index e2a7211ed19..513da086fb2 100644 --- a/src/ctapipe/image/muon/tests/test_processor.py +++ b/src/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/src/ctapipe/image/reducer.py b/src/ctapipe/image/reducer.py index 535aca63355..b607f7fe25d 100644 --- a/src/ctapipe/image/reducer.py +++ b/src/ctapipe/image/reducer.py @@ -6,7 +6,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, @@ -178,15 +178,13 @@ def __init__( self.image_extractor_type = [("type", "*", name)] self.image_extractors[name] = image_extractor - def select_pixels(self, waveforms, tel_id=None, selected_gain_channel=None): + def select_pixels(self, waveforms, tel_id, selected_gain_channel=None): camera_geom = self.subarray.tel[tel_id].camera.geometry # Pulse-integrate waveforms extractor = self.image_extractors[self.image_extractor_type.tel[tel_id]] # do not treat broken pixels in data volume reduction - broken_pixels = np.zeros( - (waveforms.shape[-3], camera_geom.n_pixels), dtype=bool - ) - dl1: DL1CameraContainer = extractor( + broken_pixels = np.zeros((waveforms.shape[0], camera_geom.n_pixels), dtype=bool) + dl1: DL1TelescopeContainer = extractor( waveforms, tel_id=tel_id, selected_gain_channel=selected_gain_channel, diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 067a07801aa..c9ade19e71e 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -902,18 +902,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 - - waveforms = event.r1.tel[tel_id].waveform - selected_gain_channel = np.zeros(waveforms.shape[-2], dtype=np.int64) - broken_pixels = event.mon.tel[ - tel_id - ].pixel_status.hardware_failing_pixels + for tel_id, tel_event in filter(is_flashcam, event.tel.items()): + true_charge = tel_event.simulation.true_image + + waveforms = tel_event.r1.waveform + selected_gain_channel = np.zeros(waveforms.shape[0], dtype=np.int64) + broken_pixels = tel_event.mon.pixel_status.hardware_failing_pixels dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid diff --git a/src/ctapipe/image/tests/test_image_processor.py b/src/ctapipe/image/tests/test_image_processor.py index 327b17aebf6..5000821aa38 100644 --- a/src/ctapipe/image/tests/test_image_processor.py +++ b/src/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/src/ctapipe/instrument/tests/test_trigger.py b/src/ctapipe/instrument/tests/test_trigger.py index f188b32a5db..87a68832c34 100644 --- a/src/ctapipe/instrument/tests/test_trigger.py +++ b/src/ctapipe/instrument/tests/test_trigger.py @@ -5,23 +5,26 @@ import pytest from numpy.testing import assert_equal -from ctapipe.containers import ArrayEventContainer +from ctapipe.containers import ( + DL0SubarrayContainer, + SubarrayEventContainer, + SubarrayTriggerContainer, + TelescopeEventContainer, +) from ctapipe.instrument import SubarrayDescription 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)) @@ -39,41 +42,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 not trigger(event) - 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 not trigger(event) - 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) - 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) - 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) - 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) - 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)))) @@ -116,10 +113,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): @@ -239,6 +234,14 @@ def test_software_trigger_simtel_process(tmp_path): assert len(events_no_trigger) > len(events_trigger) +def _dummy_event(tel_ids): + event = SubarrayEventContainer() + event.dl0.trigger.tels_with_trigger = tel_ids + for tel_id in tel_ids: + event.tel[tel_id] = TelescopeEventContainer() + return event + + def test_different_telescope_type_same_string(subarray_prod5_paranal): """ Test that for a subarray that contains two different telescope types @@ -278,27 +281,21 @@ def test_different_telescope_type_same_string(subarray_prod5_paranal): ) # all four LSTs, nothing to change - event = ArrayEventContainer() - event.trigger.tels_with_trigger = [1, 2, 3, 4] + event = _dummy_event([1, 2, 3, 4]) assert trigger(event) - assert_equal(event.trigger.tels_with_trigger, [1, 2, 3, 4]) + assert_equal(event.dl0.trigger.tels_with_trigger, [1, 2, 3, 4]) # two LSTs, nothing to change - event = ArrayEventContainer() - event.trigger.tels_with_trigger = [1, 2] + event = _dummy_event([1, 2]) assert trigger(event) - assert_equal(event.trigger.tels_with_trigger, [1, 2]) + assert_equal(event.dl0.trigger.tels_with_trigger, [1, 2]) # two LSTs, nothing to change - event = ArrayEventContainer() - event.trigger.tels_with_trigger = [1, 3] + event = _dummy_event([1, 3]) assert trigger(event) - assert_equal(event.trigger.tels_with_trigger, [1, 3]) + assert_equal(event.dl0.trigger.tels_with_trigger, [1, 3]) # one LST, remove - event = ArrayEventContainer() - event.trigger.tels_with_trigger = [ - 1, - ] + event = _dummy_event([1]) assert not trigger(event) - assert_equal(event.trigger.tels_with_trigger, []) + assert_equal(event.dl0.trigger.tels_with_trigger, []) diff --git a/src/ctapipe/instrument/trigger.py b/src/ctapipe/instrument/trigger.py index f5c519ee356..f2269c87a5f 100644 --- a/src/ctapipe/instrument/trigger.py +++ b/src/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 @@ -70,7 +70,7 @@ def __init__(self, subarray, *args, **kwargs): self._ids_by_type[tel_str] = set() self._ids_by_type[tel_str].update(self.subarray.get_tel_ids_for_type(tel)) - 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 @@ -92,7 +92,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) tels_in_event = tels_with_trigger.intersection(tel_ids) if len(tels_in_event) < min_tels: @@ -107,24 +107,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/src/ctapipe/io/datawriter.py b/src/ctapipe/io/datawriter.py index d0f93905f66..5042c257cf4 100644 --- a/src/ctapipe/io/datawriter.py +++ b/src/ctapipe/io/datawriter.py @@ -13,10 +13,10 @@ from traitlets import Dict, Instance from ..containers import ( - ArrayEventContainer, SimulatedShowerDistribution, + SubarrayEventContainer, TelescopeConfigurationIndexContainer, - TelEventIndexContainer, + 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 @@ -297,50 +289,180 @@ 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_trigger(event) - # write fixed pointing only for simulation, observed data will have monitoring - if self._is_simulation: - self._write_constant_pointing(event) + self._write_simulation_subarray(event) + self._write_dl0_subarray(event) + self._write_dl2_subarray(event) + + for tel_event in event.tel.values(): + 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_r1_waveforms: - self._write_r1_telescope_events(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) + self._write_constant_pointing(tel_event) if self.write_r0_waveforms: - self._write_r0_telescope_events(event) + self._write_r0_telescope(tel_event) + + if self.write_r1_waveforms: + self._write_r1_telescope(tel_event) - # write telescope event data - self._write_dl1_telescope_events(event) + # always write dl0, contains needed trigger info etc. + self._write_dl0_telescope(tel_event) + + if self.write_dl1_images or self.write_dl1_parameters: + self._write_dl1_telescope(tel_event) - # write DL2 info if requested if self.write_dl2: - 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_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_dl1_parameters: + self._writer.write( + table_name=f"dl1/event/telescope/parameters/{table_name}", + containers=[tel_index, *tel_event.dl1.parameters.values()], + ) + + if self.write_dl1_images: + if tel_event.dl1.image is None: + raise ValueError( + "DataWriter.write_dl1_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_dl1_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 _write_constant_pointing(self, event): + def _write_constant_pointing(self, tel_event): """ Write pointing configuration from event data assuming fixed pointing over the OB. @@ -353,19 +475,22 @@ def _write_constant_pointing(self, event): So we write the first pointing information for each telescope into the configuration table. """ - obs_id = event.index.obs_id - - for tel_id, pointing in event.pointing.tel.items(): - if tel_id in self._constant_telescope_pointing_written[obs_id]: - continue - index = TelescopeConfigurationIndexContainer( - obs_id=obs_id, - tel_id=tel_id, - ) - self._writer.write( - f"configuration/telescope/pointing/tel_{tel_id:03d}", (index, pointing) - ) - self._constant_telescope_pointing_written[obs_id].add(tel_id) + obs_id = tel_event.index.obs_id + tel_id = tel_event.index.tel_id + + if tel_id in self._constant_telescope_pointing_written[obs_id]: + return + + index = TelescopeConfigurationIndexContainer( + obs_id=obs_id, + tel_id=tel_id, + ) + table_name = self.table_name(tel_id) + self._writer.write( + f"configuration/telescope/pointing/{table_name}", + (index, tel_event.pointing), + ) + self._constant_telescope_pointing_written[obs_id].add(tel_id) def finish(self): """called after all events are done""" @@ -464,11 +589,12 @@ 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, ) + writer.exclude("/dl0/event/telescope/trigger", "trigger_pixels") writer.exclude("/dl1/event/telescope/trigger", "trigger_pixels") writer.exclude("/dl1/event/telescope/images/.*", "parameters") writer.exclude("/simulation/event/telescope/images/.*", "true_parameters") @@ -624,145 +750,9 @@ def fill_from_simtel( ) def table_name(self, tel_id): - """construct dataset table names depending on chosen split method""" + """Create table name for a telescope-wise table for given ``tel_id``""" 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 - """ - - 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_dl1_parameters: - self._writer.write( - table_name=f"dl1/event/telescope/parameters/{table_name}", - containers=[tel_index, *dl1_camera.parameters.values()], - ) - - if self.write_dl1_images: - if dl1_camera.image is None: - raise ValueError( - "DataWriter.write_dl1_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_dl1_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/src/ctapipe/io/eventsource.py b/src/ctapipe/io/eventsource.py index b6edb9859fe..09d508bec96 100644 --- a/src/ctapipe/io/eventsource.py +++ b/src/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 @@ -74,7 +74,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. @@ -302,7 +302,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/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index afbe2153d2c..ec2d0dfc282 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -11,39 +11,44 @@ from ctapipe.instrument.optics import FocalLengthKind from ..containers import ( - ArrayEventContainer, CameraHillasParametersContainer, CameraTimingParametersContainer, ConcentrationContainer, CoordinateFrameType, DispContainer, - 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 @@ -141,8 +146,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 ---------- @@ -210,13 +215,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("."))) @@ -343,7 +347,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 @@ -374,7 +378,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( @@ -406,14 +410,15 @@ 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_) + waveform_readers = None 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 } @@ -428,7 +433,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 @@ -573,13 +578,13 @@ 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"}, ) @@ -592,88 +597,95 @@ def _generator(self): 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, ) - # 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 - for tel_id in full_tels_with_trigger: - tel_index, tel_trigger = next(telescope_trigger_reader) - - 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) + event.simulation.shower = next(mc_shower_reader) for kind, readers in dl2_readers.items(): - c = getattr(data.dl2.stereo, kind) + c = getattr(event.dl2, 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) - self._fill_telescope_pointing(data, pointing_interpolator) + self._fill_array_pointing(event) - for tel_id in data.trigger.tel.keys(): - key = f"tel_{tel_id:03d}" + # 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 + for tel_id in full_tels_with_trigger: + tel_index, tel_trigger = next(telescope_trigger_reader) if self.allowed_tels and tel_id not in self.allowed_tels: continue - if key in true_impact_readers: - data.simulation.tel[tel_id].impact = next(true_impact_readers[key]) + key = f"tel_{tel_id:03d}" + + containers = { + "index": tel_index, + "dl0": DL0TelescopeContainer(trigger=tel_trigger), + } - if DataLevel.R1 in self.datalevels: - data.r1.tel[tel_id] = next(waveform_readers[key]) + containers["pointing"] = self._fill_telescope_pointing( + tel_index.obs_id, + tel_id, + time=tel_trigger.time, + tel_pointing_interpolator=pointing_interpolator, + ) + if waveform_readers is not None: + r1 = next(waveform_readers[key]) # TODO: Delete if we drop support for datamodel versions < v6.0.0 - r1_waveform = data.r1.tel[tel_id].waveform - if r1_waveform.ndim == 2: - data.r1.tel[tel_id].waveform = r1_waveform[np.newaxis, ...] + if r1.waveform is not None and r1.waveform.ndim == 2: + r1.waveform = r1.waveform[np.newaxis, ...] + containers["r1"] = r1 - if self.has_simulated_dl1: - simulated = data.simulation.tel[tel_id] + if self.is_simulation: + containers["simulation"] = SimulationTelescopeContainer() + + 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 probably 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], @@ -684,7 +696,8 @@ def _generator(self): ) if self.has_simulated_dl1: simulated_params = next(simulated_param_readers[key]) - simulated.true_parameters = ImageParametersContainer( + sim = containers["simulation"] + sim.true_parameters = ImageParametersContainer( hillas=simulated_params[0], leakage=simulated_params[1], concentration=simulated_params[2], @@ -694,14 +707,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]) @@ -710,7 +724,15 @@ 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 def _fill_array_pointing(self, event): @@ -721,34 +743,33 @@ def _fill_array_pointing(self, event): ob = self.observation_blocks[obs_id] frame = ob.subarray_pointing_frame if frame is CoordinateFrameType.ALTAZ: - event.pointing.array_azimuth = ob.subarray_pointing_lon - event.pointing.array_altitude = ob.subarray_pointing_lat + event.pointing.azimuth = ob.subarray_pointing_lon + event.pointing.altitude = ob.subarray_pointing_lat elif frame is CoordinateFrameType.ICRS: - event.pointing.array_ra = ob.subarray_pointing_lon - event.pointing.array_dec = ob.subarray_pointing_lat + event.pointing.ra = ob.subarray_pointing_lon + event.pointing.dec = ob.subarray_pointing_lat else: raise ValueError(f"Unsupported pointing frame: {frame}") - def _fill_telescope_pointing(self, event, tel_pointing_interpolator=None): + def _fill_telescope_pointing( + self, obs_id, tel_id, time, tel_pointing_interpolator=None + ): """ Fill the telescope pointing information of a given event """ if tel_pointing_interpolator is not None: - for tel_id, trigger in event.trigger.tel.items(): - alt, az = tel_pointing_interpolator(tel_id, trigger.time) - event.pointing.tel[tel_id] = TelescopePointingContainer( - altitude=alt, - azimuth=az, - ) - return + alt, az = tel_pointing_interpolator(tel_id, time) + return TelescopePointingContainer( + altitude=alt, + azimuth=az, + ) - for tel_id in event.trigger.tels_with_trigger: - tel_pointing = self._constant_telescope_pointing.get(tel_id) - if tel_pointing is None: - continue + tel_pointing = self._constant_telescope_pointing.get(tel_id) + if tel_pointing is None: + return TelescopePointingContainer() - current = tel_pointing.loc[event.index.obs_id] - event.pointing.tel[tel_id] = TelescopePointingContainer( - altitude=current["telescope_pointing_altitude"], - azimuth=current["telescope_pointing_azimuth"], - ) + current = tel_pointing.loc[obs_id] + return TelescopePointingContainer( + altitude=current["telescope_pointing_altitude"], + azimuth=current["telescope_pointing_azimuth"], + ) diff --git a/src/ctapipe/io/hdf5merger.py b/src/ctapipe/io/hdf5merger.py index 4161dba0878..28f532212eb 100644 --- a/src/ctapipe/io/hdf5merger.py +++ b/src/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, } @@ -323,13 +327,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: @@ -357,13 +362,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/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py index e4ff831dacb..9eb838fb7ea 100644 --- a/src/ctapipe/io/simteleventsource.py +++ b/src/ctapipe/io/simteleventsource.py @@ -21,32 +21,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 ( @@ -762,25 +768,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: @@ -790,28 +793,29 @@ 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 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] @@ -824,19 +828,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"] @@ -844,14 +846,14 @@ 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) select_gain = self.select_gain is True or ( self.select_gain is None - and trigger.event_type is EventType.SUBARRAY + and array_trigger.event_type is EventType.SUBARRAY ) if select_gain: gain_selector = self.gain_selector @@ -866,24 +868,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"] - dl1_calib = data.calibration.tel[tel_id].dl1 - dl1_calib.time_shift = time_calib + calibration = TelescopeCalibrationContainer( + dl1=DL1CameraCalibrationContainer( + time_shift=time_calib, + ) + ) - 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] @@ -970,7 +991,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"] ): @@ -1002,25 +1029,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/src/ctapipe/io/tableloader.py b/src/ctapipe/io/tableloader.py index 36a954fb7de..f5ec3f58170 100644 --- a/src/ctapipe/io/tableloader.py +++ b/src/ctapipe/io/tableloader.py @@ -21,7 +21,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" @@ -269,7 +272,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 _check_args(self, **kwargs): """Checking args: @@ -317,7 +332,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"]] @@ -415,7 +430,13 @@ def read_subarray_events( simulated = updated_args["simulated"] observation_info = updated_args["observation_info"] - 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) @@ -532,7 +553,7 @@ def _read_telescope_events_for_id( 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, @@ -850,7 +871,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/src/ctapipe/io/tests/test_datawriter.py b/src/ctapipe/io/tests/test_datawriter.py index 72450ea1093..7acdbfe8a27 100644 --- a/src/ctapipe/io/tests/test_datawriter.py +++ b/src/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, @@ -198,13 +195,13 @@ def test_roundtrip(tmpdir: Path): # make sure it is readable by the event source and matches the images 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/src/ctapipe/io/tests/test_event_source.py b/src/ctapipe/io/tests/test_event_source.py index c4df13b41e0..f977fa6068c 100644 --- a/src/ctapipe/io/tests/test_event_source.py +++ b/src/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/src/ctapipe/io/tests/test_hdf5.py b/src/ctapipe/io/tests/test_hdf5.py index 681567068a6..e5372f9598e 100644 --- a/src/ctapipe/io/tests/test_hdf5.py +++ b/src/ctapipe/io/tests/test_hdf5.py @@ -12,9 +12,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 @@ -27,7 +27,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 @@ -93,12 +93,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/src/ctapipe/io/tests/test_hdf5eventsource.py b/src/ctapipe/io/tests/test_hdf5eventsource.py index eef8d2b0fa2..d5aa44d24d5 100644 --- a/src/ctapipe/io/tests/test_hdf5eventsource.py +++ b/src/ctapipe/io/tests/test_hdf5eventsource.py @@ -79,14 +79,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)) @@ -96,8 +96,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): @@ -106,12 +106,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)) @@ -122,14 +122,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)) @@ -137,12 +138,11 @@ 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_pointing_divergent(dl1_divergent_file): @@ -157,34 +157,34 @@ def test_pointing_divergent(dl1_divergent_file): for event, simtel_event in zip_longest(source, simtel_source): assert event.index.event_id == simtel_event.index.event_id assert u.isclose( - event.pointing.array_azimuth, - simtel_event.pointing.array_azimuth, + event.pointing.azimuth, + simtel_event.pointing.azimuth, ) assert u.isclose( - event.pointing.array_altitude, - simtel_event.pointing.array_altitude, + event.pointing.altitude, + simtel_event.pointing.altitude, ) - assert event.pointing.tel.keys() == simtel_event.pointing.tel.keys() - for tel in event.pointing.tel: + assert event.tel.keys() == simtel_event.tel.keys() + for tel_id in event.tel: assert u.isclose( - event.pointing.tel[tel].azimuth, - simtel_event.pointing.tel[tel].azimuth, + event.tel[tel_id].pointing.azimuth, + simtel_event.tel[tel_id].pointing.azimuth, ) assert u.isclose( - event.pointing.tel[tel].altitude, - simtel_event.pointing.tel[tel].altitude, + event.tel[tel_id].pointing.altitude, + simtel_event.tel[tel_id].pointing.altitude, ) 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 @@ -198,8 +198,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 @@ -215,18 +215,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 @@ -235,7 +235,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 ) @@ -244,7 +245,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 ) @@ -258,7 +259,6 @@ def test_interpolate_pointing(dl1_mon_pointing_file): with HDF5EventSource(dl1_mon_pointing_file) as source: for e in source: - assert set(e.pointing.tel.keys()) == set(e.trigger.tels_with_trigger) - for pointing in e.pointing.tel.values(): - assert not np.isnan(pointing.azimuth) - assert not np.isnan(pointing.altitude) + for tel_event in e.tel.values(): + assert not np.isnan(tel_event.pointing.azimuth) + assert not np.isnan(tel_event.pointing.altitude) diff --git a/src/ctapipe/io/tests/test_simteleventsource.py b/src/ctapipe/io/tests/test_simteleventsource.py index c614bab2145..4280a8c96f7 100644 --- a/src/ctapipe/io/tests/test_simteleventsource.py +++ b/src/ctapipe/io/tests/test_simteleventsource.py @@ -140,9 +140,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 @@ -176,16 +176,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) + 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 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 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(): @@ -199,11 +203,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(): @@ -232,7 +233,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 @@ -244,11 +245,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(): @@ -257,8 +264,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(): @@ -390,7 +397,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 @@ -399,10 +406,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, ) @@ -413,10 +420,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, ) @@ -428,7 +435,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, @@ -437,11 +444,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 @@ -451,8 +458,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/src/ctapipe/io/tests/test_table_loader.py b/src/ctapipe/io/tests/test_table_loader.py index 398b58d4277..32984ca5b22 100644 --- a/src/ctapipe/io/tests/test_table_loader.py +++ b/src/ctapipe/io/tests/test_table_loader.py @@ -177,7 +177,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) @@ -262,7 +262,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 @@ -400,6 +400,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(path) as loader: diff --git a/src/ctapipe/io/tests/test_toysource.py b/src/ctapipe/io/tests/test_toysource.py index f8e58c4d439..f029938f003 100644 --- a/src/ctapipe/io/tests/test_toysource.py +++ b/src/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/src/ctapipe/io/toymodel.py b/src/ctapipe/io/toymodel.py index 16270c4eef8..4ac451d9adb 100644 --- a/src/ctapipe/io/toymodel.py +++ b/src/ctapipe/io/toymodel.py @@ -8,11 +8,11 @@ import numpy as np from ..containers import ( - ArrayEventContainer, - DL1CameraContainer, - EventIndexContainer, + DL1TelescopeContainer, ObservationBlockContainer, SchedulingBlockContainer, + SubarrayEventContainer, + SubarrayEventIndexContainer, ) from ..core import TelescopeComponent, traits from .datalevels import DataLevel @@ -97,15 +97,12 @@ def _generator(self): def generate_event(self): from ..image import toymodel - 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(): @@ -149,6 +146,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/src/ctapipe/reco/hillas_intersection.py b/src/ctapipe/reco/hillas_intersection.py index 603293ba883..4caebce00d0 100644 --- a/src/ctapipe/reco/hillas_intersection.py +++ b/src/ctapipe/reco/hillas_intersection.py @@ -114,7 +114,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. @@ -123,20 +123,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/src/ctapipe/reco/hillas_reconstructor.py b/src/ctapipe/reco/hillas_reconstructor.py index f3252993445..458e9dee998 100644 --- a/src/ctapipe/reco/hillas_reconstructor.py +++ b/src/ctapipe/reco/hillas_reconstructor.py @@ -139,7 +139,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. @@ -149,7 +149,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 @@ -170,7 +170,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 ) @@ -194,9 +194,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, @@ -238,10 +236,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, ) @@ -262,7 +259,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/src/ctapipe/reco/impact.py b/src/ctapipe/reco/impact.py index b35c0af9a95..f85c549a3ab 100644 --- a/src/ctapipe/reco/impact.py +++ b/src/ctapipe/reco/impact.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python3 """ Implementation of the ImPACT reconstruction algorithm """ @@ -178,21 +177,21 @@ def __call__(self, event): Parameters ---------- event : container - `ctapipe.containers.ArrayEventContainer` + `ctapipe.containers.SubarrayEventContainer` """ try: hillas_dict = self._create_hillas_dict(event) except (TooFewTelescopesException, InvalidWidthException): - event.dl2.stereo.geometry[self.__class__.__name__] = INVALID_GEOMETRY - event.dl2.stereo.energy[self.__class__.__name__] = INVALID_ENERGY + event.dl2.geometry[self.__class__.__name__] = INVALID_GEOMETRY + event.dl2.energy[self.__class__.__name__] = INVALID_ENERGY 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(), ) @@ -202,10 +201,11 @@ def __call__(self, event): # Finally get the telescope images and and the selection masks mask_dict, image_dict, time_dict = {}, {}, {} for tel_id in hillas_dict.keys(): - image = event.dl1.tel[tel_id].image + dl1 = event.tel[tel_id].dl1 + image = dl1.image image_dict[tel_id] = image - time_dict[tel_id] = event.dl1.tel[tel_id].peak_time - mask = event.dl1.tel[tel_id].image_mask + time_dict[tel_id] = dl1.peak_time + mask = dl1.image_mask # Dilate the images around the original cleaning to help the fit for _ in range(3): @@ -216,15 +216,14 @@ def __call__(self, event): # Next, we look for geometry and energy seeds from previously applied reconstructors. # Both need to be present at elast once for ImPACT to run. - reco_geom_pred = event.dl2.stereo.geometry - + reco_geom_pred = event.dl2.geometry valid_geometry_seed = False for geom_pred in reco_geom_pred.values(): if geom_pred.is_valid: valid_geometry_seed = True break - reco_energy_pred = event.dl2.stereo.energy + reco_energy_pred = event.dl2.energy valid_energy_seed = False for E_pred in reco_energy_pred.values(): @@ -233,8 +232,8 @@ def __call__(self, event): break if valid_geometry_seed is False or valid_energy_seed is False: - event.dl2.stereo.geometry[self.__class__.__name__] = INVALID_GEOMETRY - event.dl2.stereo.energy[self.__class__.__name__] = INVALID_ENERGY + event.dl2.geometry[self.__class__.__name__] = INVALID_GEOMETRY + event.dl2.energy[self.__class__.__name__] = INVALID_ENERGY self._store_impact_parameter(event) return @@ -249,8 +248,8 @@ def __call__(self, event): mask_dict=mask_dict, time_dict=time_dict, ) - event.dl2.stereo.geometry[self.__class__.__name__] = shower_result - event.dl2.stereo.energy[self.__class__.__name__] = energy_result + event.dl2.geometry[self.__class__.__name__] = shower_result + event.dl2.energy[self.__class__.__name__] = energy_result self._store_impact_parameter(event) diff --git a/src/ctapipe/reco/preprocessing.py b/src/ctapipe/reco/preprocessing.py index eb887ea3906..e707c594461 100644 --- a/src/ctapipe/reco/preprocessing.py +++ b/src/ctapipe/reco/preprocessing.py @@ -16,7 +16,7 @@ from ctapipe.coordinates import MissingFrameAttributeWarning, TelescopeFrame -from ..containers import ArrayEventContainer +from ..containers import SubarrayEventContainer LOG = logging.getLogger(__name__) @@ -68,13 +68,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 @@ -87,10 +87,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, @@ -98,7 +99,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, @@ -107,7 +108,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/src/ctapipe/reco/reconstructor.py b/src/ctapipe/reco/reconstructor.py index 24b708e9d4a..e67e0cda0c3 100644 --- a/src/ctapipe/reco/reconstructor.py +++ b/src/ctapipe/reco/reconstructor.py @@ -7,7 +7,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 Integer, List @@ -90,7 +90,7 @@ def __init__(self, subarray, atmosphere_profile=None, **kwargs): self.atmosphere_profile = atmosphere_profile @abstractmethod - def __call__(self, event: ArrayEventContainer): + def __call__(self, event: SubarrayEventContainer): """ Perform stereo reconstruction on event. @@ -99,7 +99,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. @@ -161,9 +161,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: @@ -186,16 +186,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( @@ -207,12 +208,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/src/ctapipe/reco/shower_processor.py b/src/ctapipe/reco/shower_processor.py index 1512080a0ed..416cf70ae44 100644 --- a/src/ctapipe/reco/shower_processor.py +++ b/src/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 @@ -83,13 +83,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/src/ctapipe/reco/sklearn.py b/src/ctapipe/reco/sklearn.py index 6b6c9cc133f..1e8ecec92e5 100644 --- a/src/ctapipe/reco/sklearn.py +++ b/src/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 ( @@ -169,15 +169,14 @@ def __init__( self.prefix = self.model_cls @abstractmethod - def __call__(self, event: ArrayEventContainer) -> None: - """ - Event-wise prediction for the EventSource-Loop. + 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 @@ -388,8 +387,8 @@ class EnergyRegressor(SKLearnRegressionReconstructor): target = "true_energy" property = ReconstructionProperty.ENERGY - 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) @@ -412,7 +411,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) @@ -453,8 +452,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] @@ -475,7 +474,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) @@ -707,18 +706,17 @@ def _predict(self, key, table): return prediction, score, valid - def __call__(self, event: ArrayEventContainer) -> None: - """ - Event-wise prediction for the EventSource-Loop. + def __call__(self, event: SubarrayEventContainer) -> None: + """Event-wise prediction for the EventSource-Loop. Fills the event.dl2.tel[tel_id].disp[prefix] container and event.dl2.tel[tel_id].geometry[prefix] container. 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) @@ -734,8 +732,7 @@ def __call__(self, event: ArrayEventContainer) -> None: parameter=disp[0], sign_score=sign_score[0], ) - - 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) @@ -744,8 +741,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()) @@ -774,8 +771,8 @@ def __call__(self, event: ArrayEventContainer) -> None: disp_container.prefix = f"{self.prefix}_tel" 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/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index 7c48a971fe5..645e762d507 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/src/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 @@ -78,7 +78,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. """ @@ -155,15 +155,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: @@ -185,7 +184,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, @@ -198,12 +197,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: @@ -216,7 +217,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 = [] @@ -224,13 +225,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) @@ -258,7 +261,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), @@ -267,7 +270,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/src/ctapipe/reco/tests/test_HillasReconstructor.py b/src/ctapipe/reco/tests/test_HillasReconstructor.py index d0f1126d47c..54283efc9e9 100644 --- a/src/ctapipe/reco/tests/test_HillasReconstructor.py +++ b/src/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/src/ctapipe/reco/tests/test_ImPACT.py b/src/ctapipe/reco/tests/test_ImPACT.py index 35e0fc7c20b..94e45c0390c 100644 --- a/src/ctapipe/reco/tests/test_ImPACT.py +++ b/src/ctapipe/reco/tests/test_ImPACT.py @@ -213,10 +213,10 @@ def test_selected_subarray( energy_test.energy = 500 * u.GeV energy_test.is_valid = True - event.dl2.stereo.geometry["test"] = shower_test - event.dl2.stereo.energy["test_energy"] = energy_test + event.dl2.geometry["test"] = shower_test + event.dl2.energy["test_energy"] = energy_test reconstructor = ImPACTReconstructor(subarray, table_profile) reconstructor.root_dir = str(tmp_path) reconstructor(event) - assert event.dl2.stereo.geometry["ImPACTReconstructor"].is_valid - assert event.dl2.stereo.energy["ImPACTReconstructor"].is_valid + assert event.dl2.geometry["ImPACTReconstructor"].is_valid + assert event.dl2.energy["ImPACTReconstructor"].is_valid diff --git a/src/ctapipe/reco/tests/test_hillas_intersection.py b/src/ctapipe/reco/tests/test_hillas_intersection.py index bfef423cc78..8a15efcc52a 100644 --- a/src/ctapipe/reco/tests/test_hillas_intersection.py +++ b/src/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/src/ctapipe/reco/tests/test_preprocessing.py b/src/ctapipe/reco/tests/test_preprocessing.py index 861fb33404e..ef8b4c3651e 100644 --- a/src/ctapipe/reco/tests/test_preprocessing.py +++ b/src/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/src/ctapipe/reco/tests/test_reconstruction_methods.py b/src/ctapipe/reco/tests/test_reconstruction_methods.py index 4088ed94f0e..8ad133151d6 100644 --- a/src/ctapipe/reco/tests/test_reconstruction_methods.py +++ b/src/ctapipe/reco/tests/test_reconstruction_methods.py @@ -47,7 +47,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/src/ctapipe/reco/tests/test_shower_processor.py b/src/ctapipe/reco/tests/test_shower_processor.py index 1560bf837e4..6145103e87c 100644 --- a/src/ctapipe/reco/tests/test_shower_processor.py +++ b/src/ctapipe/reco/tests/test_shower_processor.py @@ -67,11 +67,8 @@ def test_shower_processor_geometry( process_shower(example_event_copy) - print(reconstructor_types) for reco_type in reconstructor_types: - DL2a = example_event_copy.dl2.stereo.geometry[reco_type] - print(reco_type) - print(DL2a) + DL2a = example_event_copy.dl2.geometry[reco_type] assert isfinite(DL2a.alt) assert isfinite(DL2a.az) assert isfinite(DL2a.core_x) @@ -97,7 +94,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/src/ctapipe/reco/tests/test_stereo_combination.py b/src/ctapipe/reco/tests/test_stereo_combination.py index 7187262af27..3b05cdbc333 100644 --- a/src/ctapipe/reco/tests/test_stereo_combination.py +++ b/src/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/src/ctapipe/tools/apply_models.py b/src/ctapipe/tools/apply_models.py index 6ec237006e0..2a107fd1b48 100644 --- a/src/ctapipe/tools/apply_models.py +++ b/src/ctapipe/tools/apply_models.py @@ -234,9 +234,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/src/ctapipe/tools/display_dl1.py b/src/ctapipe/tools/display_dl1.py index 7a27305b60b..30baab8dc0a 100644 --- a/src/ctapipe/tools/display_dl1.py +++ b/src/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/src/ctapipe/tools/tests/test_apply_models.py b/src/ctapipe/tools/tests/test_apply_models.py index 6958a3b5812..1bb33629ddc 100644 --- a/src/ctapipe/tools/tests/test_apply_models.py +++ b/src/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}" @@ -73,7 +76,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) @@ -158,6 +161,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 = ( @@ -170,6 +174,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/src/ctapipe/tools/tests/test_process.py b/src/ctapipe/tools/tests/test_process.py index 769574c71f4..c77c0428c6c 100644 --- a/src/ctapipe/tools/tests/test_process.py +++ b/src/ctapipe/tools/tests/test_process.py @@ -87,9 +87,11 @@ 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 @@ -445,7 +447,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/src/ctapipe/tools/tests/test_process_ml.py b/src/ctapipe/tools/tests/test_process_ml.py index 6a25f503aa5..4d5aaad405c 100644 --- a/src/ctapipe/tools/tests/test_process_ml.py +++ b/src/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 @@ -169,7 +169,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/src/ctapipe/utils/event_type_filter.py b/src/ctapipe/utils/event_type_filter.py index 6369608a4d4..81db188a292 100644 --- a/src/ctapipe/utils/event_type_filter.py +++ b/src/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/src/ctapipe/utils/tests/test_event_filter.py b/src/ctapipe/utils/tests/test_event_filter.py index 6a0fc5990e8..a70e93551d2 100644 --- a/src/ctapipe/utils/tests/test_event_filter.py +++ b/src/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/src/ctapipe/visualization/mpl_array.py b/src/ctapipe/visualization/mpl_array.py index 3bc0064f1b6..c1ed5d4fba7 100644 --- a/src/ctapipe/visualization/mpl_array.py +++ b/src/ctapipe/visualization/mpl_array.py @@ -278,7 +278,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/src/ctapipe/visualization/tests/test_bokeh.py b/src/ctapipe/visualization/tests/test_bokeh.py index a6ee6906865..907c9a0759c 100644 --- a/src/ctapipe/visualization/tests/test_bokeh.py +++ b/src/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) @@ -66,8 +66,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/src/ctapipe/visualization/tests/test_mpl.py b/src/ctapipe/visualization/tests/test_mpl.py index 70776678d11..8c05f2b110c 100644 --- a/src/ctapipe/visualization/tests/test_mpl.py +++ b/src/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/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()