Merge pull request #2343 from cta-observatory/starting_grammage

Remove SimulationConfigContainer.injection_height, add SimulatedShowercontainer.starting_grammage

ctapipe/containers.py

@@ -586,8 +586,11 @@ class SimulatedShowerContainer(Container):
     core_x = Field(nan * u.m, "Simulated core position (x)", unit=u.m)
     core_y = Field(nan * u.m, "Simulated core position (y)", unit=u.m)
     h_first_int = Field(nan * u.m, "Height of first interaction", unit=u.m)
-    x_max = Field(
-        nan * u.g / (u.cm**2), "Simulated Xmax value", unit=u.g / (u.cm**2)
+    x_max = Field(nan * u.g / u.cm**2, "Simulated Xmax value", unit=u.g / u.cm**2)
+    starting_grammage = Field(
+        nan * u.g / u.cm**2,
+        "Grammage (mass overburden) where the particle was injected into the atmosphere",
+        unit=u.g / u.cm**2,
     )
     shower_primary_id = Field(
         np.int16(np.iinfo(np.int16).max),
@@ -692,9 +695,6 @@ class SimulationConfigContainer(Container):
     core_pos_mode = Field(
         nan, description="Core Position Mode (0=Circular, 1=Rectangular)"
     )
-    injection_height = Field(
-        nan * u.m, description="Height of particle injection", unit=u.m
-    )
     atmosphere = Field(nan * u.m, description="Atmospheric model number")
     corsika_iact_options = Field(
         nan, description="CORSIKA simulation options for IACTs"

ctapipe/io/simteleventsource.py

@@ -88,6 +88,7 @@
 
 _half_pi = 0.5 * np.pi
 _half_pi_maxval = (1 + 1e-6) * _half_pi
+_float32_nan = np.float32(np.nan)
 
 
 def _clip_altitude_if_close(altitude):
@@ -117,7 +118,7 @@ class MirrorClass(enum.Enum):
     DUAL_MIRROR = 2
 
 
-X_MAX_UNIT = u.g / (u.cm**2)
+GRAMMAGE_UNIT = u.g / (u.cm**2)
 NANOSECONDS_PER_DAY = (1 * u.day).to_value(u.ns)
 
 
@@ -990,7 +991,6 @@ def _parse_simulation_header(self):
             max_scatter_range=mc_run_head["core_range"][1] * u.m,
             min_scatter_range=mc_run_head["core_range"][0] * u.m,
             core_pos_mode=mc_run_head["core_pos_mode"],
-            injection_height=mc_run_head["injection_height"] * u.m,
             atmosphere=mc_run_head["atmosphere"],
             corsika_iact_options=mc_run_head["corsika_iact_options"],
             corsika_low_E_model=mc_run_head["corsika_low_E_model"],
@@ -1053,6 +1053,9 @@ def _fill_simulated_event_information(array_event):
             core_x=u.Quantity(mc_event["xcore"], u.m),
             core_y=u.Quantity(mc_event["ycore"], u.m),
             h_first_int=u.Quantity(mc_shower["h_first_int"], u.m),
-            x_max=u.Quantity(mc_shower["xmax"], X_MAX_UNIT),
+            x_max=u.Quantity(mc_shower["xmax"], GRAMMAGE_UNIT),
             shower_primary_id=mc_shower["primary_id"],
+            starting_grammage=u.Quantity(
+                mc_shower.get("depth_start", _float32_nan), GRAMMAGE_UNIT
+            ),
         )

ctapipe/io/tests/test_simteleventsource.py

@@ -591,3 +591,11 @@ def test_float32_pihalf(sign):
     assert shower.alt.value == sign * np.pi / 2
     # check we cana create a Latitude:
     Latitude(shower.alt.value, u.rad)
+
+
+def test_starting_grammage():
+    path = "dataset://lst_muons.simtel.zst"
+
+    with SimTelEventSource(path, focal_length_choice="EQUIVALENT") as source:
+        e = next(iter(source))
+        assert e.simulation.shower.starting_grammage == 580 * u.g / u.cm**2

docs/changes/2343.feature.rst

@@ -0,0 +1,5 @@
+Remove the ``injection_height`` field from the ``SimulationConfigContainer``,
+this field was always empty and is never filled by ``sim_telarray``.
+
+Add the corresponding ``starting_grammage`` field to the ``SimulatedShowerContainer``,
+where it is actually available.