Fix h_max definition

ctapipe/containers.py

@@ -920,12 +920,17 @@ 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,
         (
-            "direction validity flag. True if the shower direction"
+            "Geometry validity flag. True if the shower geometry"
             "was properly reconstructed by the algorithm"
         ),
     )

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):
@@ -260,7 +261,8 @@ 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(
+
+        h_max = self.reconstruct_h_max(
             nom.fov_lon,
             nom.fov_lat,
             tilt.x,
@@ -287,8 +289,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__,
         )
@@ -309,7 +311,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)
@@ -381,7 +383,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()
@@ -437,7 +440,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
     ):
         """
@@ -495,25 +498,19 @@ def reconstruct_xmax(
             np.array(ty),
         )
         weight = np.array(amp)
-        mean_height = 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
-        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 += 2100  # TODO: replace with instrument info
-
-        if mean_height > 100000 or np.isnan(mean_height):
-            mean_height = 100000
+        mean_height += self.subarray.reference_location.geodetic.height.to_value(u.m)
 
-        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)
+        if mean_height > H_MAX_UPPER_LIMIT_M:
+            mean_height = np.nan
 
-        return mean_height
+        return u.Quantity(mean_height, u.m)
 
     @staticmethod
     def intersect_lines(xp1, yp1, phi1, xp2, yp2, phi2):

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,10 @@ 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)
+        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)
@@ -412,18 +416,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
+        astropy.unit.Quantity:
+            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

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)

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_xmax(
+    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):

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.