From 4028688f774a3c7da5849f4edc6c7abfc85d6730 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Tue, 2 Feb 2021 19:40:36 +0100 Subject: [PATCH 01/39] Cache files by full url/path --- ctapipe/utils/datasets.py | 4 ++-- ctapipe/utils/download.py | 19 ++++++++++--------- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/ctapipe/utils/datasets.py b/ctapipe/utils/datasets.py index b5af3d6a15c..12b632b86df 100644 --- a/ctapipe/utils/datasets.py +++ b/ctapipe/utils/datasets.py @@ -38,7 +38,7 @@ def get_searchpath_dirs(searchpath=os.getenv("CTAPIPE_SVC_PATH")): else: searchpaths = [Path(p) for p in os.path.expandvars(searchpath).split(":")] - searchpaths.append(get_cache_path("")) + searchpaths.append(get_cache_path(DEFAULT_URL)) return searchpaths @@ -175,7 +175,7 @@ def try_filetypes(basename, role, file_types, **kwargs): # look first in cache so we don't have to try non-existing downloads for ext, reader in file_types.items(): filename = basename + ext - cache_path = get_cache_path(filename) + cache_path = get_cache_path(os.path.join(DEFAULT_URL, filename)) if cache_path.exists(): path = cache_path break diff --git a/ctapipe/utils/download.py b/ctapipe/utils/download.py index b6b3f9eb18f..3500aa89176 100644 --- a/ctapipe/utils/download.py +++ b/ctapipe/utils/download.py @@ -61,14 +61,15 @@ def download_file(url, path, auth=None, chunk_size=10240, progress=False): part_file.rename(path) -def get_cache_path(path, cache_name="ctapipe", env_override="CTAPIPE_CACHE"): +def get_cache_path(url, cache_name="ctapipe", env_override="CTAPIPE_CACHE"): if os.getenv(env_override): base = Path(os.environ["CTAPIPE_CACHE"]) else: base = Path(os.environ["HOME"]) / ".cache" / cache_name - # need to make it relative - path = str(path).lstrip("/") + url = urlparse(url) + + path = os.path.join(url.netloc.rstrip("/"), url.path.lstrip("/")) path = base / path path.parent.mkdir(parents=True, exist_ok=True) return path @@ -109,18 +110,18 @@ def download_file_cached( path: pathlib.Path the full path to the downloaded data. """ - path = get_cache_path(name, cache_name=cache_name) + log.debug(f"File {name} is not available in cache, downloading.") + + base_url = os.environ.get(env_prefix + "URL", default_url).rstrip("/") + url = base_url + "/" + str(name).lstrip("/") + + path = get_cache_path(url, cache_name=cache_name) # if we already dowloaded the file, just use it if path.is_file(): log.debug(f"File {name} is available in cache.") return path - log.debug(f"File {name} is not available in cache, downloading.") - - base_url = os.environ.get(env_prefix + "URL", default_url).rstrip("/") - url = base_url + "/" + str(name).lstrip("/") - if auth is True: try: auth = ( From 7d839871f5eec97e8aceacc34065676e86e6981f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Tue, 2 Feb 2021 20:00:10 +0100 Subject: [PATCH 02/39] Fix part file not deleted on ctrl c --- ctapipe/utils/download.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/ctapipe/utils/download.py b/ctapipe/utils/download.py index 3500aa89176..1e87636b7a2 100644 --- a/ctapipe/utils/download.py +++ b/ctapipe/utils/download.py @@ -4,6 +4,7 @@ import logging from tqdm import tqdm from urllib.parse import urlparse +import time log = logging.getLogger(__name__) @@ -51,7 +52,7 @@ def download_file(url, path, auth=None, chunk_size=10240, progress=False): for chunk in r.iter_content(chunk_size=chunk_size): f.write(chunk) pbar.update(len(chunk)) - except Exception: + except: # we really want to catch everythin here # cleanup part file if something goes wrong if part_file.is_file(): part_file.unlink() @@ -116,6 +117,12 @@ def download_file_cached( url = base_url + "/" + str(name).lstrip("/") path = get_cache_path(url, cache_name=cache_name) + part_file = path.with_suffix(path.suffix + ".part") + + if part_file.is_file(): + log.warning("Another download for this file is already running, waiting.") + while part_file.is_file(): + time.sleep(1) # if we already dowloaded the file, just use it if path.is_file(): From b52bf0c2bf8e8cef873ffd6ce8d424a6283c3f1b Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Wed, 3 Feb 2021 19:43:16 +0100 Subject: [PATCH 03/39] Add container for nominal reconstruction of hillas parameters - x,y replaced by lon and lat - units changed to deg from m --- ctapipe/containers.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index b6e94e091f8..595aa460c31 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -122,6 +122,30 @@ class HillasParametersContainer(Container): kurtosis = Field(nan, "measure of the tailedness") +class NominalHillasParametersContainer(Container): + container_prefix = "hillas_nominal" + + intensity = Field(nan, "total intensity (size)") + + lon = Field(nan * u.deg, "centroid longitude in the nominal frame", unit=u.deg) + lat = Field(nan * u.deg, "centroid latitude in the nominal frame", unit=u.deg) + r = Field( + nan * u.deg, "radial coordinate of centroid in the nominal frame", unit=u.deg + ) + phi = Field( + nan * u.deg, "polar coordinate of centroid in the nominal frame", unit=u.deg + ) + + length = Field(nan * u.deg, "standard deviation along the major-axis", unit=u.deg) + length_uncertainty = Field(nan * u.deg, "uncertainty of length", unit=u.deg) + width = Field(nan * u.deg, "standard spread along the minor-axis", unit=u.deg) + width_uncertainty = Field(nan * u.deg, "uncertainty of width", unit=u.deg) + psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) + + skewness = Field(nan, "measure of the asymmetry") + kurtosis = Field(nan, "measure of the tailedness") + + class LeakageContainer(Container): """ Fraction of signal in 1 or 2-pixel width border from the edge of the From 4000a083e81af24eaab45c2ce02b36e4cdf144d7 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Wed, 3 Feb 2021 19:44:05 +0100 Subject: [PATCH 04/39] Add method to calculate hillas parameters in the nominal frame - Requires pixel positions in the nominal frame instead of a geom object - Most of the code is a duplicate of the reconstruction in the camera frame for now --- ctapipe/image/__init__.py | 1 + ctapipe/image/hillas.py | 117 ++++++++++++++++++++++++-- ctapipe/image/tests/test_hillas.py | 127 +++++++++++++++++++++++------ 3 files changed, 213 insertions(+), 32 deletions(-) diff --git a/ctapipe/image/__init__.py b/ctapipe/image/__init__.py index c2d9a3f7dd5..35a74cf4b41 100644 --- a/ctapipe/image/__init__.py +++ b/ctapipe/image/__init__.py @@ -1,5 +1,6 @@ from .hillas import ( hillas_parameters, + nominal_hillas_parameters, HillasParameterizationError, camera_to_shower_coordinates, ) diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 8cd87afaa08..ef2c7a70234 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -8,16 +8,13 @@ import numpy as np from astropy.coordinates import Angle from astropy.units import Quantity -from ..containers import HillasParametersContainer +from ..containers import HillasParametersContainer, NominalHillasParametersContainer HILLAS_ATOL = np.finfo(np.float64).eps -__all__ = [ - "hillas_parameters", - "HillasParameterizationError", -] +__all__ = ["hillas_parameters", "HillasParameterizationError"] def camera_to_shower_coordinates(x, y, cog_x, cog_y, psi): @@ -211,3 +208,113 @@ def hillas_parameters(geom, image): skewness=skewness_long, kurtosis=kurtosis_long, ) + + +def nominal_hillas_parameters(pixel_positions, image): + """ + Reconstruct hillas parameters given the pixel positions in the nominal frame. + """ + unit = pixel_positions.fov_lon.unit + pix_lon = Quantity(np.asanyarray(pixel_positions.fov_lon, dtype=np.float64)).value + pix_lat = Quantity(np.asanyarray(pixel_positions.fov_lat, dtype=np.float64)).value + image = np.asanyarray(image, dtype=np.float64) + image = np.ma.filled(image, 0) + msg = "Image and pixel shape do not match" + assert pix_lon.shape == pix_lat.shape == image.shape, msg + + size = np.sum(image) + + if size == 0.0: + raise HillasParameterizationError("size=0, cannot calculate HillasParameters") + + # calculate the cog as the mean of the coordinates weighted with the image + cog_lon = np.average(pix_lon, weights=image) + cog_lat = np.average(pix_lat, weights=image) + + # polar coordinates of the cog + cog_r = np.linalg.norm([cog_lon, cog_lat]) + cog_phi = np.arctan2(cog_lon, cog_lat) + + # do the PCA for the hillas parameters + delta_lon = pix_lon - cog_lon + delta_lat = pix_lat - cog_lat + + # The ddof=0 makes this comparable to the other methods, + # but ddof=1 should be more correct, mostly affects small showers + # on a percent level + cov = np.cov(delta_lon, delta_lat, aweights=image, ddof=0) + eig_vals, eig_vecs = np.linalg.eigh(cov) + + # round eig_vals to get rid of nans when eig val is something like -8.47032947e-22 + near_zero = np.isclose(eig_vals, 0, atol=HILLAS_ATOL) + eig_vals[near_zero] = 0 + + # width and length are eigen values of the PCA + width, length = np.sqrt(eig_vals) + + # psi is the angle of the eigenvector to length to the x-axis + vx, vy = eig_vecs[0, 1], eig_vecs[1, 1] + + # avoid divide by 0 warnings + if length == 0: + psi = skewness_long = kurtosis_long = np.nan + else: + if vx != 0: + psi = np.arctan(vy / vx) + else: + psi = np.pi / 2 + + # calculate higher order moments along shower axes + longitudinal = delta_lon * np.cos(psi) + delta_lat * np.sin(psi) + + m3_long = np.average(longitudinal ** 3, weights=image) + skewness_long = m3_long / length ** 3 + + m4_long = np.average(longitudinal ** 4, weights=image) + kurtosis_long = m4_long / length ** 4 + + # Compute of the Hillas parameters uncertainties. + # Implementation described in [hillas_uncertainties]_ This is an internal MAGIC document + # not generally accessible. + + # intermediate variables + cos_2psi = np.cos(2 * psi) + a = (1 + cos_2psi) / 2 + b = (1 - cos_2psi) / 2 + c = np.sin(2 * psi) + + A = ((delta_lon ** 2.0) - cov[0][0]) / size + B = ((delta_lat ** 2.0) - cov[1][1]) / size + C = ((delta_lon * delta_lat) - cov[0][1]) / size + + # Hillas's uncertainties + + # avoid divide by 0 warnings + if length == 0: + length_uncertainty = np.nan + else: + length_uncertainty = np.sqrt( + np.sum(((((a * A) + (b * B) + (c * C))) ** 2.0) * image) + ) / (2 * length) + + if width == 0: + width_uncertainty = np.nan + else: + width_uncertainty = np.sqrt( + np.sum(((((b * A) + (a * B) + (-c * C))) ** 2.0) * image) + ) / (2 * width) + + return NominalHillasParametersContainer( + lon=u.Quantity(cog_lon, unit), + lat=u.Quantity(cog_lat, unit), + r=u.Quantity(cog_r, unit), + phi=Angle(cog_phi, unit=u.rad), + intensity=size, + length=u.Quantity(length, unit), + length_uncertainty=u.Quantity(length_uncertainty, unit), + width=u.Quantity(width, unit), + width_uncertainty=u.Quantity(width_uncertainty, unit), + psi=Angle(psi, unit=u.rad), + skewness=skewness_long, + kurtosis=kurtosis_long, + ) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index 4e557d68b2f..f1648835066 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -1,8 +1,13 @@ from ctapipe.instrument import CameraGeometry from ctapipe.image import tailcuts_clean, toymodel -from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError +from ctapipe.image.hillas import ( + hillas_parameters, + nominal_hillas_parameters, + HillasParameterizationError, +) from ctapipe.containers import HillasParametersContainer from astropy.coordinates import Angle +from ctapipe.coordinates import NominalFrame, CameraFrame from astropy import units as u import numpy as np from numpy import isclose, zeros_like @@ -10,6 +15,9 @@ from pytest import approx import itertools import pytest +from astropy.time import Time +from astropy.coordinates import EarthLocation +from astropy.coordinates import SkyCoord, AltAz def create_sample_image( @@ -28,11 +36,7 @@ def create_sample_image( model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) # generate toymodel image in camera for this shower model. - image, _, _ = model.generate_image( - geom, - intensity=1500, - nsb_level_pe=3, - ) + image, _, _ = model.generate_image(geom, intensity=1500, nsb_level_pe=3) # calculate pixels likely containing signal clean_mask = tailcuts_clean(geom, image, 10, 5) @@ -131,18 +135,10 @@ def test_with_toy(): for psi in psis: # make a toymodel shower model - model = toymodel.Gaussian( - x=x, - y=y, - width=width, - length=length, - psi=psi, - ) + model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) image, signal, noise = model.generate_image( - geom, - intensity=intensity, - nsb_level_pe=5, + geom, intensity=intensity, nsb_level_pe=5 ) result = hillas_parameters(geom, signal) @@ -178,19 +174,10 @@ def test_skewness(): for x, y, psi, skew in itertools.product(xs, ys, psis, skews): # make a toymodel shower model model = toymodel.SkewedGaussian( - x=x, - y=y, - width=width, - length=length, - psi=psi, - skewness=skew, + x=x, y=y, width=width, length=length, psi=psi, skewness=skew ) - _, signal, _ = model.generate_image( - geom, - intensity=intensity, - nsb_level_pe=5, - ) + _, signal, _ = model.generate_image(geom, intensity=intensity, nsb_level_pe=5) result = hillas_parameters(geom, signal) @@ -268,3 +255,89 @@ def test_single_pixel(): assert hillas.length.value == 0 assert hillas.width.value == 0 assert np.isnan(hillas.psi) + + +def test_reconstruction_in_nominal_frame(): + np.random.seed(42) + + obstime = Time("2013-11-01T03:00") + location = EarthLocation.of_site("Roque de los Muchachos") + focal_length = 28 * u.m + horizon_frame = AltAz(location=location, obstime=obstime) + pointing = SkyCoord(alt=70 * u.deg, az=180 * u.deg, frame=horizon_frame) + nominal_frame = NominalFrame(origin=pointing, obstime=obstime, location=location) + camera_frame = CameraFrame( + telescope_pointing=pointing, + focal_length=focal_length, + obstime=obstime, + location=location, + ) + + geom = CameraGeometry.from_name("LSTCam") + nominal_pixel_positions = SkyCoord( + x=geom.pix_x, y=geom.pix_y, frame=camera_frame + ).transform_to(nominal_frame) + + width = 0.03 * u.m + length = 0.15 * u.m + intensity = 500 + + xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m) + ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m) + psis = Angle([-90, -45, 0, 45, 90], unit="deg") + + def distance(coord): + return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2 + + for x, y in zip(xs, ys): + for psi in psis: + # make a toymodel shower model + model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) + + image, signal, noise = model.generate_image( + geom, intensity=intensity, nsb_level_pe=5 + ) + + nominal_result = nominal_hillas_parameters(nominal_pixel_positions, signal) + assert u.isclose(np.abs(nominal_result.lon), 1 * u.deg, rtol=0.1) + assert u.isclose(np.abs(nominal_result.lat), 1 * u.deg, rtol=0.1) + assert u.isclose(nominal_result.width, 0.06 * u.deg, rtol=0.1) + assert u.isclose(nominal_result.width_uncertainty, 0.002 * u.deg, rtol=0.4) + assert u.isclose(nominal_result.length, 0.3 * u.deg, rtol=0.1) + assert u.isclose(nominal_result.length_uncertainty, 0.01 * u.deg, rtol=0.4) + assert signal.sum() == nominal_result.intensity + + # Compare results with calculation in the camera frame + camera_result = hillas_parameters(geom, signal) + + transformed_cog = SkyCoord( + fov_lon=nominal_result.lon, + fov_lat=nominal_result.lat, + frame=nominal_frame, + ).transform_to(camera_frame) + assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01) + assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01) + + main_edges = u.Quantity([-nominal_result.length, nominal_result.length]) + main_lon = main_edges * np.cos(nominal_result.psi) + nominal_result.lon + main_lat = main_edges * np.sin(nominal_result.psi) + nominal_result.lat + cam_main_axis = SkyCoord( + fov_lon=main_lon, fov_lat=main_lat, frame=nominal_frame + ).transform_to(camera_frame) + transformed_length = distance(cam_main_axis) + assert u.isclose(transformed_length, camera_result.length, rtol=0.01) + + secondary_edges = u.Quantity( + [-nominal_result.length, nominal_result.length] + ) + secondary_lon = ( + secondary_edges * np.cos(nominal_result.psi) + nominal_result.lon + ) + secondary_lat = ( + secondary_edges * np.sin(nominal_result.psi) + nominal_result.lat + ) + cam_secondary_edges = SkyCoord( + fov_lon=secondary_lon, fov_lat=secondary_lat, frame=nominal_frame + ).transform_to(camera_frame) + transformed_width = distance(cam_secondary_edges) + assert u.isclose(transformed_width, camera_result.length, rtol=0.01) From 6aa2eb2482f5c605d29bd042e421c5783884b740 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 4 Feb 2021 16:28:36 +0100 Subject: [PATCH 05/39] Remove second hillas parameters function - Decide which container to return based on pixel position unit --- ctapipe/image/__init__.py | 1 - ctapipe/image/hillas.py | 152 ++++++----------------------- ctapipe/image/tests/test_hillas.py | 16 ++- 3 files changed, 35 insertions(+), 134 deletions(-) diff --git a/ctapipe/image/__init__.py b/ctapipe/image/__init__.py index 35a74cf4b41..c2d9a3f7dd5 100644 --- a/ctapipe/image/__init__.py +++ b/ctapipe/image/__init__.py @@ -1,6 +1,5 @@ from .hillas import ( hillas_parameters, - nominal_hillas_parameters, HillasParameterizationError, camera_to_shower_coordinates, ) diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index ef2c7a70234..7fe9ed8e17c 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -194,127 +194,33 @@ def hillas_parameters(geom, image): np.sum(((((b * A) + (a * B) + (-c * C))) ** 2.0) * image) ) / (2 * width) - return HillasParametersContainer( - x=u.Quantity(cog_x, unit), - y=u.Quantity(cog_y, unit), - r=u.Quantity(cog_r, unit), - phi=Angle(cog_phi, unit=u.rad), - intensity=size, - length=u.Quantity(length, unit), - length_uncertainty=u.Quantity(length_uncertainty, unit), - width=u.Quantity(width, unit), - width_uncertainty=u.Quantity(width_uncertainty, unit), - psi=Angle(psi, unit=u.rad), - skewness=skewness_long, - kurtosis=kurtosis_long, - ) - - -def nominal_hillas_parameters(pixel_positions, image): - """ - Reconstruct hillas parameters given the pixel positions in the nominal frame. - """ - unit = pixel_positions.fov_lon.unit - pix_lon = Quantity(np.asanyarray(pixel_positions.fov_lon, dtype=np.float64)).value - pix_lat = Quantity(np.asanyarray(pixel_positions.fov_lat, dtype=np.float64)).value - image = np.asanyarray(image, dtype=np.float64) - image = np.ma.filled(image, 0) - msg = "Image and pixel shape do not match" - assert pix_lon.shape == pix_lat.shape == image.shape, msg - - size = np.sum(image) - - if size == 0.0: - raise HillasParameterizationError("size=0, cannot calculate HillasParameters") - - # calculate the cog as the mean of the coordinates weighted with the image - cog_lon = np.average(pix_lon, weights=image) - cog_lat = np.average(pix_lat, weights=image) - - # polar coordinates of the cog - cog_r = np.linalg.norm([cog_lon, cog_lat]) - cog_phi = np.arctan2(cog_lon, cog_lat) - - # do the PCA for the hillas parameters - delta_lon = pix_lon - cog_lon - delta_lat = pix_lat - cog_lat - - # The ddof=0 makes this comparable to the other methods, - # but ddof=1 should be more correct, mostly affects small showers - # on a percent level - cov = np.cov(delta_lon, delta_lat, aweights=image, ddof=0) - eig_vals, eig_vecs = np.linalg.eigh(cov) - - # round eig_vals to get rid of nans when eig val is something like -8.47032947e-22 - near_zero = np.isclose(eig_vals, 0, atol=HILLAS_ATOL) - eig_vals[near_zero] = 0 - - # width and length are eigen values of the PCA - width, length = np.sqrt(eig_vals) - - # psi is the angle of the eigenvector to length to the x-axis - vx, vy = eig_vecs[0, 1], eig_vecs[1, 1] - - # avoid divide by 0 warnings - if length == 0: - psi = skewness_long = kurtosis_long = np.nan - else: - if vx != 0: - psi = np.arctan(vy / vx) - else: - psi = np.pi / 2 - - # calculate higher order moments along shower axes - longitudinal = delta_lon * np.cos(psi) + delta_lat * np.sin(psi) - - m3_long = np.average(longitudinal ** 3, weights=image) - skewness_long = m3_long / length ** 3 - - m4_long = np.average(longitudinal ** 4, weights=image) - kurtosis_long = m4_long / length ** 4 - - # Compute of the Hillas parameters uncertainties. - # Implementation described in [hillas_uncertainties]_ This is an internal MAGIC document - # not generally accessible. - - # intermediate variables - cos_2psi = np.cos(2 * psi) - a = (1 + cos_2psi) / 2 - b = (1 - cos_2psi) / 2 - c = np.sin(2 * psi) - - A = ((delta_lon ** 2.0) - cov[0][0]) / size - B = ((delta_lat ** 2.0) - cov[1][1]) / size - C = ((delta_lon * delta_lat) - cov[0][1]) / size - - # Hillas's uncertainties - - # avoid divide by 0 warnings - if length == 0: - length_uncertainty = np.nan + if unit.is_equivalent(u.m): + return HillasParametersContainer( + x=u.Quantity(cog_x, unit), + y=u.Quantity(cog_y, unit), + r=u.Quantity(cog_r, unit), + phi=Angle(cog_phi, unit=u.rad), + intensity=size, + length=u.Quantity(length, unit), + length_uncertainty=u.Quantity(length_uncertainty, unit), + width=u.Quantity(width, unit), + width_uncertainty=u.Quantity(width_uncertainty, unit), + psi=Angle(psi, unit=u.rad), + skewness=skewness_long, + kurtosis=kurtosis_long, + ) else: - length_uncertainty = np.sqrt( - np.sum(((((a * A) + (b * B) + (c * C))) ** 2.0) * image) - ) / (2 * length) - - if width == 0: - width_uncertainty = np.nan - else: - width_uncertainty = np.sqrt( - np.sum(((((b * A) + (a * B) + (-c * C))) ** 2.0) * image) - ) / (2 * width) - - return NominalHillasParametersContainer( - lon=u.Quantity(cog_lon, unit), - lat=u.Quantity(cog_lat, unit), - r=u.Quantity(cog_r, unit), - phi=Angle(cog_phi, unit=u.rad), - intensity=size, - length=u.Quantity(length, unit), - length_uncertainty=u.Quantity(length_uncertainty, unit), - width=u.Quantity(width, unit), - width_uncertainty=u.Quantity(width_uncertainty, unit), - psi=Angle(psi, unit=u.rad), - skewness=skewness_long, - kurtosis=kurtosis_long, - ) + return NominalHillasParametersContainer( + lon=u.Quantity(cog_x, unit), + lat=u.Quantity(cog_y, unit), + r=u.Quantity(cog_r, unit), + phi=Angle(cog_phi, unit=u.rad), + intensity=size, + length=u.Quantity(length, unit), + length_uncertainty=u.Quantity(length_uncertainty, unit), + width=u.Quantity(width, unit), + width_uncertainty=u.Quantity(width_uncertainty, unit), + psi=Angle(psi, unit=u.rad), + skewness=skewness_long, + kurtosis=kurtosis_long, + ) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index f1648835066..fce6b91b726 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -1,10 +1,6 @@ from ctapipe.instrument import CameraGeometry from ctapipe.image import tailcuts_clean, toymodel -from ctapipe.image.hillas import ( - hillas_parameters, - nominal_hillas_parameters, - HillasParameterizationError, -) +from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError from ctapipe.containers import HillasParametersContainer from astropy.coordinates import Angle from ctapipe.coordinates import NominalFrame, CameraFrame @@ -274,9 +270,6 @@ def test_reconstruction_in_nominal_frame(): ) geom = CameraGeometry.from_name("LSTCam") - nominal_pixel_positions = SkyCoord( - x=geom.pix_x, y=geom.pix_y, frame=camera_frame - ).transform_to(nominal_frame) width = 0.03 * u.m length = 0.15 * u.m @@ -289,7 +282,7 @@ def test_reconstruction_in_nominal_frame(): def distance(coord): return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2 - for x, y in zip(xs, ys): + for i, (x, y) in enumerate(zip(xs, ys)): for psi in psis: # make a toymodel shower model model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) @@ -298,7 +291,10 @@ def distance(coord): geom, intensity=intensity, nsb_level_pe=5 ) - nominal_result = nominal_hillas_parameters(nominal_pixel_positions, signal) + geom.frame = camera_frame + geom_nom = geom.transform_to(nominal_frame) + + nominal_result = hillas_parameters(geom_nom, signal) assert u.isclose(np.abs(nominal_result.lon), 1 * u.deg, rtol=0.1) assert u.isclose(np.abs(nominal_result.lat), 1 * u.deg, rtol=0.1) assert u.isclose(nominal_result.width, 0.06 * u.deg, rtol=0.1) From f554d8fd3d9a9d7ccba5f681513b50b7ec3c1fa5 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 4 Feb 2021 21:00:04 +0100 Subject: [PATCH 06/39] Add containers for nominal frame and fix hillas-related parameters --- ctapipe/containers.py | 43 ++++++++++++++++++++++++++++++++-- ctapipe/image/concentration.py | 21 ++++++++++++----- ctapipe/image/timing.py | 31 ++++++++++++++++++------ 3 files changed, 80 insertions(+), 15 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 595aa460c31..7674cdb8c6c 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -29,6 +29,8 @@ "MonitoringCameraContainer", "MonitoringContainer", "MorphologyContainer", + "NominalHillasParametersContainer", + "NominalTimingParametersContainer", "ParticleClassificationContainer", "PedestalContainer", "PixelStatusContainer", @@ -204,6 +206,24 @@ class TimingParametersContainer(Container): ) +class NominalTimingParametersContainer(Container): + """ + Slope and Intercept of a linear regression of the arrival times + along the shower main axis in the nominal frame + """ + + container_prefix = "nominal_timing" + slope = Field( + nan / u.deg, "Slope of arrival times along main shower axis", unit=1 / u.deg + ) + intercept = Field(nan, "intercept of arrival times along main shower axis") + deviation = Field( + nan, + "Root-mean-square deviation of the pulse times " + "with respect to the predicted time", + ) + + class MorphologyContainer(Container): """ Parameters related to pixels surviving image cleaning """ @@ -250,6 +270,23 @@ class ImageParametersContainer(Container): ) +class NominalImageParametersContainer(Container): + """ Collection of image parameters """ + + container_prefix = "params" + hillas = Field(NominalHillasParametersContainer(), "Hillas Parameters") + timing = Field(NominalTimingParametersContainer(), "Timing Parameters") + leakage = Field(LeakageContainer(), "Leakage Parameters") + concentration = Field(ConcentrationContainer(), "Concentration Parameters") + morphology = Field(MorphologyContainer(), "Image Morphology Parameters") + intensity_statistics = Field( + IntensityStatisticsContainer(), "Intensity image statistics" + ) + peak_time_statistics = Field( + PeakTimeStatisticsContainer(), "Peak time image statistics" + ) + + class DL1CameraContainer(Container): """ Storage of output of camera calibration e.g the final calibrated @@ -278,7 +315,7 @@ class DL1CameraContainer(Container): ndim=1, ) - parameters = Field(None, "Image parameters", type=ImageParametersContainer) + parameters = Field(None, "Image parameters", type=NominalImageParametersContainer) class DL1Container(Container): @@ -431,7 +468,9 @@ class SimulatedCameraContainer(Container): ) true_parameters = Field( - None, "Parameters derived from the true_image", type=ImageParametersContainer + None, + "Parameters derived from the true_image", + type=NominalImageParametersContainer, ) diff --git a/ctapipe/image/concentration.py b/ctapipe/image/concentration.py index 956dc2a94a8..8aa5c0da651 100644 --- a/ctapipe/image/concentration.py +++ b/ctapipe/image/concentration.py @@ -1,7 +1,11 @@ import numpy as np import astropy.units as u -from ..containers import ConcentrationContainer +from ..containers import ( + ConcentrationContainer, + HillasParametersContainer, + NominalHillasParametersContainer, +) from .hillas import camera_to_shower_coordinates from ..utils.quantities import all_to_value @@ -19,11 +23,16 @@ def concentration_parameters(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 - ) + if isinstance(h, HillasParametersContainer): + 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 + ) + elif isinstance(h, NominalHillasParametersContainer): + unit = h.lon.unit + pix_x, pix_y, x, y, length, width = all_to_value( + geom.pix_x, geom.pix_y, h.lon, h.lat, h.length, h.width, unit=unit + ) delta_x = pix_x - x delta_y = pix_y - y diff --git a/ctapipe/image/timing.py b/ctapipe/image/timing.py index e6fa99c8942..87b1404ba0b 100644 --- a/ctapipe/image/timing.py +++ b/ctapipe/image/timing.py @@ -4,7 +4,12 @@ import numpy as np import astropy.units as u -from ..containers import TimingParametersContainer +from ..containers import ( + TimingParametersContainer, + NominalTimingParametersContainer, + HillasParametersContainer, + NominalHillasParametersContainer, +) from .hillas import camera_to_shower_coordinates from ..utils.quantities import all_to_value from ..fitting import lts_linear_regression @@ -58,9 +63,16 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N raise ValueError("The non-masked pixels must verify signal >= 0") 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 - ) + if isinstance(h, HillasParametersContainer): + 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 + ) + elif isinstance(h, NominalHillasParametersContainer): + unit = h.lon.unit + pix_x, pix_y, x, y, length, width = all_to_value( + geom.pix_x, geom.pix_y, h.lon, h.lat, h.length, h.width, unit=unit + ) longi, _ = camera_to_shower_coordinates( pix_x, pix_y, x, y, hillas_parameters.psi.to_value(u.rad) @@ -73,6 +85,11 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N # recalculate for all points deviation = rmse(longi * beta[0] + beta[1], peak_time) - return TimingParametersContainer( - slope=beta[0] / unit, intercept=beta[1], deviation=deviation - ) + if unit.is_equivalent(u.m): + return TimingParametersContainer( + slope=beta[0] / unit, intercept=beta[1], deviation=deviation + ) + else: + return NominalTimingParametersContainer( + slope=beta[0] / unit, intercept=beta[1], deviation=deviation + ) From 4f61b739394fe96f99e23d5c8fb6581d6788bede Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 4 Feb 2021 21:02:30 +0100 Subject: [PATCH 07/39] Use image parameters in nominal frame, bump data level version - Save image parameters calculated in the nominal frame - This affects the hillas parameters and the timing slope - The coordinate transformations happen in the image processor --- ctapipe/image/image_processor.py | 50 +++++++++++++---- ctapipe/io/dl1eventsource.py | 75 ++++++++++++++++++------- ctapipe/io/dl1writer.py | 3 +- ctapipe/io/tests/test_dl1eventsource.py | 2 +- ctapipe/tools/tests/test_tools.py | 4 +- 5 files changed, 99 insertions(+), 35 deletions(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index b461469341f..a5aa55724b0 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -1,14 +1,17 @@ """ + High level image processing (ImageProcessor Component) """ - +from astropy.coordinates import SkyCoord, AltAz +from astropy.time import Time +from ctapipe.coordinates import NominalFrame, CameraFrame, EarthLocation from ..containers import ( ArrayEventContainer, - ImageParametersContainer, IntensityStatisticsContainer, + NominalImageParametersContainer, + NominalTimingParametersContainer, PeakTimeStatisticsContainer, - TimingParametersContainer, ) from ..core import QualityQuery, TelescopeComponent from ..core.traits import List, create_class_enum_trait @@ -24,8 +27,8 @@ ) -DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer() -DEFAULT_TIMING_PARAMETERS = TimingParametersContainer() +DEFAULT_IMAGE_PARAMETERS = NominalImageParametersContainer() +DEFAULT_TIMING_PARAMETERS = NominalTimingParametersContainer() DEFAULT_PEAKTIME_STATISTICS = PeakTimeStatisticsContainer() @@ -88,8 +91,8 @@ def __call__(self, event: ArrayEventContainer): self._process_telescope_event(event) def _parameterize_image( - self, tel_id, image, signal_pixels, peak_time=None - ) -> ImageParametersContainer: + self, tel_id, image, signal_pixels, geometry, peak_time=None + ) -> NominalImageParametersContainer: """Apply image cleaning and calculate image features Parameters @@ -105,12 +108,10 @@ def _parameterize_image( Returns ------- - ImageParametersContainer: + NominalImageParametersContainer: cleaning mask, parameters """ - tel = self.subarray.tel[tel_id] - geometry = tel.camera.geometry image_selected = image[signal_pixels] # check if image can be parameterized: @@ -151,7 +152,7 @@ def _parameterize_image( timing = DEFAULT_TIMING_PARAMETERS peak_time_statistics = DEFAULT_PEAKTIME_STATISTICS - return ImageParametersContainer( + return NominalImageParametersContainer( hillas=hillas, timing=timing, leakage=leakage, @@ -169,8 +170,33 @@ def _process_telescope_event(self, event): """ Loop over telescopes and process the calibrated images into parameters """ + location = EarthLocation.of_site("Roque de los Muchachos") + obstime = Time(event.trigger.time, format="mjd") + altaz = AltAz(location=location, obstime=obstime) + array_pointing = SkyCoord( + alt=event.pointing.array_azimuth, + az=event.pointing.array_altitude, + frame=altaz, + ) + nominal_frame = NominalFrame( + origin=array_pointing, obstime=obstime, location=location + ) for tel_id, dl1_camera in event.dl1.tel.items(): + tel_pointings = SkyCoord( + alt=event.pointing.tel[tel_id].azimuth, + az=event.pointing.tel[tel_id].altitude, + frame=altaz, + ) + camera_frame = CameraFrame( + telescope_pointing=tel_pointings, + focal_length=self.subarray.tel[tel_id].optics.equivalent_focal_length, + obstime=obstime, + location=location, + ) + geometry = self.subarray.tel[tel_id].camera.geometry + geometry.frame = camera_frame + geometry_nominal = geometry.transform_to(nominal_frame) # compute image parameters only if requested to write them dl1_camera.image_mask = self.clean( @@ -184,6 +210,7 @@ def _process_telescope_event(self, event): image=dl1_camera.image, signal_pixels=dl1_camera.image_mask, peak_time=dl1_camera.peak_time, + geometry=geometry_nominal, ) self.log.debug("params: %s", dl1_camera.parameters.as_dict(recursive=True)) @@ -197,6 +224,7 @@ def _process_telescope_event(self, event): tel_id, image=sim_camera.true_image, signal_pixels=sim_camera.true_image > 0, + geometry=geometry_nominal, peak_time=None, # true image from simulation has no peak time ) self.log.debug( diff --git a/ctapipe/io/dl1eventsource.py b/ctapipe/io/dl1eventsource.py index 544c83415dc..2e4c3c4d029 100644 --- a/ctapipe/io/dl1eventsource.py +++ b/ctapipe/io/dl1eventsource.py @@ -12,6 +12,7 @@ ArrayEventContainer, EventIndexContainer, HillasParametersContainer, + NominalHillasParametersContainer, IntensityStatisticsContainer, LeakageContainer, MorphologyContainer, @@ -19,8 +20,10 @@ SimulatedShowerContainer, PeakTimeStatisticsContainer, TimingParametersContainer, + NominalTimingParametersContainer, TriggerContainer, ImageParametersContainer, + NominalImageParametersContainer, ) from ctapipe.utils import IndexFinder @@ -28,7 +31,7 @@ logger = logging.getLogger(__name__) -COMPATIBLE_DL1_VERSIONS = ["v1.0.0", "v1.0.1", "v1.0.2", "v1.0.3"] +COMPATIBLE_DL1_VERSIONS = ["v1.0.0", "v1.0.1", "v1.0.2", "v1.0.3", "v1.0.4"] class DL1EventSource(EventSource): @@ -228,8 +231,16 @@ def _generate_events(self): tel.name: HDF5TableReader(self.file_).read( f"/dl1/event/telescope/parameters/{tel.name}", containers=[ - HillasParametersContainer(), - TimingParametersContainer(), + ( + NominalHillasParametersContainer() + if (self.datamodel_version >= "1.0.4") + else HillasParametersContainer() + ), + ( + NominalTimingParametersContainer() + if (self.datamodel_version >= "1.0.4") + else TimingParametersContainer() + ), LeakageContainer(), ConcentrationContainer(), MorphologyContainer(), @@ -245,7 +256,11 @@ def _generate_events(self): tel.name: HDF5TableReader(self.file_).read( f"/simulation/event/telescope/parameters/{tel.name}", containers=[ - HillasParametersContainer(), + ( + NominalHillasParametersContainer() + if (self.datamodel_version >= "1.0.4") + else HillasParametersContainer() + ), LeakageContainer(), ConcentrationContainer(), MorphologyContainer(), @@ -353,15 +368,26 @@ def _generate_events(self): # Best would probbaly be if we could directly read # into the ImageParametersContainer params = next(param_readers[f"tel_{tel:03d}"]) - dl1.parameters = ImageParametersContainer( - hillas=params[0], - timing=params[1], - leakage=params[2], - concentration=params[3], - morphology=params[4], - intensity_statistics=params[5], - peak_time_statistics=params[6], - ) + if self.datamodel_version >= "1.0.4": + dl1.parameters = NominalImageParametersContainer( + hillas=params[0], + timing=params[1], + leakage=params[2], + concentration=params[3], + morphology=params[4], + intensity_statistics=params[5], + peak_time_statistics=params[6], + ) + else: + dl1.parameters = ImageParametersContainer( + hillas=params[0], + timing=params[1], + leakage=params[2], + concentration=params[3], + morphology=params[4], + intensity_statistics=params[5], + peak_time_statistics=params[6], + ) if self.has_simulated_dl1: if f"tel_{tel:03d}" not in param_readers.keys(): @@ -374,13 +400,22 @@ def _generate_events(self): simulated_params = next( simulated_param_readers[f"tel_{tel:03d}"] ) - simulated.true_parameters = ImageParametersContainer( - hillas=simulated_params[0], - leakage=simulated_params[1], - concentration=simulated_params[2], - morphology=simulated_params[3], - intensity_statistics=simulated_params[4], - ) + if self.datamodel_version >= "1.0.4": + simulated.true_parameters = NominalImageParametersContainer( + hillas=simulated_params[0], + leakage=simulated_params[1], + concentration=simulated_params[2], + morphology=simulated_params[3], + intensity_statistics=simulated_params[4], + ) + else: + simulated.true_parameters = ImageParametersContainer( + hillas=simulated_params[0], + leakage=simulated_params[1], + concentration=simulated_params[2], + morphology=simulated_params[3], + intensity_statistics=simulated_params[4], + ) yield data diff --git a/ctapipe/io/dl1writer.py b/ctapipe/io/dl1writer.py index 0afd67a02c7..b18a2ce2480 100644 --- a/ctapipe/io/dl1writer.py +++ b/ctapipe/io/dl1writer.py @@ -28,9 +28,10 @@ tables.parameters.NODE_CACHE_SLOTS = 3000 # fixes problem with too many datasets -DL1_DATA_MODEL_VERSION = "v1.0.3" +DL1_DATA_MODEL_VERSION = "v1.0.4" DL1_DATA_MODEL_CHANGE_HISTORY = """ - v1.0.3: true_image dtype changed from float32 to int32 +- v1.0.4: hillas and timing parameters saved in nominal frame (degree) as opposed to camera frame (m) """ PROV = Provenance() diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index 2ac2239de0b..830e8af243d 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -79,7 +79,7 @@ def test_dl1_data(dl1_file): for event in source: for tel in event.dl1.tel: assert event.dl1.tel[tel].image.any() - assert event.dl1.tel[tel].parameters.hillas.x != np.nan + assert event.dl1.tel[tel].parameters.hillas.lon != np.nan def test_pointing(dl1_file): diff --git a/ctapipe/tools/tests/test_tools.py b/ctapipe/tools/tests/test_tools.py index 255c663f6e2..f4ddf8295ec 100644 --- a/ctapipe/tools/tests/test_tools.py +++ b/ctapipe/tools/tests/test_tools.py @@ -283,7 +283,7 @@ def test_merge(tmpdir): "obs_id", "event_id", "tel_id", - "hillas_intensity", + "hillas_nominal_intensity", "concentration_cog", "leakage_pixels_width_1", ) @@ -368,7 +368,7 @@ def test_stage_1_dl1(tmpdir, dl1_image_file, dl1_parameters_file): "obs_id", "event_id", "tel_id", - "hillas_intensity", + "hillas_nominal_intensity", "concentration_cog", "leakage_pixels_width_1", ) From 3863addd1fe5afe7d6535c663535855dfc56917f Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 5 Feb 2021 11:02:16 +0100 Subject: [PATCH 08/39] Fix import --- ctapipe/image/image_processor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index a5aa55724b0..a08f15a5988 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -4,7 +4,8 @@ """ from astropy.coordinates import SkyCoord, AltAz from astropy.time import Time -from ctapipe.coordinates import NominalFrame, CameraFrame, EarthLocation +from ctapipe.coordinates import NominalFrame, CameraFrame +from astropy.coordinates import EarthLocation from ..containers import ( ArrayEventContainer, From 34f9d5f8438b59f8aa8b4bf5dfe5509d8c9fb442 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 5 Feb 2021 14:38:40 +0100 Subject: [PATCH 09/39] Clean up container structure, use telescope frame for image parameters - Stage 1 and the connected tools now calculate the hillas parameters in the telescope frame - Hillas and TimingParametersContainer now store values as deg. Old containers persist as CameraHillas/TimingParametersContainer - In order to keep a single ImageParametersContainer, base containers have been added for hillas and timing parameters - Failing tests adapted: x/y -> lon/lat, reconstruction algorithms use the CameraHillasParametersContainer --- ctapipe/containers.py | 108 ++++++++---------- ctapipe/image/concentration.py | 6 +- ctapipe/image/hillas.py | 6 +- ctapipe/image/image_processor.py | 50 ++------ ctapipe/image/tests/test_hillas.py | 80 +++++++------ ctapipe/image/tests/test_timing_parameters.py | 8 +- ctapipe/image/timing.py | 12 +- ctapipe/io/dl1eventsource.py | 70 ++++-------- ctapipe/io/dl1writer.py | 2 +- ctapipe/io/tests/test_hdf5.py | 26 +++-- ctapipe/reco/tests/test_ImPACT.py | 4 +- .../reco/tests/test_hillas_intersection.py | 24 ++-- ctapipe/tools/tests/test_tools.py | 4 +- ctapipe/visualization/tests/test_mpl.py | 8 +- 14 files changed, 174 insertions(+), 234 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 7674cdb8c6c..62fcadeda5f 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -29,8 +29,8 @@ "MonitoringCameraContainer", "MonitoringContainer", "MorphologyContainer", - "NominalHillasParametersContainer", - "NominalTimingParametersContainer", + "CameraHillasParametersContainer", + "CameraTimingParametersContainer", "ParticleClassificationContainer", "PedestalContainer", "PixelStatusContainer", @@ -104,11 +104,14 @@ class TelEventIndexContainer(Container): tel_id = Field(0, "telescope identifier") -class HillasParametersContainer(Container): - container_prefix = "hillas" - +class BaseHillasParametersContainer(Container): intensity = Field(nan, "total intensity (size)") + skewness = Field(nan, "measure of the asymmetry") + kurtosis = Field(nan, "measure of the tailedness") + +class CameraHillasParametersContainer(BaseHillasParametersContainer): + container_prefix = "camera_frame_hillas" x = Field(nan * u.m, "centroid x coordinate", unit=u.m) y = Field(nan * u.m, "centroid x coordinate", unit=u.m) r = Field(nan * u.m, "radial coordinate of centroid", unit=u.m) @@ -120,22 +123,28 @@ class HillasParametersContainer(Container): width_uncertainty = Field(nan * u.m, "uncertainty of width", unit=u.m) psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) - skewness = Field(nan, "measure of the asymmetry") - kurtosis = Field(nan, "measure of the tailedness") - - -class NominalHillasParametersContainer(Container): - container_prefix = "hillas_nominal" - - intensity = Field(nan, "total intensity (size)") - lon = Field(nan * u.deg, "centroid longitude in the nominal frame", unit=u.deg) - lat = Field(nan * u.deg, "centroid latitude in the nominal frame", unit=u.deg) +class HillasParametersContainer(BaseHillasParametersContainer): + container_prefix = "hillas" + lon = Field( + nan * u.deg, + "centroid longitude in a sky frame e.g. the telescope frame", + unit=u.deg, + ) + lat = Field( + nan * u.deg, + "centroid latitude in a sky frame e.g. the telescope frame", + unit=u.deg, + ) r = Field( - nan * u.deg, "radial coordinate of centroid in the nominal frame", unit=u.deg + nan * u.deg, + "radial coordinate of centroid in a sky frame e.g. the telescope frame", + unit=u.deg, ) phi = Field( - nan * u.deg, "polar coordinate of centroid in the nominal frame", unit=u.deg + nan * u.deg, + "polar coordinate of centroid in a sky frame e.g. the telescope frame", + unit=u.deg, ) length = Field(nan * u.deg, "standard deviation along the major-axis", unit=u.deg) @@ -144,9 +153,6 @@ class NominalHillasParametersContainer(Container): width_uncertainty = Field(nan * u.deg, "uncertainty of width", unit=u.deg) psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) - skewness = Field(nan, "measure of the asymmetry") - kurtosis = Field(nan, "measure of the tailedness") - class LeakageContainer(Container): """ @@ -188,40 +194,37 @@ class ConcentrationContainer(Container): pixel = Field(nan, "Percentage of photo-electrons in the brightest pixel") -class TimingParametersContainer(Container): +class BaseTimingParametersContainer(Container): + intercept = Field(nan, "intercept of arrival times along main shower axis") + deviation = Field( + nan, + "Root-mean-square deviation of the pulse times " + "with respect to the predicted time", + ) + + +class CameraTimingParametersContainer(BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times - along the shower main axis + along the shower main axis in the camera frame. """ - container_prefix = "timing" + container_prefix = "camera_frame_timing" slope = Field( nan / u.m, "Slope of arrival times along main shower axis", unit=1 / u.m ) - intercept = Field(nan, "intercept of arrival times along main shower axis") - deviation = Field( - nan, - "Root-mean-square deviation of the pulse times " - "with respect to the predicted time", - ) -class NominalTimingParametersContainer(Container): +class TimingParametersContainer(BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times - along the shower main axis in the nominal frame + along the shower main axis in a sky frame e.g. the telescope frame. """ - container_prefix = "nominal_timing" + container_prefix = "timing" slope = Field( nan / u.deg, "Slope of arrival times along main shower axis", unit=1 / u.deg ) - intercept = Field(nan, "intercept of arrival times along main shower axis") - deviation = Field( - nan, - "Root-mean-square deviation of the pulse times " - "with respect to the predicted time", - ) class MorphologyContainer(Container): @@ -257,25 +260,8 @@ class ImageParametersContainer(Container): """ Collection of image parameters """ container_prefix = "params" - hillas = Field(HillasParametersContainer(), "Hillas Parameters") - timing = Field(TimingParametersContainer(), "Timing Parameters") - leakage = Field(LeakageContainer(), "Leakage Parameters") - concentration = Field(ConcentrationContainer(), "Concentration Parameters") - morphology = Field(MorphologyContainer(), "Image Morphology Parameters") - intensity_statistics = Field( - IntensityStatisticsContainer(), "Intensity image statistics" - ) - peak_time_statistics = Field( - PeakTimeStatisticsContainer(), "Peak time image statistics" - ) - - -class NominalImageParametersContainer(Container): - """ Collection of image parameters """ - - container_prefix = "params" - hillas = Field(NominalHillasParametersContainer(), "Hillas Parameters") - timing = Field(NominalTimingParametersContainer(), "Timing Parameters") + hillas = Field(BaseHillasParametersContainer(), "Hillas Parameters") + timing = Field(BaseTimingParametersContainer(), "Timing Parameters") leakage = Field(LeakageContainer(), "Leakage Parameters") concentration = Field(ConcentrationContainer(), "Concentration Parameters") morphology = Field(MorphologyContainer(), "Image Morphology Parameters") @@ -306,7 +292,6 @@ class DL1CameraContainer(Container): dtype=np.float32, ndim=1, ) - image_mask = Field( None, "Boolean numpy array where True means the pixel has passed cleaning." @@ -314,8 +299,7 @@ class DL1CameraContainer(Container): dtype=np.bool, ndim=1, ) - - parameters = Field(None, "Image parameters", type=NominalImageParametersContainer) + parameters = Field(None, "Image parameters", type=ImageParametersContainer) class DL1Container(Container): @@ -468,9 +452,7 @@ class SimulatedCameraContainer(Container): ) true_parameters = Field( - None, - "Parameters derived from the true_image", - type=NominalImageParametersContainer, + None, "Parameters derived from the true_image", type=ImageParametersContainer ) diff --git a/ctapipe/image/concentration.py b/ctapipe/image/concentration.py index 8aa5c0da651..43258bd4ea9 100644 --- a/ctapipe/image/concentration.py +++ b/ctapipe/image/concentration.py @@ -4,7 +4,7 @@ from ..containers import ( ConcentrationContainer, HillasParametersContainer, - NominalHillasParametersContainer, + CameraHillasParametersContainer, ) from .hillas import camera_to_shower_coordinates from ..utils.quantities import all_to_value @@ -23,12 +23,12 @@ def concentration_parameters(geom, image, hillas_parameters): """ h = hillas_parameters - if isinstance(h, HillasParametersContainer): + if isinstance(h, CameraHillasParametersContainer): 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 ) - elif isinstance(h, NominalHillasParametersContainer): + elif isinstance(h, HillasParametersContainer): unit = h.lon.unit pix_x, pix_y, x, y, length, width = all_to_value( geom.pix_x, geom.pix_y, h.lon, h.lat, h.length, h.width, unit=unit diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 7fe9ed8e17c..13bed79c0a9 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -8,7 +8,7 @@ import numpy as np from astropy.coordinates import Angle from astropy.units import Quantity -from ..containers import HillasParametersContainer, NominalHillasParametersContainer +from ..containers import CameraHillasParametersContainer, HillasParametersContainer HILLAS_ATOL = np.finfo(np.float64).eps @@ -195,7 +195,7 @@ def hillas_parameters(geom, image): ) / (2 * width) if unit.is_equivalent(u.m): - return HillasParametersContainer( + return CameraHillasParametersContainer( x=u.Quantity(cog_x, unit), y=u.Quantity(cog_y, unit), r=u.Quantity(cog_r, unit), @@ -210,7 +210,7 @@ def hillas_parameters(geom, image): kurtosis=kurtosis_long, ) else: - return NominalHillasParametersContainer( + return HillasParametersContainer( lon=u.Quantity(cog_x, unit), lat=u.Quantity(cog_y, unit), r=u.Quantity(cog_r, unit), diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index a08f15a5988..803845933b1 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -2,16 +2,13 @@ High level image processing (ImageProcessor Component) """ -from astropy.coordinates import SkyCoord, AltAz -from astropy.time import Time -from ctapipe.coordinates import NominalFrame, CameraFrame -from astropy.coordinates import EarthLocation +from ctapipe.coordinates import TelescopeFrame from ..containers import ( ArrayEventContainer, IntensityStatisticsContainer, - NominalImageParametersContainer, - NominalTimingParametersContainer, + ImageParametersContainer, + TimingParametersContainer, PeakTimeStatisticsContainer, ) from ..core import QualityQuery, TelescopeComponent @@ -28,8 +25,8 @@ ) -DEFAULT_IMAGE_PARAMETERS = NominalImageParametersContainer() -DEFAULT_TIMING_PARAMETERS = NominalTimingParametersContainer() +DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer() +DEFAULT_TIMING_PARAMETERS = TimingParametersContainer() DEFAULT_PEAKTIME_STATISTICS = PeakTimeStatisticsContainer() @@ -87,13 +84,14 @@ def __init__( ) self.check_image = ImageQualityQuery(parent=self) self._is_simulation = is_simulation + self.telescope_frame = TelescopeFrame() def __call__(self, event: ArrayEventContainer): self._process_telescope_event(event) def _parameterize_image( self, tel_id, image, signal_pixels, geometry, peak_time=None - ) -> NominalImageParametersContainer: + ) -> ImageParametersContainer: """Apply image cleaning and calculate image features Parameters @@ -109,7 +107,7 @@ def _parameterize_image( Returns ------- - NominalImageParametersContainer: + ImageParametersContainer: cleaning mask, parameters """ @@ -153,7 +151,7 @@ def _parameterize_image( timing = DEFAULT_TIMING_PARAMETERS peak_time_statistics = DEFAULT_PEAKTIME_STATISTICS - return NominalImageParametersContainer( + return ImageParametersContainer( hillas=hillas, timing=timing, leakage=leakage, @@ -171,33 +169,9 @@ def _process_telescope_event(self, event): """ Loop over telescopes and process the calibrated images into parameters """ - location = EarthLocation.of_site("Roque de los Muchachos") - obstime = Time(event.trigger.time, format="mjd") - altaz = AltAz(location=location, obstime=obstime) - array_pointing = SkyCoord( - alt=event.pointing.array_azimuth, - az=event.pointing.array_altitude, - frame=altaz, - ) - nominal_frame = NominalFrame( - origin=array_pointing, obstime=obstime, location=location - ) - for tel_id, dl1_camera in event.dl1.tel.items(): - tel_pointings = SkyCoord( - alt=event.pointing.tel[tel_id].azimuth, - az=event.pointing.tel[tel_id].altitude, - frame=altaz, - ) - camera_frame = CameraFrame( - telescope_pointing=tel_pointings, - focal_length=self.subarray.tel[tel_id].optics.equivalent_focal_length, - obstime=obstime, - location=location, - ) geometry = self.subarray.tel[tel_id].camera.geometry - geometry.frame = camera_frame - geometry_nominal = geometry.transform_to(nominal_frame) + geometry_telescope = geometry.transform_to(self.telescope_frame) # compute image parameters only if requested to write them dl1_camera.image_mask = self.clean( @@ -211,7 +185,7 @@ def _process_telescope_event(self, event): image=dl1_camera.image, signal_pixels=dl1_camera.image_mask, peak_time=dl1_camera.peak_time, - geometry=geometry_nominal, + geometry=geometry_telescope, ) self.log.debug("params: %s", dl1_camera.parameters.as_dict(recursive=True)) @@ -225,7 +199,7 @@ def _process_telescope_event(self, event): tel_id, image=sim_camera.true_image, signal_pixels=sim_camera.true_image > 0, - geometry=geometry_nominal, + geometry=geometry_telescope, peak_time=None, # true image from simulation has no peak time ) self.log.debug( diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index fce6b91b726..07e21c79db8 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -1,9 +1,12 @@ from ctapipe.instrument import CameraGeometry from ctapipe.image import tailcuts_clean, toymodel from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError -from ctapipe.containers import HillasParametersContainer -from astropy.coordinates import Angle -from ctapipe.coordinates import NominalFrame, CameraFrame +from ctapipe.containers import ( + HillasParametersContainer, + CameraHillasParametersContainer, +) +from astropy.coordinates import Angle, SkyCoord +from ctapipe.coordinates import TelescopeFrame, CameraFrame from astropy import units as u import numpy as np from numpy import isclose, zeros_like @@ -11,9 +14,6 @@ from pytest import approx import itertools import pytest -from astropy.time import Time -from astropy.coordinates import EarthLocation -from astropy.coordinates import SkyCoord, AltAz def create_sample_image( @@ -109,6 +109,11 @@ def test_hillas_container(): geom, image = create_sample_image_zeros(psi="0d") params = hillas_parameters(geom, image) + assert isinstance(params, CameraHillasParametersContainer) + + geom.frame = CameraFrame(focal_length=28 * u.m) + geom_telescope_frame = geom.transform_to(TelescopeFrame()) + params = hillas_parameters(geom_telescope_frame, image) assert isinstance(params, HillasParametersContainer) @@ -253,23 +258,15 @@ def test_single_pixel(): assert np.isnan(hillas.psi) -def test_reconstruction_in_nominal_frame(): +def test_reconstruction_in_telescope_frame(): np.random.seed(42) - obstime = Time("2013-11-01T03:00") - location = EarthLocation.of_site("Roque de los Muchachos") - focal_length = 28 * u.m - horizon_frame = AltAz(location=location, obstime=obstime) - pointing = SkyCoord(alt=70 * u.deg, az=180 * u.deg, frame=horizon_frame) - nominal_frame = NominalFrame(origin=pointing, obstime=obstime, location=location) - camera_frame = CameraFrame( - telescope_pointing=pointing, - focal_length=focal_length, - obstime=obstime, - location=location, - ) + telescope_frame = TelescopeFrame() + camera_frame = CameraFrame(focal_length=28 * u.m) geom = CameraGeometry.from_name("LSTCam") + geom.frame = camera_frame + geom_nom = geom.transform_to(telescope_frame) width = 0.03 * u.m length = 0.15 * u.m @@ -291,49 +288,50 @@ def distance(coord): geom, intensity=intensity, nsb_level_pe=5 ) - geom.frame = camera_frame - geom_nom = geom.transform_to(nominal_frame) - - nominal_result = hillas_parameters(geom_nom, signal) - assert u.isclose(np.abs(nominal_result.lon), 1 * u.deg, rtol=0.1) - assert u.isclose(np.abs(nominal_result.lat), 1 * u.deg, rtol=0.1) - assert u.isclose(nominal_result.width, 0.06 * u.deg, rtol=0.1) - assert u.isclose(nominal_result.width_uncertainty, 0.002 * u.deg, rtol=0.4) - assert u.isclose(nominal_result.length, 0.3 * u.deg, rtol=0.1) - assert u.isclose(nominal_result.length_uncertainty, 0.01 * u.deg, rtol=0.4) - assert signal.sum() == nominal_result.intensity + telescope_result = hillas_parameters(geom_nom, signal) + assert u.isclose(np.abs(telescope_result.lon), 1 * u.deg, rtol=0.1) + assert u.isclose(np.abs(telescope_result.lat), 1 * u.deg, rtol=0.1) + assert u.isclose(telescope_result.width, 0.06 * u.deg, rtol=0.1) + assert u.isclose( + telescope_result.width_uncertainty, 0.002 * u.deg, rtol=0.4 + ) + assert u.isclose(telescope_result.length, 0.3 * u.deg, rtol=0.1) + assert u.isclose( + telescope_result.length_uncertainty, 0.01 * u.deg, rtol=0.4 + ) + assert signal.sum() == telescope_result.intensity # Compare results with calculation in the camera frame camera_result = hillas_parameters(geom, signal) transformed_cog = SkyCoord( - fov_lon=nominal_result.lon, - fov_lat=nominal_result.lat, - frame=nominal_frame, + fov_lon=telescope_result.lon, + fov_lat=telescope_result.lat, + frame=telescope_frame, ).transform_to(camera_frame) assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01) assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01) - main_edges = u.Quantity([-nominal_result.length, nominal_result.length]) - main_lon = main_edges * np.cos(nominal_result.psi) + nominal_result.lon - main_lat = main_edges * np.sin(nominal_result.psi) + nominal_result.lat + main_edges = u.Quantity([-telescope_result.length, telescope_result.length]) + main_lon = main_edges * np.cos(telescope_result.psi) + telescope_result.lon + main_lat = main_edges * np.sin(telescope_result.psi) + telescope_result.lat cam_main_axis = SkyCoord( - fov_lon=main_lon, fov_lat=main_lat, frame=nominal_frame + fov_lon=main_lon, fov_lat=main_lat, frame=telescope_frame ).transform_to(camera_frame) transformed_length = distance(cam_main_axis) assert u.isclose(transformed_length, camera_result.length, rtol=0.01) secondary_edges = u.Quantity( - [-nominal_result.length, nominal_result.length] + [-telescope_result.length, telescope_result.length] ) secondary_lon = ( - secondary_edges * np.cos(nominal_result.psi) + nominal_result.lon + secondary_edges * np.cos(telescope_result.psi) + telescope_result.lon ) secondary_lat = ( - secondary_edges * np.sin(nominal_result.psi) + nominal_result.lat + secondary_edges * np.sin(telescope_result.psi) + telescope_result.lat ) cam_secondary_edges = SkyCoord( - fov_lon=secondary_lon, fov_lat=secondary_lat, frame=nominal_frame + fov_lon=secondary_lon, fov_lat=secondary_lat, frame=telescope_frame ).transform_to(camera_frame) transformed_width = distance(cam_secondary_edges) assert u.isclose(transformed_width, camera_result.length, rtol=0.01) diff --git a/ctapipe/image/tests/test_timing_parameters.py b/ctapipe/image/tests/test_timing_parameters.py index 268a80666f7..132e5ab2225 100644 --- a/ctapipe/image/tests/test_timing_parameters.py +++ b/ctapipe/image/tests/test_timing_parameters.py @@ -2,7 +2,7 @@ import astropy.units as u from numpy.testing import assert_allclose from ctapipe.instrument.camera import CameraGeometry -from ctapipe.containers import HillasParametersContainer +from ctapipe.containers import CameraHillasParametersContainer def test_psi_0(): @@ -17,7 +17,7 @@ def test_psi_0(): deviation = 0.1 geom = CameraGeometry.from_name("LSTCam") - hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) + hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) random = np.random.RandomState(1) peak_time = intercept + grad * geom.pix_x.value @@ -47,7 +47,7 @@ def test_psi_20(): geom = CameraGeometry.from_name("LSTCam") psi = 20 * u.deg - hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=psi) + hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=psi) random = np.random.RandomState(1) peak_time = intercept + grad * ( @@ -77,7 +77,7 @@ def test_ignore_negative(): deviation = 0.1 geom = CameraGeometry.from_name("LSTCam") - hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) + hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) random = np.random.RandomState(1) peak_time = intercept + grad * geom.pix_x.value diff --git a/ctapipe/image/timing.py b/ctapipe/image/timing.py index 87b1404ba0b..c5f128d4749 100644 --- a/ctapipe/image/timing.py +++ b/ctapipe/image/timing.py @@ -5,10 +5,10 @@ import numpy as np import astropy.units as u from ..containers import ( + CameraTimingParametersContainer, TimingParametersContainer, - NominalTimingParametersContainer, + CameraHillasParametersContainer, HillasParametersContainer, - NominalHillasParametersContainer, ) from .hillas import camera_to_shower_coordinates from ..utils.quantities import all_to_value @@ -63,12 +63,12 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N raise ValueError("The non-masked pixels must verify signal >= 0") h = hillas_parameters - if isinstance(h, HillasParametersContainer): + if isinstance(h, CameraHillasParametersContainer): 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 ) - elif isinstance(h, NominalHillasParametersContainer): + elif isinstance(h, HillasParametersContainer): unit = h.lon.unit pix_x, pix_y, x, y, length, width = all_to_value( geom.pix_x, geom.pix_y, h.lon, h.lat, h.length, h.width, unit=unit @@ -86,10 +86,10 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N deviation = rmse(longi * beta[0] + beta[1], peak_time) if unit.is_equivalent(u.m): - return TimingParametersContainer( + return CameraTimingParametersContainer( slope=beta[0] / unit, intercept=beta[1], deviation=deviation ) else: - return NominalTimingParametersContainer( + return TimingParametersContainer( slope=beta[0] / unit, intercept=beta[1], deviation=deviation ) diff --git a/ctapipe/io/dl1eventsource.py b/ctapipe/io/dl1eventsource.py index 2e4c3c4d029..c65e11b6e63 100644 --- a/ctapipe/io/dl1eventsource.py +++ b/ctapipe/io/dl1eventsource.py @@ -11,19 +11,18 @@ ConcentrationContainer, ArrayEventContainer, EventIndexContainer, + CameraHillasParametersContainer, HillasParametersContainer, - NominalHillasParametersContainer, IntensityStatisticsContainer, LeakageContainer, MorphologyContainer, SimulationConfigContainer, SimulatedShowerContainer, PeakTimeStatisticsContainer, + CameraTimingParametersContainer, TimingParametersContainer, - NominalTimingParametersContainer, TriggerContainer, ImageParametersContainer, - NominalImageParametersContainer, ) from ctapipe.utils import IndexFinder @@ -232,14 +231,14 @@ def _generate_events(self): f"/dl1/event/telescope/parameters/{tel.name}", containers=[ ( - NominalHillasParametersContainer() + HillasParametersContainer() if (self.datamodel_version >= "1.0.4") - else HillasParametersContainer() + else CameraHillasParametersContainer() ), ( - NominalTimingParametersContainer() + TimingParametersContainer() if (self.datamodel_version >= "1.0.4") - else TimingParametersContainer() + else CameraTimingParametersContainer() ), LeakageContainer(), ConcentrationContainer(), @@ -257,9 +256,9 @@ def _generate_events(self): f"/simulation/event/telescope/parameters/{tel.name}", containers=[ ( - NominalHillasParametersContainer() + HillasParametersContainer() if (self.datamodel_version >= "1.0.4") - else HillasParametersContainer() + else CameraHillasParametersContainer() ), LeakageContainer(), ConcentrationContainer(), @@ -368,27 +367,15 @@ def _generate_events(self): # Best would probbaly be if we could directly read # into the ImageParametersContainer params = next(param_readers[f"tel_{tel:03d}"]) - if self.datamodel_version >= "1.0.4": - dl1.parameters = NominalImageParametersContainer( - hillas=params[0], - timing=params[1], - leakage=params[2], - concentration=params[3], - morphology=params[4], - intensity_statistics=params[5], - peak_time_statistics=params[6], - ) - else: - dl1.parameters = ImageParametersContainer( - hillas=params[0], - timing=params[1], - leakage=params[2], - concentration=params[3], - morphology=params[4], - intensity_statistics=params[5], - peak_time_statistics=params[6], - ) - + dl1.parameters = ImageParametersContainer( + hillas=params[0], + timing=params[1], + leakage=params[2], + concentration=params[3], + morphology=params[4], + intensity_statistics=params[5], + peak_time_statistics=params[6], + ) if self.has_simulated_dl1: if f"tel_{tel:03d}" not in param_readers.keys(): logger.debug( @@ -400,22 +387,13 @@ def _generate_events(self): simulated_params = next( simulated_param_readers[f"tel_{tel:03d}"] ) - if self.datamodel_version >= "1.0.4": - simulated.true_parameters = NominalImageParametersContainer( - hillas=simulated_params[0], - leakage=simulated_params[1], - concentration=simulated_params[2], - morphology=simulated_params[3], - intensity_statistics=simulated_params[4], - ) - else: - simulated.true_parameters = ImageParametersContainer( - hillas=simulated_params[0], - leakage=simulated_params[1], - concentration=simulated_params[2], - morphology=simulated_params[3], - intensity_statistics=simulated_params[4], - ) + simulated.true_parameters = ImageParametersContainer( + hillas=simulated_params[0], + leakage=simulated_params[1], + concentration=simulated_params[2], + morphology=simulated_params[3], + intensity_statistics=simulated_params[4], + ) yield data diff --git a/ctapipe/io/dl1writer.py b/ctapipe/io/dl1writer.py index b18a2ce2480..5443d3fc00e 100644 --- a/ctapipe/io/dl1writer.py +++ b/ctapipe/io/dl1writer.py @@ -31,7 +31,7 @@ DL1_DATA_MODEL_VERSION = "v1.0.4" DL1_DATA_MODEL_CHANGE_HISTORY = """ - v1.0.3: true_image dtype changed from float32 to int32 -- v1.0.4: hillas and timing parameters saved in nominal frame (degree) as opposed to camera frame (m) +- v1.0.4: hillas and timing parameters saved in telescope frame (degree) as opposed to camera frame (m) """ PROV = Provenance() diff --git a/ctapipe/io/tests/test_hdf5.py b/ctapipe/io/tests/test_hdf5.py index a6485d02347..976db925888 100644 --- a/ctapipe/io/tests/test_hdf5.py +++ b/ctapipe/io/tests/test_hdf5.py @@ -52,7 +52,7 @@ def test_write_container(temp_h5_file): def test_read_multiple_containers(): hillas_parameter_container = HillasParametersContainer( - x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m + lon=1 * u.deg, lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) leakage_container = LeakageContainer( @@ -66,7 +66,7 @@ def test_read_multiple_containers(): writer.write("params", [hillas_parameter_container, leakage_container]) df = pd.read_hdf(f.name, key="/dl1/params") - assert "hillas_x" in df.columns + assert "hillas_lon" in df.columns assert "leakage_pixels_width_1" in df.columns # test reading both containers separately @@ -110,7 +110,7 @@ def test_read_multiple_containers(): def test_read_without_prefixes(): hillas_parameter_container = HillasParametersContainer( - x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m + lon=1 * u.deg, lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) leakage_container = LeakageContainer( @@ -124,7 +124,7 @@ def test_read_without_prefixes(): writer.write("params", [hillas_parameter_container, leakage_container]) df = pd.read_hdf(f.name, key="/dl1/params") - assert "x" in df.columns + assert "lon" in df.columns assert "pixels_width_1" in df.columns # call with prefixes=False @@ -168,10 +168,18 @@ def test_read_without_prefixes(): def test_read_duplicated_container_types(): hillas_config_1 = HillasParametersContainer( - x=1 * u.m, y=2 * u.m, length=3 * u.m, width=4 * u.m, prefix="hillas_1" + lon=1 * u.deg, + lat=2 * u.deg, + length=3 * u.deg, + width=4 * u.deg, + prefix="hillas_1", ) hillas_config_2 = HillasParametersContainer( - x=2 * u.m, y=3 * u.m, length=4 * u.m, width=5 * u.m, prefix="hillas_2" + lon=2 * u.deg, + lat=3 * u.deg, + length=4 * u.deg, + width=5 * u.deg, + prefix="hillas_2", ) with tempfile.NamedTemporaryFile() as f: @@ -179,8 +187,8 @@ def test_read_duplicated_container_types(): writer.write("params", [hillas_config_1, hillas_config_2]) df = pd.read_hdf(f.name, key="/dl1/params") - assert "hillas_1_x" in df.columns - assert "hillas_2_x" in df.columns + assert "hillas_1_lon" in df.columns + assert "hillas_2_lat" in df.columns with HDF5TableReader(f.name) as reader: generator = reader.read( @@ -203,7 +211,7 @@ def test_read_duplicated_container_types(): def test_custom_prefix(): container = HillasParametersContainer( - x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m + lon=1 * u.deg, lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) container.prefix = "custom" with tempfile.NamedTemporaryFile() as f: diff --git a/ctapipe/reco/tests/test_ImPACT.py b/ctapipe/reco/tests/test_ImPACT.py index 4990d4f50a5..90d8749a39d 100644 --- a/ctapipe/reco/tests/test_ImPACT.py +++ b/ctapipe/reco/tests/test_ImPACT.py @@ -8,7 +8,7 @@ ReconstructedEnergyContainer, ) from ctapipe.reco.ImPACT import ImPACTReconstructor -from ctapipe.containers import HillasParametersContainer +from ctapipe.containers import CameraHillasParametersContainer from astropy.coordinates import Angle, AltAz, SkyCoord @@ -18,7 +18,7 @@ def setup_class(self): self.impact_reco = ImPACTReconstructor(root_dir=".") self.horizon_frame = AltAz() - self.h1 = HillasParametersContainer( + self.h1 = CameraHillasParametersContainer( x=1 * u.deg, y=1 * u.deg, r=1 * u.deg, diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py index 38c3fd027a5..dad559005f2 100644 --- a/ctapipe/reco/tests/test_hillas_intersection.py +++ b/ctapipe/reco/tests/test_hillas_intersection.py @@ -4,7 +4,7 @@ import numpy as np from astropy.coordinates import SkyCoord from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame -from ctapipe.containers import HillasParametersContainer +from ctapipe.containers import CameraHillasParametersContainer from ctapipe.io import EventSource @@ -81,12 +81,12 @@ def test_intersection_xmax_reco(): focal_length = 28 * u.m hillas_dict = { - 1: HillasParametersContainer( + 1: CameraHillasParametersContainer( x=-(delta / focal_length) * u.rad, y=((0 * u.m) / focal_length) * u.rad, intensity=1, ), - 2: HillasParametersContainer( + 2: CameraHillasParametersContainer( x=((0 * u.m) / focal_length) * u.rad, y=-(delta / focal_length) * u.rad, intensity=1, @@ -119,9 +119,9 @@ def test_intersection_reco_impact_point_tilted(): tel_y_dict = {1: delta, 2: delta, 3: -delta} hillas_dict = { - 1: HillasParametersContainer(intensity=100, psi=-90 * u.deg), - 2: HillasParametersContainer(intensity=100, psi=-45 * u.deg), - 3: HillasParametersContainer(intensity=100, psi=0 * u.deg), + 1: CameraHillasParametersContainer(intensity=100, psi=-90 * u.deg), + 2: CameraHillasParametersContainer(intensity=100, psi=-45 * u.deg), + 3: CameraHillasParametersContainer(intensity=100, psi=0 * u.deg), } reco_konrad = hill_inter.reconstruct_tilted( @@ -144,9 +144,9 @@ def test_intersection_weighting_spoiled_parameters(): # telescope 2 have a spoiled reconstruction (45 instead of -45) hillas_dict = { - 1: HillasParametersContainer(intensity=10000, psi=-90 * u.deg), - 2: HillasParametersContainer(intensity=1, psi=45 * u.deg), - 3: HillasParametersContainer(intensity=10000, psi=0 * u.deg), + 1: CameraHillasParametersContainer(intensity=10000, psi=-90 * u.deg), + 2: CameraHillasParametersContainer(intensity=1, psi=45 * u.deg), + 3: CameraHillasParametersContainer(intensity=10000, psi=0 * u.deg), } reco_konrad_spoiled = hill_inter.reconstruct_tilted( @@ -189,21 +189,21 @@ def test_intersection_nominal_reconstruction(): cog_coords_nom_3 = cog_coords_camera_3.transform_to(nominal_frame) # x-axis is along the altitude and y-axis is along the azimuth - hillas_1 = HillasParametersContainer( + hillas_1 = CameraHillasParametersContainer( x=cog_coords_nom_1.fov_lat, y=cog_coords_nom_1.fov_lon, intensity=100, psi=0 * u.deg, ) - hillas_2 = HillasParametersContainer( + hillas_2 = CameraHillasParametersContainer( x=cog_coords_nom_2.fov_lat, y=cog_coords_nom_2.fov_lon, intensity=100, psi=45 * u.deg, ) - hillas_3 = HillasParametersContainer( + hillas_3 = CameraHillasParametersContainer( x=cog_coords_nom_3.fov_lat, y=cog_coords_nom_3.fov_lon, intensity=100, diff --git a/ctapipe/tools/tests/test_tools.py b/ctapipe/tools/tests/test_tools.py index f4ddf8295ec..255c663f6e2 100644 --- a/ctapipe/tools/tests/test_tools.py +++ b/ctapipe/tools/tests/test_tools.py @@ -283,7 +283,7 @@ def test_merge(tmpdir): "obs_id", "event_id", "tel_id", - "hillas_nominal_intensity", + "hillas_intensity", "concentration_cog", "leakage_pixels_width_1", ) @@ -368,7 +368,7 @@ def test_stage_1_dl1(tmpdir, dl1_image_file, dl1_parameters_file): "obs_id", "event_id", "tel_id", - "hillas_nominal_intensity", + "hillas_intensity", "concentration_cog", "leakage_pixels_width_1", ) diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 89a1b60ca0c..2b0f21d2211 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -10,7 +10,7 @@ TelescopeDescription, PixelShape, ) -from ctapipe.containers import HillasParametersContainer +from ctapipe.containers import CameraHillasParametersContainer import numpy as np from astropy import units as u @@ -103,8 +103,8 @@ def test_array_display(): # test using hillas params: hillas_dict = { - 1: HillasParametersContainer(length=100.0 * u.m, psi=90 * u.deg), - 2: HillasParametersContainer(length=20000 * u.cm, psi="95deg"), + 1: CameraHillasParametersContainer(length=100.0 * u.m, psi=90 * u.deg), + 2: CameraHillasParametersContainer(length=20000 * u.cm, psi="95deg"), } grad = 2 @@ -112,7 +112,7 @@ def test_array_display(): geom = CameraGeometry.from_name("LSTCam") rot_angle = 20 * u.deg - hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=rot_angle) + hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=rot_angle) timing_rot20 = timing_parameters( geom, From 9fb048158df72e06aca81ba2671c4a63bb1fd9e4 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 5 Feb 2021 16:40:32 +0100 Subject: [PATCH 10/39] Precompute camera geometries --- ctapipe/image/image_processor.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index 803845933b1..9800e323fb8 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -84,7 +84,13 @@ def __init__( ) self.check_image = ImageQualityQuery(parent=self) self._is_simulation = is_simulation - self.telescope_frame = TelescopeFrame() + telescope_frame = TelescopeFrame() + self.telescope_frame_geometries = { + tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to( + telescope_frame + ) + for tel_id in self.subarray.tel + } def __call__(self, event: ArrayEventContainer): self._process_telescope_event(event) @@ -170,9 +176,8 @@ def _process_telescope_event(self, event): Loop over telescopes and process the calibrated images into parameters """ for tel_id, dl1_camera in event.dl1.tel.items(): - geometry = self.subarray.tel[tel_id].camera.geometry - geometry_telescope = geometry.transform_to(self.telescope_frame) + geometry_telescope = self.telescope_frame_geometries[tel_id] # compute image parameters only if requested to write them dl1_camera.image_mask = self.clean( tel_id=tel_id, From 8f9a81a8704327e454f85d17407c95d60c560c0e Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 5 Feb 2021 18:34:21 +0100 Subject: [PATCH 11/39] Fix unit test - Value changed slightly due to differences in the camera readout definitions (#1588) --- ctapipe/image/tests/test_hillas.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index 07e21c79db8..0be1770b6a0 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -291,7 +291,7 @@ def distance(coord): telescope_result = hillas_parameters(geom_nom, signal) assert u.isclose(np.abs(telescope_result.lon), 1 * u.deg, rtol=0.1) assert u.isclose(np.abs(telescope_result.lat), 1 * u.deg, rtol=0.1) - assert u.isclose(telescope_result.width, 0.06 * u.deg, rtol=0.1) + assert u.isclose(telescope_result.width, 0.065 * u.deg, rtol=0.1) assert u.isclose( telescope_result.width_uncertainty, 0.002 * u.deg, rtol=0.4 ) From 3bd826b778e6c7e3c9ef8f7cabf79b8a72ba7458 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Sat, 6 Feb 2021 19:01:15 +0100 Subject: [PATCH 12/39] Rename lon/lat to fov_lon/fov_lat --- ctapipe/containers.py | 10 +++++----- ctapipe/image/concentration.py | 4 ++-- ctapipe/image/hillas.py | 4 ++-- ctapipe/image/tests/test_hillas.py | 22 ++++++++++++++-------- ctapipe/image/timing.py | 4 ++-- ctapipe/io/tests/test_dl1eventsource.py | 2 +- ctapipe/io/tests/test_hdf5.py | 22 +++++++++++----------- ctapipe/visualization/tests/test_mpl.py | 2 +- 8 files changed, 38 insertions(+), 32 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 62fcadeda5f..f8bd83c2bf6 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -126,19 +126,19 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer): class HillasParametersContainer(BaseHillasParametersContainer): container_prefix = "hillas" - lon = Field( + fov_lon = Field( nan * u.deg, - "centroid longitude in a sky frame e.g. the telescope frame", + "longitude angle in a spherical system centered on the pointing position (TelescopeFrame)", unit=u.deg, ) - lat = Field( + fov_lat = Field( nan * u.deg, - "centroid latitude in a sky frame e.g. the telescope frame", + "latitude angle in a spherical system centered on the pointing position (TelescopeFrame)", unit=u.deg, ) r = Field( nan * u.deg, - "radial coordinate of centroid in a sky frame e.g. the telescope frame", + "radial coordinate in a spherical system centered on the pointing position (TelescopeFrame)", unit=u.deg, ) phi = Field( diff --git a/ctapipe/image/concentration.py b/ctapipe/image/concentration.py index 43258bd4ea9..d5f06bbbb3f 100644 --- a/ctapipe/image/concentration.py +++ b/ctapipe/image/concentration.py @@ -29,9 +29,9 @@ def concentration_parameters(geom, image, hillas_parameters): geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, unit=unit ) elif isinstance(h, HillasParametersContainer): - unit = h.lon.unit + unit = h.fov_lon.unit pix_x, pix_y, x, y, length, width = all_to_value( - geom.pix_x, geom.pix_y, h.lon, h.lat, h.length, h.width, unit=unit + geom.pix_x, geom.pix_y, h.fov_lon, h.fov_lat, h.length, h.width, unit=unit ) delta_x = pix_x - x diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 13bed79c0a9..962e265616f 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -211,8 +211,8 @@ def hillas_parameters(geom, image): ) else: return HillasParametersContainer( - lon=u.Quantity(cog_x, unit), - lat=u.Quantity(cog_y, unit), + fov_lon=u.Quantity(cog_x, unit), + fov_lat=u.Quantity(cog_y, unit), r=u.Quantity(cog_r, unit), phi=Angle(cog_phi, unit=u.rad), intensity=size, diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index 0be1770b6a0..6a03ab4c80d 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -289,8 +289,8 @@ def distance(coord): ) telescope_result = hillas_parameters(geom_nom, signal) - assert u.isclose(np.abs(telescope_result.lon), 1 * u.deg, rtol=0.1) - assert u.isclose(np.abs(telescope_result.lat), 1 * u.deg, rtol=0.1) + assert u.isclose(np.abs(telescope_result.fov_lon), 1 * u.deg, rtol=0.1) + assert u.isclose(np.abs(telescope_result.fov_lat), 1 * u.deg, rtol=0.1) assert u.isclose(telescope_result.width, 0.065 * u.deg, rtol=0.1) assert u.isclose( telescope_result.width_uncertainty, 0.002 * u.deg, rtol=0.4 @@ -305,16 +305,20 @@ def distance(coord): camera_result = hillas_parameters(geom, signal) transformed_cog = SkyCoord( - fov_lon=telescope_result.lon, - fov_lat=telescope_result.lat, + fov_lon=telescope_result.fov_lon, + fov_lat=telescope_result.fov_lat, frame=telescope_frame, ).transform_to(camera_frame) assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01) assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01) main_edges = u.Quantity([-telescope_result.length, telescope_result.length]) - main_lon = main_edges * np.cos(telescope_result.psi) + telescope_result.lon - main_lat = main_edges * np.sin(telescope_result.psi) + telescope_result.lat + main_lon = ( + main_edges * np.cos(telescope_result.psi) + telescope_result.fov_lon + ) + main_lat = ( + main_edges * np.sin(telescope_result.psi) + telescope_result.fov_lat + ) cam_main_axis = SkyCoord( fov_lon=main_lon, fov_lat=main_lat, frame=telescope_frame ).transform_to(camera_frame) @@ -325,10 +329,12 @@ def distance(coord): [-telescope_result.length, telescope_result.length] ) secondary_lon = ( - secondary_edges * np.cos(telescope_result.psi) + telescope_result.lon + secondary_edges * np.cos(telescope_result.psi) + + telescope_result.fov_lon ) secondary_lat = ( - secondary_edges * np.sin(telescope_result.psi) + telescope_result.lat + secondary_edges * np.sin(telescope_result.psi) + + telescope_result.fov_lat ) cam_secondary_edges = SkyCoord( fov_lon=secondary_lon, fov_lat=secondary_lat, frame=telescope_frame diff --git a/ctapipe/image/timing.py b/ctapipe/image/timing.py index c5f128d4749..a069f512567 100644 --- a/ctapipe/image/timing.py +++ b/ctapipe/image/timing.py @@ -69,9 +69,9 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, unit=unit ) elif isinstance(h, HillasParametersContainer): - unit = h.lon.unit + unit = h.fov_lon.unit pix_x, pix_y, x, y, length, width = all_to_value( - geom.pix_x, geom.pix_y, h.lon, h.lat, h.length, h.width, unit=unit + geom.pix_x, geom.pix_y, h.fov_lon, h.fov_lat, h.length, h.width, unit=unit ) longi, _ = camera_to_shower_coordinates( diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index 830e8af243d..f8108cdd8f0 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -79,7 +79,7 @@ def test_dl1_data(dl1_file): for event in source: for tel in event.dl1.tel: assert event.dl1.tel[tel].image.any() - assert event.dl1.tel[tel].parameters.hillas.lon != np.nan + assert event.dl1.tel[tel].parameters.hillas.fov_lon != np.nan def test_pointing(dl1_file): diff --git a/ctapipe/io/tests/test_hdf5.py b/ctapipe/io/tests/test_hdf5.py index 976db925888..98a88b79776 100644 --- a/ctapipe/io/tests/test_hdf5.py +++ b/ctapipe/io/tests/test_hdf5.py @@ -52,7 +52,7 @@ def test_write_container(temp_h5_file): def test_read_multiple_containers(): hillas_parameter_container = HillasParametersContainer( - lon=1 * u.deg, lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg + fov_lon=1 * u.deg, fov_lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) leakage_container = LeakageContainer( @@ -66,7 +66,7 @@ def test_read_multiple_containers(): writer.write("params", [hillas_parameter_container, leakage_container]) df = pd.read_hdf(f.name, key="/dl1/params") - assert "hillas_lon" in df.columns + assert "hillas_fov_lon" in df.columns assert "leakage_pixels_width_1" in df.columns # test reading both containers separately @@ -110,7 +110,7 @@ def test_read_multiple_containers(): def test_read_without_prefixes(): hillas_parameter_container = HillasParametersContainer( - lon=1 * u.deg, lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg + fov_lon=1 * u.deg, fov_lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) leakage_container = LeakageContainer( @@ -124,7 +124,7 @@ def test_read_without_prefixes(): writer.write("params", [hillas_parameter_container, leakage_container]) df = pd.read_hdf(f.name, key="/dl1/params") - assert "lon" in df.columns + assert "fov_lon" in df.columns assert "pixels_width_1" in df.columns # call with prefixes=False @@ -168,15 +168,15 @@ def test_read_without_prefixes(): def test_read_duplicated_container_types(): hillas_config_1 = HillasParametersContainer( - lon=1 * u.deg, - lat=2 * u.deg, + fov_lon=1 * u.deg, + fov_lat=2 * u.deg, length=3 * u.deg, width=4 * u.deg, prefix="hillas_1", ) hillas_config_2 = HillasParametersContainer( - lon=2 * u.deg, - lat=3 * u.deg, + fov_lon=2 * u.deg, + fov_lat=3 * u.deg, length=4 * u.deg, width=5 * u.deg, prefix="hillas_2", @@ -187,8 +187,8 @@ def test_read_duplicated_container_types(): writer.write("params", [hillas_config_1, hillas_config_2]) df = pd.read_hdf(f.name, key="/dl1/params") - assert "hillas_1_lon" in df.columns - assert "hillas_2_lat" in df.columns + assert "hillas_1_fov_lon" in df.columns + assert "hillas_2_fov_lat" in df.columns with HDF5TableReader(f.name) as reader: generator = reader.read( @@ -211,7 +211,7 @@ def test_read_duplicated_container_types(): def test_custom_prefix(): container = HillasParametersContainer( - lon=1 * u.deg, lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg + fov_lon=1 * u.deg, fov_lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) container.prefix = "custom" with tempfile.NamedTemporaryFile() as f: diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index fdba638ed02..6babf7496d3 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -49,7 +49,7 @@ def test_hillas_overlay(): from ctapipe.visualization import CameraDisplay disp = CameraDisplay(CameraGeometry.from_name("LSTCam")) - hillas = HillasParametersContainer( + hillas = CameraHillasParametersContainer( x=0.1 * u.m, y=-0.1 * u.m, length=0.5 * u.m, width=0.2 * u.m, psi=90 * u.deg ) From 8ded4e7a8144f57ff41d453d0d00674d3d3d5b60 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Sun, 7 Feb 2021 16:05:50 +0100 Subject: [PATCH 13/39] Fix use of ImageParametersContainer - Avoid base containers and thus lacking columns in the output file --- ctapipe/image/image_processor.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index 9800e323fb8..66c5013d022 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -10,6 +10,7 @@ ImageParametersContainer, TimingParametersContainer, PeakTimeStatisticsContainer, + HillasParametersContainer, ) from ..core import QualityQuery, TelescopeComponent from ..core.traits import List, create_class_enum_trait @@ -25,7 +26,10 @@ ) -DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer() +# avoid use of base containers for unparameterized images +DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer( + hillas=HillasParametersContainer(), timing=TimingParametersContainer() +) DEFAULT_TIMING_PARAMETERS = TimingParametersContainer() DEFAULT_PEAKTIME_STATISTICS = PeakTimeStatisticsContainer() From 2647f8aaa05b73358d2e1f561d07f2969920b74a Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 8 Feb 2021 11:49:50 +0100 Subject: [PATCH 14/39] Address codacy complaints --- ctapipe/containers.py | 42 +++++++++++++++++++++--------- ctapipe/image/hillas.py | 29 ++++++++++----------- ctapipe/image/tests/test_hillas.py | 18 ++++++------- ctapipe/image/timing.py | 7 +++-- 4 files changed, 55 insertions(+), 41 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index f8bd83c2bf6..4c4644b6e9a 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -105,12 +105,23 @@ class TelEventIndexContainer(Container): class BaseHillasParametersContainer(Container): + """ + Base container for hillas parameters to + allow the CameraHillasParametersContainer to + be assigned to an ImageParametersContainer as well. + """ + intensity = Field(nan, "total intensity (size)") skewness = Field(nan, "measure of the asymmetry") kurtosis = Field(nan, "measure of the tailedness") class CameraHillasParametersContainer(BaseHillasParametersContainer): + """ + Hillas Parameters in the camera frame. The cog position + is given in meter from the camera center. + """ + container_prefix = "camera_frame_hillas" x = Field(nan * u.m, "centroid x coordinate", unit=u.m) y = Field(nan * u.m, "centroid x coordinate", unit=u.m) @@ -125,27 +136,25 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer): class HillasParametersContainer(BaseHillasParametersContainer): + """ + Hillas Parameters in a spherical system centered on the pointing position + (TelescopeFrame). The cog position is given as offset in + longitude and latitude in degree. + """ + container_prefix = "hillas" fov_lon = Field( nan * u.deg, - "longitude angle in a spherical system centered on the pointing position (TelescopeFrame)", + "longitude angle in a spherical system centered on the pointing position", unit=u.deg, ) fov_lat = Field( nan * u.deg, - "latitude angle in a spherical system centered on the pointing position (TelescopeFrame)", - unit=u.deg, - ) - r = Field( - nan * u.deg, - "radial coordinate in a spherical system centered on the pointing position (TelescopeFrame)", - unit=u.deg, - ) - phi = Field( - nan * u.deg, - "polar coordinate of centroid in a sky frame e.g. the telescope frame", + "latitude angle in a spherical system centered on the pointing position", unit=u.deg, ) + r = Field(nan * u.deg, "radial coordinate of centroid", unit=u.deg) + phi = Field(nan * u.deg, "polar coordinate of centroid", unit=u.deg) length = Field(nan * u.deg, "standard deviation along the major-axis", unit=u.deg) length_uncertainty = Field(nan * u.deg, "uncertainty of length", unit=u.deg) @@ -195,6 +204,12 @@ class ConcentrationContainer(Container): class BaseTimingParametersContainer(Container): + """ + Base container for timing parameters to + allow the CameraTimingParametersContainer to + be assigned to an ImageParametersContainer as well. + """ + intercept = Field(nan, "intercept of arrival times along main shower axis") deviation = Field( nan, @@ -218,7 +233,8 @@ class CameraTimingParametersContainer(BaseTimingParametersContainer): class TimingParametersContainer(BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times - along the shower main axis in a sky frame e.g. the telescope frame. + along the shower main axis in a + spherical system centered on the pointing position (TelescopeFrame) """ container_prefix = "timing" diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index 962e265616f..24c3c57c629 100644 --- a/ctapipe/image/hillas.py +++ b/ctapipe/image/hillas.py @@ -209,18 +209,17 @@ def hillas_parameters(geom, image): skewness=skewness_long, kurtosis=kurtosis_long, ) - else: - return HillasParametersContainer( - fov_lon=u.Quantity(cog_x, unit), - fov_lat=u.Quantity(cog_y, unit), - r=u.Quantity(cog_r, unit), - phi=Angle(cog_phi, unit=u.rad), - intensity=size, - length=u.Quantity(length, unit), - length_uncertainty=u.Quantity(length_uncertainty, unit), - width=u.Quantity(width, unit), - width_uncertainty=u.Quantity(width_uncertainty, unit), - psi=Angle(psi, unit=u.rad), - skewness=skewness_long, - kurtosis=kurtosis_long, - ) + return HillasParametersContainer( + fov_lon=u.Quantity(cog_x, unit), + fov_lat=u.Quantity(cog_y, unit), + r=u.Quantity(cog_r, unit), + phi=Angle(cog_phi, unit=u.rad), + intensity=size, + length=u.Quantity(length, unit), + length_uncertainty=u.Quantity(length_uncertainty, unit), + width=u.Quantity(width, unit), + width_uncertainty=u.Quantity(width_uncertainty, unit), + psi=Angle(psi, unit=u.rad), + skewness=skewness_long, + kurtosis=kurtosis_long, + ) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index 6a03ab4c80d..93802be368a 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -1,12 +1,5 @@ -from ctapipe.instrument import CameraGeometry -from ctapipe.image import tailcuts_clean, toymodel -from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError -from ctapipe.containers import ( - HillasParametersContainer, - CameraHillasParametersContainer, -) -from astropy.coordinates import Angle, SkyCoord from ctapipe.coordinates import TelescopeFrame, CameraFrame +from astropy.coordinates import Angle, SkyCoord from astropy import units as u import numpy as np from numpy import isclose, zeros_like @@ -14,6 +7,13 @@ from pytest import approx import itertools import pytest +from ctapipe.instrument import CameraGeometry +from ctapipe.image import tailcuts_clean, toymodel +from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError +from ctapipe.containers import ( + HillasParametersContainer, + CameraHillasParametersContainer, +) def create_sample_image( @@ -279,7 +279,7 @@ def test_reconstruction_in_telescope_frame(): def distance(coord): return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2 - for i, (x, y) in enumerate(zip(xs, ys)): + for x, y in zip(xs, ys): for psi in psis: # make a toymodel shower model model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) diff --git a/ctapipe/image/timing.py b/ctapipe/image/timing.py index a069f512567..b7eb81c2f96 100644 --- a/ctapipe/image/timing.py +++ b/ctapipe/image/timing.py @@ -89,7 +89,6 @@ def timing_parameters(geom, image, peak_time, hillas_parameters, cleaning_mask=N return CameraTimingParametersContainer( slope=beta[0] / unit, intercept=beta[1], deviation=deviation ) - else: - return TimingParametersContainer( - slope=beta[0] / unit, intercept=beta[1], deviation=deviation - ) + return TimingParametersContainer( + slope=beta[0] / unit, intercept=beta[1], deviation=deviation + ) From dbb0f23d7db8fff371a00f9717d77cc257ff273e Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 8 Feb 2021 12:07:32 +0100 Subject: [PATCH 15/39] Fix default ImageParametersContainer - Sets the Hillas/Timing parameters container in the telescope frame as default values as opposed to the base containers --- ctapipe/containers.py | 12 ++++++++++-- ctapipe/image/image_processor.py | 5 +---- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 4c4644b6e9a..9145284a167 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -276,8 +276,16 @@ class ImageParametersContainer(Container): """ Collection of image parameters """ container_prefix = "params" - hillas = Field(BaseHillasParametersContainer(), "Hillas Parameters") - timing = Field(BaseTimingParametersContainer(), "Timing Parameters") + hillas = Field( + HillasParametersContainer(), + "Hillas Parameters", + type=BaseHillasParametersContainer, + ) + timing = Field( + TimingParametersContainer(), + "Timing Parameters", + type=BaseTimingParametersContainer, + ) leakage = Field(LeakageContainer(), "Leakage Parameters") concentration = Field(ConcentrationContainer(), "Concentration Parameters") morphology = Field(MorphologyContainer(), "Image Morphology Parameters") diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index 66c5013d022..60a8475589c 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -10,7 +10,6 @@ ImageParametersContainer, TimingParametersContainer, PeakTimeStatisticsContainer, - HillasParametersContainer, ) from ..core import QualityQuery, TelescopeComponent from ..core.traits import List, create_class_enum_trait @@ -27,9 +26,7 @@ # avoid use of base containers for unparameterized images -DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer( - hillas=HillasParametersContainer(), timing=TimingParametersContainer() -) +DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer() DEFAULT_TIMING_PARAMETERS = TimingParametersContainer() DEFAULT_PEAKTIME_STATISTICS = PeakTimeStatisticsContainer() From b413f85a03a9289d3a60f0a9c8012575033fbfcf Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 8 Feb 2021 12:18:43 +0100 Subject: [PATCH 16/39] Change datamodel version - Replace v1.0.4 with v1.1.0 since column names have changed on top of the units --- ctapipe/io/dl1eventsource.py | 8 ++++---- ctapipe/io/dl1writer.py | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/ctapipe/io/dl1eventsource.py b/ctapipe/io/dl1eventsource.py index c65e11b6e63..ad375819137 100644 --- a/ctapipe/io/dl1eventsource.py +++ b/ctapipe/io/dl1eventsource.py @@ -30,7 +30,7 @@ logger = logging.getLogger(__name__) -COMPATIBLE_DL1_VERSIONS = ["v1.0.0", "v1.0.1", "v1.0.2", "v1.0.3", "v1.0.4"] +COMPATIBLE_DL1_VERSIONS = ["v1.0.0", "v1.0.1", "v1.0.2", "v1.0.3", "v1.1.0"] class DL1EventSource(EventSource): @@ -232,12 +232,12 @@ def _generate_events(self): containers=[ ( HillasParametersContainer() - if (self.datamodel_version >= "1.0.4") + if (self.datamodel_version >= "1.1.0") else CameraHillasParametersContainer() ), ( TimingParametersContainer() - if (self.datamodel_version >= "1.0.4") + if (self.datamodel_version >= "1.1.0") else CameraTimingParametersContainer() ), LeakageContainer(), @@ -257,7 +257,7 @@ def _generate_events(self): containers=[ ( HillasParametersContainer() - if (self.datamodel_version >= "1.0.4") + if (self.datamodel_version >= "1.1.0") else CameraHillasParametersContainer() ), LeakageContainer(), diff --git a/ctapipe/io/dl1writer.py b/ctapipe/io/dl1writer.py index 5443d3fc00e..44a734df1a7 100644 --- a/ctapipe/io/dl1writer.py +++ b/ctapipe/io/dl1writer.py @@ -28,10 +28,10 @@ tables.parameters.NODE_CACHE_SLOTS = 3000 # fixes problem with too many datasets -DL1_DATA_MODEL_VERSION = "v1.0.4" +DL1_DATA_MODEL_VERSION = "v1.1.0" DL1_DATA_MODEL_CHANGE_HISTORY = """ - v1.0.3: true_image dtype changed from float32 to int32 -- v1.0.4: hillas and timing parameters saved in telescope frame (degree) as opposed to camera frame (m) +- v1.1.0: hillas and timing parameters saved in telescope frame (degree) as opposed to camera frame (m) """ PROV = Provenance() From 4dc22431d5abdda0271267350adb0de81b6bf3ca Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 9 Jul 2021 19:23:42 +0200 Subject: [PATCH 17/39] Use fov_lon/lat instead of x/y --- ctapipe/visualization/tests/test_mpl.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 8ae0b4d6dae..8e0de18578b 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -178,10 +178,10 @@ def test_array_display(): # test using hillas params for divergent pointing in telescopeframe: hillas_dict = { 1: HillasParametersContainer( - x=1.0 * u.deg, y=1.0 * u.deg, length=1.0 * u.deg, psi=90 * u.deg + fov_lon=1.0 * u.deg, fov_lat=1.0 * u.deg, length=1.0 * u.deg, psi=90 * u.deg ), 2: HillasParametersContainer( - x=1.0 * u.deg, y=1.0 * u.deg, length=1.0 * u.deg, psi=95 * u.deg + fov_lon=1.0 * u.deg, fov_lat=1.0 * u.deg, length=1.0 * u.deg, psi=95 * u.deg ), } ad.set_vector_hillas( @@ -196,10 +196,10 @@ def test_array_display(): # test using hillas params for parallel pointing in telescopeframe: hillas_dict = { 1: HillasParametersContainer( - x=1.0 * u.deg, y=1.0 * u.deg, length=1.0 * u.deg, psi=90 * u.deg + fov_lon=1.0 * u.deg, fov_lat=1.0 * u.deg, length=1.0 * u.deg, psi=90 * u.deg ), 2: HillasParametersContainer( - x=1.0 * u.deg, y=1.0 * u.deg, length=1.0 * u.deg, psi=95 * u.deg + fov_lon=1.0 * u.deg, fov_lat=1.0 * u.deg, length=1.0 * u.deg, psi=95 * u.deg ), } ad.set_vector_hillas( From 02754011d8b5d1801e27624019b85be03d2bf705 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 9 Jul 2021 19:25:30 +0200 Subject: [PATCH 18/39] Distinguish between CameraHillasParametersContainer and HillasParametersContainer --- ctapipe/reco/hillas_reconstructor.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py index e76a913f13b..52e49d1e768 100644 --- a/ctapipe/reco/hillas_reconstructor.py +++ b/ctapipe/reco/hillas_reconstructor.py @@ -8,7 +8,10 @@ InvalidWidthException, TooFewTelescopesException, ) -from ctapipe.containers import ReconstructedGeometryContainer +from ctapipe.containers import ( + ReconstructedGeometryContainer, + CameraHillasParametersContainer, +) from itertools import combinations from ctapipe.coordinates import ( @@ -349,7 +352,9 @@ def initialize_hillas_planes( focal_length = subarray.tel[tel_id].optics.equivalent_focal_length - if moments.x.unit.is_equivalent(u.m): # Image parameters are in CameraFrame + if isinstance( + moments, CameraHillasParametersContainer + ): # Image parameters are in CameraFrame # we just need any point on the main shower axis a bit away from the cog p2_x = moments.x + 0.1 * self._cam_radius_m[camera] * np.cos( @@ -369,17 +374,19 @@ def initialize_hillas_planes( else: # Image parameters are already in TelescopeFrame # we just need any point on the main shower axis a bit away from the cog - p2_delta_alt = moments.y + 0.1 * self._cam_radius_deg[camera] * np.sin( - moments.psi - ) - p2_delta_az = moments.x + 0.1 * self._cam_radius_deg[camera] * np.cos( - moments.psi - ) + p2_delta_alt = moments.fov_lat + 0.1 * self._cam_radius_deg[ + camera + ] * np.sin(moments.psi) + p2_delta_az = moments.fov_lon + 0.1 * self._cam_radius_deg[ + camera + ] * np.cos(moments.psi) telescope_frame = TelescopeFrame(telescope_pointing=pointing) cog_coord = SkyCoord( - fov_lon=moments.x, fov_lat=moments.y, frame=telescope_frame + fov_lon=moments.fov_lon, + fov_lat=moments.fov_lat, + frame=telescope_frame, ) p2_coord = SkyCoord( fov_lon=p2_delta_az, fov_lat=p2_delta_alt, frame=telescope_frame From 2b108b07be8111ff83542bafeec207e59fdb415f Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 22 Jul 2021 13:21:19 +0200 Subject: [PATCH 19/39] Adapt hillas intersection to work with hillas parameters in the telescope and camera frame --- ctapipe/reco/hillas_intersection.py | 78 ++++++++++++------- .../reco/tests/test_hillas_intersection.py | 70 +++++++++-------- 2 files changed, 86 insertions(+), 62 deletions(-) diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py index 362d9bb0597..d037258a7bd 100644 --- a/ctapipe/reco/hillas_intersection.py +++ b/ctapipe/reco/hillas_intersection.py @@ -16,13 +16,18 @@ InvalidWidthException, TooFewTelescopesException, ) -from ctapipe.containers import ReconstructedGeometryContainer +from ctapipe.containers import ( + ReconstructedGeometryContainer, + CameraHillasParametersContainer, + HillasParametersContainer, +) from ctapipe.instrument import get_atmosphere_profile_functions from astropy.coordinates import SkyCoord from ctapipe.coordinates import ( NominalFrame, CameraFrame, + TelescopeFrame, TiltedGroundFrame, project_to_ground, MissingFrameAttributeWarning, @@ -144,33 +149,49 @@ def predict(self, hillas_dict, subarray, array_pointing, telescopes_pointings=No nom_frame = NominalFrame(origin=array_pointing) - hillas_dict_mod = copy.deepcopy(hillas_dict) - - for tel_id, hillas in hillas_dict_mod.items(): - # prevent from using rads instead of meters as inputs - assert hillas.x.to(u.m).unit == u.Unit("m") - - focal_length = subarray.tel[tel_id].optics.equivalent_focal_length - - camera_frame = CameraFrame( - telescope_pointing=telescopes_pointings[tel_id], - focal_length=focal_length, + hillas_dict_mod = {} + + for tel_id, hillas in hillas_dict.items(): + if isinstance(hillas, CameraHillasParametersContainer): + focal_length = subarray.tel[tel_id].optics.equivalent_focal_length + camera_frame = CameraFrame( + telescope_pointing=telescopes_pointings[tel_id], + focal_length=focal_length, + ) + cog_coords = SkyCoord(x=hillas.x, y=hillas.y, frame=camera_frame) + cog_coords_nom = cog_coords.transform_to(nom_frame) + else: + telescope_frame = TelescopeFrame( + telescope_pointing=telescopes_pointings[tel_id] + ) + cog_coords = SkyCoord( + fov_lon=hillas.fov_lon, + fov_lat=hillas.fov_lat, + frame=telescope_frame, + ) + cog_coords_nom = cog_coords.transform_to(nom_frame) + hillas_dict_mod[tel_id] = HillasParametersContainer( + fov_lon=cog_coords_nom.fov_lon, + fov_lat=cog_coords_nom.fov_lat, + psi=hillas.psi, + width=hillas.width, + length=hillas.length, + intensity=hillas.intensity, ) - cog_coords = SkyCoord(x=hillas.x, y=hillas.y, frame=camera_frame) - cog_coords_nom = cog_coords.transform_to(nom_frame) - hillas.x = cog_coords_nom.fov_lat - hillas.y = cog_coords_nom.fov_lon - src_x, src_y, err_x, err_y = self.reconstruct_nominal(hillas_dict_mod) + src_fov_lon, src_fov_lat, err_fov_lon, err_fov_lat = self.reconstruct_nominal( + hillas_dict_mod + ) core_x, core_y, core_err_x, core_err_y = self.reconstruct_tilted( hillas_dict_mod, tel_x, tel_y ) - err_x *= u.rad - err_y *= u.rad + err_fov_lon *= u.rad + err_fov_lat *= u.rad - nom = SkyCoord(fov_lat=src_x * u.rad, fov_lon=src_y * u.rad, frame=nom_frame) - # nom = sky_pos.transform_to(nom_frame) + nom = SkyCoord( + fov_lon=src_fov_lon * u.rad, fov_lat=src_fov_lat * u.rad, frame=nom_frame + ) sky_pos = nom.transform_to(array_pointing.frame) tilt = SkyCoord(x=core_x * u.m, y=core_y * u.m, frame=tilted_frame) grd = project_to_ground(tilt) @@ -185,7 +206,7 @@ def predict(self, hillas_dict, subarray, array_pointing, telescopes_pointings=No 90 * u.deg - array_pointing.alt, ) - src_error = np.sqrt(err_x ** 2 + err_y ** 2) + src_error = np.sqrt(err_fov_lon ** 2 + err_fov_lat ** 2) result = ReconstructedGeometryContainer( alt=sky_pos.altaz.alt.to(u.rad), @@ -202,7 +223,6 @@ def predict(self, hillas_dict, subarray, array_pointing, telescopes_pointings=No h_max_uncert=u.Quantity(np.nan * x_max.unit), goodness_of_fit=np.nan, ) - return result def reconstruct_nominal(self, hillas_parameters): @@ -232,8 +252,8 @@ def reconstruct_nominal(self, hillas_parameters): map( lambda h: [ h[0].psi.to_value(u.rad), - h[0].x.to_value(u.rad), - h[0].y.to_value(u.rad), + h[0].fov_lon.to_value(u.rad), + h[0].fov_lat.to_value(u.rad), h[0].intensity, ], hillas_pairs, @@ -246,8 +266,8 @@ def reconstruct_nominal(self, hillas_parameters): map( lambda h: [ h[1].psi.to_value(u.rad), - h[1].x.to_value(u.rad), - h[1].y.to_value(u.rad), + h[1].fov_lon.to_value(u.rad), + h[1].fov_lat.to_value(u.rad), h[1].intensity, ], hillas_pairs, @@ -389,8 +409,8 @@ def reconstruct_xmax( # Loops over telescopes in event for tel in hillas_parameters.keys(): - cog_x.append(hillas_parameters[tel].x.to_value(u.rad)) - cog_y.append(hillas_parameters[tel].y.to_value(u.rad)) + cog_x.append(hillas_parameters[tel].fov_lon.to_value(u.rad)) + cog_y.append(hillas_parameters[tel].fov_lat.to_value(u.rad)) amp.append(hillas_parameters[tel].intensity) tx.append(tel_x[tel].to_value(u.m)) diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py index 7398862c7ac..66f450c80a4 100644 --- a/ctapipe/reco/tests/test_hillas_intersection.py +++ b/ctapipe/reco/tests/test_hillas_intersection.py @@ -3,8 +3,11 @@ from numpy.testing import assert_allclose import numpy as np from astropy.coordinates import SkyCoord -from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame -from ctapipe.containers import CameraHillasParametersContainer +from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame, TelescopeFrame +from ctapipe.containers import ( + CameraHillasParametersContainer, + HillasParametersContainer, +) from ctapipe.io import EventSource @@ -81,14 +84,14 @@ def test_intersection_xmax_reco(): focal_length = 28 * u.m hillas_dict = { - 1: CameraHillasParametersContainer( - x=-(delta / focal_length) * u.rad, - y=((0 * u.m) / focal_length) * u.rad, + 1: HillasParametersContainer( + fov_lon=-(delta / focal_length) * u.rad, + fov_lat=((0 * u.m) / focal_length) * u.rad, intensity=1, ), - 2: CameraHillasParametersContainer( - x=((0 * u.m) / focal_length) * u.rad, - y=-(delta / focal_length) * u.rad, + 2: HillasParametersContainer( + fov_lon=((0 * u.m) / focal_length) * u.rad, + fov_lat=-(delta / focal_length) * u.rad, intensity=1, ), } @@ -119,9 +122,9 @@ def test_intersection_reco_impact_point_tilted(): tel_y_dict = {1: delta, 2: delta, 3: -delta} hillas_dict = { - 1: CameraHillasParametersContainer(intensity=100, psi=-90 * u.deg), - 2: CameraHillasParametersContainer(intensity=100, psi=-45 * u.deg), - 3: CameraHillasParametersContainer(intensity=100, psi=0 * u.deg), + 1: HillasParametersContainer(intensity=100, psi=-90 * u.deg), + 2: HillasParametersContainer(intensity=100, psi=-45 * u.deg), + 3: HillasParametersContainer(intensity=100, psi=0 * u.deg), } reco_konrad = hill_inter.reconstruct_tilted( @@ -144,9 +147,9 @@ def test_intersection_weighting_spoiled_parameters(): # telescope 2 have a spoiled reconstruction (45 instead of -45) hillas_dict = { - 1: CameraHillasParametersContainer(intensity=10000, psi=-90 * u.deg), - 2: CameraHillasParametersContainer(intensity=1, psi=45 * u.deg), - 3: CameraHillasParametersContainer(intensity=10000, psi=0 * u.deg), + 1: HillasParametersContainer(intensity=10000, psi=-90 * u.deg), + 2: HillasParametersContainer(intensity=1, psi=45 * u.deg), + 3: HillasParametersContainer(intensity=10000, psi=0 * u.deg), } reco_konrad_spoiled = hill_inter.reconstruct_tilted( @@ -180,32 +183,31 @@ def test_intersection_nominal_reconstruction(): focal_length=focal_length, telescope_pointing=array_direction ) - cog_coords_camera_1 = SkyCoord(x=delta, y=0 * u.m, frame=camera_frame) - cog_coords_camera_2 = SkyCoord(x=delta / 0.7, y=delta / 0.7, frame=camera_frame) - cog_coords_camera_3 = SkyCoord(x=0 * u.m, y=delta, frame=camera_frame) + cog_coords_camera_1 = SkyCoord(y=delta, x=0 * u.m, frame=camera_frame) + cog_coords_camera_2 = SkyCoord(y=delta / 0.7, x=delta / 0.7, frame=camera_frame) + cog_coords_camera_3 = SkyCoord(y=0 * u.m, x=delta, frame=camera_frame) cog_coords_nom_1 = cog_coords_camera_1.transform_to(nominal_frame) cog_coords_nom_2 = cog_coords_camera_2.transform_to(nominal_frame) cog_coords_nom_3 = cog_coords_camera_3.transform_to(nominal_frame) - # x-axis is along the altitude and y-axis is along the azimuth - hillas_1 = CameraHillasParametersContainer( - x=cog_coords_nom_1.fov_lat, - y=cog_coords_nom_1.fov_lon, + hillas_1 = HillasParametersContainer( + fov_lat=cog_coords_nom_1.fov_lat, + fov_lon=cog_coords_nom_1.fov_lon, intensity=100, psi=0 * u.deg, ) - hillas_2 = CameraHillasParametersContainer( - x=cog_coords_nom_2.fov_lat, - y=cog_coords_nom_2.fov_lon, + hillas_2 = HillasParametersContainer( + fov_lat=cog_coords_nom_2.fov_lat, + fov_lon=cog_coords_nom_2.fov_lon, intensity=100, psi=45 * u.deg, ) - hillas_3 = CameraHillasParametersContainer( - x=cog_coords_nom_3.fov_lat, - y=cog_coords_nom_3.fov_lon, + hillas_3 = HillasParametersContainer( + fov_lat=cog_coords_nom_3.fov_lat, + fov_lon=cog_coords_nom_3.fov_lon, intensity=100, psi=90 * u.deg, ) @@ -257,10 +259,11 @@ def test_reconstruction(): hillas_dict = {} telescope_pointings = {} - + telescope_frame = TelescopeFrame() for tel_id, dl1 in event.dl1.tel.items(): geom = source.subarray.tel[tel_id].camera.geometry + geom_tel_frame = geom.transform_to(telescope_frame) telescope_pointings[tel_id] = SkyCoord( alt=event.pointing.tel[tel_id].altitude, @@ -273,7 +276,7 @@ def test_reconstruction(): ) try: - moments = hillas_parameters(geom[mask], dl1.image[mask]) + moments = hillas_parameters(geom_tel_frame[mask], dl1.image[mask]) hillas_dict[tel_id] = moments except HillasParameterizationError as e: print(e) @@ -320,15 +323,16 @@ def test_reconstruction_works(subarray_and_event_gamma_off_axis_500_gev): for tel_id, pointing in event.pointing.tel.items() if tel_id in hillas_dict } + true_coord = SkyCoord( + alt=event.simulation.shower.alt, az=event.simulation.shower.az, frame=AltAz() + ) + nom_frame = NominalFrame(origin=array_pointing) + true_nom = true_coord.transform_to(nom_frame) result = reconstructor.predict( hillas_dict, subarray, array_pointing, telescope_pointings ) - reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz()) - true_coord = SkyCoord( - alt=event.simulation.shower.alt, az=event.simulation.shower.az, frame=AltAz() - ) assert reco_coord.separation(true_coord) < 0.1 * u.deg From c03e5c07e33c7f3db4b4459d6a0558841ef4f402 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 22 Jul 2021 13:48:54 +0200 Subject: [PATCH 20/39] Remove duplicated tests of the results - Comparing values between both reconstructions should handle all cases - Comparing results in the telescope frame to fixed values was error prone, because the image is generated in the camera frame requiring explicit conversion every time --- ctapipe/image/tests/test_hillas.py | 72 +++++++++++++----------------- 1 file changed, 32 insertions(+), 40 deletions(-) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index 7c800b9c488..fc938c258d8 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -280,31 +280,43 @@ def test_reconstruction_in_telescope_frame(): def distance(coord): return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2 + def get_transformed_length(telescope_hillas, telescope_frame, camera_frame): + main_edges = u.Quantity([-telescope_hillas.length, telescope_hillas.length]) + main_lon = main_edges * np.cos(telescope_hillas.psi) + telescope_hillas.fov_lon + main_lat = main_edges * np.sin(telescope_hillas.psi) + telescope_hillas.fov_lat + cam_main_axis = SkyCoord( + fov_lon=main_lon, fov_lat=main_lat, frame=telescope_frame + ).transform_to(camera_frame) + transformed_length = distance(cam_main_axis) + return transformed_length + + def get_transformed_width(telescope_hillas, telescope_frame, camera_frame): + secondary_edges = u.Quantity([-telescope_hillas.width, telescope_hillas.width]) + secondary_lon = ( + secondary_edges * np.cos(telescope_hillas.psi) + telescope_result.fov_lon + ) + secondary_lat = ( + secondary_edges * np.sin(telescope_hillas.psi) + telescope_result.fov_lat + ) + cam_secondary_edges = SkyCoord( + fov_lon=secondary_lon, fov_lat=secondary_lat, frame=telescope_frame + ).transform_to(camera_frame) + transformed_width = distance(cam_secondary_edges) + return transformed_width + for x, y in zip(xs, ys): for psi in psis: - # make a toymodel shower model + # generate a toy image model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi) - image, signal, noise = model.generate_image( geom, intensity=intensity, nsb_level_pe=5 ) telescope_result = hillas_parameters(geom_nom, signal) - assert u.isclose(np.abs(telescope_result.fov_lon), 1 * u.deg, rtol=0.1) - assert u.isclose(np.abs(telescope_result.fov_lat), 1 * u.deg, rtol=0.1) - assert u.isclose(telescope_result.width, 0.065 * u.deg, rtol=0.1) - assert u.isclose( - telescope_result.width_uncertainty, 0.002 * u.deg, rtol=0.4 - ) - assert u.isclose(telescope_result.length, 0.3 * u.deg, rtol=0.1) - assert u.isclose( - telescope_result.length_uncertainty, 0.01 * u.deg, rtol=0.4 - ) - assert signal.sum() == telescope_result.intensity - - # Compare results with calculation in the camera frame camera_result = hillas_parameters(geom, signal) + assert camera_result.intensity == telescope_result.intensity + # Compare results in both frames transformed_cog = SkyCoord( fov_lon=telescope_result.fov_lon, fov_lat=telescope_result.fov_lat, @@ -313,32 +325,12 @@ def distance(coord): assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01) assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01) - main_edges = u.Quantity([-telescope_result.length, telescope_result.length]) - main_lon = ( - main_edges * np.cos(telescope_result.psi) + telescope_result.fov_lon - ) - main_lat = ( - main_edges * np.sin(telescope_result.psi) + telescope_result.fov_lat + transformed_length = get_transformed_length( + telescope_result, telescope_frame, camera_frame ) - cam_main_axis = SkyCoord( - fov_lon=main_lon, fov_lat=main_lat, frame=telescope_frame - ).transform_to(camera_frame) - transformed_length = distance(cam_main_axis) assert u.isclose(transformed_length, camera_result.length, rtol=0.01) - secondary_edges = u.Quantity( - [-telescope_result.length, telescope_result.length] + transformed_width = get_transformed_width( + telescope_result, telescope_frame, camera_frame ) - secondary_lon = ( - secondary_edges * np.cos(telescope_result.psi) - + telescope_result.fov_lon - ) - secondary_lat = ( - secondary_edges * np.sin(telescope_result.psi) - + telescope_result.fov_lat - ) - cam_secondary_edges = SkyCoord( - fov_lon=secondary_lon, fov_lat=secondary_lat, frame=telescope_frame - ).transform_to(camera_frame) - transformed_width = distance(cam_secondary_edges) - assert u.isclose(transformed_width, camera_result.length, rtol=0.01) + assert u.isclose(transformed_width, camera_result.width, rtol=0.01) From d6641e2b0bf9fcbc50bd0ff136cf804b8921fb55 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Thu, 22 Jul 2021 19:22:27 +0200 Subject: [PATCH 21/39] Remove unused declarations and imports --- ctapipe/reco/hillas_intersection.py | 1 - ctapipe/reco/tests/test_hillas_intersection.py | 7 +------ 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py index d037258a7bd..89f781f777f 100644 --- a/ctapipe/reco/hillas_intersection.py +++ b/ctapipe/reco/hillas_intersection.py @@ -32,7 +32,6 @@ project_to_ground, MissingFrameAttributeWarning, ) -import copy import warnings from ctapipe.core import traits diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py index 66f450c80a4..e05af2d370c 100644 --- a/ctapipe/reco/tests/test_hillas_intersection.py +++ b/ctapipe/reco/tests/test_hillas_intersection.py @@ -4,10 +4,7 @@ import numpy as np from astropy.coordinates import SkyCoord from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame, TelescopeFrame -from ctapipe.containers import ( - CameraHillasParametersContainer, - HillasParametersContainer, -) +from ctapipe.containers import HillasParametersContainer from ctapipe.io import EventSource @@ -327,8 +324,6 @@ def test_reconstruction_works(subarray_and_event_gamma_off_axis_500_gev): alt=event.simulation.shower.alt, az=event.simulation.shower.az, frame=AltAz() ) - nom_frame = NominalFrame(origin=array_pointing) - true_nom = true_coord.transform_to(nom_frame) result = reconstructor.predict( hillas_dict, subarray, array_pointing, telescope_pointings ) From bfa7505c94b50ddb69044c44ec87417a6164292c Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 10:36:56 +0200 Subject: [PATCH 22/39] Fix datamodel version --- ctapipe/io/datawriter.py | 2 +- ctapipe/io/dl1eventsource.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py index f45061e0f75..39931ff39c8 100644 --- a/ctapipe/io/datawriter.py +++ b/ctapipe/io/datawriter.py @@ -32,10 +32,10 @@ DATA_MODEL_VERSION = "v1.2.0" DATA_MODEL_CHANGE_HISTORY = """ +- v1.3.0: hillas and timing parameters saved in telescope frame (degree) as opposed to camera frame (m) - v1.2.0: change to more general data model, including also DL2 (DL1 unchanged) - v1.1.0: images and peak_times can be stored as scaled integers - v1.0.3: true_image dtype changed from float32 to int32 -- v1.1.0: hillas and timing parameters saved in telescope frame (degree) as opposed to camera frame (m) """ PROV = Provenance() diff --git a/ctapipe/io/dl1eventsource.py b/ctapipe/io/dl1eventsource.py index f1f55297a13..b9209ac0303 100644 --- a/ctapipe/io/dl1eventsource.py +++ b/ctapipe/io/dl1eventsource.py @@ -261,12 +261,12 @@ def _generate_events(self): containers=[ ( HillasParametersContainer() - if (self.datamodel_version >= "1.1.0") + if (self.datamodel_version >= "1.3.0") else CameraHillasParametersContainer() ), ( TimingParametersContainer() - if (self.datamodel_version >= "1.1.0") + if (self.datamodel_version >= "1.3.0") else CameraTimingParametersContainer() ), LeakageContainer(), @@ -286,7 +286,7 @@ def _generate_events(self): containers=[ ( HillasParametersContainer() - if (self.datamodel_version >= "1.1.0") + if (self.datamodel_version >= "1.3.0") else CameraHillasParametersContainer() ), LeakageContainer(), From ff648fb4fd36e6ba0384c5d4daa50f1c3213dd63 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 10:51:36 +0200 Subject: [PATCH 23/39] Fix unit test - The test file gets generated by calling stage1 on a simtel file, so this has to take the current datamodel into account --- ctapipe/io/tests/test_dl1eventsource.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index f84e5b64bb2..d9bb0c140c3 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -68,7 +68,9 @@ def test_simulation_info(dl1_file): for tel in event.simulation.tel: assert tel in event.simulation.tel assert event.simulation.tel[tel].true_image is not None - assert event.simulation.tel[tel].true_parameters.hillas.x != np.nan + assert ( + event.simulation.tel[tel].true_parameters.hillas.fov_lon != np.nan + ) def test_dl1_a_only_data(dl1_image_file): @@ -82,7 +84,7 @@ def test_dl1_b_only_data(dl1_parameters_file): with DL1EventSource(input_url=dl1_parameters_file) as source: for event in source: for tel in event.dl1.tel: - assert event.dl1.tel[tel].parameters.hillas.x != np.nan + assert event.dl1.tel[tel].parameters.hillas.fov_lon != np.nan def test_dl1_data(dl1_file): From 46ae8922da320fd52a5898d33af84e5d7cccba4b Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 11:04:49 +0200 Subject: [PATCH 24/39] Fix comparison to nan in unit tests --- ctapipe/io/tests/test_dl1eventsource.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index d9bb0c140c3..fa72f4e8a30 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -68,8 +68,8 @@ def test_simulation_info(dl1_file): for tel in event.simulation.tel: assert tel in event.simulation.tel assert event.simulation.tel[tel].true_image is not None - assert ( - event.simulation.tel[tel].true_parameters.hillas.fov_lon != np.nan + assert not np.isnan( + event.simulation.tel[tel].true_parameters.hillas.fov_lon ) @@ -84,7 +84,7 @@ def test_dl1_b_only_data(dl1_parameters_file): with DL1EventSource(input_url=dl1_parameters_file) as source: for event in source: for tel in event.dl1.tel: - assert event.dl1.tel[tel].parameters.hillas.fov_lon != np.nan + assert not np.isnan(event.dl1.tel[tel].parameters.hillas.fov_lon) def test_dl1_data(dl1_file): @@ -92,7 +92,7 @@ def test_dl1_data(dl1_file): for event in source: for tel in event.dl1.tel: assert event.dl1.tel[tel].image.any() - assert event.dl1.tel[tel].parameters.hillas.fov_lon != np.nan + assert not np.isnan(event.dl1.tel[tel].parameters.hillas.fov_lon) def test_pointing(dl1_file): From 98b77949c6a91f4450ab3524d0b70bc4e07beb2b Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 12:09:29 +0200 Subject: [PATCH 25/39] Fix reading of dl1 files - Add v1.3 to the list of supported datamodel versions - The renaming of the parameter containers changed the default container prefixes (HillasParametersContainer -> CameraHillasParametersContainer leads to the camera_hillas prefix), breaking backwards compatibility - This is fixed by specifying explicit prefixes when reading = "1.3.0") - else CameraHillasParametersContainer() + if (self.datamodel_version >= "v1.3.0") + else CameraHillasParametersContainer(prefix="hillas") ), ( TimingParametersContainer() - if (self.datamodel_version >= "1.3.0") - else CameraTimingParametersContainer() + if (self.datamodel_version >= "v1.3.0") + else CameraTimingParametersContainer(prefix="timing") ), LeakageContainer(), ConcentrationContainer(), @@ -286,8 +294,8 @@ def _generate_events(self): containers=[ ( HillasParametersContainer() - if (self.datamodel_version >= "1.3.0") - else CameraHillasParametersContainer() + if (self.datamodel_version >= "v1.3.0") + else CameraHillasParametersContainer(prefix="hillas") ), LeakageContainer(), ConcentrationContainer(), From 6a0f5d49985dfc06aeff59e4f9b18f886e601973 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 12:13:48 +0200 Subject: [PATCH 26/39] Fix tests for the dl1eventsource - Improper comparison to np.nan lead to tests effectively being skipped - The assertion would have failed actually, because the test file contains nan events (which is expected) - The new test checks if all entries are nan instead. Due to the event loop in the event source, this requires collecting the values in a list first --- ctapipe/io/tests/test_dl1eventsource.py | 48 +++++++++++++++++-------- 1 file changed, 34 insertions(+), 14 deletions(-) diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index fa72f4e8a30..8a71b33314a 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -4,15 +4,6 @@ from ctapipe.io import EventSource import astropy.units as u import numpy as np -import tempfile -import pytest - -d = tempfile.TemporaryDirectory() - - -@pytest.fixture(scope="module") -def dl1_dir(tmp_path_factory): - return tmp_path_factory.mktemp("dl1") def test_is_compatible(dl1_file): @@ -60,17 +51,28 @@ def test_allowed_tels(dl1_file): def test_simulation_info(dl1_file): + """ + Test that the simulated event information is plausible. + In particular this means simulated event information is finite + for all events and parameters calculated on the true images + are not all nan with the same number of nans in different columns. + """ + reco_lons = [] + reco_concentrations = [] with DL1EventSource(input_url=dl1_file) as source: for event in source: assert np.isfinite(event.simulation.shower.energy) - # the currently used file does not include true dl1 information - # this is skipped for that reason for tel in event.simulation.tel: assert tel in event.simulation.tel assert event.simulation.tel[tel].true_image is not None - assert not np.isnan( + reco_lons.append( event.simulation.tel[tel].true_parameters.hillas.fov_lon ) + reco_concentrations.append( + event.simulation.tel[tel].true_parameters.concentration.core + ) + assert not np.isnan(reco_lons).all() + assert sum(np.isnan(reco_lons)) == sum(np.isnan(reco_concentrations)) def test_dl1_a_only_data(dl1_image_file): @@ -81,18 +83,36 @@ def test_dl1_a_only_data(dl1_image_file): def test_dl1_b_only_data(dl1_parameters_file): + reco_lons = [] + reco_concentrations = [] with DL1EventSource(input_url=dl1_parameters_file) as source: for event in source: for tel in event.dl1.tel: - assert not np.isnan(event.dl1.tel[tel].parameters.hillas.fov_lon) + reco_lons.append( + event.simulation.tel[tel].true_parameters.hillas.fov_lon + ) + reco_concentrations.append( + event.simulation.tel[tel].true_parameters.concentration.core + ) + assert not np.isnan(reco_lons).all() + assert sum(np.isnan(reco_lons)) == sum(np.isnan(reco_concentrations)) def test_dl1_data(dl1_file): + reco_lons = [] + reco_concentrations = [] with DL1EventSource(input_url=dl1_file) as source: for event in source: for tel in event.dl1.tel: assert event.dl1.tel[tel].image.any() - assert not np.isnan(event.dl1.tel[tel].parameters.hillas.fov_lon) + reco_lons.append( + event.simulation.tel[tel].true_parameters.hillas.fov_lon + ) + reco_concentrations.append( + event.simulation.tel[tel].true_parameters.concentration.core + ) + assert not np.isnan(reco_lons).all() + assert sum(np.isnan(reco_lons)) == sum(np.isnan(reco_concentrations)) def test_pointing(dl1_file): From b281696e44121df6bc3a9f7a36e535cf8a225139 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 12:47:24 +0200 Subject: [PATCH 27/39] Write out correct datamodel version --- ctapipe/io/datawriter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py index 39931ff39c8..4919301a50c 100644 --- a/ctapipe/io/datawriter.py +++ b/ctapipe/io/datawriter.py @@ -30,7 +30,7 @@ tables.parameters.NODE_CACHE_SLOTS = 3000 # fixes problem with too many datasets -DATA_MODEL_VERSION = "v1.2.0" +DATA_MODEL_VERSION = "v1.3.0" DATA_MODEL_CHANGE_HISTORY = """ - v1.3.0: hillas and timing parameters saved in telescope frame (degree) as opposed to camera frame (m) - v1.2.0: change to more general data model, including also DL2 (DL1 unchanged) From ffb8de13f935f56d7ad3583d6de0cad20c68da0b Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 13:09:55 +0200 Subject: [PATCH 28/39] Fix unit test --- ctapipe/io/tests/test_dl1eventsource.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index 8a71b33314a..364aaffe3b7 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -66,7 +66,7 @@ def test_simulation_info(dl1_file): assert tel in event.simulation.tel assert event.simulation.tel[tel].true_image is not None reco_lons.append( - event.simulation.tel[tel].true_parameters.hillas.fov_lon + event.simulation.tel[tel].true_parameters.hillas.fov_lon.value ) reco_concentrations.append( event.simulation.tel[tel].true_parameters.concentration.core @@ -89,7 +89,7 @@ def test_dl1_b_only_data(dl1_parameters_file): for event in source: for tel in event.dl1.tel: reco_lons.append( - event.simulation.tel[tel].true_parameters.hillas.fov_lon + event.simulation.tel[tel].true_parameters.hillas.fov_lon.value ) reco_concentrations.append( event.simulation.tel[tel].true_parameters.concentration.core @@ -106,7 +106,7 @@ def test_dl1_data(dl1_file): for tel in event.dl1.tel: assert event.dl1.tel[tel].image.any() reco_lons.append( - event.simulation.tel[tel].true_parameters.hillas.fov_lon + event.simulation.tel[tel].true_parameters.hillas.fov_lon.value ) reco_concentrations.append( event.simulation.tel[tel].true_parameters.concentration.core From 54f2e066d6cbf700457ffd84c1db248d40bbd590 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 15:14:14 +0200 Subject: [PATCH 29/39] Try to fix the docs by making the base containers private --- ctapipe/containers.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index d75c1e5b0d2..257867bb520 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -109,7 +109,7 @@ class TelEventIndexContainer(Container): tel_id = Field(0, "telescope identifier") -class BaseHillasParametersContainer(Container): +class _BaseHillasParametersContainer(Container): """ Base container for hillas parameters to allow the CameraHillasParametersContainer to @@ -121,7 +121,7 @@ class BaseHillasParametersContainer(Container): kurtosis = Field(nan, "measure of the tailedness") -class CameraHillasParametersContainer(BaseHillasParametersContainer): +class CameraHillasParametersContainer(_BaseHillasParametersContainer): """ Hillas Parameters in the camera frame. The cog position is given in meter from the camera center. @@ -140,7 +140,7 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer): psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) -class HillasParametersContainer(BaseHillasParametersContainer): +class HillasParametersContainer(_BaseHillasParametersContainer): """ Hillas Parameters in a spherical system centered on the pointing position (TelescopeFrame). The cog position is given as offset in @@ -208,7 +208,7 @@ class ConcentrationContainer(Container): pixel = Field(nan, "Percentage of photo-electrons in the brightest pixel") -class BaseTimingParametersContainer(Container): +class _BaseTimingParametersContainer(Container): """ Base container for timing parameters to allow the CameraTimingParametersContainer to @@ -223,7 +223,7 @@ class BaseTimingParametersContainer(Container): ) -class CameraTimingParametersContainer(BaseTimingParametersContainer): +class CameraTimingParametersContainer(_BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times along the shower main axis in the camera frame. @@ -235,7 +235,7 @@ class CameraTimingParametersContainer(BaseTimingParametersContainer): ) -class TimingParametersContainer(BaseTimingParametersContainer): +class TimingParametersContainer(_BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times along the shower main axis in a @@ -291,12 +291,12 @@ class ImageParametersContainer(Container): hillas = Field( HillasParametersContainer(), "Hillas Parameters", - type=BaseHillasParametersContainer, + type=_BaseHillasParametersContainer, ) timing = Field( TimingParametersContainer(), "Timing Parameters", - type=BaseTimingParametersContainer, + type=_BaseTimingParametersContainer, ) leakage = Field(LeakageContainer(), "Leakage Parameters") concentration = Field(ConcentrationContainer(), "Concentration Parameters") From 19a30cc3616831a3272572ee117e82cc52e56dd5 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 15:14:42 +0200 Subject: [PATCH 30/39] Fix codacy complaints --- ctapipe/image/image_processor.py | 2 +- ctapipe/image/tests/test_hillas.py | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index d2156e11928..85e0f323b74 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -227,4 +227,4 @@ def _process_telescope_event(self, event): event.simulation.tel[tel_id].true_parameters.as_dict( recursive=True ), - ) \ No newline at end of file + ) diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index fc938c258d8..16fff8c0c27 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -1,4 +1,3 @@ -from ctapipe.coordinates import TelescopeFrame, CameraFrame from astropy.coordinates import Angle, SkyCoord from astropy import units as u import numpy as np @@ -6,6 +5,7 @@ from pytest import approx import itertools import pytest +from ctapipe.coordinates import TelescopeFrame, CameraFrame from ctapipe.instrument import CameraGeometry from ctapipe.image import tailcuts_clean, toymodel from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError @@ -260,6 +260,10 @@ def test_single_pixel(): def test_reconstruction_in_telescope_frame(): + """ + Compare the reconstruction in the telescope + and camera frame. + """ np.random.seed(42) telescope_frame = TelescopeFrame() From 52a46b6b31e8eca2fbe0a2b8b4cd93113387ef1b Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Fri, 23 Jul 2021 15:27:39 +0200 Subject: [PATCH 31/39] Add base containers to __all__ - This seems to be necessary for the sphinx docs --- ctapipe/containers.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index 257867bb520..baaa64d5b6e 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -31,6 +31,7 @@ "MonitoringCameraContainer", "MonitoringContainer", "MorphologyContainer", + "BaseHillasParametersContainer", "CameraHillasParametersContainer", "CameraTimingParametersContainer", "ParticleClassificationContainer", @@ -48,6 +49,7 @@ "SimulatedShowerDistribution", "SimulationConfigContainer", "TelEventIndexContainer", + "BaseTimingParametersContainer", "TimingParametersContainer", "TriggerContainer", "WaveformCalibrationContainer", @@ -109,7 +111,7 @@ class TelEventIndexContainer(Container): tel_id = Field(0, "telescope identifier") -class _BaseHillasParametersContainer(Container): +class BaseHillasParametersContainer(Container): """ Base container for hillas parameters to allow the CameraHillasParametersContainer to @@ -121,7 +123,7 @@ class _BaseHillasParametersContainer(Container): kurtosis = Field(nan, "measure of the tailedness") -class CameraHillasParametersContainer(_BaseHillasParametersContainer): +class CameraHillasParametersContainer(BaseHillasParametersContainer): """ Hillas Parameters in the camera frame. The cog position is given in meter from the camera center. @@ -140,7 +142,7 @@ class CameraHillasParametersContainer(_BaseHillasParametersContainer): psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) -class HillasParametersContainer(_BaseHillasParametersContainer): +class HillasParametersContainer(BaseHillasParametersContainer): """ Hillas Parameters in a spherical system centered on the pointing position (TelescopeFrame). The cog position is given as offset in @@ -208,7 +210,7 @@ class ConcentrationContainer(Container): pixel = Field(nan, "Percentage of photo-electrons in the brightest pixel") -class _BaseTimingParametersContainer(Container): +class BaseTimingParametersContainer(Container): """ Base container for timing parameters to allow the CameraTimingParametersContainer to @@ -223,7 +225,7 @@ class _BaseTimingParametersContainer(Container): ) -class CameraTimingParametersContainer(_BaseTimingParametersContainer): +class CameraTimingParametersContainer(BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times along the shower main axis in the camera frame. @@ -235,7 +237,7 @@ class CameraTimingParametersContainer(_BaseTimingParametersContainer): ) -class TimingParametersContainer(_BaseTimingParametersContainer): +class TimingParametersContainer(BaseTimingParametersContainer): """ Slope and Intercept of a linear regression of the arrival times along the shower main axis in a @@ -291,12 +293,12 @@ class ImageParametersContainer(Container): hillas = Field( HillasParametersContainer(), "Hillas Parameters", - type=_BaseHillasParametersContainer, + type=BaseHillasParametersContainer, ) timing = Field( TimingParametersContainer(), "Timing Parameters", - type=_BaseTimingParametersContainer, + type=BaseTimingParametersContainer, ) leakage = Field(LeakageContainer(), "Leakage Parameters") concentration = Field(ConcentrationContainer(), "Concentration Parameters") From 91fdebac8232e9f3879c5307ea7f2266ac31e408 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 2 Aug 2021 17:26:24 +0200 Subject: [PATCH 32/39] Add option to choose camera frame instead of telescope frame --- ctapipe/image/image_processor.py | 27 +++++++++++------- ctapipe/image/tests/test_image_processor.py | 31 +++++++++++++++++++++ 2 files changed, 48 insertions(+), 10 deletions(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index 85e0f323b74..ce8eae86616 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -65,6 +65,7 @@ def __init__( self, subarray: SubarrayDescription, is_simulation, + use_telescope_frame=True, config=None, parent=None, **kwargs, @@ -94,13 +95,15 @@ def __init__( ) self.check_image = ImageQualityQuery(parent=self) self._is_simulation = is_simulation - telescope_frame = TelescopeFrame() - self.telescope_frame_geometries = { - tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to( - telescope_frame - ) - for tel_id in self.subarray.tel - } + self._use_telescope_frame = use_telescope_frame + if self.use_telescope_frame: + telescope_frame = TelescopeFrame() + self.telescope_frame_geometries = { + tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to( + telescope_frame + ) + for tel_id in self.subarray.tel + } def __call__(self, event: ArrayEventContainer): self._process_telescope_event(event) @@ -191,7 +194,11 @@ def _process_telescope_event(self, event): """ for tel_id, dl1_camera in event.dl1.tel.items(): - geometry_telescope = self.telescope_frame_geometries[tel_id] + if self.use_telescope_frame: + # Use the transformed geometries + geometry = self.telescope_frame_geometries[tel_id] + else: + geometry = dl1_camera.geometry # compute image parameters only if requested to write them dl1_camera.image_mask = self.clean( tel_id=tel_id, @@ -204,7 +211,7 @@ def _process_telescope_event(self, event): image=dl1_camera.image, signal_pixels=dl1_camera.image_mask, peak_time=dl1_camera.peak_time, - geometry=geometry_telescope, + geometry=geometry, ) self.log.debug("params: %s", dl1_camera.parameters.as_dict(recursive=True)) @@ -218,7 +225,7 @@ def _process_telescope_event(self, event): tel_id, image=sim_camera.true_image, signal_pixels=sim_camera.true_image > 0, - geometry=geometry_telescope, + geometry=geometry, peak_time=None, # true image from simulation has no peak time default=DEFAULT_TRUE_IMAGE_PARAMETERS, ) diff --git a/ctapipe/image/tests/test_image_processor.py b/ctapipe/image/tests/test_image_processor.py index 59b788d29ba..603f0b75a49 100644 --- a/ctapipe/image/tests/test_image_processor.py +++ b/ctapipe/image/tests/test_image_processor.py @@ -26,6 +26,37 @@ def test_image_processor(example_event, example_subarray): for dl1tel in example_event.dl1.tel.values(): assert isfinite(dl1tel.image_mask.sum()) assert isfinite(dl1tel.parameters.hillas.length.value) + dl1tel.parameters.hillas.length.to("meter") + assert isfinite(dl1tel.parameters.timing.slope.value) + assert isfinite(dl1tel.parameters.leakage.pixels_width_1) + assert isfinite(dl1tel.parameters.concentration.cog) + assert isfinite(dl1tel.parameters.morphology.num_pixels) + assert isfinite(dl1tel.parameters.intensity_statistics.max) + assert isfinite(dl1tel.parameters.peak_time_statistics.max) + + process_images.check_image.to_table() + + +def test_image_processor_camera_frame(example_event, example_subarray): + """ ensure we get parameters in the camera frame if explicitly specified """ + + calibrate = CameraCalibrator(subarray=example_subarray) + process_images = ImageProcessor( + subarray=example_subarray, + is_simulation=True, + use_telescope_frame=False, + image_cleaner_type="MARSImageCleaner", + ) + + assert isinstance(process_images.clean, MARSImageCleaner) + + calibrate(example_event) + process_images(example_event) + + for dl1tel in example_event.dl1.tel.values(): + assert isfinite(dl1tel.image_mask.sum()) + assert isfinite(dl1tel.parameters.hillas.length.value) + dl1tel.parameters.hillas.length.to("deg") assert isfinite(dl1tel.parameters.timing.slope.value) assert isfinite(dl1tel.parameters.leakage.pixels_width_1) assert isfinite(dl1tel.parameters.concentration.cog) From 42ce511540fbfacfbae03c19fdce752d6de87a31 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 2 Aug 2021 17:29:44 +0200 Subject: [PATCH 33/39] Remove unused import --- ctapipe/io/tests/test_dl1eventsource.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index 47e00d8765b..236bbcc682b 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -1,6 +1,5 @@ import astropy.units as u import numpy as np -import pytest from ctapipe.io import DataLevel, EventSource from ctapipe.io.dl1eventsource import DL1EventSource from ctapipe.utils import get_dataset_path From 25cca1a5a2a4560b9b2a817a3a3f23e63a8856b7 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 2 Aug 2021 17:40:21 +0200 Subject: [PATCH 34/39] Add missing underscore --- ctapipe/image/image_processor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index ce8eae86616..fef48e40db0 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -96,7 +96,7 @@ def __init__( self.check_image = ImageQualityQuery(parent=self) self._is_simulation = is_simulation self._use_telescope_frame = use_telescope_frame - if self.use_telescope_frame: + if self._use_telescope_frame: telescope_frame = TelescopeFrame() self.telescope_frame_geometries = { tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to( @@ -194,7 +194,7 @@ def _process_telescope_event(self, event): """ for tel_id, dl1_camera in event.dl1.tel.items(): - if self.use_telescope_frame: + if self._use_telescope_frame: # Use the transformed geometries geometry = self.telescope_frame_geometries[tel_id] else: From a62e95fa2caeebc02397d14fa552d579d59aba11 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 2 Aug 2021 17:54:27 +0200 Subject: [PATCH 35/39] Fix selection of camera frame geometry --- ctapipe/image/image_processor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index fef48e40db0..c6424c74a99 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -198,7 +198,7 @@ def _process_telescope_event(self, event): # Use the transformed geometries geometry = self.telescope_frame_geometries[tel_id] else: - geometry = dl1_camera.geometry + self.subarray.tel[tel_id].camera.geometry # compute image parameters only if requested to write them dl1_camera.image_mask = self.clean( tel_id=tel_id, From 85e065970069843fb93c047e6bd9645b97f689f8 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 2 Aug 2021 17:55:14 +0200 Subject: [PATCH 36/39] Fix test converting to wrong unit --- ctapipe/image/tests/test_image_processor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ctapipe/image/tests/test_image_processor.py b/ctapipe/image/tests/test_image_processor.py index 603f0b75a49..76215e98c84 100644 --- a/ctapipe/image/tests/test_image_processor.py +++ b/ctapipe/image/tests/test_image_processor.py @@ -26,7 +26,7 @@ def test_image_processor(example_event, example_subarray): for dl1tel in example_event.dl1.tel.values(): assert isfinite(dl1tel.image_mask.sum()) assert isfinite(dl1tel.parameters.hillas.length.value) - dl1tel.parameters.hillas.length.to("meter") + dl1tel.parameters.hillas.length.to("deg") assert isfinite(dl1tel.parameters.timing.slope.value) assert isfinite(dl1tel.parameters.leakage.pixels_width_1) assert isfinite(dl1tel.parameters.concentration.cog) @@ -56,7 +56,7 @@ def test_image_processor_camera_frame(example_event, example_subarray): for dl1tel in example_event.dl1.tel.values(): assert isfinite(dl1tel.image_mask.sum()) assert isfinite(dl1tel.parameters.hillas.length.value) - dl1tel.parameters.hillas.length.to("deg") + dl1tel.parameters.hillas.length.to("meter") assert isfinite(dl1tel.parameters.timing.slope.value) assert isfinite(dl1tel.parameters.leakage.pixels_width_1) assert isfinite(dl1tel.parameters.concentration.cog) From 6f7c5bc9e1be6e8dd9bfcdbf7c245401214c5784 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Mon, 2 Aug 2021 18:05:09 +0200 Subject: [PATCH 37/39] Add missing variable name --- ctapipe/image/image_processor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index c6424c74a99..0787de70407 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -198,7 +198,7 @@ def _process_telescope_event(self, event): # Use the transformed geometries geometry = self.telescope_frame_geometries[tel_id] else: - self.subarray.tel[tel_id].camera.geometry + geometry = self.subarray.tel[tel_id].camera.geometry # compute image parameters only if requested to write them dl1_camera.image_mask = self.clean( tel_id=tel_id, From dee99b0f426f1ca4a6ba82db7106a8e1b28ad184 Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Tue, 3 Aug 2021 19:41:12 +0200 Subject: [PATCH 38/39] Use traits for the selection of the frame --- ctapipe/image/image_processor.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py index 0787de70407..2b6f2de65c9 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -12,7 +12,7 @@ PeakTimeStatisticsContainer, ) from ..core import QualityQuery, TelescopeComponent -from ..core.traits import List, create_class_enum_trait +from ..core.traits import Bool, List, create_class_enum_trait from ..instrument import SubarrayDescription from . import ( ImageCleaner, @@ -60,12 +60,15 @@ class ImageProcessor(TelescopeComponent): image_cleaner_type = create_class_enum_trait( base_class=ImageCleaner, default_value="TailcutsImageCleaner" ) + use_telescope_frame = Bool( + default_value=True, + help="Whether to calculate parameters in the telescope or camera frame", + ).tag(config=True) def __init__( self, subarray: SubarrayDescription, is_simulation, - use_telescope_frame=True, config=None, parent=None, **kwargs, @@ -95,8 +98,7 @@ def __init__( ) self.check_image = ImageQualityQuery(parent=self) self._is_simulation = is_simulation - self._use_telescope_frame = use_telescope_frame - if self._use_telescope_frame: + if self.use_telescope_frame: telescope_frame = TelescopeFrame() self.telescope_frame_geometries = { tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to( @@ -194,7 +196,7 @@ def _process_telescope_event(self, event): """ for tel_id, dl1_camera in event.dl1.tel.items(): - if self._use_telescope_frame: + if self.use_telescope_frame: # Use the transformed geometries geometry = self.telescope_frame_geometries[tel_id] else: From 2e33e0083f6ea4d8c405d124c5e0b45192db3fda Mon Sep 17 00:00:00 2001 From: Lukas Nickel Date: Wed, 4 Aug 2021 11:45:21 +0200 Subject: [PATCH 39/39] Fix unit test --- ctapipe/image/tests/test_image_processor.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ctapipe/image/tests/test_image_processor.py b/ctapipe/image/tests/test_image_processor.py index b6e64f454be..badb32786ca 100644 --- a/ctapipe/image/tests/test_image_processor.py +++ b/ctapipe/image/tests/test_image_processor.py @@ -41,7 +41,6 @@ def test_image_processor_camera_frame(example_event, example_subarray): calibrate = CameraCalibrator(subarray=example_subarray) process_images = ImageProcessor( subarray=example_subarray, - is_simulation=True, use_telescope_frame=False, image_cleaner_type="MARSImageCleaner", )