Properly transform between camera and telescope frame in muon fitter

ctapipe/image/muon/intensity_fitter.py

@@ -12,7 +12,6 @@
 
 import numpy as np
 from astropy import units as u
-from astropy.coordinates import SkyCoord
 from iminuit import Minuit
 from numba import double, vectorize
 from scipy.constants import alpha
@@ -21,7 +20,7 @@
 from ctapipe.image.pixel_likelihood import neg_log_likelihood_approx
 
 from ...containers import MuonEfficiencyContainer
-from ...coordinates import CameraFrame, TelescopeFrame
+from ...coordinates import TelescopeFrame
 from ...core import TelescopeComponent
 from ...core.traits import FloatTelescopeParameter, IntTelescopeParameter
 
@@ -309,7 +308,8 @@ def image_prediction_no_units(
 
 def build_negative_log_likelihood(
     image,
-    telescope_description,
+    optics,
+    geometry_tel_frame,
     mask,
     oversampling,
     min_lambda,
@@ -330,30 +330,20 @@ def build_negative_log_likelihood(
     """
 
     # get all the neeed values and transform them into appropriate units
-    optics = telescope_description.optics
     mirror_area = optics.mirror_area.to_value(u.m**2)
     mirror_radius = np.sqrt(mirror_area / np.pi)
 
-    focal_length = optics.equivalent_focal_length
-
-    cam = telescope_description.camera.geometry
-    camera_frame = CameraFrame(focal_length=focal_length, rotation=cam.cam_rotation)
-    cam_coords = SkyCoord(x=cam.pix_x, y=cam.pix_y, frame=camera_frame)
-    tel_coords = cam_coords.transform_to(TelescopeFrame())
-
     # Use only a subset of pixels, indicated by mask:
-    pixel_x = tel_coords.fov_lon.to_value(u.rad)
-    pixel_y = tel_coords.fov_lat.to_value(u.rad)
+    pixel_x = geometry_tel_frame.pix_x.to_value(u.rad)
+    pixel_y = geometry_tel_frame.pix_y.to_value(u.rad)
 
     if mask is not None:
         pixel_x = pixel_x[mask]
         pixel_y = pixel_y[mask]
         image = image[mask]
         pedestal = pedestal[mask]
 
-    pixel_diameter = 2 * (
-        np.sqrt(cam.pix_area[0] / np.pi) / focal_length * u.rad
-    ).to_value(u.rad)
+    pixel_diameter = geometry_tel_frame.pixel_width[0].to_value(u.rad)
 
     min_lambda = min_lambda.to_value(u.m)
     max_lambda = max_lambda.to_value(u.m)
@@ -466,6 +456,14 @@ class MuonIntensityFitter(TelescopeComponent):
         help="Oversampling for the line integration", default_value=3
     ).tag(config=True)
 
+    def __init__(self, subarray, **kwargs):
+        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):
         """
 
@@ -498,9 +496,10 @@ def __call__(self, tel_id, center_x, center_y, radius, image, pedestal, mask=Non
             )
 
         negative_log_likelihood = build_negative_log_likelihood(
-            image,
-            telescope,
-            mask,
+            image=image,
+            optics=telescope.optics,
+            geometry_tel_frame=self._geometries_tel_frame[tel_id],
+            mask=mask,
             oversampling=self.oversampling.tel[tel_id],
             min_lambda=self.min_lambda_m.tel[tel_id] * u.m,
             max_lambda=self.max_lambda_m.tel[tel_id] * u.m,

ctapipe/image/muon/tests/test_intensity_fit.py

@@ -20,18 +20,19 @@ def test_chord_length():
 
 
 def test_muon_efficiency_fit(prod5_lst, reference_location):
-    from ctapipe.coordinates import CameraFrame, TelescopeFrame
+    from ctapipe.coordinates import TelescopeFrame
     from ctapipe.image.muon.intensity_fitter import (
         MuonIntensityFitter,
         image_prediction,
     )
     from ctapipe.instrument import SubarrayDescription
 
+    tel_id = 1
     telescope = prod5_lst
     subarray = SubarrayDescription(
         name="LSTMono",
-        tel_positions={0: [0, 0, 0] * u.m},
-        tel_descriptions={0: telescope},
+        tel_positions={tel_id: [0, 0, 0] * u.m},
+        tel_descriptions={tel_id: telescope},
         reference_location=reference_location,
     )
 
@@ -43,29 +44,18 @@ def test_muon_efficiency_fit(prod5_lst, reference_location):
     phi = 0 * u.rad
     efficiency = 0.5
 
-    focal_length = telescope.optics.equivalent_focal_length
-    geom = telescope.camera.geometry
+    geom = telescope.camera.geometry.transform_to(TelescopeFrame())
     mirror_radius = np.sqrt(telescope.optics.mirror_area / np.pi)
-    pixel_diameter = (
-        2
-        * u.rad
-        * (np.sqrt(geom.pix_area / np.pi) / focal_length).to_value(
-            u.dimensionless_unscaled
-        )
-    )
 
-    tel = CameraFrame(
-        x=geom.pix_x,
-        y=geom.pix_y,
-        focal_length=focal_length,
-        rotation=geom.cam_rotation,
-    ).transform_to(TelescopeFrame())
-    x = tel.fov_lon
-    y = tel.fov_lat
+    pixel_diameter = geom.pixel_width[0]
+    x = geom.pix_x
+    y = geom.pix_y
+
+    fitter = MuonIntensityFitter(subarray=subarray)
 
     image = image_prediction(
         mirror_radius,
-        hole_radius=0 * u.m,
+        hole_radius=fitter.hole_radius_m.tel[tel_id] * u.m,
         impact_parameter=impact_parameter,
         phi=phi,
         center_x=center_x,
@@ -74,12 +64,11 @@ def test_muon_efficiency_fit(prod5_lst, reference_location):
         ring_width=ring_width,
         pixel_x=x,
         pixel_y=y,
-        pixel_diameter=pixel_diameter[0],
+        pixel_diameter=pixel_diameter,
     )
 
-    fitter = MuonIntensityFitter(subarray=subarray)
     result = fitter(
-        tel_id=0,
+        tel_id=tel_id,
         center_x=center_x,
         center_y=center_y,
         radius=radius,

docs/changes/2464.bugfix.rst

@@ -0,0 +1,6 @@
+Properly transform pixel coordinates between ``CameraFrame``
+and ``TelescopeFrame`` in ``MuonIntensityFitter`` taking.
+Before, ``MuonIntensityFitter`` always used the equivalent focal
+length for transformations, now it is using the focal length
+attached to the ``CameraGeometry``, thus respecting the
+``focal_length_choice`` options of the event sources.