Repository: cta-observatory/ctapipe
Pull request #2018 (2022-07-15 10:12:12)

Get reference_location and optics/camera name from simtel metadata if available
This completely removes the guessing if the metadata first included in the current prod6 iterations is present.

* We take the reference location from the `ALTITUDE`, `*LATITUDE`, `*LONGITUDE` attributes.
* Optics and camera names from `OPTICS_CONFIG_NAME` and `CAMERA_CONFIG_NAME`
* subarray name from `ARRAY_CONFIG_NAME`
* number of mirrors from `MIRROR_CLASS`


Patch:
diff --git a/ctapipe/instrument/guess.py b/ctapipe/instrument/guess.py
index 838a72e4..856be851 100644
--- a/ctapipe/instrument/guess.py
+++ b/ctapipe/instrument/guess.py
@@ -9,6 +9,13 @@ from collections import namedtuple
 import astropy.units as u
 import numpy as np
 
+__all__ = [
+    "GuessingKey",
+    "GuessingResult",
+    "guess_telescope",
+    "unknown_telescope",
+    "type_from_mirror_area",
+]
 
 GuessingKey = namedtuple(
     "GuessingKey",
@@ -99,21 +106,27 @@ def guess_telescope(n_pixels, focal_length, num_mirror_tiles=None):
         raise ValueError(f"Unknown telescope: n_pixel={n_pixels}, f={focal_length}")
 
 
-def unknown_telescope(mirror_area, n_pixels, n_mirrors=-1):
-    """Create a GuessingResult for an unknown_telescope"""
-    mirror_area = u.Quantity(mirror_area, u.m**2).to_value(u.m**2)
+@u.quantity_input(mirror_area=u.m**2)
+def type_from_mirror_area(mirror_area):
+    mirror_diameter = (2 * np.sqrt(mirror_area / np.pi)).to_value(u.m)
 
-    mirror_diameter = 2 * np.sqrt(mirror_area / np.pi)
     if mirror_diameter < 8:
-        telescope_type = "SST"
-    elif mirror_diameter < 16:
-        telescope_type = "MST"
-    else:
-        telescope_type = "LST"
+        return "SST"
+
+    if mirror_diameter < 16:
+        return "MST"
 
+    return "LST"
+
+
+def unknown_telescope(mirror_area, n_pixels, n_mirrors=-1):
+    """Create a GuessingResult for an unknown_telescope"""
+    # this allows passing a plain number in square meter and any quantity
+    # with an area unit
+    mirror_area = u.Quantity(mirror_area, u.m**2)
     return GuessingResult(
-        type=telescope_type,
-        name=f"UNKNOWN-{mirror_area:.0f}M2",
+        type=type_from_mirror_area(mirror_area),
+        name=f"UNKNOWN-{mirror_area.to_value(u.m**2):.0f}M2",
         camera_name=f"UNKNOWN-{n_pixels}PX",
         n_mirrors=n_mirrors,
     )
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 196bda33..6a95d8bf 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -1,3 +1,4 @@
+import enum
 import warnings
 from gzip import GzipFile
 from io import BufferedReader
@@ -5,7 +6,7 @@ from pathlib import Path
 
 import numpy as np
 from astropy import units as u
-from astropy.coordinates import Angle
+from astropy.coordinates import Angle, EarthLocation
 from astropy.time import Time
 from eventio.file_types import is_eventio
 from eventio.simtel.simtelfile import SimTelFile
@@ -40,14 +41,16 @@ from ..instrument import (
     TelescopeDescription,
 )
 from ..instrument.camera import UnknownPixelShapeWarning
-from ..instrument.guess import guess_telescope, unknown_telescope
+from ..instrument.guess import (
+    GuessingResult,
+    guess_telescope,
+    type_from_mirror_area,
+    unknown_telescope,
+)
 from ..reco.impact_distance import shower_impact_distance
 from .datalevels import DataLevel
 from .eventsource import EventSource
 
-X_MAX_UNIT = u.g / (u.cm**2)
-
-
 __all__ = ["SimTelEventSource"]
 
 # Mapping of SimTelArray Calibration trigger types to EventType:
@@ -61,6 +64,16 @@ SIMTEL_TO_CTA_EVENT_TYPE = {
 }
 
 
+@enum.unique
+class MirrorClass(enum.Enum):
+    """Enum for the sim_telarray MIRROR_CLASS values"""
+
+    SINGLE_SEGMENTED_MIRROR = 0
+    SINGLE_SOLID_PARABOLIC_MIRROR = 1
+    DUAL_MIRROR = 2
+
+
+X_MAX_UNIT = u.g / (u.cm**2)
 NANOSECONDS_PER_DAY = (1 * u.day).to_value(u.ns)
 
 
@@ -71,6 +84,22 @@ def parse_simtel_time(simtel_time):
     )
 
 
+def _location_from_meta(global_meta):
+    """Extract reference location of subarray from metadata"""
+    lat = global_meta.get(b"*LATITUDE")
+    lon = global_meta.get(b"*LONGITUDE")
+    height = global_meta.get(b"ALTITUDE")
+
+    if lat is None or lon is None or height is None:
+        return None
+
+    return EarthLocation(
+        lon=float(lon) * u.deg,
+        lat=float(lat) * u.deg,
+        height=float(height) * u.m,
+    )
+
+
 def build_camera(cam_settings, pixel_settings, telescope, frame):
     pixel_shape = cam_settings["pixel_shape"][0]
     try:
@@ -109,6 +138,27 @@ def build_camera(cam_settings, pixel_settings, telescope, frame):
     )
 
 
+def _telescope_from_meta(telescope_meta, mirror_area):
+    optics_name = telescope_meta.get(b"OPTICS_CONFIG_NAME")
+    camera_name = telescope_meta.get(b"CAMERA_CONFIG_NAME")
+    mirror_class = telescope_meta.get(b"MIRROR_CLASS")
+
+    if optics_name is None or camera_name is None or mirror_class is None:
+        return None
+
+    telescope_type = type_from_mirror_area(mirror_area)
+
+    mirror_class = MirrorClass(int(mirror_class))
+    n_mirrors = 2 if mirror_class is MirrorClass.DUAL_MIRROR else 1
+
+    return GuessingResult(
+        type=telescope_type,
+        name=optics_name.decode("utf-8"),
+        camera_name=camera_name.decode("utf-8"),
+        n_mirrors=n_mirrors,
+    )
+
+
 def apply_simtel_r1_calibration(
     r0_waveforms, pedestal, dc_to_pe, gain_selector, calib_scale=1.0, calib_shift=0.0
 ):
@@ -167,7 +217,32 @@ def apply_simtel_r1_calibration(
 
 
 class SimTelEventSource(EventSource):
-    """Read events from a SimTelArray data file (in EventIO format)."""
+    """
+    Read events from a SimTelArray data file (in EventIO format).
+
+    ctapipe makes use of the sim_telarray metadata system to fill some
+    information not directly accessible from the data inside the files itself.
+    Make sure you set this parameters in the simulation configuration to fully
+    make use of ctapipe. In future, ctapipe might require these metadata parameters.
+
+    This includes:
+
+    * Reference Point of the telescope coordinates. Make sure to include the
+        user-defined parameters ``LONGITUDE`` and ``LATITUDE`` with the geodetic
+        coordinates of the array reference point. Also make sure the ``ALTITUDE``
+        config parameter is included in the global metadata.
+
+    * Names of the optical structures and the cameras are read from
+        ``OPTICS_CONFIG_NAME`` and ``CAMERA_CONFIG_NAME``, make sure to include
+        these in the telescope meta.
+
+    * The ``MIRROR_CLASS`` should also be included in the telescope meta
+        to correctly setup coordinate transforms.
+
+    If these parameters are not included in the input data, ctapipe will
+    fallback guesses these based on avaible data and the list of known telescopes
+    for `ctapipe.instrument.guess_telescope`.
+    """
 
     skip_calibration_events = Bool(True, help="Skip calibration events").tag(
         config=True
@@ -317,14 +392,29 @@ class SimTelEventSource(EventSource):
                 cam_settings.get("effective_focal_length", np.nan), u.m
             )
 
-            try:
-                telescope = guess_telescope(
-                    n_pixels,
-                    equivalent_focal_length,
-                    cam_settings["n_mirrors"],
+            telescope = _telescope_from_meta(
+                self.file_.telescope_meta.get(tel_id, {}),
+                mirror_area,
+            )
+
+            if telescope is None:
+                try:
+                    telescope = guess_telescope(
+                        n_pixels,
+                        equivalent_focal_length,
+                        cam_settings["n_mirrors"],
+                    )
+                except ValueError:
+                    telescope = unknown_telescope(mirror_area, n_pixels)
+
+                # TODO: switch to warning or even an exception once
+                # we start relying on this.
+                self.log.info(
+                    "Could not determine telescope from sim_telarray metadata,"
+                    " guessing using builtin lookup-table: %d: %s",
+                    tel_id,
+                    telescope,
                 )
-            except ValueError:
-                telescope = unknown_telescope(mirror_area, n_pixels)
 
             optics = OpticsDescription(
                 name=telescope.name,
@@ -369,10 +459,12 @@ class SimTelEventSource(EventSource):
             self.telescope_indices_original[tel_id] = tel_idx
             tel_positions[tel_id] = header["tel_pos"][tel_idx] * u.m
 
+        name = self.file_.global_meta.get(b"ARRAY_CONFIG_NAME", b"MonteCarloArray")
         subarray = SubarrayDescription(
-            name="MonteCarloArray",
+            name=name.decode(),
             tel_positions=tel_positions,
             tel_descriptions=tel_descriptions,
+            reference_location=_location_from_meta(self.file_.global_meta),
         )
 
         self.n_telescopes_original = len(subarray)


Test patch:
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index 6694d926..1da8b056 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -450,3 +450,61 @@ def test_extracted_calibevents():
             assert e.simulation is not None
             assert e.simulation.shower is None
         assert i == 4
+
+
+def test_simtel_metadata(monkeypatch):
+    from ctapipe.instrument import guess
+
+    # prod6 is the first prod to use the metadata system
+    path = "dataset://gamma_prod6_preliminary.simtel.zst"
+
+    with monkeypatch.context() as m:
+        # remove all guessing keys so we cannot use guessing
+        m.setattr(guess, "LOOKUP_TREE", {})
+
+        with SimTelEventSource(path) as source:
+            subarray = source.subarray
+
+    assert subarray.name == "Paranal-prod6"
+    assert subarray.tel[1].camera.camera_name == "LSTcam"
+    assert subarray.tel[1].optics.name == "LST"
+
+    assert subarray.tel[5].camera.camera_name == "FlashCam"
+    assert subarray.tel[5].optics.name == "MST"
+
+    assert subarray.tel[50].camera.camera_name == "SST-Camera"
+    assert subarray.tel[50].optics.name == "SST"
+
+
+def test_simtel_no_metadata(monkeypatch):
+    from ctapipe.instrument import guess
+
+    # prod5 was before the metadata system was introduced
+    path = "dataset://gamma_prod5.simtel.zst"
+
+    # this will use the guessing system
+    with SimTelEventSource(path) as source:
+        subarray = source.subarray
+
+    assert subarray.name == "MonteCarloArray"
+    assert subarray.tel[1].camera.camera_name == "LSTCam"
+    assert subarray.tel[1].optics.name == "LST"
+
+    assert subarray.tel[5].camera.camera_name == "FlashCam"
+    assert subarray.tel[5].optics.name == "MST"
+
+    assert subarray.tel[50].camera.camera_name == "CHEC"
+    assert subarray.tel[50].optics.name == "ASTRI"
+
+    # check we get all unknown telescopes if we remove the guessing keys
+    with monkeypatch.context() as m:
+        # remove all guessing keys so we cannot use guessing
+        m.setattr(guess, "LOOKUP_TREE", {})
+
+        with SimTelEventSource(path) as source:
+            subarray = source.subarray
+
+        assert all(
+            [t.camera.camera_name.startswith("UNKNOWN") for t in subarray.tel.values()]
+        )
+        assert all([t.optics.name.startswith("UNKNOWN") for t in subarray.tel.values()])