From ce28ea28ec6c3cfaf310972230e72581264eda63 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Wed, 20 Sep 2023 14:58:49 +0200
Subject: [PATCH 01/14] make definition of h_max consistent

---
 ctapipe/reco/hillas_intersection.py  | 40 +++++++++++++---------------
 ctapipe/reco/hillas_reconstructor.py | 32 +++++++++++++++-------
 2 files changed, 41 insertions(+), 31 deletions(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index d7288b58576..d0d20b9ef7c 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -260,16 +260,20 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
         sky_pos = nom.transform_to(array_pointing.frame)
         tilt = SkyCoord(x=core_x * u.m, y=core_y * u.m, z=0 * u.m, frame=tilted_frame)
         grd = project_to_ground(tilt)
-        x_max = self.reconstruct_xmax(
-            nom.fov_lon,
-            nom.fov_lat,
-            tilt.x,
-            tilt.y,
-            hillas_dict_mod,
-            tel_x,
-            tel_y,
-            90 * u.deg - array_pointing.alt,
-        )
+
+        if self.subarray.reference_location:
+            h_max = self.reconstruct_h_max(
+                nom.fov_lon,
+                nom.fov_lat,
+                tilt.x,
+                tilt.y,
+                hillas_dict_mod,
+                tel_x,
+                tel_y,
+                90 * u.deg - array_pointing.alt,
+            )
+        else:
+            h_max = u.Quantity(np.nan, u.m)
 
         src_error = np.sqrt(err_fov_lon**2 + err_fov_lat**2)
 
@@ -287,8 +291,8 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
             is_valid=True,
             alt_uncert=src_error.to(u.rad),
             az_uncert=src_error.to(u.rad),
-            h_max=x_max,
-            h_max_uncert=u.Quantity(np.nan * x_max.unit),
+            h_max=h_max,
+            h_max_uncert=u.Quantity(np.nan * h_max.unit),
             goodness_of_fit=np.nan,
             prefix=self.__class__.__name__,
         )
@@ -437,7 +441,7 @@ def reconstruct_tilted(self, hillas_parameters, tel_x, tel_y):
 
         return x_pos, y_pos, np.sqrt(var_x), np.sqrt(var_y)
 
-    def reconstruct_xmax(
+    def reconstruct_h_max(
         self, source_x, source_y, core_x, core_y, hillas_parameters, tel_x, tel_y, zen
     ):
         """
@@ -502,18 +506,12 @@ def reconstruct_xmax(
         mean_height *= np.cos(zen)
 
         # Add on the height of the detector above sea level
-        mean_height += 2100  # TODO: replace with instrument info
+        mean_height += self.subarray.reference_location.geodetic.height.to_value(u.m)
 
         if mean_height > 100000 or np.isnan(mean_height):
             mean_height = 100000
 
-        mean_height *= u.m
-        # Lookup this height in the depth tables, the convert Hmax to Xmax
-        # x_max = self.thickness_profile(mean_height.to(u.km))
-        # Convert to slant depth
-        # x_max /= np.cos(zen)
-
-        return mean_height
+        return u.Quantity(mean_height, u.m)
 
     @staticmethod
     def intersect_lines(xp1, yp1, phi1, xp2, yp2, phi2):
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index 1fbd305d7c2..ad45809e56b 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -19,6 +19,7 @@
     altaz_to_righthanded_cartesian,
     project_to_ground,
 )
+from ..instrument import SubarrayDescription
 from .reconstructor import (
     HillasGeometryReconstructor,
     InvalidWidthException,
@@ -106,7 +107,7 @@ class that reconstructs the direction of an atmospheric shower
 
     """
 
-    def __init__(self, subarray, **kwargs):
+    def __init__(self, subarray: SubarrayDescription, **kwargs):
         super().__init__(subarray=subarray, **kwargs)
         _cam_radius_m = {
             cam: cam.geometry.guess_radius().to_value(u.m)
@@ -181,7 +182,15 @@ def __call__(self, event):
         _, lat, lon = cartesian_to_spherical(*direction)
 
         # estimate max height of shower
-        h_max = self.estimate_h_max(cog_cartesian, telescope_positions)
+        if self.subarray.reference_location:
+            h_max = (
+                HillasReconstructor.estimate_relative_h_max(
+                    cog_cartesian, telescope_positions
+                )
+                + self.subarray.reference_location.geodetic.height
+            )
+        else:
+            h_max = u.Quantity(np.nan, u.m)
 
         # az is clockwise, lon counter-clockwise, make sure it stays in [0, 2pi)
         az = Longitude(-lon)
@@ -411,19 +420,22 @@ def estimate_core_position(array_pointing, psi, positions):
 
         return core_pos_ground, core_pos_tilted
 
-    @staticmethod
-    def estimate_h_max(cog_vectors, positions):
-        """
-        Estimate the max height by intersecting the lines of the cog directions of each telescope.
+    def estimate_relative_h_max(cog_vectors, positions):
+        """Estimate the relative (to the observatory) vertical height of
+        shower-max by intersecting the lines of the cog directions of each
+        telescope.
 
         Returns
         -------
         astropy.unit.Quantity
-            the estimated max height
+            the estimated height above observatory level (not sea level) of the
+            shower-max point
+
         """
         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 u.Quantity(
-            np.linalg.norm(line_line_intersection_3d(cog_vectors, positions)),
+        shower_max = u.Quantity(
+            line_line_intersection_3d(cog_vectors, positions),
             u.m,
         )
+
+        return shower_max[2]  # the z-coordinate only

From d292adc1995c58580b61603a69208b4495957177 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Wed, 20 Sep 2023 14:59:04 +0200
Subject: [PATCH 02/14] get tests to run

---
 ctapipe/conftest.py                            |  4 ++--
 ctapipe/reco/tests/test_HillasReconstructor.py | 12 +++---------
 ctapipe/reco/tests/test_hillas_intersection.py |  2 +-
 3 files changed, 6 insertions(+), 12 deletions(-)

diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index e40abb42d72..f4f57e4dad7 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -79,11 +79,11 @@ def subarray_prod3_paranal():
 
 
 @pytest.fixture(scope="session")
-def example_subarray(subarray_prod3_paranal):
+def example_subarray(subarray_prod5_paranal):
     """
     Subarray corresponding to the example event
     """
-    return subarray_prod3_paranal
+    return subarray_prod5_paranal
 
 
 @pytest.fixture(scope="function")
diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py
index 9fe7c57d4de..fd8ca438495 100644
--- a/ctapipe/reco/tests/test_HillasReconstructor.py
+++ b/ctapipe/reco/tests/test_HillasReconstructor.py
@@ -44,12 +44,9 @@ def test_h_max_results():
     )
 
     cog_cart = altaz_to_righthanded_cartesian(cog_alt, cog_az)
-
-    # creating the fit class and setting the the great circle member
-
-    # performing the direction fit with the minimisation algorithm
-    # and a seed that is perpendicular to the up direction
-    h_max_reco = HillasReconstructor.estimate_h_max(cog_cart, positions)
+    h_max_reco = HillasReconstructor.estimate_relative_h_max(
+        cog_vectors=cog_cart, positions=positions
+    )
 
     # the results should be close to the direction straight up
     assert u.isclose(h_max_reco, 1 * u.km)
@@ -210,7 +207,6 @@ def test_reconstruction_against_simulation_camera_frame(
     ],
 )
 def test_CameraFrame_against_TelescopeFrame(filename):
-
     input_file = get_dataset_path(filename)
     # "gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz"
     # )
@@ -243,7 +239,6 @@ def test_CameraFrame_against_TelescopeFrame(filename):
     reconstructed_events = 0
 
     for event_telescope_frame in source:
-
         calib(event_telescope_frame)
         # make a copy of the calibrated event for the camera frame case
         # later we clean and paramretrize the 2 events in the same way
@@ -281,7 +276,6 @@ def test_CameraFrame_against_TelescopeFrame(filename):
                 elif isinstance(cam, list):
                     assert cam == tel
                 else:
-
                     if cam == 0 or tel == 0:
                         kwargs["atol"] = 1e-6
                     assert np.isclose(cam, tel, **kwargs)
diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py
index 94492aba3d6..d9f2f3318ea 100644
--- a/ctapipe/reco/tests/test_hillas_intersection.py
+++ b/ctapipe/reco/tests/test_hillas_intersection.py
@@ -83,7 +83,7 @@ def test_intersection_xmax_reco(example_subarray):
         ),
     }
 
-    x_max = hill_inter.reconstruct_xmax(
+    x_max = hill_inter.reconstruct_h_max(
         source_x=nom_pos_reco.fov_lon,
         source_y=nom_pos_reco.fov_lat,
         core_x=0 * u.m,

From 258a98f2465a3b331c76af8a45a59b3652b66588 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Thu, 21 Sep 2023 16:42:34 +0200
Subject: [PATCH 03/14] revert change to conftest for now

---
 ctapipe/conftest.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index f4f57e4dad7..e40abb42d72 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -79,11 +79,11 @@ def subarray_prod3_paranal():
 
 
 @pytest.fixture(scope="session")
-def example_subarray(subarray_prod5_paranal):
+def example_subarray(subarray_prod3_paranal):
     """
     Subarray corresponding to the example event
     """
-    return subarray_prod5_paranal
+    return subarray_prod3_paranal
 
 
 @pytest.fixture(scope="function")

From 5d742bece5a7db5685b90a91dd299c4a155e9634 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Thu, 21 Sep 2023 16:47:01 +0200
Subject: [PATCH 04/14] deal with case where there is no reference_location

Will be fixed in another PR (should always have at least the observatory height)
---
 ctapipe/reco/hillas_intersection.py            | 18 +++++++++++++-----
 ctapipe/reco/hillas_reconstructor.py           |  2 +-
 ctapipe/reco/tests/test_hillas_intersection.py |  9 +++++----
 3 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index d0d20b9ef7c..f834ad85214 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -261,7 +261,7 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
         tilt = SkyCoord(x=core_x * u.m, y=core_y * u.m, z=0 * u.m, frame=tilted_frame)
         grd = project_to_ground(tilt)
 
-        if self.subarray.reference_location:
+        if self.subarray.reference_location is not None:
             h_max = self.reconstruct_h_max(
                 nom.fov_lon,
                 nom.fov_lat,
@@ -499,17 +499,25 @@ def reconstruct_h_max(
             np.array(ty),
         )
         weight = np.array(amp)
-        mean_height = np.sum(height * weight) / np.sum(weight)
+        mean_distance = np.sum(height * weight) / np.sum(weight)
 
         # This value is height above telescope in the tilted system,
         # we should convert to height above ground
-        mean_height *= np.cos(zen)
+        mean_height = mean_distance * np.cos(zen.to_value(u.rad))
 
         # Add on the height of the detector above sea level
-        mean_height += self.subarray.reference_location.geodetic.height.to_value(u.m)
+        if self.subarray.reference_location is not None:
+            mean_height += self.subarray.reference_location.geodetic.height.to_value(
+                u.m
+            )
+        else:
+            # FIXME: Can remoev this check once we ensure the reference_location is always loaded
+            warnings.warn(
+                "Computing h_max with no reference location. Height will be wrong."
+            )
 
         if mean_height > 100000 or np.isnan(mean_height):
-            mean_height = 100000
+            mean_height = np.nan
 
         return u.Quantity(mean_height, u.m)
 
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index ad45809e56b..6bbea1da7e5 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -182,7 +182,7 @@ def __call__(self, event):
         _, lat, lon = cartesian_to_spherical(*direction)
 
         # estimate max height of shower
-        if self.subarray.reference_location:
+        if self.subarray.reference_location is not None:
             h_max = (
                 HillasReconstructor.estimate_relative_h_max(
                     cog_cartesian, telescope_positions
diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py
index d9f2f3318ea..924699b74ed 100644
--- a/ctapipe/reco/tests/test_hillas_intersection.py
+++ b/ctapipe/reco/tests/test_hillas_intersection.py
@@ -30,8 +30,8 @@ def test_intersect():
 
 def test_parallel():
     """
-    Simple test to check the intersection of lines. Try to intersect positions at (0,0) and (0,1)
-    with angles parallel and check the behaviour
+    ?    Simple test to check the intersection of lines. Try to intersect positions at (0,0) and (0,1)
+        with angles parallel and check the behaviour
     """
     x1 = 0
     y1 = 0
@@ -83,7 +83,7 @@ def test_intersection_xmax_reco(example_subarray):
         ),
     }
 
-    x_max = hill_inter.reconstruct_h_max(
+    h_max = hill_inter.reconstruct_h_max(
         source_x=nom_pos_reco.fov_lon,
         source_y=nom_pos_reco.fov_lat,
         core_x=0 * u.m,
@@ -93,7 +93,8 @@ def test_intersection_xmax_reco(example_subarray):
         tel_y={1: (0 * u.m), 2: (150 * u.m)},
         zen=zen_pointing,
     )
-    print(x_max)
+
+    assert h_max > 0 * u.m
 
 
 def test_intersection_reco_impact_point_tilted(example_subarray):

From b26add845411fe576b86bbd6b9fa685857ced891 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Thu, 28 Sep 2023 13:01:04 +0200
Subject: [PATCH 05/14] removed checks for reference_location (not needed)

---
 ctapipe/reco/hillas_intersection.py  | 10 +---------
 ctapipe/reco/hillas_reconstructor.py | 14 +++++---------
 2 files changed, 6 insertions(+), 18 deletions(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index f834ad85214..2a4bfe49c9a 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -506,15 +506,7 @@ def reconstruct_h_max(
         mean_height = mean_distance * np.cos(zen.to_value(u.rad))
 
         # Add on the height of the detector above sea level
-        if self.subarray.reference_location is not None:
-            mean_height += self.subarray.reference_location.geodetic.height.to_value(
-                u.m
-            )
-        else:
-            # FIXME: Can remoev this check once we ensure the reference_location is always loaded
-            warnings.warn(
-                "Computing h_max with no reference location. Height will be wrong."
-            )
+        mean_height += self.subarray.reference_location.geodetic.height.to_value(u.m)
 
         if mean_height > 100000 or np.isnan(mean_height):
             mean_height = np.nan
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index 6bbea1da7e5..b4b7cd3fdf5 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -182,15 +182,10 @@ def __call__(self, event):
         _, lat, lon = cartesian_to_spherical(*direction)
 
         # estimate max height of shower
-        if self.subarray.reference_location is not None:
-            h_max = (
-                HillasReconstructor.estimate_relative_h_max(
-                    cog_cartesian, telescope_positions
-                )
-                + self.subarray.reference_location.geodetic.height
-            )
-        else:
-            h_max = u.Quantity(np.nan, u.m)
+        h_max = (
+            self.estimate_relative_h_max(cog_cartesian, telescope_positions)
+            + self.subarray.reference_location.geodetic.height
+        )
 
         # az is clockwise, lon counter-clockwise, make sure it stays in [0, 2pi)
         az = Longitude(-lon)
@@ -420,6 +415,7 @@ def estimate_core_position(array_pointing, psi, positions):
 
         return core_pos_ground, core_pos_tilted
 
+    @staticmethod
     def estimate_relative_h_max(cog_vectors, positions):
         """Estimate the relative (to the observatory) vertical height of
         shower-max by intersecting the lines of the cog directions of each

From 64426e4ace83ddf0178e0bab1e406687418d455a Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Thu, 28 Sep 2023 13:12:05 +0200
Subject: [PATCH 06/14] added initial changelog

---
 docs/changes/2403.bugfix.rst | 14 ++++++++++++++
 1 file changed, 14 insertions(+)
 create mode 100644 docs/changes/2403.bugfix.rst

diff --git a/docs/changes/2403.bugfix.rst b/docs/changes/2403.bugfix.rst
new file mode 100644
index 00000000000..611abd143ab
--- /dev/null
+++ b/docs/changes/2403.bugfix.rst
@@ -0,0 +1,14 @@
+Fixed the definition of ``h_max``, which was both inconsistent between
+`ctapipe.reco.HillasReconstructor` and `ctapipe.reco.HillasIntersection`
+implementations, and was also incorrect since it was measured from the
+observatory elevation rather than from sea level.
+
+The value of ``h_max`` is now defined as the height above sea level of the
+shower-max point (in meters), not the distance to that point. Therefore it is
+not corrected for the zenith angle of the shower. This is consistent with the
+options currently used for *CORSIKA*, where the *SLANT* option is set to false,
+meaning heights are actual heights not distances from the impact point, and
+``x_max`` is a *depth*, not a *slant depth*. Note that this definition may be
+inconsistent with other observatories where slant-depths are used, and also note
+that the slant depth or distance to shower max are the more useful quantities
+for shower physics.

From 755eeed695213a5d9c8797e2aa816ebca9a709a9 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Thu, 28 Sep 2023 14:57:21 +0200
Subject: [PATCH 07/14] remove one more redundant check for reference loc

---
 ctapipe/reco/hillas_intersection.py | 23 ++++++++++-------------
 1 file changed, 10 insertions(+), 13 deletions(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index 2a4bfe49c9a..2d6d28f172b 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -261,19 +261,16 @@ def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
         tilt = SkyCoord(x=core_x * u.m, y=core_y * u.m, z=0 * u.m, frame=tilted_frame)
         grd = project_to_ground(tilt)
 
-        if self.subarray.reference_location is not None:
-            h_max = self.reconstruct_h_max(
-                nom.fov_lon,
-                nom.fov_lat,
-                tilt.x,
-                tilt.y,
-                hillas_dict_mod,
-                tel_x,
-                tel_y,
-                90 * u.deg - array_pointing.alt,
-            )
-        else:
-            h_max = u.Quantity(np.nan, u.m)
+        h_max = self.reconstruct_h_max(
+            nom.fov_lon,
+            nom.fov_lat,
+            tilt.x,
+            tilt.y,
+            hillas_dict_mod,
+            tel_x,
+            tel_y,
+            90 * u.deg - array_pointing.alt,
+        )
 
         src_error = np.sqrt(err_fov_lon**2 + err_fov_lat**2)
 

From 604fe728a715a05171594beefdf81b8fc7d5ccef Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Thu, 28 Sep 2023 14:57:48 +0200
Subject: [PATCH 08/14] fix bug where return value was incorrect

Returning None here would cause an exception, since the values are
unpacked into a 4-tuple
---
 ctapipe/reco/hillas_intersection.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index 2d6d28f172b..b07afe42947 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -310,7 +310,7 @@ def reconstruct_nominal(self, hillas_parameters):
 
         """
         if len(hillas_parameters) < 2:
-            return None  # Throw away events with < 2 images
+            return (np.nan, np.nan, np.nan, np.nan)  # Throw away events with < 2 images
 
         # Find all pairs of Hillas parameters
         combos = itertools.combinations(list(hillas_parameters.values()), 2)
@@ -382,7 +382,8 @@ def reconstruct_tilted(self, hillas_parameters, tel_x, tel_y):
             core uncertainty X
         """
         if len(hillas_parameters) < 2:
-            return None  # Throw away events with < 2 images
+            return (np.nan, np.nan, np.nan, np.nan)  # Throw away events with < 2 images
+
         hill_list = list()
         tx = list()
         ty = list()

From 9675cc6e96771b6bd4a4c0092fe287bdc81fa45e Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Fri, 29 Sep 2023 15:39:40 +0200
Subject: [PATCH 09/14] better description of h_max

---
 ctapipe/containers.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index cb4282c5d4c..57b3c678ad5 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -920,8 +920,13 @@ class ReconstructedGeometryContainer(Container):
         "reconstructed core position uncertainty along tilted frame Y axis",
         unit=u.m,
     )
-    h_max = Field(nan * u.m, "reconstructed height of the shower maximum", unit=u.m)
+    h_max = Field(
+        nan * u.m,
+        "reconstructed vertical height (above sea level) of the shower maximum",
+        unit=u.m,
+    )
     h_max_uncert = Field(nan * u.m, "uncertainty of h_max", unit=u.m)
+
     is_valid = Field(
         False,
         (

From c2bef1b48bfbb43340894a618032cd02c74ab2ce Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Fri, 29 Sep 2023 15:40:03 +0200
Subject: [PATCH 10/14] fix doctring format

---
 ctapipe/reco/hillas_reconstructor.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index b4b7cd3fdf5..0f82e2bc51d 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -423,7 +423,7 @@ def estimate_relative_h_max(cog_vectors, positions):
 
         Returns
         -------
-        astropy.unit.Quantity
+        astropy.unit.Quantity:
             the estimated height above observatory level (not sea level) of the
             shower-max point
 

From 3c9bd5be6f6aa940a82c4c169069152875389298 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Fri, 29 Sep 2023 15:40:18 +0200
Subject: [PATCH 11/14] use relative references in changelog

---
 docs/changes/2403.bugfix.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/changes/2403.bugfix.rst b/docs/changes/2403.bugfix.rst
index 611abd143ab..042b14c59bd 100644
--- a/docs/changes/2403.bugfix.rst
+++ b/docs/changes/2403.bugfix.rst
@@ -1,5 +1,5 @@
 Fixed the definition of ``h_max``, which was both inconsistent between
-`ctapipe.reco.HillasReconstructor` and `ctapipe.reco.HillasIntersection`
+`~ctapipe.reco.HillasReconstructor` and `~ctapipe.reco.HillasIntersection`
 implementations, and was also incorrect since it was measured from the
 observatory elevation rather than from sea level.
 

From 72505520d585f96c71083e7a64b5e72fd8ab9fc7 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Tue, 3 Oct 2023 11:05:30 +0200
Subject: [PATCH 12/14] fixed phrasing

---
 ctapipe/containers.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index 57b3c678ad5..511be059a0d 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -922,7 +922,7 @@ class ReconstructedGeometryContainer(Container):
     )
     h_max = Field(
         nan * u.m,
-        "reconstructed vertical height (above sea level) of the shower maximum",
+        "reconstructed vertical height above sea level of the shower maximum",
         unit=u.m,
     )
     h_max_uncert = Field(nan * u.m, "uncertainty of h_max", unit=u.m)
@@ -930,7 +930,7 @@ class ReconstructedGeometryContainer(Container):
     is_valid = Field(
         False,
         (
-            "direction validity flag. True if the shower direction"
+            "Geometry validity flag. True if the shower geometry"
             "was properly reconstructed by the algorithm"
         ),
     )

From ad401fc5743e1606f4dceb30a90b946986273fb4 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Mon, 23 Oct 2023 14:31:31 +0200
Subject: [PATCH 13/14] define magic number

---
 ctapipe/reco/hillas_intersection.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index b07afe42947..0529dc7d3d9 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -44,6 +44,7 @@
 )
 
 FOV_ANGULAR_DISTANCE_LIMIT_RAD = (45 * u.deg).to_value(u.rad)
+H_MAX_UPPER_LIMIT_M = 100_000
 
 
 def _far_outside_fov(fov_lat, fov_lon):
@@ -506,7 +507,7 @@ def reconstruct_h_max(
         # Add on the height of the detector above sea level
         mean_height += self.subarray.reference_location.geodetic.height.to_value(u.m)
 
-        if mean_height > 100000 or np.isnan(mean_height):
+        if mean_height > H_MAX_UPPER_LIMIT_M:
             mean_height = np.nan
 
         return u.Quantity(mean_height, u.m)

From 9a1494ff2bc8c2262115344d19dbf67317068b98 Mon Sep 17 00:00:00 2001
From: Karl Kosack <karl.kosack@cea.fr>
Date: Mon, 23 Oct 2023 14:31:46 +0200
Subject: [PATCH 14/14] use np.average

---
 ctapipe/reco/hillas_intersection.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index 0529dc7d3d9..0906cf0ff4b 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -498,7 +498,7 @@ def reconstruct_h_max(
             np.array(ty),
         )
         weight = np.array(amp)
-        mean_distance = np.sum(height * weight) / np.sum(weight)
+        mean_distance = np.average(height, weights=weight)
 
         # This value is height above telescope in the tilted system,
         # we should convert to height above ground