Repository: cta-observatory/ctapipe
Pull request #2077 (2022-09-21 11:18:31)

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

codecov[bot]: # [Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
Base: **92.43**% // Head: **92.43**% // Increases project coverage by **`+0.00%`** :tada:
> Coverage data is based on head [(`d1c418e`)](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) compared to base [(`b0d9b33`)](https://codecov.io/gh/cta-observatory/ctapipe/commit/b0d9b3385d5704028822e703b9b10c5525bf83b5?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).
> Patch coverage: 92.85% of modified lines in pull request are covered.

> :exclamation: Current head d1c418e differs from pull request most recent head 0326e53. Consider uploading reports for the commit 0326e53 to get more accurate results

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##           master    #2077   +/-   ##
=======================================
  Coverage   92.43%   92.43%           
=======================================
  Files         197      197           
  Lines       16198    16218   +20     
=======================================
+ Hits        14972    14991   +19     
- Misses       1226     1227    +1     
```


| [Impacted Files](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | Coverage Δ | |
|---|---|---|
| [ctapipe/io/simteleventsource.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9pby9zaW10ZWxldmVudHNvdXJjZS5weQ==) | `96.51% <81.81%> (-0.23%)` | :arrow_down: |
| [ctapipe/coordinates/\_\_init\_\_.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9jb29yZGluYXRlcy9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [ctapipe/io/tests/test\_simteleventsource.py](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-Y3RhcGlwZS9pby90ZXN0cy90ZXN0X3NpbXRlbGV2ZW50c291cmNlLnB5) | `100.00% <100.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/cta-observatory/ctapipe/pull/2077?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).


Patch:
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 2b4b141e..4dd57cc2 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -73,6 +73,27 @@ SIMTEL_TO_CTA_EVENT_TYPE = {
     -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 @@ class SimTelEventSource(EventSource):
 
         # 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 @@ class SimTelEventSource(EventSource):
             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 @@ class SimTelEventSource(EventSource):
         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),


Test patch:
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index c9844ef0..80a23be0 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -5,7 +5,7 @@ from pathlib import Path
 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)