Update Timing Parameters

* use container for output
* use geometry, peakpos and hillas_parameters as input
* adapt tests

ctapipe/image/tests/test_timing_parameters.py

@@ -2,38 +2,50 @@
 import numpy as np
 import astropy.units as u
 from numpy.testing import assert_allclose
+from ctapipe.instrument.camera import CameraGeometry
+from ctapipe.io.containers import HillasParametersContainer
 
 
-def test_grad_fit():
+def test_psi_0():
     """
     Simple test that gradient fitting gives expected answers for perfect
     gradient
     """
-    grad = 2
-    intercept = 1
+    grad = 2.0
+    intercept = 1.0
+
+    geom = CameraGeometry.from_name("LSTCam")
+    hillas = HillasParametersContainer(x=0*u.m, y=0*u.m, psi=0*u.deg)
 
     timing = timing_parameters(
-        pix_x=np.arange(4) * u.deg,
-        pix_y=np.zeros(4) * u.deg,
-        image=np.ones(4),
-        peak_time=intercept * u.ns + grad * np.arange(4) * u.ns,
-        rotation_angle=0 * u.deg
+        geom, 
+        image=np.ones(geom.n_pixels),
+        peakpos=intercept + grad * geom.pix_x.value,
+        hillas_parameters = hillas,
     )
 
     # Test we get the values we put in back out again
-    assert_allclose(timing.gradient, grad * u.ns / u.deg)
-    assert_allclose(timing.intercept, intercept * u.deg)
+    assert_allclose(timing.slope, grad / geom.pix_x.unit)
+    assert_allclose(timing.intercept, intercept)
+
+
+def test_psi_20():
 
     # Then try a different rotation angle
-    rot_angle = 20 * u.deg
-    timing_rot20 = timing_parameters(
-        pix_x=np.arange(4) * u.deg,
-        pix_y=np.zeros(4) * u.deg,
-        image=np.ones(4),
-        peak_time=intercept * u.ns +
-        grad * np.arange(4) * u.ns,
-        rotation_angle=rot_angle
+    grad = 2
+    intercept = 1
+
+    geom = CameraGeometry.from_name("LSTCam")
+    psi = 20*u.deg
+    hillas = HillasParametersContainer(x=0*u.m, y=0*u.m, psi=psi)
+
+    timing = timing_parameters(
+        geom, 
+        image=np.ones(geom.n_pixels),
+        peakpos=intercept + grad * (np.cos(psi) * geom.pix_x.value + np.sin(psi) * geom.pix_y.value),
+        hillas_parameters = hillas,
     )
-    # Test the output again makes sense
-    assert_allclose(timing_rot20.gradient, timing.gradient / np.cos(rot_angle))
-    assert_allclose(timing_rot20.intercept, timing.intercept)
+
+    # Test we get the values we put in back out again
+    assert_allclose(timing.slope, grad / geom.pix_x.unit)
+    assert_allclose(timing.intercept, intercept)

ctapipe/image/timing_parameters.py

@@ -1,80 +1,40 @@
-# Licensed under a 3-clause BSD style license - see LICENSE.rst
-# -*- coding: UTF-8 -*-
 """
 Image timing-based shower image parametrization.
 """
 
 from collections import namedtuple
 import numpy as np
 from astropy.units import Quantity
+from ctapipe.io.containers import TimingParametersContainer
 
 __all__ = [
-    'TimingParameters',
     'timing_parameters'
 ]
 
-TimingParameters = namedtuple(
-    "TimingParameters",
-    "gradient, intercept"
-)
 
-
-class TimingParameterizationError(RuntimeError):
-    pass
-
-
-def rotate_translate(pixel_pos_x, pixel_pos_y, psi):
-    """
-    Function to perform rotation and translation of pixel lists
-
-    Parameters
-    ----------
-    pixel_pos_x: ndarray
-        Array of pixel x positions
-    pixel_pos_y: ndarray
-        Array of pixel x positions
-    phi: float
-        Rotation angle of pixels
-
-    Returns
-    -------
-        ndarray,ndarray: Transformed pixel x and y coordinates
-
-    """
-    # Here we are derotating, so we use a minus sign
-    psi = -psi
-    pixel_pos_rot_x = pixel_pos_x * np.cos(psi) - pixel_pos_y * np.sin(psi)
-    pixel_pos_rot_y = pixel_pos_x * np.sin(psi) + pixel_pos_y * np.cos(psi)
-    return pixel_pos_rot_x, pixel_pos_rot_y
-
-
-def timing_parameters(pix_x, pix_y, image, peak_time, rotation_angle):
+def timing_parameters(geom, image, peakpos, hillas_parameters):
     """
     Function to extract timing parameters from a cleaned image
 
     Parameters
     ----------
-    pix_x : array_like
-        Pixel x-coordinate
-    pix_y : array_like
-        Pixel y-coordinate
+    geom: ctapipe.instrument.CameraGeometry
+        Camera geometry
     image : array_like
-        Pixel values corresponding
-    peak_time : array_like
-        Pixel times corresponding
-    rotation_angle: float
-        Rotation angle for the image major axis, namely psi (in radians)
+        Pixel values
+    peakpos : array_like
+        Pixel peak positions array
+    hillas_parameters: ctapipe.io.containers.HillasParametersContainer
+        Result of hillas_parameters
 
     Returns
     -------
-    timing_parameters: TimingParameters
+    timing_parameters: TimingParametersContainer
     """
 
-    unit = Quantity(pix_x).unit
-    pix_x = Quantity(np.asanyarray(pix_x, dtype=np.float64)).value
-    pix_y = Quantity(np.asanyarray(pix_y, dtype=np.float64)).value
-    image = np.asanyarray(image, dtype=np.float64)
-    peak_time = np.asanyarray(peak_time, dtype=np.float64)
+    unit = geom.pix_x.unit
+    pix_x = geom.pix_x.value
+    pix_y = geom.pix_y.value
 
     # select only the pixels in the cleaned image that are greater than zero.
     # This is to allow to use a dilated mask (which might be better):
@@ -84,23 +44,20 @@ def timing_parameters(pix_x, pix_y, image, peak_time, rotation_angle):
     pix_x = pix_x[mask]
     pix_y = pix_y[mask]
     image = image[mask]
-    peak_time = peak_time[mask]
+    peakpos = peakpos[mask]
 
     assert pix_x.shape == image.shape
     assert pix_y.shape == image.shape
-    assert peak_time.shape == image.shape
-
-    # since the values of peak_pos are integers, sometimes the mask is constant
-    # for all the selected pixels. Asking for a mask to have at least 3 different
-    # values for the peak_pos in the selected pixels seems reasonable.
-
-    if np.unique(peak_time).size > 2:
-        pix_x_rot, pix_y_rot = rotate_translate(pix_x, pix_y, rotation_angle)
-        gradient, intercept = np.polyfit(x=pix_x_rot, y=peak_time,
-                                         deg=1, w=np.sqrt(image))
-    else:
-        gradient = 0.
-        intercept = 0.
-
-    return TimingParameters(gradient=gradient * (peak_time.unit / unit),
-                            intercept=intercept * unit)
+    assert peakpos.shape == image.shape
+
+    longi, trans = geom.get_shower_coordinates(
+        hillas_parameters.x,
+        hillas_parameters.y,
+        hillas_parameters.psi
+    )
+    slope, intercept = np.polyfit(longi, peakpos, deg=1, w=np.sqrt(image))
+
+    return TimingParametersContainer(
+        slope=slope / unit,
+        intercept=intercept,
+    )

ctapipe/io/containers.py

@@ -43,6 +43,7 @@
     'HillasParametersContainer',
     'LeakageContainer',
     'ConcentrationContainer',
+    'TimingParametersContainer',
 ]
 
 
@@ -718,3 +719,11 @@ class ConcentrationContainer(Container):
         nan,
         'Percentage of photo-electrons in the brightest pixel'
     )
+
+class TimingParametersContainer(Container):
+    """
+    Slope and Intercept of a linear regression of the arrival times
+    along the shower main axis
+    """
+    slope = Field(nan, 'Slope of arrival times along main shower axis')
+    intercept = Field(nan, 'intercept of arrival times along main shower axis')

ctapipe/visualization/tests/test_mpl.py

@@ -88,17 +88,19 @@ def test_array_display():
     grad = 2
     intercept = 1
 
+    geom = CameraGeometry.from_name("LSTCam")
     rot_angle = 20 * u.deg
+    hillas = HillasParametersContainer(x=0*u.m, y=0*u.m, psi=rot_angle)
+
     timing_rot20 = timing_parameters(
-        pix_x=arange(4) * u.deg,
-        pix_y=zeros(4) * u.deg,
-        image=ones(4),
-        peak_time=intercept * u.ns + grad * arange(4) * u.ns,
-        rotation_angle=rot_angle
+        geom,
+        image=ones(geom.n_pixels),
+        peakpos=intercept + grad * geom.pix_x.value,
+        hillas_parameters = hillas,
     )
     gradient_dict = {
-        1: timing_rot20.gradient.value,
-        2: timing_rot20.gradient.value,
+        1: timing_rot20.slope.value,
+        2: timing_rot20.slope.value,
     }
     ad.set_vector_hillas(hillas_dict=hillas_dict,
                          length=500,