From 6d55373348c1c9728b9f99215c65b13e683d2ffe Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Mon, 9 Dec 2024 15:09:46 +0100 Subject: [PATCH 01/12] Added fields to the MuonParametersContainer --- src/ctapipe/containers.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 311142311a3..b3b0b2d833e 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1189,6 +1189,22 @@ class MuonParametersContainer(Container): nan * u.deg**2, "MSE of the deviation of all pixels after cleaning from the ring fit", ) + ring_size = Field(nan, "Number of pixels in the ring") + size_outside = Field(nan, "Number of pixels outside the ring") + num_pixels_in_ring = Field(nan, "Number of pixels in the ring") + mean_pixel_outside_ring = Field(nan, "Mean pixel charge outside the ring") + standard_dev = Field( + nan * u.deg, + "Standard deviation of the light distribution along the ring radius." + ) + 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): From a1ab58af9997287e419aa05ec209246d7aa1b002 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Mon, 9 Dec 2024 15:10:41 +0100 Subject: [PATCH 02/12] Added methods to calculate new fields --- src/ctapipe/image/muon/features.py | 105 ++++++++++++++++++++++++++++- 1 file changed, 104 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index 352073ed959..39147a506fe 100644 --- a/src/ctapipe/image/muon/features.py +++ b/src/ctapipe/image/muon/features.py @@ -10,7 +10,6 @@ "ring_containment", ] - def mean_squared_error(pixel_x, pixel_y, weights, radius, center_x, center_y): """ Calculate the weighted mean squared error for a circle @@ -158,3 +157,107 @@ 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_size: float + Sum of the p.e. inside the integration area of the ring + size_outside: float + Sum of the photons outside the ring integration area that passed the cleaning + num_pixels_in_ring: int + Number of pixels inside the ring integration area that passed the cleaning + mean_pixel_outside_ring: float + Mean intensity of the pixels outside the ring, but still close to it + """ + + 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_size = np.sum(pix_ring) + size_outside = np.sum(pix_outside_ring * image_mask) + num_pixels_in_ring = np.sum(dist_mask & image_mask) + mean_pixel_outside_ring = (np.sum(pix_ring_2) / len(pix_ring_2)) + + return ring_size, size_outside, num_pixels_in_ring, mean_pixel_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 : float + Standard deviation of the light distribution along the ring radius. + skewness : float + Skewness of the light distribution along the ring radius. + excess_kurtosis : float + Excess kurtosis of the light distribution along the ring radius. + """ + + + if np.sum(image) == 0: + return np.nan * u.deg, np.nan, np.nan + + x0 = center_x.to_value(u.deg) + y0 = center_y.to_value(u.deg) + pix_x = pixel_x.to_value(u.deg) + pix_y = pixel_y.to_value(u.deg) + pixel_r = np.sqrt((pix_x - x0) ** 2 + (pix_y - y0) ** 2) + + 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. + + return standard_dev * u.deg, skewness, excess_kurtosis From 385e0397a66d685030dd0d1aef7adadf4c30cae4 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Mon, 9 Dec 2024 15:32:34 +0100 Subject: [PATCH 03/12] Added tests for new muon feature methods --- .../image/muon/tests/test_muon_features.py | 69 ++++++++++++++++++- 1 file changed, 68 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/tests/test_muon_features.py b/src/ctapipe/image/muon/tests/test_muon_features.py index 85a950ec68e..287588216e6 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 ( + ring_completeness, + ring_containment, + ring_size_parameters, + radial_light_distribution +) def test_ring_containment(): @@ -50,3 +55,65 @@ 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) From dca7ff8add1c03a25b25ab3eb1ebe27c7582b24f Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Mon, 9 Dec 2024 15:43:45 +0100 Subject: [PATCH 04/12] Added processing of the new features during muon analysis --- src/ctapipe/image/muon/processor.py | 52 ++++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index 30b382aef2c..99d745f6511 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -17,6 +17,8 @@ mean_squared_error, ring_completeness, ring_containment, + ring_size_parameters, + radial_light_distribution ) 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) @@ -173,7 +191,7 @@ def _process_telescope_event(self, event, tel_id): parameters = self._calculate_muon_parameters( tel_id, image, dl1.image_mask, ring ) - + checks = self.ring_query(parameters=parameters, ring=ring, mask=mask) if not all(checks): event.muon.tel[tel_id] = MuonTelescopeContainer( @@ -268,9 +286,41 @@ def _calculate_muon_parameters( ring.center_fov_lat, ) + ( + ring_size, + size_outside, + num_pixels_in_ring, + mean_pixel_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_size=ring_size, + size_outside=size_outside, + num_pixels_in_ring=num_pixels_in_ring, + mean_pixel_outside_ring=mean_pixel_outside_ring, + standard_dev=standard_dev, + skewness=skewness, + excess_kurtosis=excess_kurtosis, ) From 69b496a56ce08d2860e88a2780f409e818d0bae0 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Tue, 10 Dec 2024 11:10:45 +0100 Subject: [PATCH 05/12] Fixed lint issues --- src/ctapipe/containers.py | 10 ++--- src/ctapipe/image/muon/features.py | 37 +++++++++++-------- src/ctapipe/image/muon/processor.py | 20 +++++----- .../image/muon/tests/test_muon_features.py | 15 +++++--- 4 files changed, 45 insertions(+), 37 deletions(-) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index b3b0b2d833e..aaa68f8bbd8 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1195,15 +1195,11 @@ class MuonParametersContainer(Container): mean_pixel_outside_ring = Field(nan, "Mean pixel charge outside the ring") standard_dev = Field( nan * u.deg, - "Standard deviation of the light distribution along the ring radius." - ) - skewness = Field( - nan, - "Skewness of the light distribution along the ring radius." + "Standard deviation of the light distribution along the ring radius.", ) + 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." + nan, "Excess kurtosis of the light distribution along the ring radius." ) diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index 39147a506fe..f3871a8f908 100644 --- a/src/ctapipe/image/muon/features.py +++ b/src/ctapipe/image/muon/features.py @@ -10,6 +10,7 @@ "ring_containment", ] + def mean_squared_error(pixel_x, pixel_y, weights, radius, center_x, center_y): """ Calculate the weighted mean squared error for a circle @@ -160,7 +161,15 @@ def ring_containment(radius, center_x, center_y, camera_radius): def ring_size_parameters( - radius, center_x, center_y, pixel_x, pixel_y, ring_integration_width, outer_ring_width, image, image_mask + 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. @@ -202,19 +211,18 @@ def ring_size_parameters( 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) - ) + + 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_size = np.sum(pix_ring) size_outside = np.sum(pix_outside_ring * image_mask) num_pixels_in_ring = np.sum(dist_mask & image_mask) - mean_pixel_outside_ring = (np.sum(pix_ring_2) / len(pix_ring_2)) - + mean_pixel_outside_ring = np.sum(pix_ring_2) / len(pix_ring_2) + return ring_size, size_outside, num_pixels_in_ring, mean_pixel_outside_ring @@ -244,10 +252,9 @@ def radial_light_distribution(center_x, center_y, pixel_x, pixel_y, image): Excess kurtosis of the light distribution along the ring radius. """ - if np.sum(image) == 0: return np.nan * u.deg, np.nan, np.nan - + x0 = center_x.to_value(u.deg) y0 = center_y.to_value(u.deg) pix_x = pixel_x.to_value(u.deg) @@ -256,8 +263,8 @@ def radial_light_distribution(center_x, center_y, pixel_x, pixel_y, image): 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. + 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 * u.deg, skewness, excess_kurtosis diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index 99d745f6511..dc928a57496 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -15,10 +15,10 @@ from .features import ( intensity_ratio_inside_ring, mean_squared_error, + radial_light_distribution, ring_completeness, ring_containment, ring_size_parameters, - radial_light_distribution ) from .intensity_fitter import MuonIntensityFitter from .ring_fitter import MuonRingFitter @@ -98,15 +98,15 @@ class MuonProcessor(TelescopeComponent): ).tag(config=True) ring_integration_width = FloatTelescopeParameter( - default_value=0.25, + default_value=0.25, help=( - "Width of the ring in units of the ring radius, " + "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, + default_value=0.2, help=( "Width of the outer ring in units of the ring radius, " "used for computing the charge outside the ring." @@ -191,7 +191,7 @@ def _process_telescope_event(self, event, tel_id): parameters = self._calculate_muon_parameters( tel_id, image, dl1.image_mask, ring ) - + checks = self.ring_query(parameters=parameters, ring=ring, mask=mask) if not all(checks): event.muon.tel[tel_id] = MuonTelescopeContainer( @@ -288,9 +288,9 @@ def _calculate_muon_parameters( ( ring_size, - size_outside, - num_pixels_in_ring, - mean_pixel_outside_ring + size_outside, + num_pixels_in_ring, + mean_pixel_outside_ring, ) = ring_size_parameters( ring.radius, ring.center_fov_lon, @@ -322,5 +322,5 @@ def _calculate_muon_parameters( mean_pixel_outside_ring=mean_pixel_outside_ring, standard_dev=standard_dev, skewness=skewness, - excess_kurtosis=excess_kurtosis, + 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 287588216e6..c2cb5ad27fe 100644 --- a/src/ctapipe/image/muon/tests/test_muon_features.py +++ b/src/ctapipe/image/muon/tests/test_muon_features.py @@ -4,10 +4,10 @@ import numpy as np from ctapipe.image.muon.features import ( - ring_completeness, + radial_light_distribution, + ring_completeness, ring_containment, - ring_size_parameters, - radial_light_distribution + ring_size_parameters, ) @@ -68,7 +68,12 @@ def test_ring_size_parameters(): 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( + ( + ring_size, + size_outside, + num_pixels_in_ring, + mean_pixel_outside_ring, + ) = ring_size_parameters( radius, center_x, center_y, @@ -77,7 +82,7 @@ def test_ring_size_parameters(): ring_integration_width, outer_ring_width, image, - image_mask + image_mask, ) assert ring_size > 0 From a54658a0488bbfd4eb472f2887e11f3bfcfbb490 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Wed, 18 Dec 2024 16:02:22 +0100 Subject: [PATCH 06/12] Updated field naming in the muon parameter containers --- src/ctapipe/containers.py | 12 ++++++++---- src/ctapipe/image/muon/features.py | 13 +++++++++---- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index aaa68f8bbd8..b509d84bda0 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1189,10 +1189,14 @@ class MuonParametersContainer(Container): nan * u.deg**2, "MSE of the deviation of all pixels after cleaning from the ring fit", ) - ring_size = Field(nan, "Number of pixels in the ring") - size_outside = Field(nan, "Number of pixels outside the ring") - num_pixels_in_ring = Field(nan, "Number of pixels in the ring") - mean_pixel_outside_ring = Field(nan, "Mean pixel charge outside the ring") + 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(nan, "Number of pixels in the ring") + mean_intensity_outside_ring = Field( + nan, "Mean intensity of pixels outside the ring" + ) standard_dev = Field( nan * u.deg, "Standard deviation of the light distribution along the ring radius.", diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index f3871a8f908..99e7e19e79f 100644 --- a/src/ctapipe/image/muon/features.py +++ b/src/ctapipe/image/muon/features.py @@ -244,12 +244,17 @@ def radial_light_distribution(center_x, center_y, pixel_x, pixel_y, image): Returns ------- - standard_dev : float - Standard deviation of the light distribution along the ring radius. + 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 light distribution along the ring radius. + Skewness of the radial light distribution. + Measures the asymmetry of the distribution around the mean radius. + excess_kurtosis : float - Excess kurtosis of the light distribution along the ring radius. + Excess kurtosis of the radial light distribution. + Indicates the "tailedness" of the distribution compared to a normal distribution. """ if np.sum(image) == 0: From d2677d903a5324630d6dcbbdee3d4ba96744b0e5 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Wed, 18 Dec 2024 16:06:57 +0100 Subject: [PATCH 07/12] Updated doc string in muon parameters container --- src/ctapipe/containers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index b509d84bda0..7619c88c417 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1199,7 +1199,7 @@ class MuonParametersContainer(Container): ) standard_dev = Field( nan * u.deg, - "Standard deviation of the light distribution along the ring radius.", + "Standard deviation of the radial light distribution.", ) skewness = Field(nan, "Skewness of the light distribution along the ring radius.") excess_kurtosis = Field( From 87831c19dd96a57c0b0730bc90ca9dc9b007fa33 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Wed, 18 Dec 2024 16:09:14 +0100 Subject: [PATCH 08/12] Proper variables naming to be consistent with recent changes in container --- src/ctapipe/image/muon/features.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index 99e7e19e79f..1159bbdad94 100644 --- a/src/ctapipe/image/muon/features.py +++ b/src/ctapipe/image/muon/features.py @@ -197,13 +197,13 @@ def ring_size_parameters( Returns ------- - ring_size: float + ring_intensity: float Sum of the p.e. inside the integration area of the ring - size_outside: float + intensity_outside_ring: float Sum of the photons outside the ring integration area that passed the cleaning - num_pixels_in_ring: int + n_pixels_in_ring: int Number of pixels inside the ring integration area that passed the cleaning - mean_pixel_outside_ring: float + mean_intensity_outside_ring: float Mean intensity of the pixels outside the ring, but still close to it """ @@ -218,12 +218,17 @@ def ring_size_parameters( ) pix_ring_2 = image[dist_mask_2] - ring_size = np.sum(pix_ring) - size_outside = np.sum(pix_outside_ring * image_mask) - num_pixels_in_ring = np.sum(dist_mask & image_mask) - mean_pixel_outside_ring = np.sum(pix_ring_2) / len(pix_ring_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_size, size_outside, num_pixels_in_ring, mean_pixel_outside_ring + 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): From d1924ad0e63ee92a33e8abba2e31f2463a3a07b1 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Wed, 18 Dec 2024 16:11:20 +0100 Subject: [PATCH 09/12] Update of muon processor with new muon parameters field names --- src/ctapipe/image/muon/processor.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index dc928a57496..ac11a9a368b 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -287,10 +287,10 @@ def _calculate_muon_parameters( ) ( - ring_size, - size_outside, - num_pixels_in_ring, - mean_pixel_outside_ring, + ring_intensity, + intensity_outside_ring, + n_pixels_in_ring, + mean_intensity_outside_ring, ) = ring_size_parameters( ring.radius, ring.center_fov_lon, @@ -316,10 +316,10 @@ def _calculate_muon_parameters( completeness=completeness, intensity_ratio=intensity_ratio, mean_squared_error=mse, - ring_size=ring_size, - size_outside=size_outside, - num_pixels_in_ring=num_pixels_in_ring, - mean_pixel_outside_ring=mean_pixel_outside_ring, + 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, From 68d291b9cc861276702d77375499ebf46e1ba5db Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Thu, 19 Dec 2024 11:42:29 +0100 Subject: [PATCH 10/12] Doc string update for ring size parameters --- src/ctapipe/containers.py | 8 ++++++-- src/ctapipe/image/muon/features.py | 6 ++++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 7619c88c417..ce4cc51fc7c 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1193,9 +1193,13 @@ class MuonParametersContainer(Container): 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(nan, "Number of pixels in the ring") + n_pixels_in_ring = Field( + nan, + "Number of pixels inside the ring integration area that passed the cleaning", + ) mean_intensity_outside_ring = Field( - nan, "Mean intensity of pixels outside the ring" + nan, + "Mean intensity of pixels inside the region limited by ring integration width and outer ring width.", ) standard_dev = Field( nan * u.deg, diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index 1159bbdad94..3c430150191 100644 --- a/src/ctapipe/image/muon/features.py +++ b/src/ctapipe/image/muon/features.py @@ -200,11 +200,13 @@ def ring_size_parameters( ring_intensity: float Sum of the p.e. inside the integration area of the ring intensity_outside_ring: float - Sum of the photons outside the ring integration area that passed the cleaning + 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 ring, but still close to it + 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) From 81fecaaa4c957df59a5b470cd04628529c0a92e3 Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Thu, 19 Dec 2024 14:19:22 +0100 Subject: [PATCH 11/12] Better units conversion in radial light distribution method --- src/ctapipe/image/muon/features.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/image/muon/features.py b/src/ctapipe/image/muon/features.py index 3c430150191..c31e6d461eb 100644 --- a/src/ctapipe/image/muon/features.py +++ b/src/ctapipe/image/muon/features.py @@ -267,11 +267,7 @@ def radial_light_distribution(center_x, center_y, pixel_x, pixel_y, image): if np.sum(image) == 0: return np.nan * u.deg, np.nan, np.nan - x0 = center_x.to_value(u.deg) - y0 = center_y.to_value(u.deg) - pix_x = pixel_x.to_value(u.deg) - pix_y = pixel_y.to_value(u.deg) - pixel_r = np.sqrt((pix_x - x0) ** 2 + (pix_y - y0) ** 2) + pixel_r = np.hypot(pixel_x - center_x, pixel_y - center_y) mean = np.average(pixel_r, weights=image) delta_r = pixel_r - mean @@ -279,4 +275,8 @@ def radial_light_distribution(center_x, center_y, pixel_x, pixel_y, 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 * u.deg, skewness, excess_kurtosis + return ( + standard_dev, + skewness.to_value(u.dimensionless_unscaled), + excess_kurtosis.to_value(u.dimensionless_unscaled), + ) From 6140d3006fc6d0ae6600a5e93e3039164038dcad Mon Sep 17 00:00:00 2001 From: Vadym Voitsekhovskyi Date: Thu, 19 Dec 2024 14:28:36 +0100 Subject: [PATCH 12/12] Resolving problem of unit inconsistency nan vs int --- src/ctapipe/containers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index ce4cc51fc7c..831ca92ebd5 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1194,7 +1194,7 @@ class MuonParametersContainer(Container): ) intensity_outside_ring = Field(nan, "Sum of the pixel charges outside the ring") n_pixels_in_ring = Field( - nan, + -1, "Number of pixels inside the ring integration area that passed the cleaning", ) mean_intensity_outside_ring = Field(