Merge pull request #1910 from cta-observatory/effective_focal_length

Use effective focal length by default, always guess telescope names from nominal

ctapipe/conftest.py

@@ -44,7 +44,8 @@ def _global_example_event():
 
     print("******************** LOAD TEST EVENT ***********************")
 
-    with SimTelEventSource(input_url=filename) as reader:
+    # FIXME: switch to prod5b+ file that contains effective focal length
+    with SimTelEventSource(input_url=filename, focal_length_choice='nominal') as reader:
         event = next(iter(reader))
 
     return event
@@ -59,7 +60,7 @@ def example_subarray():
 
     print("******************** LOAD TEST EVENT ***********************")
 
-    with SimTelEventSource(input_url=filename) as reader:
+    with SimTelEventSource(input_url=filename, focal_length_choice='nominal') as reader:
         return reader.subarray
 
 
@@ -84,7 +85,7 @@ def _subarray_and_event_gamma_off_axis_500_gev():
 
     path = get_dataset_path("lst_prod3_calibration_and_mcphotons.simtel.zst")
 
-    with SimTelEventSource(path) as source:
+    with SimTelEventSource(path, focal_length_choice='nominal') as source:
         it = iter(source)
         # we want the second event, first event is a corner clipper
         next(it)
@@ -353,6 +354,7 @@ def dl1_muon_file(dl1_tmp_path):
             "--write-images",
             "--DataWriter.write_parameters=False",
             "--DataWriter.Contact.name=αℓℓ the äüöß",
+            "--SimTelEventSource.focal_length_choice=nominal",
         ]
         assert run_tool(ProcessorTool(), argv=argv, cwd=dl1_tmp_path) == 0
         return output

ctapipe/io/eventseeker.py

@@ -28,7 +28,7 @@ class EventSeeker(Component):
     To obtain a particular event in a simtel file:
 
     >>> from ctapipe.io import SimTelEventSource
-    >>> event_source = SimTelEventSource(input_url="dataset://gamma_test_large.simtel.gz")
+    >>> event_source = SimTelEventSource(input_url="dataset://gamma_test_large.simtel.gz", focal_length_choice="nominal")
     >>> seeker = EventSeeker(event_source=event_source)
     >>> event = seeker.get_event_index(2)
     >>> print(event.count)
@@ -37,7 +37,7 @@ class EventSeeker(Component):
     To obtain a particular event in a simtel file from its event_id:
 
     >>> from ctapipe.io import SimTelEventSource
-    >>> event_source = SimTelEventSource(input_url="dataset://gamma_test_large.simtel.gz", back_seekable=True)
+    >>> event_source = SimTelEventSource(input_url="dataset://gamma_test_large.simtel.gz", back_seekable=True, focal_length_choice="nominal")
     >>> seeker = EventSeeker(event_source=event_source)
     >>> event = seeker.get_event_id(31007)
     >>> print(event.count)

ctapipe/io/eventsource.py

@@ -36,7 +36,7 @@ class EventSource(Component):
     appropriate subclass if a compatible source is found for the given
     ``input_url``.
 
-    >>> EventSource(input_url="dataset://gamma_test_large.simtel.gz")
+    >>> EventSource(input_url="dataset://gamma_prod5.simtel.zst")
     <ctapipe.io.simteleventsource.SimTelEventSource ...>
 
     An ``EventSource`` can also be created through the configuration system,
@@ -45,7 +45,7 @@ class EventSource(Component):
     >>> self.source = EventSource(parent=self) # doctest: +SKIP
 
     To loop through the events in a file:
-    >>> source = EventSource(input_url="dataset://gamma_test_large.simtel.gz", max_events=2)
+    >>> source = EventSource(input_url="dataset://gamma_prod5.simtel.zst", max_events=2)
     >>> for event in source:
     ...     print(event.count)
     0
@@ -58,7 +58,7 @@ class EventSource(Component):
     It is encouraged to use ``EventSource`` in a context manager to ensure
     the correct cleanups are performed when you are finished with the source:
 
-    >>> with EventSource(input_url="dataset://gamma_test_large.simtel.gz", max_events=2) as source:
+    >>> with EventSource(input_url="dataset://gamma_prod5.simtel.zst", max_events=2) as source:
     ...    for event in source:
     ...        print(event.count)
     0

ctapipe/io/simteleventsource.py

@@ -180,7 +180,7 @@ class SimTelEventSource(EventSource):
 
     focal_length_choice = CaselessStrEnum(
         ["nominal", "effective"],
-        default_value="nominal",
+        default_value="effective",
         help=(
             "if both nominal and effective focal lengths are available in the "
             "SimTelArray file, which one to use when constructing the "
@@ -306,26 +306,31 @@ def prepare_subarray_info(self, telescope_descriptions, header):
             pixel_settings = telescope_description["pixel_settings"]
 
             n_pixels = cam_settings["n_pixels"]
-            focal_length = u.Quantity(cam_settings["focal_length"], u.m)
             mirror_area = u.Quantity(cam_settings["mirror_area"], u.m**2)
 
-            if self.focal_length_choice == "effective":
-                try:
-                    focal_length = u.Quantity(
-                        cam_settings["effective_focal_length"], u.m
-                    )
-                except KeyError as err:
-                    raise RuntimeError(
-                        f"the SimTelEventSource option 'focal_length_choice' was set to "
-                        f"{self.focal_length_choice}, but the effective focal length "
-                        f"was not present in the file. ({err})"
-                    )
+            nominal_focal_length = u.Quantity(cam_settings["focal_length"], u.m)
+            effective_focal_length = u.Quantity(
+                cam_settings.get("effective_focal_length", np.nan), u.m
+            )
 
             try:
-                telescope = guess_telescope(n_pixels, focal_length)
+                telescope = guess_telescope(n_pixels, nominal_focal_length)
             except ValueError:
                 telescope = unknown_telescope(mirror_area, n_pixels)
 
+            if self.focal_length_choice == "effective":
+                if np.isnan(effective_focal_length):
+                    raise RuntimeError(
+                        f"`SimTelEventSource.focal_length_choice` was set to"
+                        f" {self.focal_length_choice!r}, but the effective focal length"
+                        f" was not present in the file. "
+                        " Use nominal focal length or adapt your simulation configuration"
+                        " to include the effective focal length"
+                    )
+                focal_length = effective_focal_length
+            else:
+                focal_length = nominal_focal_length
+
             optics = OpticsDescription(
                 name=telescope.name,
                 num_mirrors=telescope.n_mirrors,

ctapipe/io/tests/conftest.py

@@ -1,6 +1,5 @@
 import pytest
 from ctapipe.io import EventSource, DataWriter
-from ctapipe.utils import get_dataset_path
 
 
 @pytest.fixture(scope="session")
@@ -9,11 +8,10 @@ def r1_path(tmp_path_factory):
 
 
 @pytest.fixture(scope="session")
-def r1_hdf5_file(r1_path):
+def r1_hdf5_file(prod5_proton_simtel_path, r1_path):
     source = EventSource(
-        get_dataset_path("gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"),
+        prod5_proton_simtel_path,
         max_events=5,
-        allowed_tels=[1, 2, 3, 4],
     )
 
     path = r1_path / "test_r1.h5"

ctapipe/io/tests/test_datawriter.py

@@ -51,6 +51,7 @@ def test_write(tmpdir: Path):
         get_dataset_path("gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"),
         max_events=20,
         allowed_tels=[1, 2, 3, 4],
+        focal_length_choice='nominal',
     )
     calibrate = CameraCalibrator(subarray=source.subarray)
 
@@ -131,6 +132,7 @@ def test_roundtrip(tmpdir: Path):
         get_dataset_path("gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"),
         max_events=20,
         allowed_tels=[1, 2, 3, 4],
+        focal_length_choice='nominal',
     )
     calibrate = CameraCalibrator(subarray=source.subarray)
 
@@ -203,7 +205,7 @@ def test_dl1writer_no_events(tmpdir: Path):
 
     output_path = Path(tmpdir / "no_events.dl1.h5")
     dataset = "lst_prod3_calibration_and_mcphotons.simtel.zst"
-    with EventSource(get_dataset_path(dataset)) as source:
+    with EventSource(get_dataset_path(dataset), focal_length_choice='nominal') as source:
         # exhaust source
         for _ in source:
             pass
@@ -233,7 +235,7 @@ def test_dl1writer_no_events(tmpdir: Path):
 def test_metadata(tmpdir: Path):
     output_path = Path(tmpdir / "metadata.dl1.h5")
 
-    dataset = "lst_prod3_calibration_and_mcphotons.simtel.zst"
+    dataset = "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
 
     config = Config(
         {

ctapipe/io/tests/test_event_source.py

@@ -6,6 +6,8 @@
 
 from ctapipe.core import Component
 
+prod5_path = "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
+
 
 def test_construct():
     # at least one of input_url / parent / config is required
@@ -43,20 +45,20 @@ def datalevels(self):
 
 
 def test_can_be_implemented():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     test_reader = DummyReader(input_url=dataset)
     assert test_reader is not None
 
 
 def test_is_iterable():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     test_reader = DummyReader(input_url=dataset)
     for _ in test_reader:
         pass
 
 
 def test_function():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     reader = EventSource(input_url=dataset)
     assert isinstance(reader, SimTelEventSource)
     assert reader.input_url == dataset
@@ -75,15 +77,15 @@ def test_function_nonexistant_file():
 
 
 def test_from_config():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     config = Config({"EventSource": {"input_url": dataset}})
     reader = EventSource(config=config)
     assert isinstance(reader, SimTelEventSource)
     assert reader.input_url == dataset
 
 
 def test_from_config_parent():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
 
     class Parent(Component):
         def __init__(self, config=None, parent=None):
@@ -109,7 +111,7 @@ def __init__(self, config=None, parent=None):
 
 def test_from_config_default():
     old_default = EventSource.input_url.default_value
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     EventSource.input_url.default_value = dataset
     config = Config()
     reader = EventSource(config=config, parent=None)
@@ -119,18 +121,18 @@ def test_from_config_default():
 
 
 def test_from_config_invalid_type():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     EventSource.input_url.default_value = dataset
     config = Config({"EventSource": {"input_url": 124}})
     with pytest.raises(TraitError):
         EventSource(config=config, parent=None)
 
 
 def test_event_source_input_url_config_override():
-    dataset1 = get_dataset_path("gamma_test_large.simtel.gz")
-    dataset2 = get_dataset_path(
+    dataset1 = get_dataset_path(
         "gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"
     )
+    dataset2 = get_dataset_path(prod5_path)
 
     config = Config({"EventSource": {"input_url": dataset1}})
     reader = EventSource(input_url=dataset2, config=config)
@@ -141,29 +143,29 @@ def test_event_source_input_url_config_override():
 
 def test_max_events():
     max_events = 10
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     reader = EventSource(input_url=dataset, max_events=max_events)
     assert reader.max_events == max_events
 
 
 def test_max_events_from_config():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     max_events = 10
     config = Config({"EventSource": {"input_url": dataset, "max_events": max_events}})
     reader = EventSource(config=config)
     assert reader.max_events == max_events
 
 
 def test_allowed_tels():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     reader = EventSource(input_url=dataset)
     assert reader.allowed_tels is None
     reader = EventSource(input_url=dataset, allowed_tels={1, 3})
-    assert len(reader.allowed_tels) == 2
+    assert reader.allowed_tels == {1, 3}
 
 
 def test_allowed_tels_from_config():
-    dataset = get_dataset_path("gamma_test_large.simtel.gz")
+    dataset = get_dataset_path(prod5_path)
     config = Config({"EventSource": {"input_url": dataset, "allowed_tels": {1, 3}}})
     reader = EventSource(config=config, parent=None)
-    assert len(reader.allowed_tels) == 2
+    assert reader.allowed_tels == {1, 3}

ctapipe/io/tests/test_eventseeker.py

@@ -8,7 +8,11 @@
 
 def test_eventseeker():
 
-    with SimTelEventSource(input_url=dataset, back_seekable=True) as reader:
+    with SimTelEventSource(
+        input_url=dataset,
+        back_seekable=True,
+        focal_length_choice='nominal',
+    ) as reader:
 
         seeker = EventSeeker(event_source=reader)
 
@@ -33,7 +37,8 @@ def test_eventseeker():
             seeker.get_event_index(dict())
 
     with SimTelEventSource(
-        input_url=dataset, max_events=5, back_seekable=True
+        input_url=dataset, max_events=5, back_seekable=True,
+        focal_length_choice='nominal',
     ) as reader:
         seeker = EventSeeker(event_source=reader)
         with pytest.raises(IndexError):
@@ -42,7 +47,10 @@ def test_eventseeker():
 
 
 def test_eventseeker_edit():
-    with SimTelEventSource(input_url=dataset, back_seekable=True) as reader:
+    with SimTelEventSource(
+        input_url=dataset, back_seekable=True,
+        focal_length_choice='nominal',
+    ) as reader:
         seeker = EventSeeker(event_source=reader)
         event = seeker.get_event_index(1)
         assert event.count == 1
@@ -54,7 +62,10 @@ def test_eventseeker_edit():
 
 def test_eventseeker_simtel():
     # Ensure the EventSeeker can forward seek even if back-seeking is not possible
-    with SimTelEventSource(input_url=dataset, back_seekable=False) as reader:
+    with SimTelEventSource(
+        input_url=dataset, back_seekable=False,
+        focal_length_choice='nominal',
+    ) as reader:
         seeker = EventSeeker(event_source=reader)
         event = seeker.get_event_index(1)
         assert event.count == 1

ctapipe/io/tests/test_hdf5eventsource.py

@@ -136,4 +136,4 @@ def test_read_r1(r1_hdf5_file):
             pass
 
         assert e is not None
-        assert e.count == 4
+        assert e.count == 3

ctapipe/io/tests/test_prod2.py

@@ -9,7 +9,10 @@
 
 def test_eventio_prod2():
     with pytest.warns(UnknownPixelShapeWarning):
-        with SimTelEventSource(input_url=dataset) as reader:
+        with SimTelEventSource(
+            input_url=dataset,
+            focal_length_choice='nominal',
+        ) as reader:
             for event in reader:
                 if event.count == 2:
                     break