Repository: cta-observatory/ctapipe
Pull request #2464 (2023-11-17 15:34:57)

Properly transform between camera and telescope frame in muon fitter
The muon code did some strange things to transform pixel coordinates from camera to telescope frame, probably from code before camera geometries could be transformed as a whole.



codecov[bot]: ## [Codecov](https://app.codecov.io/gh/cta-observatory/ctapipe/pull/2464?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`eee9b02`)](https://app.codecov.io/gh/cta-observatory/ctapipe/commit/eee9b0285e68424ac4b9580093c13607bfbeec28?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) 60.60% compared to head [(`437c75a`)](https://app.codecov.io/gh/cta-observatory/ctapipe/pull/2464?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) 60.60%.
> Report is 14 commits behind head on main.


<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main    #2464   +/-   ##
=======================================
  Coverage   60.60%   60.60%           
=======================================
  Files           3        3           
  Lines          33       33           
=======================================
  Hits           20       20           
  Misses         13       13           
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/cta-observatory/ctapipe/pull/2464?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).

maxnoe: @RuneDominik Yes, good point! I oversaw this change before.

Patch:
diff --git a/ctapipe/image/muon/intensity_fitter.py b/ctapipe/image/muon/intensity_fitter.py
index 97046316..2215b859 100644
--- a/ctapipe/image/muon/intensity_fitter.py
+++ b/ctapipe/image/muon/intensity_fitter.py
@@ -12,7 +12,6 @@ from math import erf
 
 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 scipy.ndimage import correlate1d
 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,20 +330,12 @@ 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]
@@ -351,9 +343,7 @@ def build_negative_log_likelihood(
         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 @@ class MuonIntensityFitter(TelescopeComponent):
             )
 
         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,
diff --git a/docs/changes/2464.bugfix.rst b/docs/changes/2464.bugfix.rst
new file mode 100644
index 00000000..1102ba25
--- /dev/null
+++ b/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.


Test patch:
diff --git a/ctapipe/image/muon/tests/test_intensity_fit.py b/ctapipe/image/muon/tests/test_intensity_fit.py
index 1f424178..3f255ee3 100644
--- a/ctapipe/image/muon/tests/test_intensity_fit.py
+++ b/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,