Hillas psi uncertainty ruo

docs/changes/2629.feature.rst

@@ -0,0 +1,15 @@
+This commit adds two new hillas-parameter variables:
+1. psi_uncertainty (uncertainty on the psi angle of the image)
+2. transverse_cog_uncertainty (uncertainty on the center of gravity along the transverse axis of the image)
+
+The method that estimates these uncertainties has been presented at these ASWG calls:
+https://indico.cta-observatory.org/event/5847/contributions/47481/attachments/26575/39053/ctapipe_rethink_hillas_20240718.pdf
+https://indico.cta-observatory.org/event/5900/contributions/47899/attachments/26789/39379/Psi_uncertainty_short_summary_20241022.pdf
+
+The affected files are:
+image/hillas.py
+containers.py
+
+Caution:
+These uncertainties are calculated using the the information of pixel area and the charges in the pixels.
+In order to obtain the accurate uncertainties, the charges in the image should be counted above the NSB noise fluctuation or the boundary threshold for the image cleaning.

src/ctapipe/containers.py

@@ -301,6 +301,10 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer):
     width = Field(nan * u.m, "standard spread along the minor-axis", unit=u.m)
     width_uncertainty = Field(nan * u.m, "uncertainty of width", unit=u.m)
     psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)
+    psi_uncertainty = Field(nan * u.deg, "uncertainty of psi", unit=u.deg)
+    transverse_cog_uncertainty = Field(
+        nan * u.m, "rotated centroid y coordinate uncertainty", unit=u.m
+    )
 
 
 class HillasParametersContainer(BaseHillasParametersContainer):
@@ -329,6 +333,10 @@ class HillasParametersContainer(BaseHillasParametersContainer):
     width = Field(nan * u.deg, "standard spread along the minor-axis", unit=u.deg)
     width_uncertainty = Field(nan * u.deg, "uncertainty of width", unit=u.deg)
     psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)
+    psi_uncertainty = Field(nan * u.deg, "uncertainty of psi", unit=u.deg)
+    transverse_cog_uncertainty = Field(
+        nan * u.deg, "rotated centroid y coordinate uncertainty", unit=u.deg
+    )
 
 
 class LeakageContainer(Container):

src/ctapipe/image/hillas.py

@@ -87,6 +87,7 @@ def hillas_parameters(geom, image):
     unit = geom.pix_x.unit
     pix_x = geom.pix_x.to_value(unit)
     pix_y = geom.pix_y.to_value(unit)
+    pix_area = geom.pix_area.to_value(unit**2)
     image = np.asanyarray(image, dtype=np.float64)
 
     if isinstance(image, np.ma.masked_array):
@@ -131,22 +132,39 @@ def hillas_parameters(geom, image):
     # avoid divide by 0 warnings
     # psi will be consistently defined in the range (-pi/2, pi/2)
     if length == 0:
-        psi = skewness_long = kurtosis_long = np.nan
+        psi = (
+            psi_uncert
+        ) = transverse_cog_uncert = skewness_long = kurtosis_long = np.nan
     else:
         if vx != 0:
             psi = np.arctan(vy / vx)
         else:
             psi = np.pi / 2
 
         # calculate higher order moments along shower axes
-        longitudinal = delta_x * np.cos(psi) + delta_y * np.sin(psi)
+        cos_psi = np.cos(psi)
+        sin_psi = np.sin(psi)
+        longi = delta_x * cos_psi + delta_y * sin_psi
+        trans = delta_x * -sin_psi + delta_y * cos_psi
 
-        m3_long = np.average(longitudinal**3, weights=image)
+        m3_long = np.average(longi**3, weights=image)
         skewness_long = m3_long / length**3
 
-        m4_long = np.average(longitudinal**4, weights=image)
+        m4_long = np.average(longi**4, weights=image)
         kurtosis_long = m4_long / length**4
 
+        # lsq solution to determine uncertainty on phi
+        W = np.diag(image / pix_area)
+        X = np.column_stack([longi, np.ones_like(longi)])
+        lsq_cov = np.linalg.inv(X.T @ W @ X)
+        p = lsq_cov @ X.T @ W @ trans
+        psi_uncert = np.sqrt(lsq_cov[0, 0] + p[0] * p[0]) * (
+            1.0 + pow(np.tan(width / length * np.pi / 2.0), 2)
+        )
+        transverse_cog_uncert = np.sqrt(
+            lsq_cov[1, 1] * (1.0 + np.sin(p[0]) * np.sin(p[0]))
+        )
+
     # Compute of the Hillas parameters uncertainties.
     # Implementation described in [hillas_uncertainties]_ This is an internal MAGIC document
     # not generally accessible.
@@ -190,6 +208,8 @@ def hillas_parameters(geom, image):
             width=u.Quantity(width, unit),
             width_uncertainty=u.Quantity(width_uncertainty, unit),
             psi=Angle(psi, unit=u.rad),
+            psi_uncertainty=Angle(psi_uncert, unit=u.rad),
+            transverse_cog_uncertainty=u.Quantity(transverse_cog_uncert, unit),
             skewness=skewness_long,
             kurtosis=kurtosis_long,
         )
@@ -204,6 +224,8 @@ def hillas_parameters(geom, image):
         width=u.Quantity(width, unit),
         width_uncertainty=u.Quantity(width_uncertainty, unit),
         psi=Angle(psi, unit=u.rad),
+        psi_uncertainty=Angle(psi_uncert, unit=u.rad),
+        transverse_cog_uncertainty=u.Quantity(transverse_cog_uncert, unit),
         skewness=skewness_long,
         kurtosis=kurtosis_long,
     )

src/ctapipe/visualization/mpl_camera.py

@@ -496,6 +496,24 @@ def overlay_moments(
             **kwargs,
         )
 
+        for sign in (-1, 1):
+            for direction in (-1, 1):
+                angle = params["psi_rad"] + sign * params["psi_uncert_rad"]
+                (line,) = self.axes.plot(
+                    [
+                        params["cog_x"],
+                        params["cog_x"]
+                        + direction * np.cos(angle) * 3 * params["length"],
+                    ],
+                    [
+                        params["cog_y"],
+                        params["cog_y"]
+                        + direction * np.sin(angle) * 3 * params["length"],
+                    ],
+                    **kwargs,
+                )
+                self._axes_overlays.append(line)
+
         self._axes_overlays.append(el)
 
         if with_label:

src/ctapipe/visualization/utils.py

@@ -31,11 +31,12 @@ def build_hillas_overlay(hillas, unit, with_label=True, n_sigma=1):
     psi_deg = Angle(hillas.psi).wrap_at(180 * u.deg).to_value(u.deg)
 
     ret = dict(
-        cog_x=cog_x,
-        cog_y=cog_y,
-        width=width,
-        length=length,
-        psi_rad=psi_rad,
+        cog_x=float(cog_x),
+        cog_y=float(cog_y),
+        width=float(width),
+        length=float(length),
+        psi_rad=float(psi_rad),
+        psi_uncert_rad=float(hillas.psi_uncertainty.to_value(u.rad)),
     )
 
     if not with_label:
@@ -67,9 +68,9 @@ def build_hillas_overlay(hillas, unit, with_label=True, n_sigma=1):
         label_y = cog_y - r * np.sin(psi_rad)
         rotation = psi_deg + 90
 
-    ret["rotation"] = rotation
-    ret["label_x"] = label_x
-    ret["label_y"] = label_y
+    ret["rotation"] = float(rotation)
+    ret["label_x"] = float(label_x)
+    ret["label_y"] = float(label_y)
 
     if unit == u.deg:
         sep = ""