Repository: cta-observatory/ctapipe (python)

Issue:
Changes to simulation containers
**Please describe the use case that requires this feature.**

In a discussion with @kbernloehr in the gitlab, it came up that the `injection_height` in the `mcrunheader` object, which we store in the `SimulationConfigContainer` is not used by simtel array and exists only for historical reasons, it is always -1.

The real injection height, is stored for each shower in the `MCShower` object, but currently we don't store this.

**Describe the solution you'd like**

Remove `injection_height` from the `SimulationConfigContainer` (or at least set it to nan or None in `SimtelEventSource`) and add `injection_height` to the `SimulatedShowerContainer`.




Pull request:
Remove SimulationConfigContainer.injection_height, add SimulatedShowercontainer.starting_grammage
Fixes #2335 

Patch:
diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index bcccaef7..ba8bddf8 100644
--- a/ctapipe/containers.py
+++ b/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"
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 5dd00a22..235f8edd 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -88,6 +88,7 @@ SIMTEL_TO_CTA_EVENT_TYPE = {
 
 _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 @@ class SimTelEventSource(EventSource):
             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 @@ class SimTelEventSource(EventSource):
             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
+            ),
         )
diff --git a/docs/changes/2343.feature.rst b/docs/changes/2343.feature.rst
new file mode 100644
index 00000000..6367e808
--- /dev/null
+++ b/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.


Test patch:
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index 04cf133d..e9e7ea31 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/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