Fix muon normalization for different pixel shapes (cta-observatory#2865)

* Add a function to calculate the trivial solution for the number of photons (from muon) incident on the telescope mirror.

* Add the PixelShape. For image_prediction, we fixed the normalization factor in this way.

* add Add the PixelShape

* Add test for normalisation factor

* Enhance the chord_length function by applying modulo two pi

* Change the parameter rho to absolute values instead of values relative to the mirror radius.

* Add phi0 and change phi signe.

* Add phi0.

* Documentation update. Add desctiption of phi0.

* Remove uncovered code.

* Removing unnecessary sign verification and correction.

* Add a wrapper to preserve the functional API.

* from camel case to snake case

* from camel case to snake case

* add the bugfix changelog statement

* Add the bugfix changelog statement - for chord calculation

* Get back to initail intersect_circle API

* swap phi and phi0

* remove : Fix the modulo 2 pi for chord calculation.

* white space

* Prune unrelated elements to pixel shape normalization in the fitter.

* Prune unrelated elements to pixel shape normalization in the test

* renormalise rho

* renames: 2845.bugfix.rst -> 2865.bugfix.rst

* deleted : ../docs/changes/2845.bugfix.rst

Author: burmist-git

docs/changes/2865.bugfix.rst

@@ -0,0 +1 @@
+- Fix the normalization for muon analysis to account for different PixelShapes.

src/ctapipe/image/muon/intensity_fitter.py

@@ -17,6 +17,7 @@
 from ...core.env import CTAPIPE_DISABLE_NUMBA_CACHE
 from ...core.traits import FloatTelescopeParameter, IntTelescopeParameter
 from ...exceptions import OptionalDependencyMissing
+from ...instrument.camera.geometry import PixelShape
 from ..pixel_likelihood import neg_log_likelihood_approx
 
 try:
@@ -31,6 +32,9 @@
 # ratio of the areas of the unit circle and a square of side lengths 2
 CIRCLE_SQUARE_AREA_RATIO = np.pi / 4
 
+# ratio of the areas of the unit circle and a hexagon of flat-to-flat lengths 2
+CIRCLE_HEXAGON_AREA_RATIO = np.pi / 2 / np.sqrt(3)
+
 # Sqrt of 2, as it is needed multiple times
 SQRT2 = np.sqrt(2)
 
@@ -176,6 +180,7 @@ def image_prediction(
     oversampling=3,
     min_lambda=300 * u.nm,
     max_lambda=600 * u.nm,
+    pix_type=PixelShape.HEXAGON,
 ):
     """
     Predict muon ring image from given parameters.
@@ -217,6 +222,7 @@ def image_prediction(
         oversampling=oversampling,
         min_lambda_m=min_lambda.to_value(u.m),
         max_lambda_m=max_lambda.to_value(u.m),
+        pix_type=pix_type,
     )
 
 
@@ -257,6 +263,7 @@ def image_prediction_no_units(
     oversampling=3,
     min_lambda_m=300e-9,
     max_lambda_m=600e-9,
+    pix_type=PixelShape.HEXAGON,
 ):
     """
     Unit-less version of `image_prediction`.
@@ -317,7 +324,12 @@ def image_prediction_no_units(
     # diameter. In any case, since in the end we do a data-MC comparison of the muon
     # ring analysis outputs, it is not critical that this value is exact.
 
-    pred *= CIRCLE_SQUARE_AREA_RATIO
+    if pix_type == PixelShape.HEXAGON:
+        pred *= CIRCLE_HEXAGON_AREA_RATIO
+    elif pix_type == PixelShape.SQUARE:
+        pred *= CIRCLE_SQUARE_AREA_RATIO
+    elif pix_type == PixelShape.CIRCLE:
+        return pred
 
     return pred
 
@@ -333,6 +345,7 @@ def build_negative_log_likelihood(
     spe_width,
     pedestal,
     hole_radius=0 * u.m,
+    pix_type=PixelShape.HEXAGON,
 ):
     """Create an efficient negative log_likelihood function that does
     not rely on astropy units internally by defining needed values as closures
@@ -414,6 +427,7 @@ def negative_log_likelihood(
             oversampling=oversampling,
             min_lambda_m=min_lambda,
             max_lambda_m=max_lambda,
+            pix_type=pix_type,
         )
 
         # scale prediction by optical efficiency of the telescope
@@ -535,6 +549,7 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non
             spe_width=self.spe_width.tel[tel_id],
             pedestal=pedestal,
             hole_radius=self.hole_radius_m.tel[tel_id] * u.m,
+            pix_type=telescope.camera.geometry.pix_type,
         )
         negative_log_likelihood.errordef = Minuit.LIKELIHOOD
 

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

@@ -3,6 +3,7 @@
 import astropy.units as u
 import numpy as np
 import pytest
+from scipy.constants import alpha
 
 parameter_names = [
     "radius",
@@ -24,19 +25,19 @@
         ),
         Parameters(
             radius=12,
-            rho=1.0,
+            rho=1,
             phi=90.0 * u.deg,
             expected_length=0,
         ),
         Parameters(
             radius=12,
-            rho=3.0,
+            rho=1.1,
             phi=180.0 * u.deg,
             expected_length=0,
         ),
         Parameters(
             radius=12,
-            rho=2.0,
+            rho=2,
             phi=0.0 * u.deg,
             expected_length=24,
         ),
@@ -102,6 +103,7 @@ def test_muon_efficiency_fit(prod5_lst, reference_location):
         pixel_x=x,
         pixel_y=y,
         pixel_diameter=pixel_diameter,
+        pix_type=telescope.camera.geometry.pix_type,
     )
 
     result = fitter(
@@ -145,3 +147,84 @@ def test_scts(prod5_sst, reference_location):
             image=np.zeros(telescope.camera.geometry.n_pixels),
             pedestal=np.zeros(telescope.camera.geometry.n_pixels),
         )
+
+
+def test_normalisation_factor(prod5_lst, reference_location):
+    """Test of the absolute normalization factor."""
+    from ctapipe.coordinates import TelescopeFrame
+    from ctapipe.image.muon.intensity_fitter import (
+        image_prediction,
+    )
+
+    pytest.importorskip("iminuit")
+
+    telescope = prod5_lst
+
+    geom = telescope.camera.geometry.transform_to(TelescopeFrame())
+    mirror_radius = np.sqrt(telescope.optics.mirror_area / np.pi)
+
+    pixel_diameter = geom.pixel_width[0]
+    x = geom.pix_x
+    y = geom.pix_y
+
+    image = image_prediction(
+        mirror_radius,
+        hole_radius=0 * u.m,
+        impact_parameter=0 * u.m,
+        phi=0 * u.rad,
+        center_x=0.0 * u.deg,
+        center_y=0.0 * u.deg,
+        radius=1.1 * u.deg,
+        ring_width=0.05 * u.deg,
+        pixel_x=x,
+        pixel_y=y,
+        pixel_diameter=pixel_diameter,
+        oversampling=3,
+        min_lambda=300 * u.nm,
+        max_lambda=600 * u.nm,
+        pix_type=telescope.camera.geometry.pix_type,
+    )
+
+    measured = np.sum(image)
+    expected = expected_nphot(
+        r_mirror=mirror_radius,
+        theta_cher=1.1 * u.deg,
+        lambda_min=300 * u.nm,
+        lambda_max=600 * u.nm,
+    )
+
+    assert u.isclose(measured, expected, rtol=0.02)
+
+
+def expected_nphot(r_mirror, theta_cher, lambda_min, lambda_max):
+    """
+    The trivial solution for the number of photons incident on the telescope mirror.
+
+    It is a trivial case, since we assume a muon impact at the center of the dish,
+    with no shadowing and a constant Cherenkov angle.
+    We neglect the light yield attenuation due to atmospheric absorption.
+
+    Parameters
+    ----------
+    Rmirror: quantity[length]
+        mirror radius
+    theta_cher: quantity[angle]
+        Cherenkov angle
+    lambda_min: quantity[length]
+        photon wavelength
+    lambda_max: quantity[length]
+        photon wavelength
+
+    Returns
+    -------
+    float: number of Cherenkov photons
+
+    """
+
+    return (
+        np.pi
+        * alpha
+        * r_mirror.to_value(u.m)
+        * np.sin(2 * theta_cher)
+        * (lambda_min.to_value(u.m) ** -1 - lambda_max.to_value(u.m) ** -1)
+    )