Merge pull request #2077 from cta-observatory/fix_roundng_error

Convert float32 pi/2 value to float64 pi/2 value in SimTelEventSource

ctapipe/io/simteleventsource.py

@@ -73,6 +73,27 @@
     -1: EventType.OTHER_CALIBRATION,
 }
 
+_half_pi = 0.5 * np.pi
+_half_pi_maxval = (1 + 1e-6) * _half_pi
+
+
+def _clip_altitude_if_close(altitude):
+    """
+    Round absolute values slightly larger than pi/2 in float64 to pi/2
+
+    These can come from simtel_array because float32(pi/2) > float64(pi/2)
+    and simtel using float32.
+
+    Astropy complains about these values, so we fix them here.
+    """
+    if altitude > _half_pi and altitude < _half_pi_maxval:
+        return _half_pi
+
+    if altitude < -_half_pi and altitude > -_half_pi_maxval:
+        return -_half_pi
+
+    return altitude
+
 
 @enum.unique
 class MirrorClass(enum.Enum):
@@ -700,9 +721,13 @@ def _fill_event_pointing(tracking_position):
 
         # take pointing corrected position if available
         if np.isnan(altitude_cor):
-            altitude = u.Quantity(altitude_raw, u.rad, copy=False)
+            altitude = u.Quantity(
+                _clip_altitude_if_close(altitude_raw), u.rad, copy=False
+            )
         else:
-            altitude = u.Quantity(altitude_cor, u.rad, copy=False)
+            altitude = u.Quantity(
+                _clip_altitude_if_close(altitude_cor), u.rad, copy=False
+            )
 
         return TelescopePointingContainer(azimuth=azimuth, altitude=altitude)
 
@@ -814,8 +839,8 @@ def _parse_simulation_header(self):
             detector_prog_id=mc_run_head["detector_prog_id"],
             n_showers=mc_run_head["n_showers"],
             shower_reuse=mc_run_head["n_use"],
-            max_alt=mc_run_head["alt_range"][1] * u.rad,
-            min_alt=mc_run_head["alt_range"][0] * u.rad,
+            max_alt=_clip_altitude_if_close(mc_run_head["alt_range"][1]) * u.rad,
+            min_alt=_clip_altitude_if_close(mc_run_head["alt_range"][0]) * u.rad,
             max_az=mc_run_head["az_range"][1] * u.rad,
             min_az=mc_run_head["az_range"][0] * u.rad,
             diffuse=mc_run_head["diffuse"],
@@ -878,9 +903,11 @@ def _fill_simulated_event_information(array_event):
         if mc_shower is None:
             return
 
+        alt = _clip_altitude_if_close(mc_shower["altitude"])
+
         return SimulatedShowerContainer(
             energy=u.Quantity(mc_shower["energy"], u.TeV),
-            alt=Angle(mc_shower["altitude"], u.rad),
+            alt=Angle(alt, u.rad),
             az=Angle(mc_shower["azimuth"], u.rad),
             core_x=u.Quantity(mc_event["xcore"], u.m),
             core_y=u.Quantity(mc_event["ycore"], u.m),

ctapipe/io/tests/test_simteleventsource.py

@@ -5,7 +5,7 @@
 import astropy.units as u
 import numpy as np
 import pytest
-from astropy.coordinates import Angle
+from astropy.coordinates import Angle, Latitude
 from astropy.time import Time
 from traitlets.config import Config
 
@@ -515,3 +515,33 @@ def test_simtel_no_metadata(monkeypatch):
 
         assert all([t.camera.name.startswith("UNKNOWN") for t in subarray.tel.values()])
         assert all([t.optics.name.startswith("UNKNOWN") for t in subarray.tel.values()])
+
+
+@pytest.mark.parametrize("sign", (-1, 1))
+def test_float32_pihalf(sign):
+    float32_pihalf = np.float32(sign * np.pi / 2)
+    tracking_position = {"azimuth_raw": 0, "altitude_raw": float32_pihalf}
+    pointing = SimTelEventSource._fill_event_pointing(tracking_position)
+    # check that we changed the value to float64 pi/2 to avoid astropy error
+    assert pointing.altitude.value == sign * np.pi / 2
+    # check we can create a Latitude:
+    Latitude(pointing.altitude.value, u.rad)
+
+    event = {
+        "mc_shower": {
+            "energy": 1.0,
+            "altitude": float32_pihalf,
+            "azimuth": 0.0,
+            "h_first_int": 20e3,
+            "xmax": 350,
+            "primary_id": 1,
+        },
+        "mc_event": {
+            "xcore": 0.0,
+            "ycore": 50.0,
+        },
+    }
+    shower = SimTelEventSource._fill_simulated_event_information(event)
+    assert shower.alt.value == sign * np.pi / 2
+    # check we cana create a Latitude:
+    Latitude(shower.alt.value, u.rad)