diff --git a/src/ctapipe/reco/stereo_combination.py b/src/ctapipe/reco/stereo_combination.py index 6b6586b678b..e7c2a6e6c5e 100644 --- a/src/ctapipe/reco/stereo_combination.py +++ b/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, diff --git a/src/ctapipe/reco/tests/test_stereo_combination.py b/src/ctapipe/reco/tests/test_stereo_combination.py index ccff3578f82..7f61950e5c5 100644 --- a/src/ctapipe/reco/tests/test_stereo_combination.py +++ b/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