Skip to content

Commit

Permalink
Merge pull request #2197 from cta-observatory/histogram_units
Browse files Browse the repository at this point in the history
Fix unit handling in muon ring_completeness
  • Loading branch information
maxnoe authored Jan 9, 2023
2 parents 0b58bda + bc57ba4 commit 9a83547
Show file tree
Hide file tree
Showing 2 changed files with 10 additions and 4 deletions.
11 changes: 7 additions & 4 deletions ctapipe/image/muon/features.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import astropy.units as u
import numpy as np
from ...utils.quantities import all_to_value

from ...utils.quantities import all_to_value

__all__ = [
"mean_squared_error",
Expand Down Expand Up @@ -105,8 +106,10 @@ def ring_completeness(
"""

angle = np.arctan2(pixel_y - center_y, pixel_x - center_x)
if hasattr(angle, "unit"):
angle = angle.to_value(u.rad)

hist, edges = np.histogram(angle, bins=bins, range=[-np.pi, np.pi], weights=weights)
hist, _ = np.histogram(angle, bins=bins, range=[-np.pi, np.pi], weights=weights)

bins_above_threshold = hist > threshold

Expand Down Expand Up @@ -144,7 +147,7 @@ def ring_containment(radius, center_x, center_y, camera_radius):
radius, center_x, center_y, camera_radius = all_to_value(
radius, center_x, center_y, camera_radius, unit=radius.unit
)
d = np.sqrt(center_x ** 2 + center_y ** 2)
d = np.sqrt(center_x**2 + center_y**2)

# one circle fully contained in the other
if d <= np.abs(camera_radius - radius):
Expand All @@ -154,5 +157,5 @@ def ring_containment(radius, center_x, center_y, camera_radius):
if d > (radius + camera_radius):
return 0.0

a = (radius ** 2 - camera_radius ** 2 + d ** 2) / (2 * d)
a = (radius**2 - camera_radius**2 + d**2) / (2 * d)
return np.arccos(a / radius) / np.pi
3 changes: 3 additions & 0 deletions docs/changes/2197.bugfix.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
Fix mixture of quantity and unit-less values passed to ``np.histogram``
in ``ctapipe.image.muon.ring_completeness``, which raises an error with
astropy 5.2.1.

0 comments on commit 9a83547

Please sign in to comment.