diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 311142311a3..831ca92ebd5 100644 --- a/src/ctapipe/containers.py +++ b/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): diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index 352073ed959..c31e6d461eb 100644 --- a/src/ctapipe/image/muon/features.py +++ b/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), + ) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index 30b382aef2c..ac11a9a368b 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/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, ) diff --git a/src/ctapipe/image/muon/tests/test_muon_features.py b/src/ctapipe/image/muon/tests/test_muon_features.py index 85a950ec68e..c2cb5ad27fe 100644 --- a/src/ctapipe/image/muon/tests/test_muon_features.py +++ b/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)