From 5ef63ef2a7a093455f4b0683a68653cc42e84df7 Mon Sep 17 00:00:00 2001 From: Dominik Neise Date: Thu, 13 Dec 2018 09:53:08 +0100 Subject: [PATCH] Additional metadata from pyhessio, similar to #655 (#895) --- ctapipe/io/containers.py | 15 ++++++++++ ctapipe/io/simteleventsource.py | 15 ++++++++++ ctapipe/io/tests/test_simteleventsource.py | 33 ++++++++++++++++++++++ 3 files changed, 63 insertions(+) diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py index e2c5d539a9d..d7b02ff88bd 100644 --- a/ctapipe/io/containers.py +++ b/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): diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py index 0ddfb888651..13816657f63 100644 --- a/ctapipe/io/simteleventsource.py +++ b/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 + + + diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py index 34ba28688e7..5bb9e7a6f35 100644 --- a/ctapipe/io/tests/test_simteleventsource.py +++ b/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) + )