From 28708f095f40ee5e746818a1aa11f9388d5ca115 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 21 Feb 2019 10:48:59 +0100 Subject: [PATCH] Move get_shower_coordinates to hillas, fixes #934 (#976) * Move get_shower_coordinates to hillas --- ctapipe/image/concentration.py | 3 +- ctapipe/image/hillas.py | 37 +++++++++++++++++++ ctapipe/image/timing_parameters.py | 14 ++++--- ctapipe/instrument/camera.py | 27 -------------- .../lst_analysis_bootcamp_2018.ipynb | 9 +++-- 5 files changed, 53 insertions(+), 37 deletions(-) diff --git a/ctapipe/image/concentration.py b/ctapipe/image/concentration.py index 1615b666b24..0449f1fde1a 100644 --- a/ctapipe/image/concentration.py +++ b/ctapipe/image/concentration.py @@ -1,5 +1,6 @@ import numpy as np from ..io.containers import ConcentrationContainer +from .hillas import camera_to_shower_coordinates def concentration(geom, image, hillas_parameters): @@ -21,7 +22,7 @@ def concentration(geom, image, hillas_parameters): cog_pixels = np.argsort(delta_x**2 + delta_y**2) conc_cog = np.sum(image[cog_pixels[:3]]) / h.intensity - longi, trans = geom.get_shower_coordinates(h.x, h.y, h.psi) + longi, trans = camera_to_shower_coordinates(geom.pix_x, geom.pix_y, h.x, h.y, h.psi) # get all pixels inside the hillas ellipse mask_core = (longi**2 / h.length**2) + (trans**2 / h.width**2) <= 1.0 diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 436b4c0067b..2f157af32b3 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -17,6 +17,43 @@ ] +def camera_to_shower_coordinates(x, y, cog_x, cog_y, psi): + ''' + Return longitudinal and transverse coordinates for x and y + for a given set of hillas parameters + + Parameters + ---------- + x: u.Quantity[length] + x coordinate in camera coordinates + y: u.Quantity[length] + y coordinate in camera coordinates + cog_x: u.Quantity[length] + x coordinate of center of gravity + cog_y: u.Quantity[length] + y coordinate of center of gravity + psi: Angle + orientation angle + + Returns + ------- + longitudinal: astropy.units.Quantity + longitudinal coordinates (along the shower axis) + transverse: astropy.units.Quantity + transverse coordinates (perpendicular to the shower axis) + ''' + cos_psi = np.cos(psi) + sin_psi = np.sin(psi) + + delta_x = x - cog_x + delta_y = y - cog_y + + longi = delta_x * cos_psi + delta_y * sin_psi + trans = delta_x * -sin_psi + delta_y * cos_psi + + return longi, trans + + class HillasParameterizationError(RuntimeError): pass diff --git a/ctapipe/image/timing_parameters.py b/ctapipe/image/timing_parameters.py index 1b2fcfcf677..a4cedb625fe 100644 --- a/ctapipe/image/timing_parameters.py +++ b/ctapipe/image/timing_parameters.py @@ -4,6 +4,8 @@ import numpy as np from ctapipe.io.containers import TimingParametersContainer +from .hillas import camera_to_shower_coordinates + __all__ = [ 'timing_parameters' @@ -31,26 +33,26 @@ def timing_parameters(geom, image, peakpos, hillas_parameters): """ 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. # we need to exclude possible pixels with zero signal after cleaning. mask = image > 0 - pix_x = pix_x[mask] - pix_y = pix_y[mask] + pix_x = geom.pix_x[mask] + pix_y = geom.pix_y[mask] image = image[mask] peakpos = peakpos[mask] assert pix_x.shape == image.shape, 'image shape must match geometry' assert pix_x.shape == peakpos.shape, 'peakpos shape must match geometry' - longi, trans = geom.get_shower_coordinates( + longi, trans = camera_to_shower_coordinates( + pix_x, + pix_y, hillas_parameters.x, hillas_parameters.y, hillas_parameters.psi ) - slope, intercept = np.polyfit(longi[mask], peakpos, deg=1, w=np.sqrt(image)) + slope, intercept = np.polyfit(longi.value, peakpos, deg=1, w=np.sqrt(image)) return TimingParametersContainer( slope=slope / unit, diff --git a/ctapipe/instrument/camera.py b/ctapipe/instrument/camera.py index 6bdda973d7f..1be691b45bd 100644 --- a/ctapipe/instrument/camera.py +++ b/ctapipe/instrument/camera.py @@ -533,33 +533,6 @@ def get_border_pixel_mask(self, width=1): self.border_cache[width] = mask return mask - def get_shower_coordinates(self, x, y, psi): - ''' - Return longitudinal and transverse coordinates of the pixels - for a given set of hillas parameters - - Parameters - ---------- - hillas_parameters: ctapipe.io.containers.HilllasContainer - - Returns - ------- - longitudinal: astropy.units.Quantity - longitudinal coordinates (along the shower axis) - transverse: astropy.units.Quantity - transverse coordinates (perpendicular to the shower axis) - ''' - cos_psi = np.cos(psi) - sin_psi = np.sin(psi) - - delta_x = self.pix_x - x - delta_y = self.pix_y - y - - longi = delta_x * cos_psi + delta_y * sin_psi - trans = delta_x * -sin_psi + delta_y * cos_psi - - return longi, trans - def position_to_pix_index(self, x, y): ''' Return the index of a camera pixel which contains a given position (x,y) diff --git a/docs/tutorials/lst_analysis_bootcamp_2018.ipynb b/docs/tutorials/lst_analysis_bootcamp_2018.ipynb index d4e96565f7d..4ca37d39716 100644 --- a/docs/tutorials/lst_analysis_bootcamp_2018.ipynb +++ b/docs/tutorials/lst_analysis_bootcamp_2018.ipynb @@ -465,7 +465,8 @@ "source": [ "from ctapipe.image import hillas_parameters, leakage, concentration\n", "from ctapipe.image.timing_parameters import timing_parameters\n", - "from ctapipe.image.cleaning import number_of_islands" + "from ctapipe.image.cleaning import number_of_islands\n", + "from ctapipe.image.hillas import camera_to_shower_coordinates" ] }, { @@ -519,7 +520,9 @@ "metadata": {}, "outputs": [], "source": [ - "long, trans = camera.get_shower_coordinates(hillas.x, hillas.y, hillas.psi)\n", + "long, trans = camera_to_shower_coordinates(\n", + " camera.pix_x, camera.pix_y,hillas.x, hillas.y, hillas.psi\n", + ")\n", "\n", "plt.plot(long[clean], dl1.peakpos[0][clean], 'o')\n", "plt.plot(long[clean], timing.slope * long[clean] + timing.intercept)" @@ -1112,7 +1115,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.1" + "version": "3.6.8" } }, "nbformat": 4,