Merge pull request #143 from HealthyPear/feature-impact_dist_tilted_f…

…rame Calculate impact parameter always in the tilted frame
cta-observatory · Jun 16, 2021 · 283bd23 · 283bd23
2 parents b1e1574 + 44696e2
commit 283bd23
Showing 1 changed file with 21 additions and 5 deletions.
diff --git a/protopipe/pipeline/event_preparer.py b/protopipe/pipeline/event_preparer.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 from astropy import units as u
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyCoord, AltAz
 import warnings
 from traitlets.config import Config
 from collections import namedtuple, OrderedDict
@@ -16,7 +16,10 @@
                            largest_island,
                            concentration_parameters)
 from ctapipe.utils import CutFlow
-from ctapipe.coordinates import GroundFrame, TelescopeFrame, CameraFrame
+from ctapipe.coordinates import (GroundFrame,
+                                 TelescopeFrame,
+                                 CameraFrame,
+                                 TiltedGroundFrame)
 
 # from ctapipe.image.timing_parameters import timing_parameters
 from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError
@@ -442,9 +445,12 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
             # Array pointing in AltAz frame
             az = event.pointing.array_azimuth
             alt = event.pointing.array_altitude
+            array_pointing = SkyCoord(az, alt, frame=AltAz())
 
             ground_frame = GroundFrame()
 
+            tilted_frame = TiltedGroundFrame(pointing_direction=array_pointing)
+
             for tel_id in event.r1.tel.keys():
 
                 point_azimuth_dict[tel_id] = event.pointing.tel[tel_id].azimuth
@@ -957,10 +963,20 @@ def prepare_event(self, source, return_stub=True, save_images=False, debug=False
                             frame=ground_frame,
                         )
 
-                        # Should be better handled (tilted frame)
+                        # Go back to the tilted frame
+
+                        # this should be the same...
+                        tel_tilted = tel_ground.transform_to(tilted_frame)
+
+                        # but this not
+                        core_tilted = SkyCoord(x=core_ground.x,
+                                               y=core_ground.y,
+                                               frame=tilted_frame
+                                               )
+
                         impact_dict_reco[tel_id] = np.sqrt(
-                            (core_ground.x - tel_ground.x) ** 2
-                            + (core_ground.y - tel_ground.y) ** 2
+                            (core_tilted.x - tel_tilted.x) ** 2
+                            + (core_tilted.y - tel_tilted.y) ** 2
                         )
 
             except (Exception, TooFewTelescopesException, InvalidWidthException) as e: