Muon parameters update

src/ctapipe/containers.py

@@ -1189,6 +1189,26 @@ class MuonParametersContainer(Container):
         nan * u.deg**2,
         "MSE of the deviation of all pixels after cleaning from the ring fit",
     )
+    ring_intensity = Field(
+        nan, "Sum of the pixel charges inside the integration area around a ring"
+    )
+    intensity_outside_ring = Field(nan, "Sum of the pixel charges outside the ring")
+    n_pixels_in_ring = Field(
+        -1,
+        "Number of pixels inside the ring integration area that passed the cleaning",
+    )
+    mean_intensity_outside_ring = Field(
+        nan,
+        "Mean intensity of pixels inside the region limited by ring integration width and outer ring width.",
+    )
+    standard_dev = Field(
+        nan * u.deg,
+        "Standard deviation of the radial light distribution.",
+    )
+    skewness = Field(nan, "Skewness of the light distribution along the ring radius.")
+    excess_kurtosis = Field(
+        nan, "Excess kurtosis of the light distribution along the ring radius."
+    )
 
 
 class MuonTelescopeContainer(Container):

src/ctapipe/image/muon/features.py

@@ -158,3 +158,125 @@ def ring_containment(radius, center_x, center_y, camera_radius):
 
     a = (radius**2 - camera_radius**2 + d**2) / (2 * d)
     return np.arccos(a / radius) / np.pi
+
+
+def ring_size_parameters(
+    radius,
+    center_x,
+    center_y,
+    pixel_x,
+    pixel_y,
+    ring_integration_width,
+    outer_ring_width,
+    image,
+    image_mask,
+):
+    """
+    Calculate the parameters related to the size of the ring image.
+
+    Parameters
+    ----------
+    radius: float
+        radius of the ring
+    center_x: float
+        x coordinate of the ring center
+    center_y: float
+        y coordinate of the ring center
+    pixel_x: array-like
+        x coordinates of the camera pixels
+    pixel_y: array-like
+        y coordinates of the camera pixels
+    ring_integration_width: float
+        Width of the ring in fractions of ring radius
+    outer_ring_width: float
+        Width of the outer ring in fractions of ring radius
+    image: array-like
+        Amplitude of image pixels
+    image_mask: array-like
+        mask of the camera pixels after cleaning
+
+    Returns
+    -------
+    ring_intensity: float
+        Sum of the p.e. inside the integration area of the ring
+    intensity_outside_ring: float
+        Sum of the p.e. outside the ring integration area that passed the cleaning
+    n_pixels_in_ring: int
+        Number of pixels inside the ring integration area that passed the cleaning
+    mean_intensity_outside_ring: float
+        Mean intensity of the pixels outside the integration area of the ring,
+        and restricted by the outer ring width, i.e. in the strip between
+        ring integration width and outer ring width.
+    """
+
+    dist = np.sqrt((pixel_x - center_x) ** 2 + (pixel_y - center_y) ** 2)
+    dist_mask = np.abs(dist - radius) < (radius * ring_integration_width)
+    pix_ring = image * dist_mask
+    pix_outside_ring = image * ~dist_mask
+
+    dist_mask_2 = np.logical_and(
+        ~dist_mask,
+        np.abs(dist - radius) < radius * (ring_integration_width + outer_ring_width),
+    )
+    pix_ring_2 = image[dist_mask_2]
+
+    ring_intensity = np.sum(pix_ring)
+    intensity_outside_ring = np.sum(pix_outside_ring * image_mask)
+    n_pixels_in_ring = np.sum(dist_mask & image_mask)
+    mean_intensity_outside_ring = np.sum(pix_ring_2) / len(pix_ring_2)
+
+    return (
+        ring_intensity,
+        intensity_outside_ring,
+        n_pixels_in_ring,
+        mean_intensity_outside_ring,
+    )
+
+
+def radial_light_distribution(center_x, center_y, pixel_x, pixel_y, image):
+    """
+    Calculate the radial distribution of the muon ring.
+
+    Parameters
+    ----------
+    center_x : float
+        x coordinate of the ring center.
+    center_y : float
+        y coordinate of the ring center.
+    pixel_x : array-like
+        x coordinates of the camera pixels.
+    pixel_y : array-like
+        y coordinates of the camera pixels.
+        Amplitude of image pixels.
+
+    Returns
+    -------
+    standard_dev : `astropy.units.Quantity`
+        Standard deviation of the light distribution in degrees.
+        Spread of pixel intensities around the mean radial distance from the ring center.
+
+    skewness : float
+        Skewness of the radial light distribution.
+        Measures the asymmetry of the distribution around the mean radius.
+
+    excess_kurtosis : float
+        Excess kurtosis of the radial light distribution.
+        Indicates the "tailedness" of the distribution compared to a normal distribution.
+    """
+
+    if np.sum(image) == 0:
+        return np.nan * u.deg, np.nan, np.nan
+
+    pixel_r = np.hypot(pixel_x - center_x, pixel_y - center_y)
+
+    mean = np.average(pixel_r, weights=image)
+    delta_r = pixel_r - mean
+    standard_dev = np.sqrt(np.average(delta_r**2, weights=image))
+    skewness = np.average(delta_r**3, weights=image) / standard_dev**3
+    excess_kurtosis = np.average(delta_r**4, weights=image) / standard_dev**4 - 3.0
+
+    return (
+        standard_dev,
+        skewness.to_value(u.dimensionless_unscaled),
+        excess_kurtosis.to_value(u.dimensionless_unscaled),
+    )

src/ctapipe/image/muon/processor.py

@@ -15,8 +15,10 @@
 from .features import (
     intensity_ratio_inside_ring,
     mean_squared_error,
+    radial_light_distribution,
     ring_completeness,
     ring_containment,
+    ring_size_parameters,
 )
 from .intensity_fitter import MuonIntensityFitter
 from .ring_fitter import MuonRingFitter
@@ -95,6 +97,22 @@ class MuonProcessor(TelescopeComponent):
         help="Pedestal noise rms", default_value=1.1
     ).tag(config=True)
 
+    ring_integration_width = FloatTelescopeParameter(
+        default_value=0.25,
+        help=(
+            "Width of the ring in units of the ring radius, "
+            "used for computing the ring size in charge units."
+        ),
+    ).tag(config=True)
+
+    outer_ring_width = FloatTelescopeParameter(
+        default_value=0.2,
+        help=(
+            "Width of the outer ring in units of the ring radius, "
+            "used for computing the charge outside the ring."
+        ),
+    ).tag(config=True)
+
     def __init__(self, subarray, **kwargs):
         super().__init__(subarray, **kwargs)
         self.dl1_query = ImageParameterQuery(parent=self)
@@ -268,9 +286,41 @@ def _calculate_muon_parameters(
             ring.center_fov_lat,
         )
 
+        (
+            ring_intensity,
+            intensity_outside_ring,
+            n_pixels_in_ring,
+            mean_intensity_outside_ring,
+        ) = ring_size_parameters(
+            ring.radius,
+            ring.center_fov_lon,
+            ring.center_fov_lat,
+            fov_lon,
+            fov_lat,
+            self.ring_integration_width.tel[tel_id],
+            self.outer_ring_width.tel[tel_id],
+            image,
+            clean_mask,
+        )
+
+        standard_dev, skewness, excess_kurtosis = radial_light_distribution(
+            ring.center_fov_lon,
+            ring.center_fov_lat,
+            fov_lon,
+            fov_lat,
+            image,
+        )
+
         return MuonParametersContainer(
             containment=containment,
             completeness=completeness,
             intensity_ratio=intensity_ratio,
             mean_squared_error=mse,
+            ring_intensity=ring_intensity,
+            intensity_outside_ring=intensity_outside_ring,
+            n_pixels_in_ring=n_pixels_in_ring,
+            mean_intensity_outside_ring=mean_intensity_outside_ring,
+            standard_dev=standard_dev,
+            skewness=skewness,
+            excess_kurtosis=excess_kurtosis,
         )

src/ctapipe/image/muon/tests/test_muon_features.py

@@ -3,7 +3,12 @@
 import astropy.units as u
 import numpy as np
 
-from ctapipe.image.muon.features import ring_completeness, ring_containment
+from ctapipe.image.muon.features import (
+    radial_light_distribution,
+    ring_completeness,
+    ring_containment,
+    ring_size_parameters,
+)
 
 
 def test_ring_containment():
@@ -50,3 +55,70 @@ def test_ring_completeness():
 
     assert ring_comp <= 1
     assert ring_comp >= 0
+
+
+def test_ring_size_parameters():
+    radius = 1
+    center_x = 0
+    center_y = 0
+    pixel_x = np.linspace(-2, 2, 1855)
+    pixel_y = np.linspace(-2, 2, 1855)
+    ring_integration_width = 0.25
+    outer_ring_width = 0.2
+    image = np.random.normal(loc=100, scale=10, size=1855)
+    image_mask = np.random.choice([True, False], size=1855)
+
+    (
+        ring_size,
+        size_outside,
+        num_pixels_in_ring,
+        mean_pixel_outside_ring,
+    ) = ring_size_parameters(
+        radius,
+        center_x,
+        center_y,
+        pixel_x,
+        pixel_y,
+        ring_integration_width,
+        outer_ring_width,
+        image,
+        image_mask,
+    )
+
+    assert ring_size > 0
+    assert size_outside > 0
+    assert num_pixels_in_ring > 0
+    assert mean_pixel_outside_ring > 0
+
+
+def test_radial_light_distribution():
+    center_x = 0.0 * u.deg
+    center_y = 0.0 * u.deg
+    pixel_x = np.linspace(-10, 10, 1855) * u.deg
+    pixel_y = np.linspace(-10, 10, 1855) * u.deg
+    image = np.random.normal(loc=100, scale=10, size=1855)
+
+    standard_dev, skewness, excess_kurtosis = radial_light_distribution(
+        center_x, center_y, pixel_x, pixel_y, image
+    )
+
+    assert standard_dev.unit == u.deg
+    assert np.isfinite(standard_dev.value)
+    assert np.isfinite(skewness)
+    assert np.isfinite(excess_kurtosis)
+
+
+def test_radial_light_distribution_zero_image():
+    center_x = 0.0 * u.deg
+    center_y = 0.0 * u.deg
+    pixel_x = np.linspace(-10, 10, 1855) * u.deg
+    pixel_y = np.linspace(-10, 10, 1855) * u.deg
+    image = np.zeros(1855)
+
+    standard_dev, skewness, excess_kurtosis = radial_light_distribution(
+        center_x, center_y, pixel_x, pixel_y, image
+    )
+
+    assert np.isnan(standard_dev)
+    assert np.isnan(skewness)
+    assert np.isnan(excess_kurtosis)