Hillas psi uncertainty ruo

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)
+    b_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)
+    b_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,33 @@ 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 = b_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])
+        b_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 +202,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),
+            b_uncertainty=u.Quantity(b_uncert, unit),
             skewness=skewness_long,
             kurtosis=kurtosis_long,
         )
@@ -204,6 +218,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),
+        b_uncertainty=u.Quantity(b_uncert, unit),
         skewness=skewness_long,
         kurtosis=kurtosis_long,
     )

src/ctapipe/tests/test_hillas_psi_unc_toy_model.py

@@ -0,0 +1,119 @@
+import os
+
+import astropy.units as u
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.stats import norm
+from tqdm.auto import tqdm
+
+from ctapipe.image import hillas_parameters
+from ctapipe.image.cleaning import tailcuts_clean
+from ctapipe.image.toymodel import Gaussian
+from ctapipe.instrument import CameraGeometry
+
+ctapipe_output = os.environ.get("CTAPIPE_OUTPUT_PATH")
+ctapipe_input = os.environ.get("CTAPIPE_SVC_PATH")
+
+rng = np.random.default_rng(0)
+
+# cam = CameraGeometry.from_name("LSTCam")
+cam = CameraGeometry.from_name("NectarCam")
+# cam = CameraGeometry.from_name("SCTCam")
+
+true_width = 0.05 * u.m
+true_length = 0.3 * u.m
+true_psi = 45 * u.deg
+true_x = 0.5 * u.m
+true_y = -0.2 * u.m
+
+image_intensity = 1000
+test_nsb_level_pe = 3
+# test_nsb_level_pe = 5
+
+# n_sample = 10
+n_sample = 1000
+
+model = Gaussian(true_x, true_y, true_length, true_width, true_psi)
+
+
+def sample_no_noise_no_cleaning():
+    _, signal, _ = model.generate_image(
+        cam, intensity=image_intensity, nsb_level_pe=0, rng=rng
+    )
+    h = hillas_parameters(cam, signal)
+    return h
+
+
+def sample_no_noise_with_cleaning():
+    _, signal, _ = model.generate_image(
+        cam, intensity=image_intensity, nsb_level_pe=0, rng=rng
+    )
+
+    mask = tailcuts_clean(
+        cam,
+        signal,
+        3.0 * test_nsb_level_pe,
+        2.0 * test_nsb_level_pe,
+        min_number_picture_neighbors=2,
+    )
+
+    image_clean = np.zeros_like(signal)
+    for pix in range(0, len(signal)):
+        image_clean[pix] = max(0.0, signal[pix] - 2.0 * test_nsb_level_pe)
+
+    h = hillas_parameters(cam[mask], image_clean[mask])
+    return h
+
+
+def sample_noise_with_cleaning():
+    image, _, _ = model.generate_image(
+        cam, intensity=image_intensity, nsb_level_pe=test_nsb_level_pe, rng=rng
+    )
+
+    mask = tailcuts_clean(
+        cam,
+        image,
+        3.0 * test_nsb_level_pe,
+        2.0 * test_nsb_level_pe,
+        min_number_picture_neighbors=2,
+    )
+
+    image_clean = np.zeros_like(image)
+    for pix in range(0, len(image)):
+        image_clean[pix] = max(0.0, image[pix] - 2.0 * test_nsb_level_pe)
+
+    h = hillas_parameters(cam[mask], image_clean[mask])
+    return h
+
+
+trials_no_noise_no_cleaning = [
+    sample_no_noise_no_cleaning() for _ in tqdm(range(n_sample))
+]
+trials_no_noise_with_cleaning = [
+    sample_no_noise_with_cleaning() for _ in tqdm(range(n_sample))
+]
+trials_noise_cleaning = [sample_noise_with_cleaning() for _ in tqdm(range(n_sample))]
+
+titles = [
+    "No Noise, all Pixels",
+    f"No Noise, Tailcuts({test_nsb_level_pe*3}, {test_nsb_level_pe*2})",
+    f"With Noise ({test_nsb_level_pe} p.e.), Tailcuts({test_nsb_level_pe*3}, {test_nsb_level_pe*2})",
+]
+values = [
+    trials_no_noise_no_cleaning,
+    trials_no_noise_with_cleaning,
+    trials_noise_cleaning,
+]
+
+fig, axs = plt.subplots(3, 1, constrained_layout=True, sharex=True)
+for ax, trials, title in zip(axs, values, titles):
+    pred = np.array([t.psi.to_value(u.rad) for t in trials])
+    unc = np.array([t.psi_uncertainty.to_value(u.rad) for t in trials])
+    limits = np.quantile(pred, [0.001, 0.999])
+    hist, edges, plot = ax.hist(pred, bins=51, range=limits, density=True)
+    x = np.linspace(edges[0], edges[-1], 500)
+    ax.plot(x, norm.pdf(x, pred.mean(), pred.std()))
+    ax.plot(x, norm.pdf(x, pred.mean(), unc.mean()))
+    ax.set_title(title)
+axs[2].set_xlabel("Psi / rad")
+fig.savefig(f"{ctapipe_output}/output_plots/hillas_uncertainties.png", dpi=300)

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 = ""