From 47ddd19d4c71b254a477e33fa5604083d7e1ac53 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Wed, 20 Aug 2025 18:05:29 +0200 Subject: [PATCH 01/38] add impactpoint_intensity_fitter.py --- .../muon/impactpoint_intensity_fitter.py | 165 ++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 src/ctapipe/image/muon/impactpoint_intensity_fitter.py diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py new file mode 100644 index 00000000000..a1e8dc49b4f --- /dev/null +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -0,0 +1,165 @@ +""" +Class description to be added. + +""" + +import numpy as np + +from ...containers import MuonEfficiencyContainer +from ...coordinates import TelescopeFrame +from ...core import TelescopeComponent +from ...core.traits import FloatTelescopeParameter, IntTelescopeParameter +from ...exceptions import OptionalDependencyMissing + +try: + from iminuit import Minuit +except ModuleNotFoundError: + Minuit = None + +__all__ = [ + "MuonImpactpointIntensityFitter", +] + + +def chord_length_loss_function(radius, rho, phi): + """ + Function for integrating the length of a chord across a circle (effective chord length). + + A circular mirror is used for signal, and a circular camera is used for shadowing. + + Parameters + ---------- + radius: float or ndarray + radius of circle + rho: float or ndarray + fractional distance of impact point from circle center + phi: float or ndarray in radians + rotation angles to calculate length + + Returns + ------- + float or ndarray: + effective chord length + + References + ---------- + See :cite:p:`vacanti19941`. + Equation 6: for effective chord length calculations inside/outside the ring. + Equation 7: for filtering out non-physical solutions. + + + """ + discriminant_norm = 1 - (rho**2 * np.sin(phi) ** 2) + valid = discriminant_norm >= 0 + + if not valid: + return 0 + + if rho <= 1.0: + # muon has hit the mirror + effective_chord_length = radius * ( + np.sqrt(discriminant_norm) + rho * np.cos(phi) + ) + else: + # muon did not hit the mirror + # Filtering out non-physical solutions for phi + if np.abs(phi) < np.arcsin(1.0 / rho): + effective_chord_length = 2 * radius * np.sqrt(discriminant_norm) + else: + return 0 + + return effective_chord_length + + +class MuonImpactpointIntensityFitter(TelescopeComponent): + """ + Fit muon ring images with a theoretical model to estimate optical efficiency. + + Function for producing the expected image for a given set of trial + muon parameters without using astropy units but expecting the input to + be in the correct ones. + + The image prediction function is currently modeled after :cite:p:`chalmecalvet2013`. + + For more information, also see :cite:p:`muon-review`. + """ + + spe_width = FloatTelescopeParameter( + help="Width of a single photo electron distribution", default_value=0.5 + ).tag(config=True) + + min_lambda_m = FloatTelescopeParameter( + help="Minimum wavelength for Cherenkov light in m", default_value=300e-9 + ).tag(config=True) + + max_lambda_m = FloatTelescopeParameter( + help="Minimum wavelength for Cherenkov light in m", default_value=600e-9 + ).tag(config=True) + + hole_radius_m = FloatTelescopeParameter( + help="The radius of the hole in the center of the primary mirror dish in meters." + "The hole is not circular in shape; however, it can be well approximated as a circle with the same area." + "It is defined with the flat-to-flat distance (LST: 1.51 m, MST: 1.2 m, SST: 0.78 m)." + "We approximate the hexagonal hole with a circle that has the same surface area.", + default_value=[ + ("type", "LST_*", 0.74), + ("type", "MST_*", 0.59), + ("type", "SST_1M_*", 0.38), + ], + ).tag(config=True) + + oversampling = IntTelescopeParameter( + help="Oversampling for the line integration", default_value=3 + ).tag(config=True) + + def __init__(self, subarray, **kwargs): + if Minuit is None: + raise OptionalDependencyMissing("iminuit") from None + + super().__init__(subarray=subarray, **kwargs) + self._geometries_tel_frame = { + tel_id: tel.camera.geometry.transform_to(TelescopeFrame()) + for tel_id, tel in subarray.tel.items() + } + + def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=None): + """ + + Parameters + ---------- + tel_id: int + the telescope id + center_x: Angle quantity + Initial guess for muon ring center in telescope frame + center_y: Angle quantity + Initial guess for muon ring center in telescope frame + radius: Angle quantity + Initial guess for muon ring radius in telescope frame + image: ndarray + Amplitude of image pixels + pedestal: ndarray + Pedestal standard deviation in each pixel + mask: ndarray + mask marking the pixels to be used in the likelihood fit + + Returns + ------- + MuonEfficiencyContainer + """ + telescope = self.subarray.tel[tel_id] + if telescope.optics.n_mirrors != 1: + raise NotImplementedError( + "Currently only single mirror telescopes" + f" are supported in {self.__class__.__name__}" + ) + + return MuonEfficiencyContainer( + impact=None, + impact_x=None, + impact_y=None, + width=None, + optical_efficiency=None, + is_valid=False, + parameters_at_limit=True, + likelihood_value=None, + ) From a3779e5bc49c5b124a86fd6d8d3e0c0e11649c99 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 10:41:50 +0200 Subject: [PATCH 02/38] Remove the old description. New one to be added. --- src/ctapipe/image/muon/impactpoint_intensity_fitter.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index a1e8dc49b4f..e0a873c127b 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -75,13 +75,6 @@ class MuonImpactpointIntensityFitter(TelescopeComponent): """ Fit muon ring images with a theoretical model to estimate optical efficiency. - Function for producing the expected image for a given set of trial - muon parameters without using astropy units but expecting the input to - be in the correct ones. - - The image prediction function is currently modeled after :cite:p:`chalmecalvet2013`. - - For more information, also see :cite:p:`muon-review`. """ spe_width = FloatTelescopeParameter( From d050726bae26fa10b51f90b166de954428ec3ca9 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 10:50:46 +0200 Subject: [PATCH 03/38] Removing unnecessary parameters --- src/ctapipe/image/muon/impactpoint_intensity_fitter.py | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index e0a873c127b..0b60a0d0c57 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -8,7 +8,7 @@ from ...containers import MuonEfficiencyContainer from ...coordinates import TelescopeFrame from ...core import TelescopeComponent -from ...core.traits import FloatTelescopeParameter, IntTelescopeParameter +from ...core.traits import FloatTelescopeParameter from ...exceptions import OptionalDependencyMissing try: @@ -77,10 +77,6 @@ class MuonImpactpointIntensityFitter(TelescopeComponent): """ - spe_width = FloatTelescopeParameter( - help="Width of a single photo electron distribution", default_value=0.5 - ).tag(config=True) - min_lambda_m = FloatTelescopeParameter( help="Minimum wavelength for Cherenkov light in m", default_value=300e-9 ).tag(config=True) @@ -101,10 +97,6 @@ class MuonImpactpointIntensityFitter(TelescopeComponent): ], ).tag(config=True) - oversampling = IntTelescopeParameter( - help="Oversampling for the line integration", default_value=3 - ).tag(config=True) - def __init__(self, subarray, **kwargs): if Minuit is None: raise OptionalDependencyMissing("iminuit") from None From bca80699f61be68a0ea61d0b6a891cd7de988921 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 11:52:26 +0200 Subject: [PATCH 04/38] test: add preliminary dummy test --- .../test_impactpoint_intensity_fitter.py | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py diff --git a/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py new file mode 100644 index 00000000000..84289e58fab --- /dev/null +++ b/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py @@ -0,0 +1,34 @@ +import astropy.units as u +import numpy as np +import pytest + + +def test_dummy(prod5_lst, reference_location): + from ctapipe.image.muon.impactpoint_intensity_fitter import ( + MuonImpactpointIntensityFitter, + ) + from ctapipe.instrument import SubarrayDescription + + pytest.importorskip("iminuit") + + tel_id = 1 + telescope = prod5_lst + subarray = SubarrayDescription( + name="LSTMono", + tel_positions={tel_id: [0, 0, 0] * u.m}, + tel_descriptions={tel_id: telescope}, + reference_location=reference_location, + ) + + fitter = MuonImpactpointIntensityFitter(subarray=subarray) + + fitter( + tel_id=tel_id, + center_x=0 * u.deg, + center_y=2 * u.deg, + radius=1.3 * u.deg, + image=np.zeros(telescope.camera.geometry.n_pixels), + pedestal=np.zeros(telescope.camera.geometry.n_pixels), + ) + + assert 1 From b20b256b0d9852a77fceb863cb1250383fc2eb7e Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 11:53:12 +0200 Subject: [PATCH 05/38] Add MuonImpactpointIntensityFitter to the muon processor --- src/ctapipe/image/muon/processor.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index a96a8ba4849..f47464bf8ea 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -21,7 +21,7 @@ ring_containment, ring_intensity_parameters, ) -from .intensity_fitter import MuonIntensityFitter +from .impactpoint_intensity_fitter import MuonImpactpointIntensityFitter from .ring_fitter import MuonRingFitter INVALID = MuonTelescopeContainer() @@ -131,7 +131,10 @@ def __init__(self, subarray, **kwargs): self.ring_fitter = MuonRingFitter(parent=self) - self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self) + # self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self) + self.intensity_fitter = MuonImpactpointIntensityFitter( + subarray=subarray, parent=self + ) def __call__(self, event: ArrayEventContainer): for tel_id in event.dl1.tel: From 9bfb139a71413f4177fa03deda08a8c2a3a05f79 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 11:54:47 +0200 Subject: [PATCH 06/38] Change the return value to the default MuonEfficiencyContainer for test purposes only. --- .../image/muon/impactpoint_intensity_fitter.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 0b60a0d0c57..d6ee886ab0f 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -138,13 +138,4 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non f" are supported in {self.__class__.__name__}" ) - return MuonEfficiencyContainer( - impact=None, - impact_x=None, - impact_y=None, - width=None, - optical_efficiency=None, - is_valid=False, - parameters_at_limit=True, - likelihood_value=None, - ) + return MuonEfficiencyContainer() From cebf1c24363e291b3468debc5d1e500cffb8e89f Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 11:56:38 +0200 Subject: [PATCH 07/38] change the name LSTMono -> LST --- .../image/muon/tests/test_impactpoint_intensity_fitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py index 84289e58fab..a8747eeedb7 100644 --- a/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py @@ -14,7 +14,7 @@ def test_dummy(prod5_lst, reference_location): tel_id = 1 telescope = prod5_lst subarray = SubarrayDescription( - name="LSTMono", + name="LST", tel_positions={tel_id: [0, 0, 0] * u.m}, tel_descriptions={tel_id: telescope}, reference_location=reference_location, From f4fb6687753cf7a7aa9a62ae7762af61ef2d6fa4 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 17:16:36 +0200 Subject: [PATCH 08/38] Introduce a seed muon_ring_phi_distribution_fit function for estimating the impact point. --- .../muon/impactpoint_intensity_fitter.py | 116 ++++++++++++++++++ 1 file changed, 116 insertions(+) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index d6ee886ab0f..fb90811021a 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -3,6 +3,7 @@ """ +import astropy.units as u import numpy as np from ...containers import MuonEfficiencyContainer @@ -10,6 +11,7 @@ from ...core import TelescopeComponent from ...core.traits import FloatTelescopeParameter from ...exceptions import OptionalDependencyMissing +from ...utils.quantities import all_to_value try: from iminuit import Minuit @@ -71,6 +73,107 @@ def chord_length_loss_function(radius, rho, phi): return effective_chord_length +def muon_ring_phi_distribution_fit( + x, y, mask, image, amplitude_initial=None, rho_initial=None, phi0_initial=None +): + """ + muon_ring_phi_distribution_fit. + + + Parameters + ---------- + x : array-like or astropy.units.Quantity + x-coordinates of the points. + y : array-like or astropy.units.Quantity + y-coordinates of the points. + mask : array-like of bool + Boolean mask indicating which pixels survive the cleaning process. + weights : array-like of float + Weights for the points. If not provided, all points are assigned equal weights (1). + amplitude_initial : unitless float, optional + rho_initial : astropy.units.Quantity, optional + phi0_initial : astropy.units.Quantity, optional + + Returns + ------- + amplitude : astropy.units.Quantity + Fitted amplitude. + rho : astropy.units.Quantity + Fitted rho. + phi0 : astropy.units.Quantity + Fitted phi0. + amplitude_err : astropy.units.Quantity + Fitted radius of the circle error. + rho_err : astropy.units.Quantity + Fitted x-coordinate of the circle center error. + phi0_err : astropy.units.Quantity + Fitted y-coordinate of the circle center error. + + Raises + ------ + OptionalDependencyMissing + If the iminuit package is not installed. + + Notes + ----- + The Taubin circle fit minimizes a specific loss function that balances the + squared residuals of the points from the circle with the weights. This method + is particularly useful for fitting circles to noisy data. + + References + ---------- + - Barcelona_Muons_TPA_final.pdf (slide 6) + """ + + if Minuit is None: + raise OptionalDependencyMissing("iminuit") + + original_unit = x.unit + x, y = all_to_value(x, y, unit=original_unit) + + # x_masked = x[mask] + # y_masked = y[mask] + # image_masked = image[mask] + + # minimization method + # fit = Minuit( + # make_loss_function(x_masked, y_masked, weights_masked), + # xc=xc_initial.to_value(original_unit), + # yc=yc_initial.to_value(original_unit), + # r=r_initial.to_value(original_unit), + # ) + # fit.errordef = Minuit.LEAST_SQUARES + + # set initial parameters uncertainty to a big value + # taubin_error = max_fov * 0.1 + # fit.errors["xc"] = taubin_error + # fit.errors["yc"] = taubin_error + # fit.errors["r"] = taubin_error + + # set wide rage for the minimisation + # fit.limits["xc"] = (-max_fov, max_fov) + # fit.limits["yc"] = (-max_fov, max_fov) + # fit.limits["r"] = (0, max_fov) + + # fit.migrad() + + # radius = Quantity(fit.values["r"], original_unit) + # center_x = Quantity(fit.values["xc"], original_unit) + # center_y = Quantity(fit.values["yc"], original_unit) + # radius_err = Quantity(fit.errors["r"], original_unit) + # center_x_err = Quantity(fit.errors["xc"], original_unit) + # center_y_err = Quantity(fit.errors["yc"], original_unit) + + amplitude = np.nan + rho = np.nan * original_unit + phi0 = np.nan * u.deg + amplitude_err = np.nan + rho_err = np.nan * original_unit + phi0_err = np.nan * u.deg + + return amplitude, rho, phi0, amplitude_err, rho_err, phi0_err + + class MuonImpactpointIntensityFitter(TelescopeComponent): """ Fit muon ring images with a theoretical model to estimate optical efficiency. @@ -138,4 +241,17 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non f" are supported in {self.__class__.__name__}" ) + geometry = telescope.camera.geometry.transform_to(TelescopeFrame()) + + # results_phi_dist = muon_ring_phi_distribution_fit( + muon_ring_phi_distribution_fit( + geometry.pix_x, + geometry.pix_y, + mask, + image, + amplitude_initial=None, + rho_initial=None, + phi0_initial=None, + ) + return MuonEfficiencyContainer() From f67a0c594689098b92ea622ae3e51ff4195bee86 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Thu, 21 Aug 2025 17:37:58 +0200 Subject: [PATCH 09/38] Add staticmethod to count number of calls --- .../image/muon/impactpoint_intensity_fitter.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index fb90811021a..642e0ed8721 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -180,6 +180,8 @@ class MuonImpactpointIntensityFitter(TelescopeComponent): """ + _call_counter = 0 + min_lambda_m = FloatTelescopeParameter( help="Minimum wavelength for Cherenkov light in m", default_value=300e-9 ).tag(config=True) @@ -234,6 +236,9 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non ------- MuonEfficiencyContainer """ + + MuonImpactpointIntensityFitter._call_counter += 1 + telescope = self.subarray.tel[tel_id] if telescope.optics.n_mirrors != 1: raise NotImplementedError( @@ -241,6 +246,11 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non f" are supported in {self.__class__.__name__}" ) + print( + "call_counter_increment = ", + MuonImpactpointIntensityFitter.call_counter_increment(), + ) + geometry = telescope.camera.geometry.transform_to(TelescopeFrame()) # results_phi_dist = muon_ring_phi_distribution_fit( @@ -255,3 +265,7 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non ) return MuonEfficiencyContainer() + + @staticmethod + def call_counter_increment(): + return MuonImpactpointIntensityFitter._call_counter From bbd630b61786b66a8dce7bbcf74e5fd13ff59ff1 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 22 Aug 2025 14:28:00 +0200 Subject: [PATCH 10/38] Add temporary printouts for development (to be removed) --- src/ctapipe/image/muon/processor.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index f47464bf8ea..c9e217a6aa7 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -154,6 +154,8 @@ def _process_telescope_event(self, event, tel_id): event_index = event.index event_id = event_index.event_id + print("-----> START event_id : ", event_id) + if self.subarray.tel[tel_id].optics.n_mirrors != 1: self.log.warning( f"Skipping non-single mirror telescope {tel_id}," @@ -173,10 +175,22 @@ def _process_telescope_event(self, event, tel_id): checks = self.dl1_query(dl1_params=dl1.parameters) + print( + "dl1.parameters.morphology.n_pixels = ", dl1.parameters.morphology.n_pixels + ) + print("dl1.parameters.hillas.intensity = ", dl1.parameters.hillas.intensity) + print("self.dl1_query") + rrrr = 0 if not all(checks): event.muon.tel[tel_id] = INVALID + print("self.dl1_query --> NOT OK") + rrrr = 1 + print("-----> STOP event_id : ", event_id) + print(" ") return + print("rrrr --> ", rrrr) + geometry = self.geometries[tel_id] fov_lon = geometry.pix_x fov_lat = geometry.pix_y @@ -201,6 +215,9 @@ def _process_telescope_event(self, event, tel_id): event.muon.tel[tel_id] = MuonTelescopeContainer( parameters=parameters, ring=ring ) + print("self.ring_query --> NOT OK") + print("-----> STOP event_id : ", event_id) + print(" ") return efficiency = self.intensity_fitter( @@ -223,6 +240,9 @@ def _process_telescope_event(self, event, tel_id): ring=ring, efficiency=efficiency, parameters=parameters ) + print("-----> STOP event_id : ", event_id) + print(" ") + def _calculate_muon_parameters( self, tel_id, image, clean_mask, ring ) -> MuonParametersContainer: From 1493c36d8be01a2f80bfc4ee2678d072320894f9 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sat, 23 Aug 2025 22:58:24 +0200 Subject: [PATCH 11/38] add support function for dev. only. To be removed --- .../muon/impactpoint_intensity_fitter.py | 69 ++++++++++++++++--- 1 file changed, 60 insertions(+), 9 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 642e0ed8721..6c9869c0414 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -6,6 +6,9 @@ import astropy.units as u import numpy as np +# dev to be removed +import pandas as pd + from ...containers import MuonEfficiencyContainer from ...coordinates import TelescopeFrame from ...core import TelescopeComponent @@ -73,8 +76,30 @@ def chord_length_loss_function(radius, rho, phi): return effective_chord_length +def save_histogram_to_csv(hist, csvName): + df = pd.DataFrame( + { + "x": ((np.roll(hist[1], 1) + hist[1]) / 2.0)[1:], + "y": hist[0], + } + ) + + df.to_csv(csvName, sep=" ", index=False) + + return + + def muon_ring_phi_distribution_fit( - x, y, mask, image, amplitude_initial=None, rho_initial=None, phi0_initial=None + x, + y, + mask, + image, + ring_center_x, + ring_center_y, + call_counter, + amplitude_initial=None, + rho_initial=None, + phi0_initial=None, ): """ muon_ring_phi_distribution_fit. @@ -128,12 +153,35 @@ def muon_ring_phi_distribution_fit( if Minuit is None: raise OptionalDependencyMissing("iminuit") - original_unit = x.unit - x, y = all_to_value(x, y, unit=original_unit) - - # x_masked = x[mask] - # y_masked = y[mask] - # image_masked = image[mask] + camera_unit = x.unit + ring_center_unit = ring_center_x.unit + x, y, ring_center_x, ring_center_y = all_to_value( + x, y, ring_center_x, ring_center_y, unit=camera_unit + ) + print("------------------") + print("len(x) = ", len(x)) + print("len(x_masked) = ", len(x[mask])) + print("len(y_masked) = ", len(y[mask])) + print("len(image_masked) = ", len(image[mask])) + print("ring_center_x = ", ring_center_x) + print("ring_center_y = ", ring_center_y) + print("camera_unit = ", camera_unit) + print("ring_center_unit = ", ring_center_unit) + print("++++++++++++++++++") + + hist_phi = np.histogram( + np.arctan2( + y[mask] - ring_center_y, + x[mask] - ring_center_x, + ), + weights=image[mask], + bins=np.linspace(-np.pi, np.pi, 31), + ) + hist_phi_csvName = f"hist_phi_csvName{call_counter}.csv" + save_histogram_to_csv(hist_phi, hist_phi_csvName) + + # print("np.max(phi_masked)/np.pi = ", np.max(phi_masked)/np.pi) + # print("np.min(phi_masked)/np.pi = ", np.min(phi_masked)/np.pi) # minimization method # fit = Minuit( @@ -165,10 +213,10 @@ def muon_ring_phi_distribution_fit( # center_y_err = Quantity(fit.errors["yc"], original_unit) amplitude = np.nan - rho = np.nan * original_unit + rho = np.nan * camera_unit phi0 = np.nan * u.deg amplitude_err = np.nan - rho_err = np.nan * original_unit + rho_err = np.nan * camera_unit phi0_err = np.nan * u.deg return amplitude, rho, phi0, amplitude_err, rho_err, phi0_err @@ -259,6 +307,9 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non geometry.pix_y, mask, image, + center_x, + center_y, + MuonImpactpointIntensityFitter._call_counter, amplitude_initial=None, rho_initial=None, phi0_initial=None, From 041a03a5535a8cbcc41e45911cc36fdbc07ed64a Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sat, 23 Aug 2025 23:00:35 +0200 Subject: [PATCH 12/38] add support function for dev. only. To be removed (processor.py) --- src/ctapipe/image/muon/processor.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index c9e217a6aa7..de78f0e9d6c 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -198,13 +198,21 @@ def _process_telescope_event(self, event, tel_id): # iterative ring fit. # First use cleaning pixels, then only pixels close to the ring # three iterations seems to be enough for most rings + ring_fitter_counter = 0 for _ in range(3): + print("ring_fitter_counter --> ", ring_fitter_counter) + ring_fitter_counter = ring_fitter_counter + 1 ring = self.ring_fitter(geometry, image, mask) dist = np.sqrt( (fov_lon - ring.center_fov_lon) ** 2 + (fov_lat - ring.center_fov_lat) ** 2 ) mask = np.abs(dist - ring.radius) / ring.radius < 0.4 + print("ring_fitter_counter --> np.sum(mask) = ", np.sum(mask)) + print( + "ring_fitter_counter --> np.sum(dl1.image_mask) = ", + np.sum(dl1.image_mask), + ) parameters = self._calculate_muon_parameters( tel_id, image, dl1.image_mask, ring From 552badbc1d8b195787fbe62bf01e0753de33b22e Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sat, 23 Aug 2025 23:58:09 +0200 Subject: [PATCH 13/38] Add event_id only for dev. TO be removed. --- .../muon/impactpoint_intensity_fitter.py | 30 +++++++++++++++---- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 6c9869c0414..01e07411dc6 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -3,6 +3,8 @@ """ +import os + import astropy.units as u import numpy as np @@ -76,9 +78,13 @@ def chord_length_loss_function(radius, rho, phi): return effective_chord_length -def save_histogram_to_csv(hist, csvName): +def save_histogram_to_csv(hist, csvName, event_id): + # print("event_id => ", event_id) + # print(type(event_id)) + # print(len(hist[0])) df = pd.DataFrame( { + "event_id": event_id, "x": ((np.roll(hist[1], 1) + hist[1]) / 2.0)[1:], "y": hist[0], } @@ -97,6 +103,7 @@ def muon_ring_phi_distribution_fit( ring_center_x, ring_center_y, call_counter, + event_id, amplitude_initial=None, rho_initial=None, phi0_initial=None, @@ -177,9 +184,11 @@ def muon_ring_phi_distribution_fit( weights=image[mask], bins=np.linspace(-np.pi, np.pi, 31), ) - hist_phi_csvName = f"hist_phi_csvName{call_counter}.csv" - save_histogram_to_csv(hist_phi, hist_phi_csvName) - + outdir_id = int(call_counter // 1000) + os.makedirs(f"./outdir_{outdir_id}", exist_ok=True) + # os.mkdir(f"./outdir_{outdir_id}",exist_ok=True) + hist_phi_csvName = f"./outdir_{outdir_id}/hist_phi_csvName{call_counter}.csv" + save_histogram_to_csv(hist_phi, hist_phi_csvName, event_id) # print("np.max(phi_masked)/np.pi = ", np.max(phi_masked)/np.pi) # print("np.min(phi_masked)/np.pi = ", np.min(phi_masked)/np.pi) @@ -260,7 +269,17 @@ def __init__(self, subarray, **kwargs): for tel_id, tel in subarray.tel.items() } - def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=None): + def __call__( + self, + tel_id, + center_x, + center_y, + radius, + image, + pedestal, + mask=None, + event_id=None, + ): """ Parameters @@ -310,6 +329,7 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non center_x, center_y, MuonImpactpointIntensityFitter._call_counter, + event_id, amplitude_initial=None, rho_initial=None, phi0_initial=None, From ca0bad5e1478c1d901f9677aa90b10b723bd2315 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sat, 23 Aug 2025 23:58:45 +0200 Subject: [PATCH 14/38] Add event_id only for dev. To be removed. --- src/ctapipe/image/muon/processor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index de78f0e9d6c..4672ee777d9 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -234,8 +234,9 @@ def _process_telescope_event(self, event, tel_id): ring.center_fov_lat, ring.radius, image, - mask=mask, pedestal=np.full(mask.shape, self.pedestal.tel[tel_id]), + mask=mask, + event_id=event_id, ) self.log.debug( From 64f4b8d4711a306f2e60bc5396bbf7cab6770061 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Tue, 26 Aug 2025 13:07:02 +0200 Subject: [PATCH 15/38] Add smoothing phi dist --- .../image/muon/impactpoint_intensity_fitter.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 01e07411dc6..c2fb9704bda 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -10,6 +10,7 @@ # dev to be removed import pandas as pd +from scipy.ndimage import correlate1d from ...containers import MuonEfficiencyContainer from ...coordinates import TelescopeFrame @@ -78,7 +79,7 @@ def chord_length_loss_function(radius, rho, phi): return effective_chord_length -def save_histogram_to_csv(hist, csvName, event_id): +def save_histogram_to_csv(hist, csvName, event_id, hist_phi_smooth): # print("event_id => ", event_id) # print(type(event_id)) # print(len(hist[0])) @@ -86,7 +87,7 @@ def save_histogram_to_csv(hist, csvName, event_id): { "event_id": event_id, "x": ((np.roll(hist[1], 1) + hist[1]) / 2.0)[1:], - "y": hist[0], + "y": hist_phi_smooth, } ) @@ -182,13 +183,15 @@ def muon_ring_phi_distribution_fit( x[mask] - ring_center_x, ), weights=image[mask], - bins=np.linspace(-np.pi, np.pi, 31), + bins=np.linspace(-np.pi, np.pi, 19), ) + hist_phi_smooth = (correlate1d(hist_phi[0], np.ones(2), mode="wrap", axis=0) / 2,) + outdir_id = int(call_counter // 1000) os.makedirs(f"./outdir_{outdir_id}", exist_ok=True) # os.mkdir(f"./outdir_{outdir_id}",exist_ok=True) hist_phi_csvName = f"./outdir_{outdir_id}/hist_phi_csvName{call_counter}.csv" - save_histogram_to_csv(hist_phi, hist_phi_csvName, event_id) + save_histogram_to_csv(hist_phi, hist_phi_csvName, event_id, hist_phi_smooth[0]) # print("np.max(phi_masked)/np.pi = ", np.max(phi_masked)/np.pi) # print("np.min(phi_masked)/np.pi = ", np.min(phi_masked)/np.pi) From afb162f9b7e8c896d0cb6583332df6e132848337 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 29 Aug 2025 12:28:15 +0200 Subject: [PATCH 16/38] modify chord_length with modulo --- .../muon/impactpoint_intensity_fitter.py | 99 ++++++++++++++----- 1 file changed, 73 insertions(+), 26 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index c2fb9704bda..b5a17be44e3 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -29,7 +29,7 @@ ] -def chord_length_loss_function(radius, rho, phi): +def chord_length(radius, rho, phi0, phi): """ Function for integrating the length of a chord across a circle (effective chord length). @@ -37,10 +37,10 @@ def chord_length_loss_function(radius, rho, phi): Parameters ---------- - radius: float or ndarray + radius: float radius of circle - rho: float or ndarray - fractional distance of impact point from circle center + rho: float + distance of impact point from circle center phi: float or ndarray in radians rotation angles to calculate length @@ -57,26 +57,33 @@ def chord_length_loss_function(radius, rho, phi): """ - discriminant_norm = 1 - (rho**2 * np.sin(phi) ** 2) - valid = discriminant_norm >= 0 - if not valid: - return 0 + if phi0 == 0: + if radius <= 0: + return np.zero(len(phi)) - if rho <= 1.0: - # muon has hit the mirror - effective_chord_length = radius * ( - np.sqrt(discriminant_norm) + rho * np.cos(phi) - ) - else: - # muon did not hit the mirror - # Filtering out non-physical solutions for phi - if np.abs(phi) < np.arcsin(1.0 / rho): - effective_chord_length = 2 * radius * np.sqrt(discriminant_norm) + phi_modulo = ((phi + np.pi) % (2 * np.pi) - np.pi) * np.where(phi < -1, -1, 1) + + rho_R = rho / radius + discriminant_norm = 1 - (rho_R**2 * np.sin(phi_modulo) ** 2) + + if rho_R <= 1.0: + # muon has hit the mirror + effective_chord_length = radius * ( + np.sqrt(discriminant_norm) + rho_R * np.cos(phi_modulo) + ) else: - return 0 + # muon did not hit the mirror + discriminant_norm[discriminant_norm < 0] = 0 + effective_chord_length = 2 * radius * np.sqrt(discriminant_norm) + # Filtering out non-physical solutions for phi + effective_chord_length *= np.where( + (np.abs(phi_modulo) < np.arcsin(1.0 / rho_R)), 1, 0 + ) + + return effective_chord_length - return effective_chord_length + return chord_length(radius, rho, 0, phi - phi0) def save_histogram_to_csv(hist, csvName, event_id, hist_phi_smooth): @@ -96,7 +103,7 @@ def save_histogram_to_csv(hist, csvName, event_id, hist_phi_smooth): return -def muon_ring_phi_distribution_fit( +def fit_muon_ring_phi_distribution( x, y, mask, @@ -177,27 +184,45 @@ def muon_ring_phi_distribution_fit( print("ring_center_unit = ", ring_center_unit) print("++++++++++++++++++") + n_phi_bins = 12 + n_of_smoothing_points = 1 # 1 --> no smoothing + hist_phi = np.histogram( np.arctan2( y[mask] - ring_center_y, x[mask] - ring_center_x, ), weights=image[mask], - bins=np.linspace(-np.pi, np.pi, 19), + bins=np.linspace(-np.pi, np.pi, n_phi_bins + 1), ) - hist_phi_smooth = (correlate1d(hist_phi[0], np.ones(2), mode="wrap", axis=0) / 2,) + + phi_y = ( + correlate1d(hist_phi[0], np.ones(n_of_smoothing_points), mode="wrap", axis=0) + / n_of_smoothing_points, + )[0] + phi_y_err = np.sqrt(np.abs(phi_y)) + + # weights = phi_y[phi_y > 0].astype(int) + + phi_x = ((np.roll(hist_phi[1], 1) + hist_phi[1]) / 2.0)[1:] + phi_x_err = np.ones(len(phi_x)) * np.pi / n_phi_bins + + print(phi_y) + print(phi_y_err) + print(phi_x) + print(phi_x_err) outdir_id = int(call_counter // 1000) os.makedirs(f"./outdir_{outdir_id}", exist_ok=True) # os.mkdir(f"./outdir_{outdir_id}",exist_ok=True) hist_phi_csvName = f"./outdir_{outdir_id}/hist_phi_csvName{call_counter}.csv" - save_histogram_to_csv(hist_phi, hist_phi_csvName, event_id, hist_phi_smooth[0]) + save_histogram_to_csv(hist_phi, hist_phi_csvName, event_id, phi_y) # print("np.max(phi_masked)/np.pi = ", np.max(phi_masked)/np.pi) # print("np.min(phi_masked)/np.pi = ", np.min(phi_masked)/np.pi) # minimization method # fit = Minuit( - # make_loss_function(x_masked, y_masked, weights_masked), + # phi_dist_loss_function(phi_x, phi_y, phi_x_err, phi_y_err, weights), # xc=xc_initial.to_value(original_unit), # yc=yc_initial.to_value(original_unit), # r=r_initial.to_value(original_unit), @@ -234,6 +259,28 @@ def muon_ring_phi_distribution_fit( return amplitude, rho, phi0, amplitude_err, rho_err, phi0_err +def phi_dist_loss_function(x, y, xe, ye, w): + """dist_loss_function + + x, y, xe, ye: positions of pixels surviving the cleaning + should not be quantities + w : array-like of float, weights for the points + """ + + def taubin_loss_function(xc, yc, r): + """taubin fit formula + reference : Barcelona_Muons_TPA_final.pdf (slide 6) + """ + + distance_squared = (x - xc) ** 2 + (y - yc) ** 2 + upper_term = ((w * (distance_squared - r**2)) ** 2).sum() + lower_term = (w * distance_squared).sum() + + return np.abs(upper_term) / np.abs(lower_term) + + return taubin_loss_function + + class MuonImpactpointIntensityFitter(TelescopeComponent): """ Fit muon ring images with a theoretical model to estimate optical efficiency. @@ -324,7 +371,7 @@ def __call__( geometry = telescope.camera.geometry.transform_to(TelescopeFrame()) # results_phi_dist = muon_ring_phi_distribution_fit( - muon_ring_phi_distribution_fit( + fit_muon_ring_phi_distribution( geometry.pix_x, geometry.pix_y, mask, From c19ef6e4d29bb3e6d2f0f3a6242d9b479831b229 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 29 Aug 2025 14:07:31 +0200 Subject: [PATCH 17/38] Add the loss function for the phi-distribution of the Cherenkov light. --- .../muon/impactpoint_intensity_fitter.py | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index b5a17be44e3..233305bd4a9 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -259,26 +259,23 @@ def fit_muon_ring_phi_distribution( return amplitude, rho, phi0, amplitude_err, rho_err, phi0_err -def phi_dist_loss_function(x, y, xe, ye, w): +def phi_dist_loss_function(x, y, err, w): """dist_loss_function - x, y, xe, ye: positions of pixels surviving the cleaning + x, y, err: positions of pixels surviving the cleaning should not be quantities w : array-like of float, weights for the points - """ - def taubin_loss_function(xc, yc, r): - """taubin fit formula - reference : Barcelona_Muons_TPA_final.pdf (slide 6) - """ + """ - distance_squared = (x - xc) ** 2 + (y - yc) ** 2 - upper_term = ((w * (distance_squared - r**2)) ** 2).sum() - lower_term = (w * distance_squared).sum() + def loss_function(amplitude, R_mirror, R_camera, rho, phi0): + signal = amplitude * chord_length(R_mirror, rho, phi0, x) + shadow = amplitude * chord_length(R_camera, rho, phi0, x) + diff_squared = ((signal - shadow - y) * w / err) ** 2 - return np.abs(upper_term) / np.abs(lower_term) + return diff_squared.sum() - return taubin_loss_function + return loss_function class MuonImpactpointIntensityFitter(TelescopeComponent): From 2a771f7cbe61a99801342717a0f74b8cefd4d093 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 29 Aug 2025 16:19:01 +0200 Subject: [PATCH 18/38] Add working fit with Minuit. --- .../muon/impactpoint_intensity_fitter.py | 80 ++++++++++++++----- 1 file changed, 58 insertions(+), 22 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 233305bd4a9..0b70a035b9c 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -10,6 +10,7 @@ # dev to be removed import pandas as pd +from astropy.units import Quantity from scipy.ndimage import correlate1d from ...containers import MuonEfficiencyContainer @@ -110,11 +111,10 @@ def fit_muon_ring_phi_distribution( image, ring_center_x, ring_center_y, + optics, + shadow_radius, call_counter, event_id, - amplitude_initial=None, - rho_initial=None, - phi0_initial=None, ): """ muon_ring_phi_distribution_fit. @@ -202,15 +202,22 @@ def fit_muon_ring_phi_distribution( )[0] phi_y_err = np.sqrt(np.abs(phi_y)) - # weights = phi_y[phi_y > 0].astype(int) + weights = (phi_y > 0).astype(int) phi_x = ((np.roll(hist_phi[1], 1) + hist_phi[1]) / 2.0)[1:] phi_x_err = np.ones(len(phi_x)) * np.pi / n_phi_bins + phi_err = np.sqrt(phi_x_err**2 + phi_y_err**2) + print(phi_y) print(phi_y_err) print(phi_x) print(phi_x_err) + print(phi_err) + print(weights) + print(phi_y) + print("shadow_radius = ", shadow_radius) + print("type(shadow_radius) = ", type(shadow_radius)) outdir_id = int(call_counter // 1000) os.makedirs(f"./outdir_{outdir_id}", exist_ok=True) @@ -220,27 +227,57 @@ def fit_muon_ring_phi_distribution( # print("np.max(phi_masked)/np.pi = ", np.max(phi_masked)/np.pi) # print("np.min(phi_masked)/np.pi = ", np.min(phi_masked)/np.pi) + # amplitude_initial=None, + # rho_initial=None, + # phi0_initial=None, + # minimization method - # fit = Minuit( - # phi_dist_loss_function(phi_x, phi_y, phi_x_err, phi_y_err, weights), - # xc=xc_initial.to_value(original_unit), - # yc=yc_initial.to_value(original_unit), - # r=r_initial.to_value(original_unit), - # ) - # fit.errordef = Minuit.LEAST_SQUARES + fit = Minuit( + phi_dist_loss_function(phi_x, phi_y, phi_err, weights), + amplitude=12, + R_mirror=np.sqrt(optics.mirror_area.to_value(u.m**2) / np.pi), + R_shadow=shadow_radius.to_value(u.m), + rho=5.0, + phi0=0.1, + ) + fit.errordef = Minuit.LEAST_SQUARES + + # + fit.fixed["R_mirror"] = True + fit.fixed["R_shadow"] = True # set initial parameters uncertainty to a big value # taubin_error = max_fov * 0.1 - # fit.errors["xc"] = taubin_error - # fit.errors["yc"] = taubin_error - # fit.errors["r"] = taubin_error + fit.errors["amplitude"] = 10 + fit.errors["R_mirror"] = 0.0001 + fit.errors["R_shadow"] = 0.0001 + fit.errors["rho"] = 10.0 + fit.errors["phi0"] = np.pi # set wide rage for the minimisation # fit.limits["xc"] = (-max_fov, max_fov) # fit.limits["yc"] = (-max_fov, max_fov) # fit.limits["r"] = (0, max_fov) - # fit.migrad() + fit.migrad() + + amplitude = fit.values["amplitude"] + R_mirror = Quantity(fit.values["R_mirror"], u.m) + R_shadow = Quantity(fit.values["R_shadow"], u.m) + rho = Quantity(fit.values["rho"], u.m) + phi0 = Quantity(fit.values["phi0"], u.rad) + + amplitude_err = fit.errors["amplitude"] + # R_mirror_err = Quantity(fit.errors["R_mirror"], u.m) + # R_shadow_err = Quantity(fit.errors["R_shadow"], u.m) + rho_err = Quantity(fit.errors["rho"], u.m) + phi0_err = Quantity(fit.errors["phi0"], u.rad) + + print("amplitude = ", amplitude) + print("R_mirror = ", R_mirror) + print("R_shadow = ", R_shadow) + print("rho = ", rho) + print("phi0 = ", phi0) # radius = Quantity(fit.values["r"], original_unit) # center_x = Quantity(fit.values["xc"], original_unit) @@ -268,9 +305,9 @@ def phi_dist_loss_function(x, y, err, w): """ - def loss_function(amplitude, R_mirror, R_camera, rho, phi0): + def loss_function(amplitude, R_mirror, R_shadow, rho, phi0): signal = amplitude * chord_length(R_mirror, rho, phi0, x) - shadow = amplitude * chord_length(R_camera, rho, phi0, x) + shadow = amplitude * chord_length(R_shadow, rho, phi0, x) diff_squared = ((signal - shadow - y) * w / err) ** 2 return diff_squared.sum() @@ -375,11 +412,10 @@ def __call__( image, center_x, center_y, - MuonImpactpointIntensityFitter._call_counter, - event_id, - amplitude_initial=None, - rho_initial=None, - phi0_initial=None, + optics=telescope.optics, + shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, + call_counter=MuonImpactpointIntensityFitter._call_counter, + event_id=event_id, ) return MuonEfficiencyContainer() From 3bef25e98d52fb46b9fa56719710e41f6c9f45b5 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 29 Aug 2025 16:35:43 +0200 Subject: [PATCH 19/38] Remove debugging printouts from processor.py --- src/ctapipe/image/muon/processor.py | 29 ----------------------------- 1 file changed, 29 deletions(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index 4672ee777d9..ba78bc97e02 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -154,8 +154,6 @@ def _process_telescope_event(self, event, tel_id): event_index = event.index event_id = event_index.event_id - print("-----> START event_id : ", event_id) - if self.subarray.tel[tel_id].optics.n_mirrors != 1: self.log.warning( f"Skipping non-single mirror telescope {tel_id}," @@ -175,22 +173,10 @@ def _process_telescope_event(self, event, tel_id): checks = self.dl1_query(dl1_params=dl1.parameters) - print( - "dl1.parameters.morphology.n_pixels = ", dl1.parameters.morphology.n_pixels - ) - print("dl1.parameters.hillas.intensity = ", dl1.parameters.hillas.intensity) - print("self.dl1_query") - rrrr = 0 if not all(checks): event.muon.tel[tel_id] = INVALID - print("self.dl1_query --> NOT OK") - rrrr = 1 - print("-----> STOP event_id : ", event_id) - print(" ") return - print("rrrr --> ", rrrr) - geometry = self.geometries[tel_id] fov_lon = geometry.pix_x fov_lat = geometry.pix_y @@ -198,21 +184,13 @@ def _process_telescope_event(self, event, tel_id): # iterative ring fit. # First use cleaning pixels, then only pixels close to the ring # three iterations seems to be enough for most rings - ring_fitter_counter = 0 for _ in range(3): - print("ring_fitter_counter --> ", ring_fitter_counter) - ring_fitter_counter = ring_fitter_counter + 1 ring = self.ring_fitter(geometry, image, mask) dist = np.sqrt( (fov_lon - ring.center_fov_lon) ** 2 + (fov_lat - ring.center_fov_lat) ** 2 ) mask = np.abs(dist - ring.radius) / ring.radius < 0.4 - print("ring_fitter_counter --> np.sum(mask) = ", np.sum(mask)) - print( - "ring_fitter_counter --> np.sum(dl1.image_mask) = ", - np.sum(dl1.image_mask), - ) parameters = self._calculate_muon_parameters( tel_id, image, dl1.image_mask, ring @@ -223,9 +201,6 @@ def _process_telescope_event(self, event, tel_id): event.muon.tel[tel_id] = MuonTelescopeContainer( parameters=parameters, ring=ring ) - print("self.ring_query --> NOT OK") - print("-----> STOP event_id : ", event_id) - print(" ") return efficiency = self.intensity_fitter( @@ -236,7 +211,6 @@ def _process_telescope_event(self, event, tel_id): image, pedestal=np.full(mask.shape, self.pedestal.tel[tel_id]), mask=mask, - event_id=event_id, ) self.log.debug( @@ -249,9 +223,6 @@ def _process_telescope_event(self, event, tel_id): ring=ring, efficiency=efficiency, parameters=parameters ) - print("-----> STOP event_id : ", event_id) - print(" ") - def _calculate_muon_parameters( self, tel_id, image, clean_mask, ring ) -> MuonParametersContainer: From 4dd880212e777973cc57bf3cd05519d4b9d1dcce Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 29 Aug 2025 17:06:11 +0200 Subject: [PATCH 20/38] Remove debugging printouts from impactpoint_intensity_fitter.py --- .../muon/impactpoint_intensity_fitter.py | 91 ++----------------- 1 file changed, 10 insertions(+), 81 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 0b70a035b9c..ee9fd563ab1 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -3,7 +3,6 @@ """ -import os import astropy.units as u import numpy as np @@ -67,6 +66,7 @@ def chord_length(radius, rho, phi0, phi): rho_R = rho / radius discriminant_norm = 1 - (rho_R**2 * np.sin(phi_modulo) ** 2) + discriminant_norm[discriminant_norm < 0] = 0.0 if rho_R <= 1.0: # muon has hit the mirror @@ -75,7 +75,6 @@ def chord_length(radius, rho, phi0, phi): ) else: # muon did not hit the mirror - discriminant_norm[discriminant_norm < 0] = 0 effective_chord_length = 2 * radius * np.sqrt(discriminant_norm) # Filtering out non-physical solutions for phi effective_chord_length *= np.where( @@ -113,8 +112,6 @@ def fit_muon_ring_phi_distribution( ring_center_y, optics, shadow_radius, - call_counter, - event_id, ): """ muon_ring_phi_distribution_fit. @@ -169,20 +166,10 @@ def fit_muon_ring_phi_distribution( raise OptionalDependencyMissing("iminuit") camera_unit = x.unit - ring_center_unit = ring_center_x.unit + # ring_center_unit = ring_center_x.unit x, y, ring_center_x, ring_center_y = all_to_value( x, y, ring_center_x, ring_center_y, unit=camera_unit ) - print("------------------") - print("len(x) = ", len(x)) - print("len(x_masked) = ", len(x[mask])) - print("len(y_masked) = ", len(y[mask])) - print("len(image_masked) = ", len(image[mask])) - print("ring_center_x = ", ring_center_x) - print("ring_center_y = ", ring_center_y) - print("camera_unit = ", camera_unit) - print("ring_center_unit = ", ring_center_unit) - print("++++++++++++++++++") n_phi_bins = 12 n_of_smoothing_points = 1 # 1 --> no smoothing @@ -209,28 +196,6 @@ def fit_muon_ring_phi_distribution( phi_err = np.sqrt(phi_x_err**2 + phi_y_err**2) - print(phi_y) - print(phi_y_err) - print(phi_x) - print(phi_x_err) - print(phi_err) - print(weights) - print(phi_y) - print("shadow_radius = ", shadow_radius) - print("type(shadow_radius) = ", type(shadow_radius)) - - outdir_id = int(call_counter // 1000) - os.makedirs(f"./outdir_{outdir_id}", exist_ok=True) - # os.mkdir(f"./outdir_{outdir_id}",exist_ok=True) - hist_phi_csvName = f"./outdir_{outdir_id}/hist_phi_csvName{call_counter}.csv" - save_histogram_to_csv(hist_phi, hist_phi_csvName, event_id, phi_y) - # print("np.max(phi_masked)/np.pi = ", np.max(phi_masked)/np.pi) - # print("np.min(phi_masked)/np.pi = ", np.min(phi_masked)/np.pi) - - # amplitude_initial=None, - # rho_initial=None, - # phi0_initial=None, - # minimization method fit = Minuit( phi_dist_loss_function(phi_x, phi_y, phi_err, weights), @@ -254,44 +219,24 @@ def fit_muon_ring_phi_distribution( fit.errors["rho"] = 10.0 fit.errors["phi0"] = np.pi - # set wide rage for the minimisation - # fit.limits["xc"] = (-max_fov, max_fov) - # fit.limits["yc"] = (-max_fov, max_fov) - # fit.limits["r"] = (0, max_fov) - fit.migrad() amplitude = fit.values["amplitude"] - R_mirror = Quantity(fit.values["R_mirror"], u.m) - R_shadow = Quantity(fit.values["R_shadow"], u.m) + # R_mirror = Quantity(fit.values["R_mirror"], u.m) + # R_shadow = Quantity(fit.values["R_shadow"], u.m) rho = Quantity(fit.values["rho"], u.m) phi0 = Quantity(fit.values["phi0"], u.rad) amplitude_err = fit.errors["amplitude"] - # R_mirror_err = Quantity(fit.errors["R_mirror"], u.m) - # R_shadow_err = Quantity(fit.errors["R_shadow"], u.m) rho_err = Quantity(fit.errors["rho"], u.m) phi0_err = Quantity(fit.errors["phi0"], u.rad) - print("amplitude = ", amplitude) - print("R_mirror = ", R_mirror) - print("R_shadow = ", R_shadow) - print("rho = ", rho) - print("phi0 = ", phi0) - - # radius = Quantity(fit.values["r"], original_unit) - # center_x = Quantity(fit.values["xc"], original_unit) - # center_y = Quantity(fit.values["yc"], original_unit) - # radius_err = Quantity(fit.errors["r"], original_unit) - # center_x_err = Quantity(fit.errors["xc"], original_unit) - # center_y_err = Quantity(fit.errors["yc"], original_unit) - - amplitude = np.nan - rho = np.nan * camera_unit - phi0 = np.nan * u.deg - amplitude_err = np.nan - rho_err = np.nan * camera_unit - phi0_err = np.nan * u.deg + # amplitude = np.nan + # rho = np.nan * camera_unit + # phi0 = np.nan * u.deg + # amplitude_err = np.nan + # rho_err = np.nan * camera_unit + # phi0_err = np.nan * u.deg return amplitude, rho, phi0, amplitude_err, rho_err, phi0_err @@ -321,8 +266,6 @@ class MuonImpactpointIntensityFitter(TelescopeComponent): """ - _call_counter = 0 - min_lambda_m = FloatTelescopeParameter( help="Minimum wavelength for Cherenkov light in m", default_value=300e-9 ).tag(config=True) @@ -362,7 +305,6 @@ def __call__( image, pedestal, mask=None, - event_id=None, ): """ @@ -388,8 +330,6 @@ def __call__( MuonEfficiencyContainer """ - MuonImpactpointIntensityFitter._call_counter += 1 - telescope = self.subarray.tel[tel_id] if telescope.optics.n_mirrors != 1: raise NotImplementedError( @@ -397,11 +337,6 @@ def __call__( f" are supported in {self.__class__.__name__}" ) - print( - "call_counter_increment = ", - MuonImpactpointIntensityFitter.call_counter_increment(), - ) - geometry = telescope.camera.geometry.transform_to(TelescopeFrame()) # results_phi_dist = muon_ring_phi_distribution_fit( @@ -414,12 +349,6 @@ def __call__( center_y, optics=telescope.optics, shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, - call_counter=MuonImpactpointIntensityFitter._call_counter, - event_id=event_id, ) return MuonEfficiencyContainer() - - @staticmethod - def call_counter_increment(): - return MuonImpactpointIntensityFitter._call_counter From 796140b07b3cd261b1e160648827ac37351a5e79 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 1 Sep 2025 16:23:40 +0200 Subject: [PATCH 21/38] Cleanup code, add units, and include MuonEfficiencyContainer output --- .../muon/impactpoint_intensity_fitter.py | 44 ++++++++++--------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index ee9fd563ab1..cb575bb1201 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -166,7 +166,7 @@ def fit_muon_ring_phi_distribution( raise OptionalDependencyMissing("iminuit") camera_unit = x.unit - # ring_center_unit = ring_center_x.unit + size_impact_point_unit = u.m x, y, ring_center_x, ring_center_y = all_to_value( x, y, ring_center_x, ring_center_y, unit=camera_unit ) @@ -221,24 +221,16 @@ def fit_muon_ring_phi_distribution( fit.migrad() - amplitude = fit.values["amplitude"] - # R_mirror = Quantity(fit.values["R_mirror"], u.m) - # R_shadow = Quantity(fit.values["R_shadow"], u.m) - rho = Quantity(fit.values["rho"], u.m) - phi0 = Quantity(fit.values["phi0"], u.rad) + rho = Quantity(fit.values["rho"], size_impact_point_unit) - amplitude_err = fit.errors["amplitude"] - rho_err = Quantity(fit.errors["rho"], u.m) - phi0_err = Quantity(fit.errors["phi0"], u.rad) - - # amplitude = np.nan - # rho = np.nan * camera_unit - # phi0 = np.nan * u.deg - # amplitude_err = np.nan - # rho_err = np.nan * camera_unit - # phi0_err = np.nan * u.deg - - return amplitude, rho, phi0, amplitude_err, rho_err, phi0_err + return MuonEfficiencyContainer( + impact=rho, + impact_x=rho * np.cos(Quantity(fit.values["phi0"], camera_unit)), + impact_y=rho * np.sin(Quantity(fit.values["phi0"], camera_unit)), + is_valid=fit.valid, + parameters_at_limit=fit.fmin.has_parameters_at_limit, + likelihood_value=fit.fval, + ) def phi_dist_loss_function(x, y, err, w): @@ -339,8 +331,7 @@ def __call__( geometry = telescope.camera.geometry.transform_to(TelescopeFrame()) - # results_phi_dist = muon_ring_phi_distribution_fit( - fit_muon_ring_phi_distribution( + mu_eff_container = fit_muon_ring_phi_distribution( geometry.pix_x, geometry.pix_y, mask, @@ -351,4 +342,15 @@ def __call__( shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, ) - return MuonEfficiencyContainer() + # return MuonEfficiencyContainer( + # impact=result["impact_parameter"] * u.m, + # impact_x=result["impact_parameter"] * np.cos(result["phi"]) * u.m, + # impact_y=result["impact_parameter"] * np.sin(result["phi"]) * u.m, + # width=u.Quantity(np.rad2deg(result["ring_width"]), u.deg), + # optical_efficiency=result["optical_efficiency_muon"], + # is_valid=minuit.valid, + # parameters_at_limit=minuit.fmin.has_parameters_at_limit, + # likelihood_value=minuit.fval, + # ) + + return mu_eff_container From 2362269f85e58939cc99a040649192caf4d5372f Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 12 Sep 2025 17:11:08 +0200 Subject: [PATCH 22/38] Cleaning the code and adding calculations for the initial parameters. --- .../muon/impactpoint_intensity_fitter.py | 61 ++++++++----------- 1 file changed, 24 insertions(+), 37 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index cb575bb1201..765c5fddeba 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -114,7 +114,7 @@ def fit_muon_ring_phi_distribution( shadow_radius, ): """ - muon_ring_phi_distribution_fit. + fit_muon_ring_phi_distribution Parameters @@ -125,11 +125,7 @@ def fit_muon_ring_phi_distribution( y-coordinates of the points. mask : array-like of bool Boolean mask indicating which pixels survive the cleaning process. - weights : array-like of float - Weights for the points. If not provided, all points are assigned equal weights (1). - amplitude_initial : unitless float, optional - rho_initial : astropy.units.Quantity, optional - phi0_initial : astropy.units.Quantity, optional + image : array-like of float Returns ------- @@ -146,20 +142,6 @@ def fit_muon_ring_phi_distribution( phi0_err : astropy.units.Quantity Fitted y-coordinate of the circle center error. - Raises - ------ - OptionalDependencyMissing - If the iminuit package is not installed. - - Notes - ----- - The Taubin circle fit minimizes a specific loss function that balances the - squared residuals of the points from the circle with the weights. This method - is particularly useful for fitting circles to noisy data. - - References - ---------- - - Barcelona_Muons_TPA_final.pdf (slide 6) """ if Minuit is None: @@ -171,7 +153,7 @@ def fit_muon_ring_phi_distribution( x, y, ring_center_x, ring_center_y, unit=camera_unit ) - n_phi_bins = 12 + n_phi_bins = 12 # to be added to configure file as input perameters n_of_smoothing_points = 1 # 1 --> no smoothing hist_phi = np.histogram( @@ -196,21 +178,37 @@ def fit_muon_ring_phi_distribution( phi_err = np.sqrt(phi_x_err**2 + phi_y_err**2) + total_integral = np.sum(phi_y) + + if total_integral <= 0.0: + return MuonEfficiencyContainer() + + amplitude_initial = np.nan + rho_initial = np.nan + phi0_initial = np.nan + + amplitude_initial = total_integral / 110.0 + rho_initial = 2 * (np.max(phi_y) - np.min(phi_y)) / total_integral * 110.0 + phi0_initial = phi_x[np.argmax(phi_y)] + + amplitude_initial = 12 if np.isnan(amplitude_initial) else amplitude_initial + rho_initial = 8 if np.isnan(rho_initial) is np.nan else rho_initial + phi0_initial = 0 if np.isnan(phi0_initial) is np.nan else phi0_initial + # minimization method fit = Minuit( phi_dist_loss_function(phi_x, phi_y, phi_err, weights), - amplitude=12, + amplitude=amplitude_initial, R_mirror=np.sqrt(optics.mirror_area.to_value(u.m**2) / np.pi), R_shadow=shadow_radius.to_value(u.m), - rho=5.0, - phi0=0.1, + rho=rho_initial, + phi0=phi0_initial, ) fit.errordef = Minuit.LEAST_SQUARES - # fit.fixed["R_mirror"] = True fit.fixed["R_shadow"] = True - + # # set initial parameters uncertainty to a big value # taubin_error = max_fov * 0.1 fit.errors["amplitude"] = 10 @@ -342,15 +340,4 @@ def __call__( shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, ) - # return MuonEfficiencyContainer( - # impact=result["impact_parameter"] * u.m, - # impact_x=result["impact_parameter"] * np.cos(result["phi"]) * u.m, - # impact_y=result["impact_parameter"] * np.sin(result["phi"]) * u.m, - # width=u.Quantity(np.rad2deg(result["ring_width"]), u.deg), - # optical_efficiency=result["optical_efficiency_muon"], - # is_valid=minuit.valid, - # parameters_at_limit=minuit.fmin.has_parameters_at_limit, - # likelihood_value=minuit.fval, - # ) - return mu_eff_container From b7423361e166f0076400394a2202392f95a18987 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 12 Sep 2025 17:49:56 +0200 Subject: [PATCH 23/38] add compute_muon_ring_width and comput_absolute_optical_efficiency_from_muon_ring sceletons --- .../muon/impactpoint_intensity_fitter.py | 53 +++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 765c5fddeba..277c3c9a5be 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -103,6 +103,47 @@ def save_histogram_to_csv(hist, csvName, event_id, hist_phi_smooth): return +def compute_muon_ring_width( + x, + y, + mask, + image, + ring_center_x, + ring_center_y, +): + """ + compute_muon_ring_width + + + Parameters + ---------- + + Returns + ------- + + + """ + + return 0.1 * x.unit + + +def comput_absolute_optical_efficiency_from_muon_ring(): + """ + comput_absolute_optical_efficiency_from_muon_ring + + + Parameters + ---------- + + Returns + ------- + + + """ + + return 0.2 + + def fit_muon_ring_phi_distribution( x, y, @@ -340,4 +381,16 @@ def __call__( shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, ) + mu_eff_container.width = compute_muon_ring_width( + geometry.pix_x, + geometry.pix_y, + mask, + image, + center_x, + center_y, + ) + mu_eff_container.optical_efficiency = ( + comput_absolute_optical_efficiency_from_muon_ring() + ) + return mu_eff_container From f8ee9169945fce70e5a33f0d9a1e5321294cfcf8 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 12 Sep 2025 18:20:29 +0200 Subject: [PATCH 24/38] add compute_muon_ring_width function --- .../muon/impactpoint_intensity_fitter.py | 26 ++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 277c3c9a5be..ceed3f8ab58 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -124,7 +124,31 @@ def compute_muon_ring_width( """ - return 0.1 * x.unit + radius_bin_resolution = 0.05 * u.deg + + camera_unit = x.unit + x, y, ring_center_x, ring_center_y, radius_bin_resolution = all_to_value( + x, y, ring_center_x, ring_center_y, radius_bin_resolution, unit=camera_unit + ) + + max_fov = np.abs(2 * x.max()) + + n_ring_radius_bins = int(2 * max_fov / radius_bin_resolution) + + hist_ring_radius = np.histogram( + np.sqrt((y[mask] - ring_center_y) ** 2 + (x[mask] - ring_center_x) ** 2), + weights=image[mask], + bins=np.linspace(-max_fov, max_fov, n_ring_radius_bins + 1), + ) + + bin_centers = (hist_ring_radius[1][:-1] + hist_ring_radius[1][1:]) / 2 + max_count = np.max(hist_ring_radius[0]) + half_max = max_count / 2 + indices = np.where(hist_ring_radius[0] >= half_max)[0] + fwhm = bin_centers[indices[-1]] - bin_centers[indices[0]] + print(fwhm) + + return fwhm / 2.3548 * camera_unit def comput_absolute_optical_efficiency_from_muon_ring(): From cd78b240624198170f319840e1e50c87506fc2e7 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Fri, 12 Sep 2025 18:54:39 +0200 Subject: [PATCH 25/38] update save_histogram_to_csv for tests --- .../image/muon/impactpoint_intensity_fitter.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index ceed3f8ab58..5177a36647c 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -86,18 +86,14 @@ def chord_length(radius, rho, phi0, phi): return chord_length(radius, rho, 0, phi - phi0) -def save_histogram_to_csv(hist, csvName, event_id, hist_phi_smooth): - # print("event_id => ", event_id) - # print(type(event_id)) - # print(len(hist[0])) +def save_histogram_to_csv(hist): df = pd.DataFrame( { - "event_id": event_id, "x": ((np.roll(hist[1], 1) + hist[1]) / 2.0)[1:], - "y": hist_phi_smooth, + "y": hist[0], } ) - + csvName = str("save_histogram_to_csv" + str(np.random.randint(0, 10000)) + ".csv") df.to_csv(csvName, sep=" ", index=False) return @@ -147,6 +143,7 @@ def compute_muon_ring_width( indices = np.where(hist_ring_radius[0] >= half_max)[0] fwhm = bin_centers[indices[-1]] - bin_centers[indices[0]] print(fwhm) + save_histogram_to_csv(hist_ring_radius) return fwhm / 2.3548 * camera_unit From d705e5986fdf73ed19e980c1e4c59c1f7687ca89 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 14:16:22 +0200 Subject: [PATCH 26/38] Test non-fitting algorithm: simple standard deviation calculation. Very sensitive to outliers. --- .../muon/impactpoint_intensity_fitter.py | 33 +++++++++++-------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 5177a36647c..477a2c8e359 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -128,24 +128,31 @@ def compute_muon_ring_width( ) max_fov = np.abs(2 * x.max()) - n_ring_radius_bins = int(2 * max_fov / radius_bin_resolution) - hist_ring_radius = np.histogram( - np.sqrt((y[mask] - ring_center_y) ** 2 + (x[mask] - ring_center_x) ** 2), - weights=image[mask], - bins=np.linspace(-max_fov, max_fov, n_ring_radius_bins + 1), + ring_radius = np.sqrt( + (y[mask] - ring_center_y) ** 2 + (x[mask] - ring_center_x) ** 2 + ) + weights = image[mask] + save_histogram_to_csv( + np.histogram( + ring_radius, + weights=weights, + bins=np.linspace(-max_fov, max_fov, n_ring_radius_bins + 1), + ) ) - bin_centers = (hist_ring_radius[1][:-1] + hist_ring_radius[1][1:]) / 2 - max_count = np.max(hist_ring_radius[0]) - half_max = max_count / 2 - indices = np.where(hist_ring_radius[0] >= half_max)[0] - fwhm = bin_centers[indices[-1]] - bin_centers[indices[0]] - print(fwhm) - save_histogram_to_csv(hist_ring_radius) + if weights.sum() == 0.0: + return 0.0 * u.deg + + weighted_mean = np.average( + ring_radius, + weights=weights, + ) + weighted_var = np.average((ring_radius - weighted_mean) ** 2, weights=weights) + print(np.sqrt(weighted_var) * camera_unit) - return fwhm / 2.3548 * camera_unit + return np.sqrt(weighted_var) * camera_unit def comput_absolute_optical_efficiency_from_muon_ring(): From d356921dcdb06e9248ba1849b147dea6726be988 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 15:51:19 +0200 Subject: [PATCH 27/38] add fit_muon_ring_width function --- .../muon/impactpoint_intensity_fitter.py | 85 ++++++++++++++----- 1 file changed, 66 insertions(+), 19 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 477a2c8e359..dfeab32fd6b 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -86,6 +86,10 @@ def chord_length(radius, rho, phi0, phi): return chord_length(radius, rho, 0, phi - phi0) +def gauss_pedestal(x, A, mu, sigma, pedestal): + return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + pedestal + + def save_histogram_to_csv(hist): df = pd.DataFrame( { @@ -99,13 +103,14 @@ def save_histogram_to_csv(hist): return -def compute_muon_ring_width( +def fit_muon_ring_width( x, y, mask, image, ring_center_x, ring_center_y, + radius, ): """ compute_muon_ring_width @@ -120,11 +125,20 @@ def compute_muon_ring_width( """ + if Minuit is None: + raise OptionalDependencyMissing("iminuit") + radius_bin_resolution = 0.05 * u.deg camera_unit = x.unit - x, y, ring_center_x, ring_center_y, radius_bin_resolution = all_to_value( - x, y, ring_center_x, ring_center_y, radius_bin_resolution, unit=camera_unit + x, y, ring_center_x, ring_center_y, radius, radius_bin_resolution = all_to_value( + x, + y, + ring_center_x, + ring_center_y, + radius, + radius_bin_resolution, + unit=camera_unit, ) max_fov = np.abs(2 * x.max()) @@ -134,25 +148,41 @@ def compute_muon_ring_width( (y[mask] - ring_center_y) ** 2 + (x[mask] - ring_center_x) ** 2 ) weights = image[mask] - save_histogram_to_csv( - np.histogram( - ring_radius, - weights=weights, - bins=np.linspace(-max_fov, max_fov, n_ring_radius_bins + 1), - ) - ) - if weights.sum() == 0.0: - return 0.0 * u.deg - - weighted_mean = np.average( + hist_ring_radius = np.histogram( ring_radius, weights=weights, + bins=np.linspace(-max_fov, max_fov, n_ring_radius_bins + 1), + ) + + save_histogram_to_csv(hist_ring_radius) + + ring_radius_y = hist_ring_radius[0] + ring_radius_y_err = np.sqrt(np.abs(ring_radius_y)) + + ring_radius_x = ((np.roll(hist_ring_radius[1], 1) + hist_ring_radius[1]) / 2.0)[1:] + ring_radius_x_err = np.ones(len(ring_radius_x)) * radius_bin_resolution / 2.0 + + ring_radius_err = np.sqrt(ring_radius_x_err**2 + ring_radius_y_err**2) + + # minimization method + fit = Minuit( + ring_width_loss_function(ring_radius_x, ring_radius_y, ring_radius_err), + A=np.max(ring_radius_y), + mu=radius, + sigma=radius_bin_resolution, + pedestal=0.0, ) - weighted_var = np.average((ring_radius - weighted_mean) ** 2, weights=weights) - print(np.sqrt(weighted_var) * camera_unit) + fit.errordef = Minuit.LEAST_SQUARES + # + fit.errors["A"] = 10 + fit.errors["mu"] = 0.1 + fit.errors["sigma"] = 0.1 + fit.errors["pedestal"] = 0.001 + + fit.migrad() - return np.sqrt(weighted_var) * camera_unit + return fit.values["sigma"] * camera_unit def comput_absolute_optical_efficiency_from_muon_ring(): @@ -301,7 +331,7 @@ def fit_muon_ring_phi_distribution( def phi_dist_loss_function(x, y, err, w): - """dist_loss_function + """phi_dist_loss_function x, y, err: positions of pixels surviving the cleaning should not be quantities @@ -319,6 +349,21 @@ def loss_function(amplitude, R_mirror, R_shadow, rho, phi0): return loss_function +def ring_width_loss_function(x, y, err): + """ring_width_loss_function + + x, y, err: positions of pixels surviving the cleaning + should not be quantities + w : array-like of float, weights for the points + + """ + + def loss_function_two(A, mu, sigma, pedestal): + return np.sum(((y - gauss_pedestal(x, A, mu, sigma, pedestal)) / err) ** 2) + + return loss_function_two + + class MuonImpactpointIntensityFitter(TelescopeComponent): """ Fit muon ring images with a theoretical model to estimate optical efficiency. @@ -409,14 +454,16 @@ def __call__( shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, ) - mu_eff_container.width = compute_muon_ring_width( + mu_eff_container.width = fit_muon_ring_width( geometry.pix_x, geometry.pix_y, mask, image, center_x, center_y, + radius, ) + mu_eff_container.optical_efficiency = ( comput_absolute_optical_efficiency_from_muon_ring() ) From 8868a0214e7f40c0c319d27034f66e36dd85578d Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 16:07:23 +0200 Subject: [PATCH 28/38] remove help/debug function --- src/ctapipe/image/muon/impactpoint_intensity_fitter.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index dfeab32fd6b..5f6d8d0f0f5 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -155,8 +155,6 @@ def fit_muon_ring_width( bins=np.linspace(-max_fov, max_fov, n_ring_radius_bins + 1), ) - save_histogram_to_csv(hist_ring_radius) - ring_radius_y = hist_ring_radius[0] ring_radius_y_err = np.sqrt(np.abs(ring_radius_y)) From 97184cd7ccd9e8c9c1e516edbc28b40224573e63 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 21:08:04 +0200 Subject: [PATCH 29/38] Add compute_absolute_optical_efficiency_from_muon_ring and get_measured_pe --- .../muon/impactpoint_intensity_fitter.py | 157 +++++++++++++++--- 1 file changed, 133 insertions(+), 24 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 5f6d8d0f0f5..b1dc4ee553f 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -10,6 +10,7 @@ # dev to be removed import pandas as pd from astropy.units import Quantity +from scipy.constants import alpha from scipy.ndimage import correlate1d from ...containers import MuonEfficiencyContainer @@ -64,7 +65,7 @@ def chord_length(radius, rho, phi0, phi): phi_modulo = ((phi + np.pi) % (2 * np.pi) - np.pi) * np.where(phi < -1, -1, 1) - rho_R = rho / radius + rho_R = np.abs(rho) / radius discriminant_norm = 1 - (rho_R**2 * np.sin(phi_modulo) ** 2) discriminant_norm[discriminant_norm < 0] = 0.0 @@ -103,6 +104,112 @@ def save_histogram_to_csv(hist): return +def get_measured_pe( + x, + y, + mask, + image, + ring_center_x, + ring_center_y, + ring_radius, + ring_width, + integration_window_in_simga, +): + """ + get_measured_pe + + + Parameters + ---------- + + Returns + ------- + + + """ + + camera_unit = x.unit + x, y, ring_center_x, ring_center_y, ring_radius, ring_width = all_to_value( + x, + y, + ring_center_x, + ring_center_y, + ring_radius, + ring_width, + unit=camera_unit, + ) + + ring_delta_radius = ( + np.hypot(x[mask] - ring_center_x, y[mask] - ring_center_y) - ring_radius + ) + selection_r = ring_delta_radius >= -ring_width * integration_window_in_simga + selection_l = ring_delta_radius <= ring_width * integration_window_in_simga + + # print(np.sum(image[mask][selection_r & selection_l])) + # print(np.sum(image[mask])) + + return np.sum(image[mask][selection_r & selection_l]) + + +def compute_absolute_optical_efficiency_from_muon_ring( + measured_number_pe, + radius, + min_lambda_m, + max_lambda_m, + hole_radius_m, + optics, + rho, +): + """ + compute_absolute_optical_efficiency_from_muon_ring + + + Parameters + ---------- + + Returns + ------- + + + """ + + if np.isnan(rho): + return np.nan + + # numerical integral of the chord + rho = rho.to_value(u.m) + + R_mirror = np.sqrt(optics.mirror_area.to_value(u.m**2) / np.pi) + + phi = np.linspace(-np.pi, np.pi, 10000) + + chord = chord_length(R_mirror, rho, 0.0, phi) - chord_length( + hole_radius_m.to_value(u.m), rho, 0.0, phi + ) + + chord_dPhi_integral = np.trapezoid( + chord, + phi, + ) + + # Predicted total number of Cherenkov photons falling on the mirror + pred_total_Cher_phot = ( + 0.5 + * alpha + * np.sin(2 * radius) + * (min_lambda_m.to_value(u.m) ** -1 - max_lambda_m.to_value(u.m) ** -1) + * chord_dPhi_integral + ) + + # print(measured_number_pe) + print(measured_number_pe / pred_total_Cher_phot.to_value()) + + if pred_total_Cher_phot > 0: + return measured_number_pe / pred_total_Cher_phot.to_value() + + return np.nan + + def fit_muon_ring_width( x, y, @@ -144,9 +251,8 @@ def fit_muon_ring_width( max_fov = np.abs(2 * x.max()) n_ring_radius_bins = int(2 * max_fov / radius_bin_resolution) - ring_radius = np.sqrt( - (y[mask] - ring_center_y) ** 2 + (x[mask] - ring_center_x) ** 2 - ) + ring_radius = np.hypot(x[mask] - ring_center_x, y[mask] - ring_center_y) + weights = image[mask] hist_ring_radius = np.histogram( @@ -180,24 +286,7 @@ def fit_muon_ring_width( fit.migrad() - return fit.values["sigma"] * camera_unit - - -def comput_absolute_optical_efficiency_from_muon_ring(): - """ - comput_absolute_optical_efficiency_from_muon_ring - - - Parameters - ---------- - - Returns - ------- - - - """ - - return 0.2 + return np.abs(fit.values["sigma"]) * camera_unit def fit_muon_ring_phi_distribution( @@ -316,7 +405,7 @@ def fit_muon_ring_phi_distribution( fit.migrad() - rho = Quantity(fit.values["rho"], size_impact_point_unit) + rho = Quantity(np.abs(fit.values["rho"]), size_impact_point_unit) return MuonEfficiencyContainer( impact=rho, @@ -462,8 +551,28 @@ def __call__( radius, ) + measured_number_pe = get_measured_pe( + geometry.pix_x, + geometry.pix_y, + mask, + image, + center_x, + center_y, + radius, + mu_eff_container.width, + integration_window_in_simga=3.5, + ) + mu_eff_container.optical_efficiency = ( - comput_absolute_optical_efficiency_from_muon_ring() + compute_absolute_optical_efficiency_from_muon_ring( + measured_number_pe, + radius, + self.min_lambda_m.tel[tel_id] * u.m, + self.max_lambda_m.tel[tel_id] * u.m, + self.hole_radius_m.tel[tel_id] * u.m, + telescope.optics, + mu_eff_container.impact, + ) ) return mu_eff_container From 71aca3b1fcc5481986e3528ffa9c7a97c0c352e0 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 21:11:01 +0200 Subject: [PATCH 30/38] remove prints --- src/ctapipe/image/muon/impactpoint_intensity_fitter.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index b1dc4ee553f..36942b2a2e4 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -201,9 +201,6 @@ def compute_absolute_optical_efficiency_from_muon_ring( * chord_dPhi_integral ) - # print(measured_number_pe) - print(measured_number_pe / pred_total_Cher_phot.to_value()) - if pred_total_Cher_phot > 0: return measured_number_pe / pred_total_Cher_phot.to_value() From e81f830fe7e4c6fa188099d5f367508f97344ece Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 21:12:14 +0200 Subject: [PATCH 31/38] save_histogram_to_csv debug function --- .../image/muon/impactpoint_intensity_fitter.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 36942b2a2e4..6b06a54765f 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -8,7 +8,6 @@ import numpy as np # dev to be removed -import pandas as pd from astropy.units import Quantity from scipy.constants import alpha from scipy.ndimage import correlate1d @@ -91,19 +90,6 @@ def gauss_pedestal(x, A, mu, sigma, pedestal): return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + pedestal -def save_histogram_to_csv(hist): - df = pd.DataFrame( - { - "x": ((np.roll(hist[1], 1) + hist[1]) / 2.0)[1:], - "y": hist[0], - } - ) - csvName = str("save_histogram_to_csv" + str(np.random.randint(0, 10000)) + ".csv") - df.to_csv(csvName, sep=" ", index=False) - - return - - def get_measured_pe( x, y, From 13f307754a3cecff33daa309a4d4d941c702c989 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Sun, 14 Sep 2025 21:13:18 +0200 Subject: [PATCH 32/38] remove coment --- src/ctapipe/image/muon/impactpoint_intensity_fitter.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 6b06a54765f..9e557104f24 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -6,8 +6,6 @@ import astropy.units as u import numpy as np - -# dev to be removed from astropy.units import Quantity from scipy.constants import alpha from scipy.ndimage import correlate1d From c4d27df0a686307faf29b0ea479c63ae3f86451e Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 15 Sep 2025 14:28:03 +0200 Subject: [PATCH 33/38] make the formula easier to read --- src/ctapipe/image/muon/impactpoint_intensity_fitter.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 9e557104f24..0a936147f48 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -167,9 +167,8 @@ def compute_absolute_optical_efficiency_from_muon_ring( phi = np.linspace(-np.pi, np.pi, 10000) - chord = chord_length(R_mirror, rho, 0.0, phi) - chord_length( - hole_radius_m.to_value(u.m), rho, 0.0, phi - ) + chord = chord_length(R_mirror, rho, 0.0, phi) + chord -= chord_length(hole_radius_m.to_value(u.m), rho, 0.0, phi) chord_dPhi_integral = np.trapezoid( chord, From 28c86ae4516168e1e736d42adbe1d8d962ca1466 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 22 Sep 2025 17:11:13 +0200 Subject: [PATCH 34/38] Preparation to replace the current intensity fitter with a new one - using a 1D phi distribution. --- src/ctapipe/image/muon/processor.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index ba78bc97e02..0ad5f85c4c1 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -21,7 +21,7 @@ ring_containment, ring_intensity_parameters, ) -from .impactpoint_intensity_fitter import MuonImpactpointIntensityFitter +from .intensity_fitter import MuonIntensityFitter from .ring_fitter import MuonRingFitter INVALID = MuonTelescopeContainer() @@ -131,10 +131,7 @@ def __init__(self, subarray, **kwargs): self.ring_fitter = MuonRingFitter(parent=self) - # self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self) - self.intensity_fitter = MuonImpactpointIntensityFitter( - subarray=subarray, parent=self - ) + self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self) def __call__(self, event: ArrayEventContainer): for tel_id in event.dl1.tel: From fe167fad82a541f1ebc6310d98b78f0c95670ceb Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 22 Sep 2025 17:15:15 +0200 Subject: [PATCH 35/38] put back mask argument --- src/ctapipe/image/muon/processor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py index 0ad5f85c4c1..a96a8ba4849 100644 --- a/src/ctapipe/image/muon/processor.py +++ b/src/ctapipe/image/muon/processor.py @@ -206,8 +206,8 @@ def _process_telescope_event(self, event, tel_id): ring.center_fov_lat, ring.radius, image, - pedestal=np.full(mask.shape, self.pedestal.tel[tel_id]), mask=mask, + pedestal=np.full(mask.shape, self.pedestal.tel[tel_id]), ) self.log.debug( From 677f00b412d20a146455ee4854c465337369d681 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 22 Sep 2025 17:39:55 +0200 Subject: [PATCH 36/38] Update the Class description --- src/ctapipe/image/muon/intensity_fitter.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/image/muon/intensity_fitter.py b/src/ctapipe/image/muon/intensity_fitter.py index 6a06a6b8d32..d62d7435cfd 100644 --- a/src/ctapipe/image/muon/intensity_fitter.py +++ b/src/ctapipe/image/muon/intensity_fitter.py @@ -1,10 +1,11 @@ """ -Class for performing a HESS style 2D fit of muon images +Class for the muon fit to measure: + - the impact point on the plane perpendicular to the telescope optical axis + - muon ring width + - absolute total optical throughput -To do: - - Deal with astropy untis better, currently stripped and no checks made - - unit tests - - create container class for output + The absolute optical throughput takes into account the telescope-related + efficiencies and the atmospheric transparency in the vicinity of the telescope. """ From 8d0cf9b4d79127584a99e7ad4a6801e8fd313147 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Wed, 24 Sep 2025 18:26:19 +0200 Subject: [PATCH 37/38] tmp. commit --- .../muon/impactpoint_intensity_fitter.py | 126 +++++++++++++----- 1 file changed, 94 insertions(+), 32 deletions(-) diff --git a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py index 0a936147f48..329589a1d35 100644 --- a/src/ctapipe/image/muon/impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/impactpoint_intensity_fitter.py @@ -58,7 +58,7 @@ def chord_length(radius, rho, phi0, phi): if phi0 == 0: if radius <= 0: - return np.zero(len(phi)) + return np.zeros(len(phi)) phi_modulo = ((phi + np.pi) % (2 * np.pi) - np.pi) * np.where(phi < -1, -1, 1) @@ -259,7 +259,7 @@ def fit_muon_ring_width( ) fit.errordef = Minuit.LEAST_SQUARES # - fit.errors["A"] = 10 + fit.errors["A"] = 0.01 fit.errors["mu"] = 0.1 fit.errors["sigma"] = 0.1 fit.errors["pedestal"] = 0.001 @@ -276,6 +276,7 @@ def fit_muon_ring_phi_distribution( image, ring_center_x, ring_center_y, + ring_radius, optics, shadow_radius, ): @@ -319,17 +320,21 @@ def fit_muon_ring_phi_distribution( x, y, ring_center_x, ring_center_y, unit=camera_unit ) - n_phi_bins = 12 # to be added to configure file as input perameters - n_of_smoothing_points = 1 # 1 --> no smoothing + if ring_radius.to_value(camera_unit) <= 0.0: + return MuonEfficiencyContainer() + + n_phi_bins = 30 # to be added to configure file as input perameters + n_of_smoothing_points = 1 # smoothing (1 --> no smoothing) hist_phi = np.histogram( np.arctan2( y[mask] - ring_center_y, x[mask] - ring_center_x, ), - weights=image[mask], + weights=image[mask] / np.sin(2 * ring_radius.to_value(u.rad)), bins=np.linspace(-np.pi, np.pi, n_phi_bins + 1), ) + # weights=image[mask] / np.sin(2 * ring_radius.to_value(u.rad)), phi_y = ( correlate1d(hist_phi[0], np.ones(n_of_smoothing_points), mode="wrap", axis=0) @@ -337,7 +342,8 @@ def fit_muon_ring_phi_distribution( )[0] phi_y_err = np.sqrt(np.abs(phi_y)) - weights = (phi_y > 0).astype(int) + # weights = (phi_y > 0).astype(int) + weights = np.ones(len(phi_y)) phi_x = ((np.roll(hist_phi[1], 1) + hist_phi[1]) / 2.0)[1:] phi_x_err = np.ones(len(phi_x)) * np.pi / n_phi_bins @@ -349,51 +355,105 @@ def fit_muon_ring_phi_distribution( if total_integral <= 0.0: return MuonEfficiencyContainer() - amplitude_initial = np.nan - rho_initial = np.nan - phi0_initial = np.nan + # amplitude_initial = np.nan + # rho_initial = np.nan + # phi0_initial = np.nan - amplitude_initial = total_integral / 110.0 - rho_initial = 2 * (np.max(phi_y) - np.min(phi_y)) / total_integral * 110.0 + # amplitude_initial = total_integral / 110.0 * np.sin(2 * ring_radius) + # amplitude_initial = total_integral / 110.0 + amplitude_initial = 344 + # rho_initial = 3 + rho_initial = 5 phi0_initial = phi_x[np.argmax(phi_y)] - amplitude_initial = 12 if np.isnan(amplitude_initial) else amplitude_initial - rho_initial = 8 if np.isnan(rho_initial) is np.nan else rho_initial - phi0_initial = 0 if np.isnan(phi0_initial) is np.nan else phi0_initial - # minimization method - fit = Minuit( + fit1 = Minuit( phi_dist_loss_function(phi_x, phi_y, phi_err, weights), amplitude=amplitude_initial, R_mirror=np.sqrt(optics.mirror_area.to_value(u.m**2) / np.pi), - R_shadow=shadow_radius.to_value(u.m), + R_shadow=0.0, rho=rho_initial, phi0=phi0_initial, ) - fit.errordef = Minuit.LEAST_SQUARES + fit1.errordef = Minuit.LEAST_SQUARES - fit.fixed["R_mirror"] = True - fit.fixed["R_shadow"] = True + # fit1.fixed["amplitude"] = True + fit1.fixed["R_mirror"] = True + fit1.fixed["R_shadow"] = True # # set initial parameters uncertainty to a big value # taubin_error = max_fov * 0.1 - fit.errors["amplitude"] = 10 - fit.errors["R_mirror"] = 0.0001 - fit.errors["R_shadow"] = 0.0001 - fit.errors["rho"] = 10.0 - fit.errors["phi0"] = np.pi + fit1.errors["amplitude"] = 0.001 + fit1.errors["R_mirror"] = 0.0001 + fit1.errors["R_shadow"] = 0.0001 + fit1.errors["rho"] = 0.001 + fit1.errors["phi0"] = 0.001 + fit1.migrad() - fit.migrad() + # minimization method + fit2 = Minuit( + phi_dist_loss_function(phi_x, phi_y, phi_err, weights), + amplitude=amplitude_initial, + R_mirror=np.sqrt(optics.mirror_area.to_value(u.m**2) / np.pi), + R_shadow=shadow_radius.to_value(u.m), + rho=fit1.values["rho"], + phi0=fit1.values["phi0"], + ) + fit2.errordef = Minuit.LEAST_SQUARES - rho = Quantity(np.abs(fit.values["rho"]), size_impact_point_unit) + fit2.fixed["amplitude"] = True + fit2.fixed["R_mirror"] = True + fit2.fixed["R_shadow"] = True + # + # set initial parameters uncertainty to a big value + # taubin_error = max_fov * 0.1 + fit2.errors["amplitude"] = 0.001 + fit2.errors["R_mirror"] = 0.0001 + fit2.errors["R_shadow"] = 0.0001 + fit2.errors["rho"] = 0.001 + fit2.errors["phi0"] = 0.001 + fit2.migrad() + + rho = Quantity(np.abs(fit2.values["rho"]), size_impact_point_unit) return MuonEfficiencyContainer( impact=rho, - impact_x=rho * np.cos(Quantity(fit.values["phi0"], camera_unit)), - impact_y=rho * np.sin(Quantity(fit.values["phi0"], camera_unit)), - is_valid=fit.valid, - parameters_at_limit=fit.fmin.has_parameters_at_limit, - likelihood_value=fit.fval, + impact_x=rho * np.cos(Quantity(fit2.values["phi0"], u.rad)), + impact_y=rho * np.sin(Quantity(fit2.values["phi0"], u.rad)), + chord_fit_ampl=fit2.values["amplitude"], + phi_bin_0=phi_y[0], + phi_bin_1=phi_y[1], + phi_bin_2=phi_y[2], + phi_bin_3=phi_y[3], + phi_bin_4=phi_y[4], + phi_bin_5=phi_y[5], + phi_bin_6=phi_y[6], + phi_bin_7=phi_y[7], + phi_bin_8=phi_y[8], + phi_bin_9=phi_y[9], + phi_bin_10=phi_y[10], + phi_bin_11=phi_y[11], + phi_bin_12=phi_y[12], + phi_bin_13=phi_y[13], + phi_bin_14=phi_y[14], + phi_bin_15=phi_y[15], + phi_bin_16=phi_y[16], + phi_bin_17=phi_y[17], + phi_bin_18=phi_y[18], + phi_bin_19=phi_y[19], + phi_bin_20=phi_y[20], + phi_bin_21=phi_y[21], + phi_bin_22=phi_y[22], + phi_bin_23=phi_y[23], + phi_bin_24=phi_y[24], + phi_bin_25=phi_y[25], + phi_bin_26=phi_y[26], + phi_bin_27=phi_y[27], + phi_bin_28=phi_y[28], + phi_bin_29=phi_y[29], + is_valid=fit2.valid, + parameters_at_limit=fit2.fmin.has_parameters_at_limit, + likelihood_value=fit2.fval, ) @@ -408,6 +468,7 @@ def phi_dist_loss_function(x, y, err, w): def loss_function(amplitude, R_mirror, R_shadow, rho, phi0): signal = amplitude * chord_length(R_mirror, rho, phi0, x) + # diff_squared = ((signal - y) * w / err) ** 2 shadow = amplitude * chord_length(R_shadow, rho, phi0, x) diff_squared = ((signal - shadow - y) * w / err) ** 2 @@ -517,6 +578,7 @@ def __call__( image, center_x, center_y, + radius, optics=telescope.optics, shadow_radius=self.hole_radius_m.tel[tel_id] * u.m, ) From ad13e565d8a0a0662c34334f725c6f3a762a0ee8 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Wed, 24 Sep 2025 18:27:34 +0200 Subject: [PATCH 38/38] tmp. commit --- .../tests/test_impactpoint_intensity_fitter.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py b/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py index a8747eeedb7..95cf3a6ba4d 100644 --- a/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py +++ b/src/ctapipe/image/muon/tests/test_impactpoint_intensity_fitter.py @@ -6,6 +6,7 @@ def test_dummy(prod5_lst, reference_location): from ctapipe.image.muon.impactpoint_intensity_fitter import ( MuonImpactpointIntensityFitter, + compute_absolute_optical_efficiency_from_muon_ring, ) from ctapipe.instrument import SubarrayDescription @@ -31,4 +32,17 @@ def test_dummy(prod5_lst, reference_location): pedestal=np.zeros(telescope.camera.geometry.n_pixels), ) + res0 = compute_absolute_optical_efficiency_from_muon_ring( + measured_number_pe=1.0, + radius=1.1 * u.deg, + min_lambda_m=300e-9 * u.m, + max_lambda_m=600e-9 * u.m, + hole_radius_m=0.0 * u.m, + optics=telescope.optics, + rho=0.0 * u.m, + ) + + print(" ") + print("res: ", 0.0, " ", 1.0 / res0) + assert 1