Repository: cta-observatory/ctapipe
Pull request #2036 (2022-07-21 14:36:17)

Vectorize operations in HillasReconstructor 
Timing script:

```
from ctapipe.io import EventSource
from ctapipe.reco import HillasReconstructor
import numpy as np
from timeit import repeat
from cProfile import Profile



source = EventSource("../data/gamma_20deg_180deg_run3___cta-prod6-paranal-2147m-Paranal-dark.dl1.h5", max_events=200)
events = list(source)

reco = HillasReconstructor(subarray=source.subarray)

code = '''
print("start")
for event in events:
    reco(event)
'''

with Profile() as profile:
    exec(code)

profile.dump_stats("hillas_reco.profile")


n_repeat = 5
number = 1
times = np.array(repeat(code, number=number, repeat=n_repeat, globals=globals()))
times /= number
```


On master:
```
11.73 ± 0.31 s
```

on this branch:
```
2.25 ± 0.02 s
```

Units tests are passing, but I will also process a larger DL2 run and check that results are the same.


Most time is not spend in astropy coordinate transforms and project_to_ground.

codecov[bot]: # [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
> Merging [#2036](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) (04a8a80) into [master](https://codecov.io/gh/cta-observatory/ctapipe/commit/d73bfefabfb3ed215c2d10dd794ca800cbcd4544?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) (d73bfef) will **decrease** coverage by `0.01%`.
> The diff coverage is `95.75%`.

> :exclamation: Current head 04a8a80 differs from pull request most recent head 7cf5c8f. Consider uploading reports for the commit 7cf5c8f to get more accurate results

```diff
@@            Coverage Diff             @@
##           master    #2036      +/-   ##
==========================================
- Coverage   92.43%   92.41%   -0.02%     
==========================================
  Files         196      196              
  Lines       15906    15933      +27     
==========================================
+ Hits        14703    14725      +22     
- Misses       1203     1208       +5     
```


| [Impacted Files](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | Coverage Δ | |
|---|---|---|
| [ctapipe/reco/tests/test\_reconstruction\_methods.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9yZWNvL3Rlc3RzL3Rlc3RfcmVjb25zdHJ1Y3Rpb25fbWV0aG9kcy5weQ==) | `100.00% <ø> (ø)` | |
| [ctapipe/reco/hillas\_reconstructor.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9yZWNvL2hpbGxhc19yZWNvbnN0cnVjdG9yLnB5) | `95.71% <94.28%> (-2.66%)` | :arrow_down: |
| [ctapipe/reco/tests/test\_HillasReconstructor.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9yZWNvL3Rlc3RzL3Rlc3RfSGlsbGFzUmVjb25zdHJ1Y3Rvci5weQ==) | `99.17% <97.95%> (-0.83%)` | :arrow_down: |
| [ctapipe/coordinates/\_\_init\_\_.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jb29yZGluYXRlcy9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Last update [d73bfef...7cf5c8f](https://codecov.io/gh/cta-observatory/ctapipe/pull/2036?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).

kosack: as suggested in #1922 perhaps we should also drop the CameraFrame support in HillasReconstructor (either here or in a separate PR) to simplify this code even more and remove a few branch points.
maxnoe: @kosack If you look at the code, the branching for `CameraFrame` vs. `TelescopeFrame` is now very minimal, so I'd say not really a concern here.

I am investigating a few issues with astropy concerning performance and I have identified several spots where astropy can be much improved.

For now, I would keep it in the current version I think, we already get a very noticeable speedup.
maxnoe: Will do
maxnoe: Sure, I'll wrap the angle into [-180, 180)
kosack: Astropy.coord seems to wrap to (0,360], not sure which is better (either way the user has to do a wrap transform if you are one of the usual North or South pointings).   At least we should make all Reconstructors consistent. I think HillasIntersection returns something in (0,360]. Not sure about ImPACT.
maxnoe: Ok, let's go with the astropy default then
maxnoe: I processed 100 gamma diffuse runs using 0.16.0 and this branch.

The results are not exactly identical, but they are very close for almost all events and the angular resolution is the same or even slightly better with the new code:

![angular_resolution_vectorize_hillas](https://user-images.githubusercontent.com/5488440/191551505-1db9d0b0-ce90-467e-bd03-1bd0bf83ee06.png)


Patch:
diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py
index 8b4f127f..5907d6f9 100644
--- a/ctapipe/coordinates/__init__.py
+++ b/ctapipe/coordinates/__init__.py
@@ -1,27 +1,26 @@
 """
 Coordinates.
 """
+import warnings
 
-import numpy as np
-
+import astropy.units as u
 from astropy.coordinates import (
+    CIRS,
     AltAz,
     FunctionTransformWithFiniteDifference,
-    CIRS,
     frame_transform_graph,
-    spherical_to_cartesian,
 )
-import warnings
-from .telescope_frame import TelescopeFrame
-from .nominal_frame import NominalFrame
+from erfa.ufunc import s2p as spherical_to_cartesian
+
+from .camera_frame import CameraFrame, EngineeringCameraFrame
 from .ground_frames import (
+    EastingNorthingFrame,
     GroundFrame,
     TiltedGroundFrame,
     project_to_ground,
-    EastingNorthingFrame,
 )
-from .camera_frame import CameraFrame, EngineeringCameraFrame
-
+from .nominal_frame import NominalFrame
+from .telescope_frame import TelescopeFrame
 
 __all__ = [
     "TelescopeFrame",
@@ -67,7 +66,7 @@ def altaz_to_altaz(from_coo, to_frame):
         location = to_frame.location
 
     if obstime is None or location is None:
-        return to_frame.realize_frame(from_coo.spherical)
+        return to_frame.realize_frame(from_coo.data)
 
     return from_coo.transform_to(CIRS(obstime=obstime)).transform_to(to_frame)
 
@@ -85,4 +84,8 @@ def altaz_to_righthanded_cartesian(alt, az):
     az: u.Quantity
         azimuth
     """
-    return np.array(spherical_to_cartesian(r=1, lat=alt, lon=-az))
+    if hasattr(az, "unit"):
+        az = az.to_value(u.rad)
+        alt = alt.to_value(u.rad)
+
+    return spherical_to_cartesian(-az, alt, 1.0)
diff --git a/ctapipe/coordinates/ground_frames.py b/ctapipe/coordinates/ground_frames.py
index 8b19bb4d..849fc7fa 100644
--- a/ctapipe/coordinates/ground_frames.py
+++ b/ctapipe/coordinates/ground_frames.py
@@ -14,9 +14,9 @@ TODO:
 - Tests Tests Tests!
 """
 import astropy.units as u
-from astropy.units.quantity import Quantity
 import numpy as np
 from astropy.coordinates import (
+    AffineTransform,
     AltAz,
     BaseCoordinateFrame,
     CartesianRepresentation,
@@ -24,8 +24,8 @@ from astropy.coordinates import (
     FunctionTransform,
     RepresentationMapping,
     frame_transform_graph,
-    AffineTransform,
 )
+from astropy.units.quantity import Quantity
 from numpy import cos, sin
 
 __all__ = [
@@ -216,13 +216,18 @@ def project_to_ground(tilt_system):
     z_initial = ground_system.z.value
 
     trans = get_shower_trans_matrix(
-        tilt_system.pointing_direction.az, tilt_system.pointing_direction.alt
+        tilt_system.pointing_direction.az,
+        tilt_system.pointing_direction.alt,
     )
 
     x_projected = x_initial - trans[2][0] * z_initial / trans[2][2]
     y_projected = y_initial - trans[2][1] * z_initial / trans[2][2]
 
-    return GroundFrame(x=x_projected * unit, y=y_projected * unit, z=0 * unit)
+    return GroundFrame(
+        x=u.Quantity(x_projected, unit),
+        y=u.Quantity(y_projected, unit),
+        z=u.Quantity(0, unit),
+    )
 
 
 @frame_transform_graph.transform(FunctionTransform, GroundFrame, GroundFrame)
diff --git a/ctapipe/coordinates/nominal_frame.py b/ctapipe/coordinates/nominal_frame.py
index df78124c..9e56f6f1 100644
--- a/ctapipe/coordinates/nominal_frame.py
+++ b/ctapipe/coordinates/nominal_frame.py
@@ -6,23 +6,23 @@ We are just not creating a metaclass and a factory but directly building the
 corresponding class.
 """
 import astropy.units as u
-from astropy.coordinates.matrix_utilities import (
-    rotation_matrix,
-    matrix_product,
-    matrix_transpose,
-)
 from astropy.coordinates import (
-    frame_transform_graph,
-    FunctionTransform,
-    DynamicMatrixTransform,
-    UnitSphericalRepresentation,
+    AltAz,
+    Angle,
     BaseCoordinateFrame,
     CoordinateAttribute,
-    TimeAttribute,
+    DynamicMatrixTransform,
     EarthLocationAttribute,
+    FunctionTransform,
     RepresentationMapping,
-    Angle,
-    AltAz,
+    TimeAttribute,
+    UnitSphericalRepresentation,
+    frame_transform_graph,
+)
+from astropy.coordinates.matrix_utilities import (
+    matrix_product,
+    matrix_transpose,
+    rotation_matrix,
 )
 
 __all__ = ["NominalFrame"]
@@ -73,31 +73,31 @@ class NominalFrame(BaseCoordinateFrame):
 
 
 @frame_transform_graph.transform(FunctionTransform, NominalFrame, NominalFrame)
-def skyoffset_to_skyoffset(from_telescope_coord, to_telescope_frame):
+def nominal_to_nominal(from_nominal_coord, to_nominal_frame):
     """Transform between two skyoffset frames."""
 
-    intermediate_from = from_telescope_coord.transform_to(from_telescope_coord.origin)
-    intermediate_to = intermediate_from.transform_to(to_telescope_frame.origin)
-    return intermediate_to.transform_to(to_telescope_frame)
+    intermediate_from = from_nominal_coord.transform_to(from_nominal_coord.origin)
+    intermediate_to = intermediate_from.transform_to(to_nominal_frame.origin)
+    return intermediate_to.transform_to(to_nominal_frame)
 
 
 @frame_transform_graph.transform(DynamicMatrixTransform, AltAz, NominalFrame)
-def reference_to_skyoffset(reference_frame, telescope_frame):
+def altaz_to_nominal(altaz_coord, nominal_frame):
     """Convert a reference coordinate to an sky offset frame."""
 
     # Define rotation matrices along the position angle vector, and
     # relative to the origin.
-    origin = telescope_frame.origin.spherical
+    origin = nominal_frame.origin.represent_as(UnitSphericalRepresentation)
     mat1 = rotation_matrix(-origin.lat, "y")
     mat2 = rotation_matrix(origin.lon, "z")
     return matrix_product(mat1, mat2)
 
 
 @frame_transform_graph.transform(DynamicMatrixTransform, NominalFrame, AltAz)
-def skyoffset_to_reference(skyoffset_coord, reference_frame):
+def nominal_to_altaz(nominal_coord, altaz_frame):
     """Convert an sky offset frame coordinate to the reference frame"""
 
     # use the forward transform, but just invert it
-    mat = reference_to_skyoffset(reference_frame, skyoffset_coord)
+    mat = altaz_to_nominal(altaz_frame, nominal_coord)
     # transpose is the inverse because mat is a rotation matrix
     return matrix_transpose(mat)
diff --git a/ctapipe/coordinates/telescope_frame.py b/ctapipe/coordinates/telescope_frame.py
index 88d197cf..084404d3 100644
--- a/ctapipe/coordinates/telescope_frame.py
+++ b/ctapipe/coordinates/telescope_frame.py
@@ -6,28 +6,31 @@ We are just not creating a metaclass and a factory but directly building the
 corresponding class.
 """
 import astropy.units as u
-from astropy.coordinates.matrix_utilities import (
-    rotation_matrix,
-    matrix_product,
-    matrix_transpose,
-)
 from astropy.coordinates import (
-    frame_transform_graph,
-    FunctionTransform,
-    DynamicMatrixTransform,
+    AltAz,
+    Angle,
     BaseCoordinateFrame,
     CoordinateAttribute,
-    TimeAttribute,
+    DynamicMatrixTransform,
     EarthLocationAttribute,
+    FunctionTransform,
     RepresentationMapping,
+    TimeAttribute,
     UnitSphericalRepresentation,
-    Angle,
-    AltAz,
+    frame_transform_graph,
+)
+from astropy.coordinates.matrix_utilities import (
+    matrix_product,
+    matrix_transpose,
+    rotation_matrix,
 )
 
 __all__ = ["TelescopeFrame"]
 
 
+_wrap_angle = Angle(180, unit=u.deg)
+
+
 class TelescopeFrame(BaseCoordinateFrame):
     """
     Telescope coordinate frame.
@@ -73,11 +76,11 @@ class TelescopeFrame(BaseCoordinateFrame):
 
         # make sure telescope coordinate is in range [-180°, 180°]
         if isinstance(self._data, UnitSphericalRepresentation):
-            self._data.lon.wrap_angle = Angle(180, unit=u.deg)
+            self._data.lon.wrap_angle = _wrap_angle
 
 
 @frame_transform_graph.transform(FunctionTransform, TelescopeFrame, TelescopeFrame)
-def skyoffset_to_skyoffset(from_telescope_coord, to_telescope_frame):
+def telescope_to_telescope(from_telescope_coord, to_telescope_frame):
     """Transform between two skyoffset frames."""
 
     intermediate_from = from_telescope_coord.transform_to(
@@ -90,22 +93,24 @@ def skyoffset_to_skyoffset(from_telescope_coord, to_telescope_frame):
 
 
 @frame_transform_graph.transform(DynamicMatrixTransform, AltAz, TelescopeFrame)
-def reference_to_skyoffset(reference_frame, telescope_frame):
+def altaz_to_telescope(altaz_coord, telescope_frame):
     """Convert a reference coordinate to an sky offset frame."""
 
     # Define rotation matrices along the position angle vector, and
     # relative to the telescope_pointing.
-    telescope_pointing = telescope_frame.telescope_pointing.spherical
+    telescope_pointing = telescope_frame.telescope_pointing.represent_as(
+        UnitSphericalRepresentation
+    )
     mat1 = rotation_matrix(-telescope_pointing.lat, "y")
     mat2 = rotation_matrix(telescope_pointing.lon, "z")
     return matrix_product(mat1, mat2)
 
 
 @frame_transform_graph.transform(DynamicMatrixTransform, TelescopeFrame, AltAz)
-def skyoffset_to_reference(skyoffset_coord, reference_frame):
+def telescope_to_altaz(telescope_coord, altaz_frame):
     """Convert an sky offset frame coordinate to the reference frame"""
 
     # use the forward transform, but just invert it
-    mat = reference_to_skyoffset(reference_frame, skyoffset_coord)
+    mat = altaz_to_telescope(altaz_frame, telescope_coord)
     # transpose is the inverse because mat is a rotation matrix
     return matrix_transpose(mat)
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index 9b5c20d4..6e2dc530 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -8,7 +8,7 @@ from itertools import combinations
 
 import numpy as np
 from astropy import units as u
-from astropy.coordinates import AltAz, SkyCoord, cartesian_to_spherical
+from astropy.coordinates import AltAz, Longitude, SkyCoord, cartesian_to_spherical
 
 from ctapipe.containers import (
     CameraHillasParametersContainer,
@@ -16,7 +16,6 @@ from ctapipe.containers import (
 )
 from ctapipe.coordinates import (
     CameraFrame,
-    GroundFrame,
     MissingFrameAttributeWarning,
     TelescopeFrame,
     TiltedGroundFrame,
@@ -29,7 +28,7 @@ from ctapipe.reco.reco_algorithms import (
     TooFewTelescopesException,
 )
 
-__all__ = ["HillasPlane", "HillasReconstructor"]
+__all__ = ["HillasReconstructor"]
 
 
 INVALID = ReconstructedGeometryContainer(
@@ -51,8 +50,8 @@ def angle(v1, v2):
     -------
     the angle between v1 and v2 as a dimensioned astropy quantity
     """
-    norm = np.linalg.norm(v1) * np.linalg.norm(v2)
-    return np.arccos(np.clip(v1.dot(v2) / norm, -1.0, 1.0))
+    norm = np.linalg.norm(v1, axis=-1) * np.linalg.norm(v2, axis=-1)
+    return np.arccos(np.clip(np.sum(v1 * v2, axis=-1) / norm, -1.0, 1.0))
 
 
 def normalise(vec):
@@ -67,10 +66,17 @@ def normalise(vec):
     -------
     numpy array with the same direction but length of 1
     """
-    try:
-        return vec / np.linalg.norm(vec)
-    except ZeroDivisionError:
-        return vec
+    norm = np.linalg.norm(vec, axis=-1)
+    if vec.ndim == 1:
+        try:
+            return vec / norm
+        except ZeroDivisionError:
+            return vec
+
+    result = np.zeros(vec.shape)
+    mask = norm > 0
+    result[mask] = vec[mask] / norm[mask, np.newaxis]
+    return result
 
 
 def line_line_intersection_3d(uvw_vectors, origins):
@@ -91,66 +97,6 @@ def line_line_intersection_3d(uvw_vectors, origins):
     return np.linalg.inv(S) @ C
 
 
-class HillasPlane:
-    """
-    a tiny helper class to collect some parameters for each great great
-    circle
-
-    Stores some vectors a, b, and c
-
-    These vectors are euclidean [x, y, z] where positive z values point towards the sky
-    and x and y are parallel to the ground.
-    """
-
-    def __init__(self, p1, p2, telescope_position, weight=1):
-        r"""The constructor takes two coordinates in the horizontal
-        frame (alt, az) which define a plane perpendicular
-        to the camera.
-
-        Parameters
-        ----------
-        p1: astropy.coordinates.SkyCoord
-            One of two direction vectors which define the plane.
-            This coordinate has to be defined in the ctapipe.coordinates.AltAz
-        p2: astropy.coordinates.SkyCoord
-            One of two direction vectors which define the plane.
-            This coordinate has to be defined in the ctapipe.coordinates.AltAz
-        telescope_position: np.array(3)
-            Position of the telescope on the ground
-        weight : float, optional
-            weight of this plane for later use during the reconstruction
-
-        Notes
-        -----
-        c: numpy.ndarray(3)
-            :math:`\vec c = (\vec a \times \vec b) \times \vec a`
-            :math:`\rightarrow` a and c form an orthogonal base for the
-            great circle
-            (only orthonormal if a and b are of unit-length)
-        norm: numpy.ndarray(3)
-            normal vector of the circle's plane,
-            perpendicular to a, b and c
-        """
-
-        self.pos = telescope_position
-
-        # astropy's coordinates system rotates counter clockwise. Apparently we assume it to
-        # be clockwise
-        self.a = altaz_to_righthanded_cartesian(alt=p1.alt, az=p1.az)
-        self.b = altaz_to_righthanded_cartesian(alt=p2.alt, az=p2.az)
-
-        # a and c form an orthogonal basis for the great circle
-        # not really necessary since the norm can be calculated
-        # with a and b just as well
-        self.c = np.cross(np.cross(self.a, self.b), self.a)
-        # normal vector for the plane defined by the great circle
-        self.norm = normalise(np.cross(self.a, self.c))
-        # some weight for this circle
-        # (put e.g. uncertainty on the Hillas parameters
-        # or number of PE in here)
-        self.weight = weight
-
-
 class HillasReconstructor(Reconstructor):
     """
     class that reconstructs the direction of an atmospheric shower
@@ -165,16 +111,25 @@ class HillasReconstructor(Reconstructor):
 
     def __init__(self, subarray, **kwargs):
         super().__init__(subarray=subarray, **kwargs)
-        self._cam_radius_m = {
-            cam: cam.geometry.guess_radius() for cam in subarray.camera_types
+        _cam_radius_m = {
+            cam: cam.geometry.guess_radius().to_value(u.m)
+            for cam in subarray.camera_types
         }
 
         telframe = TelescopeFrame()
-        self._cam_radius_deg = {}
-        for cam, radius_m in self._cam_radius_m.items():
-            point_cam = SkyCoord(0 * u.m, radius_m, frame=cam.geometry.frame)
+        _cam_radius_deg = {}
+        for cam, radius_m in _cam_radius_m.items():
+            point_cam = SkyCoord(0, radius_m, unit=u.m, frame=cam.geometry.frame)
             point_tel = point_cam.transform_to(telframe)
-            self._cam_radius_deg[cam] = point_tel.fov_lon
+            _cam_radius_deg[cam] = point_tel.fov_lon.to_value(u.deg)
+
+        # store for each tel_id to avoid costly hash of camera
+        self._cam_radius_m = {
+            tel_id: _cam_radius_m[t.camera] for tel_id, t in subarray.tel.items()
+        }
+        self._cam_radius_deg = {
+            tel_id: _cam_radius_deg[t.camera] for tel_id, t in subarray.tel.items()
+        }
 
     def __call__(self, event):
         """
@@ -189,92 +144,61 @@ class HillasReconstructor(Reconstructor):
         -------
         ReconstructedGeometryContainer
         """
+        warnings.filterwarnings(action="ignore", category=MissingFrameAttributeWarning)
+
         try:
             hillas_dict = self._create_hillas_dict(event)
         except (TooFewTelescopesException, InvalidWidthException):
             return INVALID
 
-        # Due to tracking the pointing of the array will never be a constant
-        array_pointing = SkyCoord(
-            az=event.pointing.array_azimuth,
-            alt=event.pointing.array_altitude,
-            frame=AltAz(),
-        )
-
-        telescope_pointings = {
-            tel_id: SkyCoord(
-                alt=event.pointing.tel[tel_id].altitude,
-                az=event.pointing.tel[tel_id].azimuth,
-                frame=AltAz(),
-            )
-            for tel_id in event.dl1.tel.keys()
-        }
-
-        event.dl2.stereo.geometry[self.__class__.__name__] = self._predict(
-            event, hillas_dict, array_pointing, telescope_pointings
-        )
-
-        self._store_impact_parameter(event)
-
-    def _predict(self, event, hillas_dict, array_pointing, telescopes_pointings):
-        """
-        The function you want to call for the reconstruction of the
-        event. It takes care of setting up the event and consecutively
-        calls the functions for the direction and core position
-        reconstruction.  Shower parameters not reconstructed by this
-        class are set to np.nan
-
-        Parameters
-        ----------
-        hillas_dict: dict
-            dictionary with telescope IDs as key and
-            HillasParametersContainer instances as values
-        inst : ctapipe.io.InstrumentContainer
-            instrumental description
-        array_pointing: SkyCoord[AltAz]
-            pointing direction of the array
-        telescopes_pointings: dict[SkyCoord[AltAz]]
-            dictionary of pointing direction per each telescope
-
-        Raises
-        ------
-        InvalidWidthException
-            if any width is np.nan or 0
-        """
-
-        # filter warnings for missing obs time. this is needed because MC data has no obs time
-        warnings.filterwarnings(action="ignore", category=MissingFrameAttributeWarning)
-
         # Here we perform some basic quality checks BEFORE applying reconstruction
         # This should be substituted by a DL1 QualityQuery specific to this
         # reconstructor
-
-        hillas_planes, psi_core_dict = self.initialize_hillas_planes(
-            hillas_dict, self.subarray, telescopes_pointings, array_pointing
-        )
+        (
+            tel_ids,
+            cog_cartesian,
+            p2_cartesian,
+            corrected_psi,
+            weights,
+            telescope_positions,
+            array_pointing,
+        ) = self.initialize_arrays(event, hillas_dict)
+
+        norm = np.cross(cog_cartesian, p2_cartesian)
+
+        # store core corrected psi values
+        for tel_id, psi in zip(tel_ids, corrected_psi):
+            event.dl1.tel[tel_id].parameters.core.psi = u.Quantity(
+                np.rad2deg(psi), u.deg
+            )
 
         # algebraic direction estimate
-        direction, err_est_dir = self.estimate_direction(hillas_planes)
+        direction, err_est_dir = self.estimate_direction(norm, weights)
 
         # array pointing is needed to define the tilted frame
         core_pos_ground, core_pos_tilted = self.estimate_core_position(
-            event, hillas_dict, array_pointing, psi_core_dict, hillas_planes
+            array_pointing, corrected_psi, telescope_positions
         )
 
         # container class for reconstructed showers
         _, lat, lon = cartesian_to_spherical(*direction)
 
         # estimate max height of shower
-        h_max = self.estimate_h_max(hillas_planes)
+        h_max = self.estimate_h_max(cog_cartesian, telescope_positions)
+
+        # az is clockwise, lon counter-clockwise, make sure it stays in [0, 2pi)
+        az = Longitude(-lon)
 
-        return ReconstructedGeometryContainer(
+        event.dl2.stereo.geometry[
+            self.__class__.__name__
+        ] = ReconstructedGeometryContainer(
             alt=lat,
-            az=-lon,  # az is clockwise, lon counter-clockwise
+            az=az,
             core_x=core_pos_ground.x,
             core_y=core_pos_ground.y,
             core_tilted_x=core_pos_tilted.x,
             core_tilted_y=core_pos_tilted.y,
-            telescopes=[tel_id for tel_id in hillas_dict.keys()],
+            telescopes=tel_ids.tolist(),
             average_intensity=np.mean([h.intensity for h in hillas_dict.values()]),
             is_valid=True,
             alt_uncert=err_est_dir,
@@ -283,12 +207,11 @@ class HillasReconstructor(Reconstructor):
             prefix=self.__class__.__name__,
         )
 
-    def initialize_hillas_planes(
-        self, hillas_dict, subarray, telescopes_pointings, array_pointing
-    ):
+        self._store_impact_parameter(event)
+
+    def initialize_arrays(self, event, hillas_dict):
         """
-        Creates a dictionary of `.HillasPlane` from a dictionary of
-        hillas parameters
+        Creates flat arrays of needed quantities from the event structure
 
         Parameters
         ----------
@@ -307,99 +230,107 @@ class HillasReconstructor(Reconstructor):
         mode and the correction to the psi angle is explained in [gasparetto]_
         section 7.1.4.
         """
+        # get one telescope id to check what frame to use
+        altaz = AltAz()
 
-        hillas_planes = {}
-
-        # dictionary to store the telescope-wise image directions
-        # to be projected on the ground and corrected in case of mispointing
-        corrected_angle_dict = {}
+        # Due to tracking the pointing of the array will never be a constant
+        array_pointing = SkyCoord(
+            az=event.pointing.array_azimuth,
+            alt=event.pointing.array_altitude,
+            frame=altaz,
+        )
 
-        k = next(iter(telescopes_pointings))
-        horizon_frame = telescopes_pointings[k].frame
-        for tel_id, moments in hillas_dict.items():
+        # create arrays / quantities of all needed things for all telescopes
+        # so we can do transformations vectorized
+        cog1 = np.empty(len(hillas_dict))
+        cog2 = np.empty(len(hillas_dict))
+        cam_radius = np.empty(len(hillas_dict))
+        psi = np.empty(len(hillas_dict))
+        weights = np.empty(len(hillas_dict))
+        focal_length = np.empty(len(hillas_dict))
+        alt = np.empty(len(hillas_dict))
+        az = np.empty(len(hillas_dict))
+        tel_ids = np.empty(len(hillas_dict), dtype=int)
+
+        hillas_in_camera_frame = False
+
+        for i, (tel_id, hillas) in enumerate(hillas_dict.items()):
+            tel_ids[i] = tel_id
+
+            pointing = event.pointing.tel[tel_id]
+            alt[i] = pointing.altitude.to_value(u.rad)
+            az[i] = pointing.azimuth.to_value(u.rad)
+
+            if isinstance(hillas, CameraHillasParametersContainer):
+                hillas_in_camera_frame = True
+                cog1[i] = hillas.x.to_value(u.m)
+                cog2[i] = hillas.y.to_value(u.m)
+                cam_radius[i] = self._cam_radius_m[tel_id]
+            else:
+                cog1[i] = hillas.fov_lon.to_value(u.deg)
+                cog2[i] = hillas.fov_lat.to_value(u.deg)
+                cam_radius[i] = self._cam_radius_deg[tel_id]
+
+            psi[i] = hillas.psi.to_value(u.rad)
+            weights[i] = hillas.intensity * hillas.length.value / hillas.width.value
 
             camera = self.subarray.tel[tel_id].camera
+            focal_length[i] = camera.geometry.frame.focal_length.to_value(u.m)
 
-            pointing = SkyCoord(
-                alt=telescopes_pointings[tel_id].alt,
-                az=telescopes_pointings[tel_id].az,
-                frame=horizon_frame,
-            )
+        indices = self.subarray.tel_index_array[tel_ids]
+        telescope_positions = self.subarray.tel_coords[indices]
+
+        telescope_pointings = SkyCoord(alt=alt, az=az, unit=u.rad, frame=altaz)
+
+        focal_length = u.Quantity(focal_length, u.m, copy=False)
+        camera_frame = CameraFrame(
+            telescope_pointing=telescope_pointings, focal_length=focal_length
+        )
+
+        telescope_frame = TelescopeFrame(telescope_pointing=telescope_pointings)
 
-            focal_length = subarray.tel[tel_id].optics.equivalent_focal_length
-
-            if isinstance(
-                moments, CameraHillasParametersContainer
-            ):  # Image parameters are in CameraFrame
-
-                # we just need any point on the main shower axis a bit away from the cog
-                p2_x = moments.x + 0.1 * self._cam_radius_m[camera] * np.cos(
-                    moments.psi
-                )
-                p2_y = moments.y + 0.1 * self._cam_radius_m[camera] * np.sin(
-                    moments.psi
-                )
-
-                camera_frame = CameraFrame(
-                    focal_length=focal_length, telescope_pointing=pointing
-                )
-
-                cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=camera_frame)
-                p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=camera_frame)
-
-            else:  # Image parameters are already in TelescopeFrame
-
-                # we just need any point on the main shower axis a bit away from the cog
-                p2_delta_alt = moments.fov_lat + 0.1 * self._cam_radius_deg[
-                    camera
-                ] * np.sin(moments.psi)
-                p2_delta_az = moments.fov_lon + 0.1 * self._cam_radius_deg[
-                    camera
-                ] * np.cos(moments.psi)
-
-                telescope_frame = TelescopeFrame(telescope_pointing=pointing)
-
-                cog_coord = SkyCoord(
-                    fov_lon=moments.fov_lon,
-                    fov_lat=moments.fov_lat,
-                    frame=telescope_frame,
-                )
-                p2_coord = SkyCoord(
-                    fov_lon=p2_delta_az, fov_lat=p2_delta_alt, frame=telescope_frame
-                )
-
-            # Coordinates in the sky
-            cog_coord = cog_coord.transform_to(horizon_frame)
-            p2_coord = p2_coord.transform_to(horizon_frame)
-
-            # re-project from sky to a "fake"-parallel-pointing telescope
-            # then recalculate the psi angle in order to be able to project
-            # it on the ground
-            # This is done to bypass divergent pointing or mispointing
-
-            camera_frame_parallel = CameraFrame(
-                focal_length=focal_length, telescope_pointing=array_pointing
+        p2_1 = cog1 + 0.1 * cam_radius * np.cos(psi)
+        p2_2 = cog2 + 0.1 * cam_radius * np.sin(psi)
+
+        if hillas_in_camera_frame:
+            cog_coord = SkyCoord(x=cog1, y=cog2, unit=u.m, frame=camera_frame)
+            p2_coord = SkyCoord(x=p2_1, y=p2_2, unit=u.m, frame=camera_frame)
+        else:
+            p2_coord = SkyCoord(
+                fov_lon=p2_1, fov_lat=p2_2, unit=u.deg, frame=telescope_frame
             )
-            cog_sky_to_parallel = cog_coord.transform_to(camera_frame_parallel)
-            p2_sky_to_parallel = p2_coord.transform_to(camera_frame_parallel)
-            angle_psi_corr = np.arctan2(
-                cog_sky_to_parallel.y - p2_sky_to_parallel.y,
-                cog_sky_to_parallel.x - p2_sky_to_parallel.x,
+            cog_coord = SkyCoord(
+                fov_lon=cog1, fov_lat=cog2, unit=u.deg, frame=telescope_frame
             )
 
-            corrected_angle_dict[tel_id] = angle_psi_corr
+        cog_coord = cog_coord.transform_to(altaz)
+        p2_coord = p2_coord.transform_to(altaz)
+        cog_cart = altaz_to_righthanded_cartesian(alt=cog_coord.alt, az=cog_coord.az)
+        p2_cart = altaz_to_righthanded_cartesian(alt=p2_coord.alt, az=p2_coord.az)
 
-            circle = HillasPlane(
-                p1=cog_coord,
-                p2=p2_coord,
-                telescope_position=subarray.positions[tel_id],
-                weight=moments.intensity * (moments.length / moments.width),
-            )
-            hillas_planes[tel_id] = circle
+        pseudo_nominal = CameraFrame(
+            telescope_pointing=array_pointing, focal_length=focal_length
+        )
+        cog_cam = cog_coord.transform_to(pseudo_nominal)
+        p2_cam = p2_coord.transform_to(pseudo_nominal)
+
+        corrected_psi = np.arctan(
+            (cog_cam.y.to_value(u.m) - p2_cam.y.to_value(u.m))
+            / (cog_cam.x.to_value(u.m) - p2_cam.x.to_value(u.m))
+        )
 
-        return hillas_planes, corrected_angle_dict
+        return (
+            tel_ids,
+            cog_cart,
+            p2_cart,
+            corrected_psi,
+            weights,
+            telescope_positions,
+            array_pointing,
+        )
 
-    def estimate_direction(self, hillas_planes):
+    @staticmethod
+    def estimate_direction(norm, weight):
         """calculates the origin of the gamma as the weighted average
         direction of the intersections of all hillas planes
 
@@ -410,35 +341,23 @@ class HillasReconstructor(Reconstructor):
         crossings : shape (n,3) list
             an error estimate
         """
+        index_a, index_b = np.array(list(combinations(range(len(norm)), 2))).T
+        crossings = np.cross(norm[index_a], norm[index_b])
+        mask = crossings[:, 2] < 0
+        crossings[mask] = -crossings[mask]
 
-        crossings = []
-        for perm in combinations(hillas_planes.values(), 2):
-            n1, n2 = perm[0].norm, perm[1].norm
-            # cross product automatically weighs in the angle between
-            # the two vectors: narrower angles have less impact,
-            # perpendicular vectors have the most
-            crossing = np.cross(n1, n2)
-
-            # two great circles cross each other twice (one would be
-            # the origin, the other one the direction of the gamma) it
-            # doesn't matter which we pick but it should at least be
-            # consistent: make sure to always take the "upper" solution
-            if crossing[2] < 0:
-                crossing *= -1
-            crossings.append(crossing * perm[0].weight * perm[1].weight)
-
-        result = normalise(np.sum(crossings, axis=0))
-        off_angles = [angle(result, cross) for cross in crossings] * u.rad
-
-        err_est_dir = np.average(
-            off_angles, weights=[len(cross) for cross in crossings]
+        result = np.average(
+            crossings, weights=weight[index_a] * weight[index_b], axis=0
         )
+        result = normalise(result)
+
+        off_angles = angle(result, crossings)
+        err_est_dir = np.mean(off_angles)
 
-        return result, err_est_dir
+        return result, u.Quantity(np.rad2deg(err_est_dir), u.deg)
 
-    def estimate_core_position(
-        self, event, hillas_dict, array_pointing, corrected_angle_dict, hillas_planes
-    ):
+    @staticmethod
+    def estimate_core_position(array_pointing, psi, positions):
         """
         Estimate the core position by intersection the major ellipse lines of each telescope.
 
@@ -469,39 +388,25 @@ class HillasReconstructor(Reconstructor):
         # of both pointing and cleaning/parametrization frame.
         # This angle will be used to visualize the telescope-wise directions of
         # the shower core the ground.
-        psi_core = corrected_angle_dict
-
-        # Record these values
-        for tel_id in hillas_dict.keys():
-            event.dl1.tel[tel_id].parameters.core.psi = psi_core[tel_id]
-
-        # Transform them for numpy
-        psi = u.Quantity(list(psi_core.values()))
 
         # Estimate the position of the shower's core
         # from the TiltedFram to the GroundFrame
 
         z = np.zeros(len(psi))
-        uvw_vectors = np.column_stack([np.cos(psi).value, np.sin(psi).value, z])
+        uvw_vectors = np.column_stack([np.cos(psi), np.sin(psi), z])
 
         tilted_frame = TiltedGroundFrame(pointing_direction=array_pointing)
-        ground_frame = GroundFrame()
-
-        positions = [
-            (
-                SkyCoord(*plane.pos, frame=ground_frame)
-                .transform_to(tilted_frame)
-                .cartesian.xyz
-            )
-            for plane in hillas_planes.values()
-        ]
-
-        core_position = line_line_intersection_3d(uvw_vectors, positions)
+        positions_tilted = positions.transform_to(tilted_frame).cartesian.xyz.T
 
+        core_position = line_line_intersection_3d(
+            uvw_vectors,
+            positions_tilted.to_value(u.m),
+        )
+        core_position = u.Quantity(core_position, u.m, copy=False)
         core_pos_tilted = SkyCoord(
-            x=core_position[0] * u.m,
-            y=core_position[1] * u.m,
-            z=0 * u.m,
+            x=core_position[0],
+            y=core_position[1],
+            z=u.Quantity(0.0, u.m),
             frame=tilted_frame,
         )
 
@@ -509,7 +414,8 @@ class HillasReconstructor(Reconstructor):
 
         return core_pos_ground, core_pos_tilted
 
-    def estimate_h_max(self, hillas_planes):
+    @staticmethod
+    def estimate_h_max(cog_vectors, positions):
         """
         Estimate the max height by intersecting the lines of the cog directions of each telescope.
 
@@ -518,8 +424,9 @@ class HillasReconstructor(Reconstructor):
         astropy.unit.Quantity
             the estimated max height
         """
-        uvw_vectors = np.array([plane.a for plane in hillas_planes.values()])
-        positions = [plane.pos for plane in hillas_planes.values()]
-
+        positions = positions.cartesian.xyz.T.to_value(u.m)
         # not sure if its better to return the length of the vector of the z component
-        return np.linalg.norm(line_line_intersection_3d(uvw_vectors, positions)) * u.m
+        return u.Quantity(
+            np.linalg.norm(line_line_intersection_3d(cog_vectors, positions)),
+            u.m,
+        )


Test patch:
diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py
index f34ea709..e76e3d42 100644
--- a/ctapipe/reco/tests/test_HillasReconstructor.py
+++ b/ctapipe/reco/tests/test_HillasReconstructor.py
@@ -8,93 +8,51 @@ from traitlets.config import Config
 
 from ctapipe.calib import CameraCalibrator
 from ctapipe.containers import HillasParametersContainer
+from ctapipe.coordinates import GroundFrame, altaz_to_righthanded_cartesian
 from ctapipe.image.image_processor import ImageProcessor
-from ctapipe.instrument import SubarrayDescription
 from ctapipe.io import SimTelEventSource
-from ctapipe.reco.hillas_reconstructor import HillasPlane, HillasReconstructor
+from ctapipe.reco.hillas_reconstructor import HillasReconstructor
 from ctapipe.utils import get_dataset_path
 
 
-def test_estimator_results(prod5_sst):
-    """
-    creating some planes pointing in different directions (two
-    north-south, two east-west) and that have a slight position errors (+-
-    0.1 m in one of the four cardinal directions"""
-    horizon_frame = AltAz()
-
-    p1 = SkyCoord(alt=43 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=47 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    circle1 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, 1, 0] * u.m)
-
-    p1 = SkyCoord(alt=44 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=46 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    circle2 = HillasPlane(p1=p1, p2=p2, telescope_position=[1, 0, 0] * u.m)
-
-    p1 = SkyCoord(alt=44.5 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=46.5 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    circle3 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, -1, 0] * u.m)
-
-    p1 = SkyCoord(alt=43.5 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=45.5 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    circle4 = HillasPlane(p1=p1, p2=p2, telescope_position=[-1, 0, 0] * u.m)
-
-    # Create a dummy subarray
-    # (not used here, but required to initialize the reconstructor)
-    subarray = SubarrayDescription(
-        "test array",
-        tel_positions={1: np.zeros(3) * u.m},
-        tel_descriptions={1: prod5_sst},
-    )
+def test_estimator_results():
+    """Test direction results from cog and point 2 (along the main axis) coordinates"""
 
-    # creating the fit class and setting the the great circle member
-    fit = HillasReconstructor(subarray)
-    hillas_planes = {1: circle1, 2: circle2, 3: circle3, 4: circle4}
+    weight = np.ones(4)
+    cog_alt = [43, 44, 44.5, 43.5] * u.deg
+    cog_az = [45, 90, 45, 90] * u.deg
 
-    # performing the direction fit with the minimisation algorithm
-    # and a seed that is perpendicular to the up direction
-    dir_fit_minimise, _ = fit.estimate_direction(hillas_planes)
-    print("direction fit test minimise:", dir_fit_minimise)
+    p2_alt = [47, 45, 46.5, 45.5] * u.deg
+    p2_az = [45, 90, 45, 90] * u.deg
 
+    cog_cart = altaz_to_righthanded_cartesian(cog_alt, cog_az)
+    p2_cart = altaz_to_righthanded_cartesian(p2_alt, p2_az)
+    norm = np.cross(cog_cart, p2_cart)
 
-def test_h_max_results(example_subarray):
-    """
-    creating some planes pointing in different directions (two
-    north-south, two east-west) and that have a slight position errors (+-
-    0.1 m in one of the four cardinal directions"""
-    horizon_frame = AltAz()
+    direction, _ = HillasReconstructor.estimate_direction(norm, weight)
+    assert np.allclose(direction, [0, 0, 1])
 
-    p1 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    circle1 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, 1, 0] * u.m)
 
-    p1 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    circle2 = HillasPlane(p1=p1, p2=p2, telescope_position=[1, 0, 0] * u.m)
+def test_h_max_results():
+    """test h_max estimation from cog coordinates"""
+    cog_alt = [45.0, 45.0, 45.0, 45.0] * u.deg
+    cog_az = [180, 90, 0, -90.0] * u.deg
 
-    p1 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=0 * u.deg, az=45 * u.deg, frame=horizon_frame)
-    circle3 = HillasPlane(p1=p1, p2=p2, telescope_position=[0, -1, 0] * u.m)
-
-    p1 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    p2 = SkyCoord(alt=0 * u.deg, az=90 * u.deg, frame=horizon_frame)
-    circle4 = HillasPlane(p1=p1, p2=p2, telescope_position=[-1, 0, 0] * u.m)
+    positions = [[1, 0, 0], [0, 1, 0], [-1, 0, 0], [0, -1, 0]] * u.km
+    positions = SkyCoord(
+        x=positions[:, 0], y=positions[:, 1], z=positions[:, 2], frame=GroundFrame()
+    )
 
-    # Create a dummy subarray
-    # (not used here, but required to initialize the reconstructor)
-    subarray = example_subarray
+    cog_cart = altaz_to_righthanded_cartesian(cog_alt, cog_az)
 
     # creating the fit class and setting the the great circle member
-    fit = HillasReconstructor(subarray)
-    hillas_planes = {1: circle1, 2: circle2, 3: circle3, 4: circle4}
 
     # performing the direction fit with the minimisation algorithm
     # and a seed that is perpendicular to the up direction
-    h_max_reco = fit.estimate_h_max(hillas_planes)
-    print("h max fit test minimise:", h_max_reco)
+    h_max_reco = HillasReconstructor.estimate_h_max(cog_cart, positions)
 
     # the results should be close to the direction straight up
-    np.testing.assert_allclose(h_max_reco.value, 0, atol=1e-8)
-    # np.testing.assert_allclose(fitted_core_position.value, [0, 0], atol=1e-3)
+    assert u.isclose(h_max_reco, 1 * u.km)
 
 
 def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev):
@@ -162,7 +120,9 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev):
     assert not result.is_valid
 
 
-def test_reconstruction_against_simulation(subarray_and_event_gamma_off_axis_500_gev):
+def test_reconstruction_against_simulation_telescope_frame(
+    subarray_and_event_gamma_off_axis_500_gev,
+):
     """Reconstruction is here done only in the TelescopeFrame,
     since the previous tests test already for the compatibility between
     frames"""
@@ -192,6 +152,45 @@ def test_reconstruction_against_simulation(subarray_and_event_gamma_off_axis_500
     # check that we are not more far than 0.1 degrees
     assert reco_coord.separation(true_coord) < 0.1 * u.deg
 
+    assert u.isclose(result.core_x, event.simulation.shower.core_x, atol=25 * u.m)
+    assert u.isclose(result.core_y, event.simulation.shower.core_y, atol=25 * u.m)
+
+
+def test_reconstruction_against_simulation_camera_frame(
+    subarray_and_event_gamma_off_axis_500_gev,
+):
+    """Reconstruction is here done only in the CameraFrame,
+    since the previous tests test already for the compatibility between
+    frames"""
+
+    # 4-LST bright event already calibrated
+    # we'll clean it and parametrize it again in the TelescopeFrame
+    subarray, event = subarray_and_event_gamma_off_axis_500_gev
+
+    # define reconstructor
+    calib = CameraCalibrator(subarray)
+    image_processor = ImageProcessor(subarray, use_telescope_frame=False)
+    reconstructor = HillasReconstructor(subarray)
+
+    # Get shower geometry
+    calib(event)
+    image_processor(event)
+    reconstructor(event)
+    result = event.dl2.stereo.geometry[reconstructor.__class__.__name__]
+
+    # get the reconstructed coordinates in the sky
+    reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz())
+    # get the simulated coordinates in the sky
+    true_coord = SkyCoord(
+        alt=event.simulation.shower.alt, az=event.simulation.shower.az, frame=AltAz()
+    )
+
+    # check that we are not more far than 0.1 degrees
+    assert reco_coord.separation(true_coord) < 0.1 * u.deg
+
+    assert u.isclose(result.core_x, event.simulation.shower.core_x, atol=25 * u.m)
+    assert u.isclose(result.core_y, event.simulation.shower.core_y, atol=25 * u.m)
+
 
 @pytest.mark.parametrize(
     "filename",
@@ -261,13 +260,32 @@ def test_CameraFrame_against_TelescopeFrame(filename):
             for field, cam in result_camera_frame.items():
                 tel = getattr(result_telescope_frame, field)
 
+                kwargs = dict(rtol=6e-3, equal_nan=True)
+
                 if hasattr(cam, "unit"):
+                    if cam.value == 0 or tel.value == 0:
+                        kwargs["atol"] = 1e-6 * cam.unit
                     assert u.isclose(
-                        cam, tel, rtol=1e-3, atol=1e-3 * tel.unit, equal_nan=True
-                    )
+                        cam, tel, **kwargs
+                    ), f"attr {field} not matching, camera: {result_camera_frame!s} telescope: {result_telescope_frame!s}"
                 elif isinstance(cam, list):
                     assert cam == tel
                 else:
-                    assert np.isclose(cam, tel, rtol=1e-3, atol=1e-3, equal_nan=True)
+
+                    if cam == 0 or tel == 0:
+                        kwargs["atol"] = 1e-6
+                    assert np.isclose(cam, tel, **kwargs)
 
     assert reconstructed_events > 0  # check that we reconstruct at least 1 event
+
+
+def test_angle():
+    from ctapipe.reco.hillas_reconstructor import angle
+
+    # test it works with single vectors
+    assert np.isclose(angle(np.array([0, 0, 1]), np.array([1, 0, 0])), np.pi / 2)
+
+    # and with an array of vectors
+    a = np.array([[1, 0, 0], [1, 0, 0]])
+    b = np.array([[1, 0, 0], [0, 1, 0]])
+    assert np.allclose(angle(a, b), [0, np.pi / 2])
diff --git a/ctapipe/reco/tests/test_reconstruction_methods.py b/ctapipe/reco/tests/test_reconstruction_methods.py
index 6acea98a..90cb7278 100644
--- a/ctapipe/reco/tests/test_reconstruction_methods.py
+++ b/ctapipe/reco/tests/test_reconstruction_methods.py
@@ -18,7 +18,6 @@ def test_reconstructors(reconstructors):
     a test of the complete fit procedure on one event including:
     • tailcut cleaning
     • hillas parametrisation
-    • HillasPlane creation
     • direction fit
     • position fit