Merge pull request #2265 from gschwefer/hillasintersection_fix

Fix HillasIntersection error for badly reconstructed events

ctapipe/reco/hillas_intersection.py

@@ -43,6 +43,29 @@
     prefix="HillasIntersection",
 )
 
+FOV_ANGULAR_DISTANCE_LIMIT_RAD = (45 * u.deg).to_value(u.rad)
+
+
+def _far_outside_fov(fov_lat, fov_lon):
+    """
+    Check if a given latitude or longiude in the FoV is further away from
+    the FoV center than `FOV_ANGULAR_DISTANCE_LIMIT`
+
+    Parameters
+    ----------
+    fov_lat : u.Quantity[angle]
+        Latitude in TelescopeFrame or NominalFrame
+    fov_lon : u.Quantity[angle]
+        Longitude in TelescopeFrame or NominalFrame
+
+    Returns
+    -------
+    bool
+    """
+    lat_outside_fov = np.abs(fov_lat) > FOV_ANGULAR_DISTANCE_LIMIT_RAD
+    lon_outside_fov = np.abs(fov_lon) > FOV_ANGULAR_DISTANCE_LIMIT_RAD
+    return lat_outside_fov or lon_outside_fov
+
 
 class HillasIntersection(HillasGeometryReconstructor):
     """
@@ -218,6 +241,12 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
         src_fov_lon, src_fov_lat, err_fov_lon, err_fov_lat = self.reconstruct_nominal(
             hillas_dict_mod
         )
+
+        # Catch events reconstructed at great angular distance from camera center
+        # and retrun INVALID container to prevent SkyCoord error below.
+        if _far_outside_fov(src_fov_lat, src_fov_lon):
+            return INVALID
+
         core_x, core_y, core_err_x, core_err_y = self.reconstruct_tilted(
             hillas_dict_mod, tel_x, tel_y
         )

ctapipe/reco/tests/test_hillas_intersection.py

@@ -1,10 +1,11 @@
 import astropy.units as u
 import numpy as np
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import AltAz, EarthLocation, SkyCoord
 from numpy.testing import assert_allclose
 
 from ctapipe.containers import HillasParametersContainer
-from ctapipe.coordinates import AltAz, CameraFrame, NominalFrame
+from ctapipe.coordinates import CameraFrame, NominalFrame
+from ctapipe.instrument import SubarrayDescription
 from ctapipe.reco.hillas_intersection import HillasIntersection
 
 
@@ -216,6 +217,52 @@ def test_intersection_nominal_reconstruction(example_subarray):
     )
 
 
+def test_badly_reconstructed_event(prod5_mst_flashcam):
+    """
+    Test that events reconstructed at large angular distance
+    from FoV center return INVALID. Event and array loosely follow an
+    actual simulation event.
+    """
+    tel_pos = {1: np.array([150, 75, 0]) * u.m, 2: np.array([150, -240, 0]) * u.m}
+    tel_desc = {1: prod5_mst_flashcam, 2: prod5_mst_flashcam}
+    reference_location = EarthLocation(
+        lon=-70.32 * u.deg,
+        lat=-24.68 * u.deg,
+        height=2147 * u.m,
+    )
+    two_tel_subarray = SubarrayDescription(
+        "two_tel_subarray",
+        tel_positions=tel_pos,
+        tel_descriptions=tel_desc,
+        reference_location=reference_location,
+    )
+    hill_inter = HillasIntersection(two_tel_subarray)
+
+    hillas_tel_1 = HillasParametersContainer(
+        fov_lon=0.6 * u.deg,
+        fov_lat=-0.3 * u.deg,
+        intensity=95,
+        psi=85.3 * u.deg,
+        length=0.09 * u.deg,
+        width=0.05 * u.deg,
+    )
+
+    hillas_tel_2 = HillasParametersContainer(
+        fov_lon=-0.1 * u.deg,
+        fov_lat=-0.4 * u.deg,
+        intensity=119,
+        psi=85.0 * u.deg,
+        length=0.15 * u.deg,
+        width=0.04 * u.deg,
+    )
+
+    hillas_dir = {1: hillas_tel_1, 2: hillas_tel_2}
+    pointing = AltAz(alt=70 * u.deg, az=0 * u.deg)
+    reco_event = hill_inter._predict(hillas_dict=hillas_dir, array_pointing=pointing)
+
+    assert not reco_event.is_valid
+
+
 def test_reconstruction_works(subarray_and_event_gamma_off_axis_500_gev):
     subarray, event = subarray_and_event_gamma_off_axis_500_gev
     reconstructor = HillasIntersection(subarray)

docs/changes/2265.bugfix.rst

@@ -0,0 +1,2 @@
+The ```HillasIntersection``` method used to fail when individual events were reconstructed to originate from a FoV offset of more than 90 degrees.
+This is now fixed by returning an INVALID container for a reconstructed offset of larger than 45 degrees.