StereoDispCombiner: Add Transformation to Nominalframe so divergent pointing is supported

src/ctapipe/reco/stereo_combination.py

@@ -2,11 +2,17 @@
 
 import astropy.units as u
 import numpy as np
-from astropy.coordinates import AltAz, CartesianRepresentation, SphericalRepresentation
+from astropy.coordinates import (
+    AltAz,
+    CartesianRepresentation,
+    SkyCoord,
+    SphericalRepresentation,
+)
 from astropy.table import Table
 from traitlets import UseEnum
 
 from ctapipe.containers import ImageParametersContainer
+from ctapipe.coordinates import NominalFrame, TelescopeFrame
 from ctapipe.core import Component, Container
 from ctapipe.core.traits import (
     Bool,
@@ -734,8 +740,8 @@ def __call__(self, event: ArrayEventContainer) -> None:
 
         if valid:
             alt, az = telescope_to_horizontal(
-                lon=stereo_fov_lon * u.deg,
-                lat=stereo_fov_lat * u.deg,
+                lon=u.Quantity(stereo_fov_lon, u.deg, copy=COPY_IF_NEEDED),
+                lat=u.Quantity(stereo_fov_lat, u.deg, copy=COPY_IF_NEEDED),
                 pointing_alt=event.monitoring.pointing.array_altitude,
                 pointing_az=event.monitoring.pointing.array_azimuth,
             )
@@ -757,6 +763,10 @@ def _collect_valid_tel_data(self, event: ArrayEventContainer):
 
         signs = np.array([-1, 1])
 
+        # Gather the subarray pointing information for the transformation to the nominal frame
+        subarray_pointing_alt = event.monitoring.pointing.array_altitude
+        subarray_pointing_az = event.monitoring.pointing.array_azimuth
+
         for tel_id, dl2 in event.dl2.tel.items():
             if not dl2.geometry[self.prefix].is_valid:
                 continue
@@ -776,12 +786,34 @@ def _collect_valid_tel_data(self, event: ArrayEventContainer):
 
             hillas_fov_lon = dl1.hillas.fov_lon.to_value(u.deg)
             hillas_fov_lat = dl1.hillas.fov_lat.to_value(u.deg)
-            hillas_psi = dl1.hillas.psi
+            hillas_psi = dl1.hillas.psi.to_value(u.rad)
             disp = disp_reco.parameter.to_value(u.deg)
 
             fov_lons = hillas_fov_lon + signs * disp * np.cos(hillas_psi)
             fov_lats = hillas_fov_lat + signs * disp * np.sin(hillas_psi)
 
+            # Convert to quantity to ensure the helper functions can handle the inputs correctly
+            fov_lons = u.Quantity(fov_lons, u.deg, copy=COPY_IF_NEEDED)
+            fov_lats = u.Quantity(fov_lats, u.deg, copy=COPY_IF_NEEDED)
+
+            # Gather the telescope pointing information for the transformation to the nominal frame
+            tel_pointing_alt = event.monitoring.tel[tel_id].pointing.altitude
+            tel_pointing_az = event.monitoring.tel[tel_id].pointing.azimuth
+
+            fov_lons, fov_lats = self._transform_to_nominal(
+                tel_pointing_alt,
+                tel_pointing_az,
+                subarray_pointing_alt,
+                subarray_pointing_az,
+                fov_lons,
+                fov_lats,
+            )
+            hillas_psi = u.Quantity(
+                np.arctan2(np.diff(fov_lats), np.diff(fov_lons))[0],
+                u.rad,
+                copy=COPY_IF_NEEDED,
+            )
+
             fov_lon_values.append(fov_lons)
             fov_lat_values.append(fov_lats)
             weights.append(self._calculate_weights(dl1) if dl1 else 1)
@@ -850,15 +882,51 @@ def predict_table(self, mono_predictions: Table) -> Table:
         valid = mono_predictions[f"{prefix_tel}_is_valid"].copy()
         self._require_disp_column(mono_predictions, prefix_tel)
 
+        # Returns values as to_value(u.deg)
+        fov_lon_values, fov_lat_values = calc_fov_lon_lat(mono_predictions, prefix_tel)
+
+        # Convert to radians for the angular difference calculation
+        fov_lon_values = u.Quantity(fov_lon_values, u.deg, copy=COPY_IF_NEEDED)
+        fov_lat_values = u.Quantity(fov_lat_values, u.deg, copy=COPY_IF_NEEDED)
+        tel_pointing_alt = mono_predictions["telescope_pointing_altitude"].quantity[
+            :, None
+        ]
+        tel_pointing_az = mono_predictions["telescope_pointing_azimuth"].quantity[
+            :, None
+        ]
+        subarray_pointing_alt = mono_predictions["subarray_pointing_lat"].quantity[
+            :, None
+        ]
+        subarray_pointing_az = mono_predictions["subarray_pointing_lon"].quantity[
+            :, None
+        ]
+
+        fov_lon_values, fov_lat_values = self._transform_to_nominal(
+            tel_pointing_alt,
+            tel_pointing_az,
+            subarray_pointing_alt,
+            subarray_pointing_az,
+            fov_lon_values,
+            fov_lat_values,
+        )
+        hillas_psis = u.Quantity(
+            np.arctan2(
+                np.diff(fov_lat_values, axis=1),
+                np.diff(fov_lon_values, axis=1),
+            )[:, 0],
+            u.rad,
+            copy=COPY_IF_NEEDED,
+        )
+
         obs_ids, event_ids, _, tel_to_array_indices = get_subarray_index(
             mono_predictions
         )
         n_array_events = len(obs_ids)
 
         valid = self._apply_min_ang_diff_cut(
-            mono_predictions=mono_predictions,
             valid=valid,
             tel_to_array_indices=tel_to_array_indices,
+            hillas_psis=hillas_psis,
         )
         valid = self._apply_n_best_tels_cut(
             mono_predictions=mono_predictions,
@@ -873,6 +941,8 @@ def predict_table(self, mono_predictions: Table) -> Table:
             tel_to_array_indices=tel_to_array_indices,
             n_array_events=n_array_events,
             prefix_tel=prefix_tel,
+            fov_lon_values=fov_lon_values[valid],
+            fov_lat_values=fov_lat_values[valid],
         )
 
         stereo_table[f"{self.prefix}_alt"] = alt
@@ -907,9 +977,9 @@ def _require_disp_column(self, mono_predictions: Table, prefix_tel: str) -> None
 
     def _apply_min_ang_diff_cut(
         self,
-        mono_predictions: Table,
         valid: np.ndarray,
         tel_to_array_indices: np.ndarray,
+        hillas_psis: u.Quantity,
     ) -> np.ndarray:
         if self.min_ang_diff is None:
             return valid
@@ -921,7 +991,7 @@ def _apply_min_ang_diff_cut(
 
         valid_idx = np.flatnonzero(valid)
         pairs_in_valid = np.flatnonzero(mask_multi2_tels).reshape(-1, 2)
-        valid_psis = mono_predictions["hillas_psi"][valid]
+        valid_psis = hillas_psis[valid]
 
         keep_pairs = check_ang_diff(
             self.min_ang_diff,
@@ -984,6 +1054,8 @@ def _compute_altaz_for_valid(
         tel_to_array_indices: np.ndarray,
         n_array_events: int,
         prefix_tel: str,
+        fov_lon_values: np.ndarray,
+        fov_lat_values: np.ndarray,
     ):
         if np.count_nonzero(valid) == 0:
             nan = u.Quantity(
@@ -994,10 +1066,6 @@ def _compute_altaz_for_valid(
         weights = self._calculate_weights(mono_predictions[valid])
         _, _, valid_multiplicity, _ = get_subarray_index(mono_predictions[valid])
 
-        fov_lon_values, fov_lat_values = calc_fov_lon_lat(
-            mono_predictions[valid], prefix_tel
-        )
-
         combs_array, combs_to_multi_indices = create_combs_array(
             valid_multiplicity.max(), self.n_tel_combinations
         )
@@ -1083,6 +1151,54 @@ def _compute_altaz_for_valid(
 
         return alt, az
 
+    @staticmethod
+    def _transform_to_nominal(
+        tel_pointing_alt,
+        tel_pointing_az,
+        subarray_pointing_alt,
+        subarray_pointing_az,
+        fov_lons,
+        fov_lats,
+    ):
+        """
+        Transform DISP candidates from telescope frames to a shared nominal frame.
+
+        Parameters
+        ----------
+        tel_pointing_alt, tel_pointing_az : astropy.units.Quantity
+            Scalar or array-valued telescope pointings defining the input frames.
+        subarray_pointing_alt, subarray_pointing_az : astropy.units.Quantity
+            Scalar or array-valued subarray pointings defining the nominal frames.
+        fov_lons, fov_lats : astropy.units.Quantity
+            Array-valued DISP candidate coordinates with angular units. The shape
+            is ``(2,)`` for one telescope event or ``(n_events, 2)`` for a table.
+
+        Returns
+        -------
+        nominal_lons, nominal_lats : tuple[numpy.ndarray, numpy.ndarray]
+            Candidate coordinates in the nominal frames, in degrees and with the
+            same shape as ``fov_lons`` and ``fov_lats``.
+        """
+        telescope_pointing = SkyCoord(
+            alt=tel_pointing_alt,
+            az=tel_pointing_az,
+            frame=AltAz(),
+        )
+        array_pointing = SkyCoord(
+            alt=subarray_pointing_alt,
+            az=subarray_pointing_az,
+            frame=AltAz(),
+        )
+        candidates = SkyCoord(
+            fov_lon=fov_lons,
+            fov_lat=fov_lats,
+            frame=TelescopeFrame(telescope_pointing=telescope_pointing),
+        ).transform_to(NominalFrame(origin=array_pointing))
+        return (
+            candidates.fov_lon.to_value(u.deg),
+            candidates.fov_lat.to_value(u.deg),
+        )
+
     def _collect_telescopes_per_event(
         self,
         mono_predictions: Table,

src/ctapipe/reco/tests/test_stereo_combination.py

@@ -110,6 +110,8 @@ def mono_table():
             ],
             "subarray_pointing_lat": 10 * [70] * u.deg,
             "subarray_pointing_lon": 10 * [0] * u.deg,
+            "telescope_pointing_altitude": 10 * [70] * u.deg,
+            "telescope_pointing_azimuth": 10 * [0] * u.deg,
         }
     )
 
@@ -392,6 +394,9 @@ def _make_disp_event(event_dict, prefix="dummy"):
                 )
             },
         )
+        pointing = event.monitoring.tel[event_dict["tel_id"][i]].pointing
+        pointing.azimuth = 0 * u.deg
+        pointing.altitude = 70 * u.deg
 
     event.monitoring.pointing = ArrayPointingContainer(
         array_azimuth=0 * u.deg, array_altitude=70 * u.deg