diff --git a/ctapipe/image/tests/test_timing_parameters.py b/ctapipe/image/tests/test_timing_parameters.py index fc40446a7ed..0b708d6acb1 100644 --- a/ctapipe/image/tests/test_timing_parameters.py +++ b/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) diff --git a/ctapipe/image/timing_parameters.py b/ctapipe/image/timing_parameters.py index d03cf4ee6ab..3be78199d55 100644 --- a/ctapipe/image/timing_parameters.py +++ b/ctapipe/image/timing_parameters.py @@ -1,5 +1,3 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -# -*- coding: UTF-8 -*- """ Image timing-based shower image parametrization. """ @@ -7,74 +5,36 @@ 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, + ) \ No newline at end of file diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py index 3cb0191dfd2..78ee62d59f3 100644 --- a/ctapipe/io/containers.py +++ b/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') \ No newline at end of file diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 195fab4598f..4f68f92e9e9 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/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,