Implement a shower processor (#1675)

* add shower processor

* add test

* should be ok...

* make the test check right things

* add class to init

* remove geometry as an option (alwayd done if processor is called)

* attempt to fix the test

* Fix basic test bug and expand class API

* Expand test to all possible usage cases

* Improve quality criteria

* remove is_simulation

* change _estimate_classification to __classify_particle_type

* Update docstrings

* Improve hillas dict creation and shower criteria, general fixes

* Update test

* Make reconstruct_energy and classify traitlets and update documentation.

* Log shower criteria without f-string

* Shower criteria table formatted string representation if log-level DEBUG

* Fix imports

* Missing list of public objects in shower_processor sub-module

* Fix docstring of ctapipe.reco.shower_processor.ShowerProcessor.__call__

* Add shower_processor submodule to automodapi

* fix docs

* Fix docs warnings

* Remove parts not implemented and make single functions public API

* Fix remaining bugs and update to new DL2 data model

* Fix docstring of ctapipe.reco.shower_processor.ShowerProcessor.reconstruct_shower

* Import reconstructor from correct submodule

* Fix docstring of ShowerProcessor

* remove stale code

* log shower geometry correctly

* remove not implemented options from test

* rename reconstruct_shower to reconstruct_geometry

Co-authored-by: Maximilian Nöthe <maximilian.noethe@tu-dortmund.de>

ctapipe/core/qualityquery.py

@@ -113,6 +113,10 @@ def _repr_html_(self):
         """display nicely in Jupyter notebooks"""
         return self.to_table()._repr_html_()
 
+    def __str__(self):
+        """Print a formatted string representation of the entire table."""
+        return self.to_table().pprint_all(show_unit=True, show_dtype=True)
+
     def __call__(self, value) -> np.ndarray:
         """
         Test that value passes all cuts

ctapipe/reco/__init__.py

@@ -1,13 +1,13 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-from .reco_algorithms import Reconstructor
 from .hillas_reconstructor import HillasReconstructor
-from .impact import ImPACTReconstructor
+from .reco_algorithms import Reconstructor
+from .shower_processor import ShowerProcessor
 from .hillas_intersection import HillasIntersection
+from .impact import ImPACTReconstructor
 
 
-__all__ = [
-    "Reconstructor",
-    "HillasReconstructor",
-    "HillasIntersection",
-    "ImPACTReconstructor",
-]
+__all__ = ["Reconstructor",
+           "ShowerProcessor",
+           "HillasReconstructor",
+           "ImPACTReconstructor",
+           "HillasIntersection"]

ctapipe/reco/shower_processor.py

@@ -0,0 +1,166 @@
+"""
+High level processing of showers.
+
+This processor will be able to process a shower/event in 3 steps:
+- shower geometry
+- estimation of energy (optional, currently unavailable)
+- estimation of classification (optional, currently unavailable)
+
+"""
+from astropy.coordinates import SkyCoord, AltAz
+
+from ctapipe.core import Component, QualityQuery
+from ctapipe.core.traits import List
+from ctapipe.containers import ArrayEventContainer, ReconstructedGeometryContainer
+from ctapipe.instrument import SubarrayDescription
+from ctapipe.reco import HillasReconstructor
+
+
+DEFAULT_SHOWER_PARAMETERS = ReconstructedGeometryContainer()
+
+
+class ShowerQualityQuery(QualityQuery):
+    """Configuring shower-wise data checks."""
+
+    quality_criteria = List(
+        default_value=[
+            ("> 50 phe", "lambda p: p.hillas.intensity > 50"),
+            ("Positive width", "lambda p: p.hillas.width.value > 0"),
+            ("> 3 pixels", "lambda p: p.morphology.num_pixels > 3"),
+        ],
+        help=QualityQuery.quality_criteria.help,
+    ).tag(config=True)
+
+
+class ShowerProcessor(Component):
+    """
+    Needs DL1_PARAMETERS as input.
+    Should be run after ImageProcessor which produces such information.
+
+    For the moment it only supports the reconstruction of the shower geometry
+    using ctapipe.reco.HillasReconstructor.
+
+    It is planned to support also energy reconstruction and particle type
+    classification.
+    """
+
+    def __init__(
+        self,
+        subarray: SubarrayDescription,
+        config=None,
+        parent=None,
+        **kwargs
+    ):
+        """
+        Parameters
+        ----------
+        subarray: SubarrayDescription
+            Description of the subarray. Provides information about the
+            camera which are useful in calibration. Also required for
+            configuring the TelescopeParameter traitlets.
+        config: traitlets.loader.Config
+            Configuration specified by config file or cmdline arguments.
+            Used to set traitlet values.
+            This is mutually exclusive with passing a ``parent``.
+        parent: ctapipe.core.Component or ctapipe.core.Tool
+            Parent of this component in the configuration hierarchy,
+            this is mutually exclusive with passing ``config``
+        """
+
+        super().__init__(config=config, parent=parent, **kwargs)
+        self.subarray = subarray
+        self.check_shower = ShowerQualityQuery(parent=self)
+        self.reconstructor = HillasReconstructor()
+
+    def reconstruct_geometry(
+        self,
+        event,
+        default=DEFAULT_SHOWER_PARAMETERS,
+    ) -> ReconstructedGeometryContainer:
+        """Perform shower reconstruction.
+
+        Parameters
+        ----------
+        event : container
+            A ``ctapipe`` event container
+        default: container
+            The default 'ReconstructedGeometryContainer' which is
+            filled with NaNs.
+        Returns
+        -------
+        ReconstructedGeometryContainer:
+            direction in the sky with uncertainty,
+            core position on the ground with uncertainty,
+            h_max with uncertainty,
+            is_valid boolean for successfull reconstruction,
+            average intensity of the intensities used for reconstruction,
+            measure of algorithm success (if fit),
+            list of tel_ids used if stereo, or None if Mono
+        """
+
+        # Select only images which pass the shower quality criteria
+        hillas_dict = {
+            tel_id: dl1.parameters.hillas
+            for tel_id, dl1 in event.dl1.tel.items()
+            if all(self.check_shower(dl1.parameters))
+        }
+        self.log.debug(
+            "shower_criteria:\n %s", self.check_shower
+        )
+
+        # Reconstruct the shower only if all shower criteria are met
+        if len(hillas_dict) > 2:
+
+            array_pointing = SkyCoord(
+                az=event.pointing.array_azimuth,
+                alt=event.pointing.array_altitude,
+                frame=AltAz(),
+            )
+
+            telescopes_pointings = {
+                tel_id: SkyCoord(
+                    alt=event.pointing.tel[tel_id].altitude,
+                    az=event.pointing.tel[tel_id].azimuth,
+                    frame=AltAz(),
+                )
+                for tel_id in hillas_dict
+            }
+
+            result = self.reconstructor.predict(hillas_dict,
+                                                self.subarray,
+                                                array_pointing,
+                                                telescopes_pointings)
+
+            return result
+
+        else:
+            self.log.debug(
+                """Less than 2 images passed the quality cuts.
+                Returning default ReconstructedGeometryContainer container"""
+            )
+            return default
+
+    def process_shower_geometry(self, event: ArrayEventContainer):
+        """Record the reconstructed shower geometry into the ArrayEventContainer."""
+
+        shower_geometry = self.reconstruct_geometry(event)
+
+        self.log.debug("shower geometry:\n %s", shower_geometry)
+
+        event.dl2.stereo.geometry["HillasReconstructor"] = shower_geometry
+
+    def __call__(self, event: ArrayEventContainer):
+        """
+        Perform the full shower geometry reconstruction on the input event.
+
+        Afterwards, optionally perform energy estimation and/or particle
+        classification (currently these two operations are not yet supported).
+
+        Parameters
+        ----------
+        event : ctapipe.containers.ArrayEventContainer
+            Top-level container for all event information.
+        """
+
+        # This is always done when calling the ShowerProcessor
+        self.process_shower_geometry(event)

ctapipe/reco/tests/test_shower_processor.py

@@ -0,0 +1,66 @@
+"""
+Tests for ShowerProcessor functionalities.
+"""
+import pytest
+from numpy import isfinite
+
+from traitlets.config.loader import Config
+
+from ctapipe.calib import CameraCalibrator
+from ctapipe.image import ImageProcessor
+from ctapipe.reco import ShowerProcessor
+
+
+def test_shower_processor_geometry(example_event, example_subarray):
+    """Ensure we get shower geometry when we input an event with parametrized images."""
+
+    calibrate = CameraCalibrator(subarray=example_subarray)
+
+    config = Config()
+
+    process_images = ImageProcessor(
+        subarray=example_subarray,
+        is_simulation=True,
+        image_cleaner_type="MARSImageCleaner",
+    )
+
+    process_shower = ShowerProcessor(
+        subarray=example_subarray
+    )
+
+    calibrate(example_event)
+    process_images(example_event)
+
+    process_shower(example_event)
+    print(process_shower.check_shower.to_table())
+
+    DL2a = example_event.dl2.stereo.geometry["HillasReconstructor"]
+    print(DL2a)
+    assert isfinite(DL2a.alt)
+    assert isfinite(DL2a.az)
+    assert isfinite(DL2a.core_x)
+    assert isfinite(DL2a.core_x)
+    assert isfinite(DL2a.core_y)
+    assert DL2a.is_valid
+    assert isfinite(DL2a.average_intensity)
+
+    # Increase some quality cuts and check that we get defaults
+    config.ShowerQualityQuery.quality_criteria = [("> 500 phes", "lambda p: p.hillas.intensity > 500")]
+
+    process_shower = ShowerProcessor(
+        config=config,
+        subarray=example_subarray
+    )
+
+    process_shower(example_event)
+    print(process_shower.check_shower.to_table())
+
+    DL2a = example_event.dl2.stereo.geometry["HillasReconstructor"]
+    print(DL2a)
+    assert not isfinite(DL2a.alt)
+    assert not isfinite(DL2a.az)
+    assert not isfinite(DL2a.core_x)
+    assert not isfinite(DL2a.core_x)
+    assert not isfinite(DL2a.core_y)
+    assert not DL2a.is_valid
+    assert not isfinite(DL2a.average_intensity)

docs/ctapipe_api/reco/index.rst

@@ -43,9 +43,12 @@ Reference/API
 =============
 
 .. automodapi:: ctapipe.reco
-
+    :no-inheritance-diagram:
 .. automodapi:: ctapipe.reco.hillas_intersection
+    :no-inheritance-diagram:
 .. automodapi:: ctapipe.reco.hillas_reconstructor
+    :no-inheritance-diagram:
 .. automodapi:: ctapipe.reco.impact
+    :no-inheritance-diagram:
 .. automodapi:: ctapipe.reco.reco_algorithms
-
+    :no-inheritance-diagram: