Additional metadata from pyhessio, similar to #655 (#895)

ctapipe/io/containers.py

@@ -270,6 +270,21 @@ class MCHeaderContainer(Container):
         "OR "
         "[0]=R.A., [1]=Declination in mode 1."
     ))
+    corsika_version = Field(0.0, "CORSIKA version * 1000")
+    simtel_version = Field(0.0, "sim_telarray version * 1000")
+    energy_range_min = Field(0.0, "Lower limit of energy range "
+                                  "of primary particle", unit=u.TeV)
+    energy_range_max = Field(0.0, "Upper limit of energy range "
+                                  "of primary particle", unit=u.TeV)
+    prod_site_B_total = Field(0.0, "total geomagnetic field", unit=u.uT)
+    prod_site_B_declination = Field(0.0, "magnetic declination", unit=u.rad)
+    prod_site_B_inclination = Field(0.0, "magnetic inclination", unit=u.rad)
+    prod_site_alt = Field(0.0, "height of observation level", unit=u.m)
+    prod_site_array = Field("None", "site array")
+    prod_site_coord = Field("None", "site (long., lat.) coordinates")
+    prod_site_subarray = Field("None", "site subarray")
+    spectral_index = Field(0.0, "Power-law spectral index of spectrum")
+
 
 
 class CentralTriggerContainer(Container):

ctapipe/io/simteleventsource.py

@@ -156,6 +156,18 @@ def __generator(self):
             data.mcheader.run_array_direction = Angle(
                 self.file_.header['direction'] * u.rad
             )
+            mc_run_head = self.file_.mc_run_headers[-1]
+            data.mcheader.corsika_version = mc_run_head['shower_prog_vers']
+            data.mcheader.simtel_version = mc_run_head['detector_prog_vers']
+            data.mcheader.energy_range_min = mc_run_head['E_range'][0] * u.TeV
+            data.mcheader.energy_range_max = mc_run_head['E_range'][1] * u.TeV
+            data.mcheader.prod_site_B_total = mc_run_head['B_total'] * u.uT
+            data.mcheader.prod_site_B_declination = Angle(
+                mc_run_head['B_declination'] * u.rad)
+            data.mcheader.prod_site_B_inclination = Angle(
+                mc_run_head['B_inclination'] * u.rad)
+            data.mcheader.prod_site_alt = mc_run_head['obsheight'] * u.m
+            data.mcheader.spectral_index = mc_run_head['spectral_index']
 
             # this should be done in a nicer way to not re-allocate the
             # data each time (right now it's just deleted and garbage
@@ -206,3 +218,6 @@ def __generator(self):
                 data.mc.tel[tel_id].azimuth_cor = tracking_position.get('azimuth_cor', 0)
                 data.mc.tel[tel_id].altitude_cor = tracking_position.get('altitude_cor', 0)
             yield data
+
+
+

ctapipe/io/tests/test_simteleventsource.py

@@ -118,3 +118,36 @@ def test_that_event_is_not_modified_after_loop():
         #      assert last_event == event
         # So for the moment we just compare event ids
         assert event.r0.event_id == last_event.r0.event_id
+
+
+def test_additional_meta_data_from_mc_header():
+    kwargs = dict(config=None, tool=None, input_url=gamma_test_path)
+
+    with SimTelEventSource(**kwargs) as reader:
+        data = next(iter(reader))
+
+    # for expectation values
+    from astropy import units as u
+    from astropy.coordinates import Angle
+
+    assert data.mcheader.corsika_version == 6990
+    assert data.mcheader.simtel_version == 1404919891
+    assert data.mcheader.spectral_index == -2.0
+
+    name_expectation = {
+        'energy_range_min': u.Quantity(3.0e-03, u.TeV),
+        'energy_range_max': u.Quantity(3.3e+02, u.TeV),
+        'prod_site_B_total': u.Quantity(27.181243896484375, u.uT),
+        'prod_site_B_declination': Angle(0.0 * u.rad),
+        'prod_site_B_inclination': Angle(-1.1581752300262451 * u.rad),
+        'prod_site_alt': 1640.0 * u.m,
+    }
+
+    for name, expectation in name_expectation.items():
+        value = getattr(data.mcheader, name)
+
+        assert value.unit == expectation.unit
+        assert np.isclose(
+            value.to_value(expectation.unit),
+            expectation.to_value(expectation.unit)
+        )