Improve performance of concentration by 2x by stripping units interna…

…lly (#1350)

* Improve performance of concentration by 2x by stripping units internally

* Also strip units in timing_parameters

ctapipe/image/concentration.py

@@ -1,6 +1,9 @@
 import numpy as np
+import astropy.units as u
+
 from ..containers import ConcentrationContainer
 from .hillas import camera_to_shower_coordinates
+from ..utils.quantities import all_to_value
 
 
 def concentration(geom, image, hillas_parameters):
@@ -14,9 +17,14 @@ def concentration(geom, image, hillas_parameters):
  """
 
  h = hillas_parameters
+ unit = h.x.unit
+
+ pix_x, pix_y, x, y, length, width = all_to_value(
+ geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, unit=unit
+ )
 
- delta_x = geom.pix_x - h.x
- delta_y = geom.pix_y - h.y
+ delta_x = pix_x - x
+ delta_y = pix_y - y
 
  # sort pixels by distance to cog
  cog_pixels = np.argsort(delta_x ** 2 + delta_y ** 2)
@@ -25,15 +33,15 @@ def concentration(geom, image, hillas_parameters):
  if hillas_parameters.width != 0:
  # get all pixels inside the hillas ellipse
  longi, trans = camera_to_shower_coordinates(
- geom.pix_x, geom.pix_y, h.x, h.y, h.psi
+ pix_x, pix_y, x, y, h.psi.to_value(u.rad)
  )
- mask_core = (longi ** 2 / h.length ** 2) + (trans ** 2 / h.width ** 2) <= 1.0
+ mask_core = (longi ** 2 / length ** 2) + (trans ** 2 / width ** 2) <= 1.0
  conc_core = image[mask_core].sum() / h.intensity
  else:
  conc_core = 0.0
 
  concentration_pixel = image.max() / h.intensity
 
  return ConcentrationContainer(
- cog=conc_cog, core=conc_core, pixel=concentration_pixel,
+ cog=conc_cog, core=conc_core, pixel=concentration_pixel
  )

ctapipe/image/timing.py

@@ -3,9 +3,11 @@
 """
 
 import numpy as np
+import astropy.units as u
 from numpy.polynomial.polynomial import polyval
-from ctapipe.containers import TimingParametersContainer
+from ..containers import TimingParametersContainer
 from .hillas import camera_to_shower_coordinates
+from ..utils.quantities import all_to_value
 
 
 __all__ = ["timing_parameters"]
@@ -44,17 +46,19 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N
  if (image < 0).any():
  raise ValueError("The non-masked pixels must verify signal >= 0")
 
- pix_x = geom.pix_x
- pix_y = geom.pix_y
+ h = hillas_parameters
+ pix_x, pix_y, x, y, length, width = all_to_value(
+ geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, unit=unit
+ )
 
- longi, trans = camera_to_shower_coordinates(
- pix_x, pix_y, hillas_parameters.x, hillas_parameters.y, hillas_parameters.psi
+ longi, _ = camera_to_shower_coordinates(
+ pix_x, pix_y, x, y, hillas_parameters.psi.to_value(u.rad)
  )
  (slope, intercept), cov = np.polyfit(
- longi.value, peak_time, deg=1, w=np.sqrt(image), cov=True,
+ longi, peak_time, deg=1, w=np.sqrt(image), cov=True
  )
  slope_err, intercept_err = np.sqrt(np.diag(cov))
- predicted_time = polyval(longi.value, (intercept, slope))
+ predicted_time = polyval(longi, (intercept, slope))
  deviation = np.sqrt(np.sum((peak_time - predicted_time) ** 2) / peak_time.size)
 
  return TimingParametersContainer(