Merge pull request #2629 from cta-observatory/hillas_psi_uncertainty_ruo

Hillas psi uncertainty ruo

Author: kosack

docs/changes/2629.feature.rst

@@ -0,0 +1,3 @@
+Add two new fields in the Hillas parameters computations:
+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)

src/ctapipe/containers.py

@@ -310,6 +310,12 @@ 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,
+        "uncertainty on the center of gravity along the transverse axis of the image",
+        unit=u.m,
+    )
 
 
 class HillasParametersContainer(BaseHillasParametersContainer):
@@ -338,6 +344,12 @@ 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,
+        "uncertainty on the center of gravity along the transverse axis of the image",
+        unit=u.deg,
+    )
 
 
 class LeakageContainer(Container):

src/ctapipe/image/hillas.py

@@ -4,6 +4,7 @@
 Hillas-style moment-based shower image parametrization.
 """
 
+
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import Angle
@@ -87,6 +88,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 +133,41 @@ 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
+        # p[0] is the psi angle in the rotated frame, which should be zero.
+        # Now we add the non-zero residual psi angle in the rotated frame to psi uncertainty
+        # We also add additional uncertainty to account for elongation of the image (i.e. width / length)
+        psi_uncert = np.sqrt(lsq_cov[0, 0] + p[0] * p[0]) * (
+            1.0 + pow(np.tan(width / length * np.pi / 2.0), 2)
+        )
+        sin_p0 = np.sin(p[0])
+        transverse_cog_uncert = np.sqrt(lsq_cov[1, 1] * (1.0 + sin_p0 * sin_p0))
+
     # Compute of the Hillas parameters uncertainties.
     # Implementation described in [hillas_uncertainties]_ This is an internal MAGIC document
     # not generally accessible.
@@ -190,6 +211,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 +227,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/image/tests/test_hillas.py

@@ -253,6 +253,69 @@ def test_single_pixel():
     assert hillas.length.value == 0
     assert hillas.width.value == 0
     assert np.isnan(hillas.psi)
+    assert np.isnan(hillas.psi_uncertainty)
+
+
+def test_psi_uncertainty():
+    x = y = np.arange(3)
+    x, y = np.meshgrid(x, y)
+
+    geom = CameraGeometry(
+        name="testcam",
+        pix_id=np.arange(9),
+        pix_x=x.ravel() * u.cm,
+        pix_y=y.ravel() * u.cm,
+        pix_type="rectangular",
+        pix_area=np.full(9, 1.0) * u.cm**2,
+    )
+
+    image = np.zeros((3, 3))
+    image[0, 0] = 2
+    image[0, 1] = 0
+    image[0, 2] = 0
+    image[1, 0] = 0
+    image[1, 1] = 5
+    image[1, 2] = 0
+    image[2, 0] = 0
+    image[2, 1] = 0
+    image[2, 2] = 2
+    image = image.ravel()
+    hillas = hillas_parameters(geom, image)
+    assert np.isclose(
+        hillas.psi_uncertainty.to_value(u.deg), 20.25711711353489, rtol=1e-05
+    )
+
+    image = np.zeros((3, 3))
+    image[0, 0] = 5
+    image[0, 1] = 0
+    image[0, 2] = 0
+    image[1, 0] = 0
+    image[1, 1] = 10
+    image[1, 2] = 0
+    image[2, 0] = 0
+    image[2, 1] = 0
+    image[2, 2] = 5
+    image = image.ravel()
+    hillas = hillas_parameters(geom, image)
+    assert np.isclose(
+        hillas.psi_uncertainty.to_value(u.deg), 12.811725781509189, rtol=1e-05
+    )
+
+    image = np.zeros((3, 3))
+    image[0, 0] = 5
+    image[0, 1] = 0
+    image[0, 2] = 2
+    image[1, 0] = 0
+    image[1, 1] = 10
+    image[1, 2] = 0
+    image[2, 0] = 0
+    image[2, 1] = 2
+    image[2, 2] = 5
+    image = image.ravel()
+    hillas = hillas_parameters(geom, image)
+    assert np.isclose(
+        hillas.psi_uncertainty.to_value(u.deg), 23.560635267712282, rtol=1e-05
+    )
 
 
 def test_reconstruction_in_telescope_frame(prod5_lst):

src/ctapipe/visualization/mpl_camera.py

@@ -511,6 +511,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 = ""