Merge pull request #2082 from cta-observatory/fix_invalid

Fix invalid case in HillasReconstructor and HillasIntersection

Author: maxnoe

ctapipe/reco/hillas_intersection.py

@@ -96,14 +96,18 @@ def __call__(self, event):
 
         Parameters
         ----------
-        event: `ctapipe.containers.ArrayEventContainer`
-            The event, needs to have dl1 parameters
+        event : `~ctapipe.containers.ArrayEventContainer`
+            The event, needs to have dl1 parameters.
+            Will be filled with the corresponding dl2 containers,
+            reconstructed stereo geometry and telescope-wise impact position.
         """
 
         try:
             hillas_dict = self._create_hillas_dict(event)
         except (TooFewTelescopesException, InvalidWidthException):
-            return INVALID
+            event.dl2.stereo.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(

ctapipe/reco/hillas_reconstructor.py

@@ -133,23 +133,23 @@ def __init__(self, subarray, **kwargs):
 
     def __call__(self, event):
         """
-        Perform the full shower geometry reconstruction on the input event.
+        Perform the full shower geometry reconstruction.
 
         Parameters
         ----------
-        event: `ctapipe.containers.ArrayEventContainer`
-            The event, needs to have dl1 parameters
-
-        Returns
-        -------
-        ReconstructedGeometryContainer
+        event : `~ctapipe.containers.ArrayEventContainer`
+            The event, needs to have dl1 parameters.
+            Will be filled with the corresponding dl2 containers,
+            reconstructed stereo geometry and telescope-wise impact position.
         """
         warnings.filterwarnings(action="ignore", category=MissingFrameAttributeWarning)
 
         try:
             hillas_dict = self._create_hillas_dict(event)
         except (TooFewTelescopesException, InvalidWidthException):
-            return INVALID
+            event.dl2.stereo.geometry[self.__class__.__name__] = INVALID
+            self._store_impact_parameter(event)
+            return
 
         # Here we perform some basic quality checks BEFORE applying reconstruction
         # This should be substituted by a DL1 QualityQuery specific to this

ctapipe/reco/reco_algorithms.py

@@ -1,5 +1,6 @@
 from abc import abstractmethod
 
+import astropy.units as u
 import numpy as np
 from astropy.coordinates import AltAz, SkyCoord
 
@@ -47,18 +48,18 @@ def __init__(self, subarray, **kwargs):
 
     @abstractmethod
     def __call__(self, event: ArrayEventContainer):
-        """overwrite this method with your favourite direction reconstruction
-        algorithm
+        """
+        Perform stereo reconstruction on event.
+
+        This method must fill the result of the reconstruction into the
+        dl2 structure of the event.
 
         Parameters
         ----------
-        tels_dict : dict
-            general dictionary containing all triggered telescopes data
-        Returns
-        -------
-        None
-        Container are directly filled into the event structure
-
+        event : `ctapipe.containers.ArrayEventContainer`
+            The event, needs to have dl1 parameters.
+            Will be filled with the corresponding dl2 containers,
+            reconstructed stereo geometry and telescope-wise impact position.
         """
 
     def _create_hillas_dict(self, event):
@@ -97,10 +98,17 @@ def _get_telescope_pointings(event):
 
     def _store_impact_parameter(self, event):
         """Compute and store the impact parameter for each reconstruction."""
-        impact_distances = shower_impact_distance(
-            shower_geom=event.dl2.stereo.geometry[self.__class__.__name__],
-            subarray=self.subarray,
-        )
+        geometry = event.dl2.stereo.geometry[self.__class__.__name__]
+
+        if geometry.is_valid:
+            impact_distances = shower_impact_distance(
+                shower_geom=geometry,
+                subarray=self.subarray,
+            )
+        else:
+            n_tels = len(self.subarray)
+            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:

ctapipe/reco/tests/test_HillasReconstructor.py

@@ -90,6 +90,8 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev):
     original_event = deepcopy(event)
 
     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 result.is_valid
 
@@ -101,6 +103,8 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev):
         invalid_event.dl1.tel[tel_id].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 not result.is_valid
 
@@ -109,13 +113,17 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev):
     invalid_event = deepcopy(original_event)
     invalid_event.dl1.tel[tel_id].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 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
     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 not result.is_valid
 
@@ -140,6 +148,8 @@ def test_reconstruction_against_simulation_telescope_frame(
     calib(event)
     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"]
 
     # get the reconstructed coordinates in the sky

ctapipe/reco/tests/test_reconstruction_methods.py

@@ -42,6 +42,8 @@ def test_reconstructors(reconstructors):
             reconstructor(event)
 
             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)