diff --git a/ctapipe/containers.py b/ctapipe/containers.py index fd46c50ba59..baaa64d5b6e 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -31,6 +31,9 @@ "MonitoringCameraContainer", "MonitoringContainer", "MorphologyContainer", + "BaseHillasParametersContainer", + "CameraHillasParametersContainer", + "CameraTimingParametersContainer", "ParticleClassificationContainer", "PedestalContainer", "PixelStatusContainer", @@ -46,6 +49,7 @@ "SimulatedShowerDistribution", "SimulationConfigContainer", "TelEventIndexContainer", + "BaseTimingParametersContainer", "TimingParametersContainer", "TriggerContainer", "WaveformCalibrationContainer", @@ -107,11 +111,25 @@ class TelEventIndexContainer(Container): tel_id = Field(0, "telescope identifier") -class HillasParametersContainer(Container): - container_prefix = "hillas" +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) r = Field(nan * u.m, "radial coordinate of centroid", unit=u.m) @@ -123,8 +141,33 @@ 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 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", + unit=u.deg, + ) + fov_lat = Field( + nan * u.deg, + "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) + 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) class LeakageContainer(Container): @@ -167,16 +210,13 @@ class ConcentrationContainer(Container): pixel = Field(nan, "Percentage of photo-electrons in the brightest pixel") -class TimingParametersContainer(Container): +class BaseTimingParametersContainer(Container): """ - Slope and Intercept of a linear regression of the arrival times - along the shower main axis + Base container for timing parameters to + allow the CameraTimingParametersContainer to + be assigned to an ImageParametersContainer as well. """ - container_prefix = "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, @@ -185,6 +225,31 @@ class TimingParametersContainer(Container): ) +class CameraTimingParametersContainer(BaseTimingParametersContainer): + """ + Slope and Intercept of a linear regression of the arrival times + along the shower main axis in the camera frame. + """ + + container_prefix = "camera_frame_timing" + slope = Field( + nan / u.m, "Slope of arrival times along main shower axis", unit=1 / u.m + ) + + +class TimingParametersContainer(BaseTimingParametersContainer): + """ + Slope and Intercept of a linear regression of the arrival times + along the shower main axis in a + spherical system centered on the pointing position (TelescopeFrame) + """ + + container_prefix = "timing" + slope = Field( + nan / u.deg, "Slope of arrival times along main shower axis", unit=1 / u.deg + ) + + class MorphologyContainer(Container): """ Parameters related to pixels surviving image cleaning """ @@ -225,8 +290,16 @@ class ImageParametersContainer(Container): """ Collection of image parameters """ container_prefix = "params" - hillas = Field(HillasParametersContainer(), "Hillas Parameters") - timing = Field(TimingParametersContainer(), "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") @@ -260,7 +333,6 @@ class DL1CameraContainer(Container): dtype=np.float32, ndim=1, ) - image_mask = Field( None, "Boolean numpy array where True means the pixel has passed cleaning." @@ -268,7 +340,6 @@ class DL1CameraContainer(Container): dtype=np.bool_, ndim=1, ) - parameters = Field(None, "Image parameters", type=ImageParametersContainer) diff --git a/ctapipe/image/concentration.py b/ctapipe/image/concentration.py index 472d9435213..fcd78cc5132 100644 --- a/ctapipe/image/concentration.py +++ b/ctapipe/image/concentration.py @@ -1,9 +1,13 @@ import numpy as np import astropy.units as u -from ..instrument import CameraGeometry -from ..containers import ConcentrationContainer +from ..containers import ( + ConcentrationContainer, + HillasParametersContainer, + CameraHillasParametersContainer, +) from .hillas import camera_to_shower_coordinates +from ..instrument import CameraGeometry from ..utils.quantities import all_to_value __all__ = ["concentration_parameters"] @@ -20,11 +24,30 @@ def concentration_parameters(geom: CameraGeometry, image, hillas_parameters): """ h = hillas_parameters - unit = h.x.unit - - pix_x, pix_y, x, y, length, width, pixel_width = all_to_value( - geom.pix_x, geom.pix_y, h.x, h.y, h.length, h.width, geom.pixel_width, unit=unit - ) + if isinstance(h, CameraHillasParametersContainer): + unit = h.x.unit + pix_x, pix_y, x, y, length, width, pixel_width = all_to_value( + geom.pix_x, + geom.pix_y, + h.x, + h.y, + h.length, + h.width, + geom.pixel_width, + unit=unit, + ) + elif isinstance(h, HillasParametersContainer): + unit = h.fov_lon.unit + pix_x, pix_y, x, y, length, width, pixel_width = all_to_value( + geom.pix_x, + geom.pix_y, + h.fov_lon, + h.fov_lat, + h.length, + h.width, + geom.pixel_width, + unit=unit, + ) delta_x = pix_x - x delta_y = pix_y - y diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py index e8312b1701d..8fdd405a857 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 CameraHillasParametersContainer, HillasParametersContainer 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): @@ -198,9 +195,24 @@ def hillas_parameters(geom, image): np.sum(((((b * A) + (a * B) + (-c * C))) ** 2.0) * image) ) / (2 * width) + if unit.is_equivalent(u.m): + return CameraHillasParametersContainer( + 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, + ) return HillasParametersContainer( - x=u.Quantity(cog_x, unit), - y=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/image_processor.py b/ctapipe/image/image_processor.py index 26dca72fb61..c39563e22ad 100644 --- a/ctapipe/image/image_processor.py +++ b/ctapipe/image/image_processor.py @@ -1,18 +1,18 @@ """ High level image processing (ImageProcessor Component) """ +from ctapipe.coordinates import TelescopeFrame import numpy as np - from ..containers import ( ArrayEventContainer, - ImageParametersContainer, IntensityStatisticsContainer, - PeakTimeStatisticsContainer, + ImageParametersContainer, TimingParametersContainer, + 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, @@ -25,6 +25,7 @@ ) +# avoid use of base containers for unparameterized images DEFAULT_IMAGE_PARAMETERS = ImageParametersContainer() DEFAULT_TRUE_IMAGE_PARAMETERS = ImageParametersContainer() DEFAULT_TRUE_IMAGE_PARAMETERS.intensity_statistics = IntensityStatisticsContainer( @@ -59,6 +60,10 @@ 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, config=None, parent=None, **kwargs @@ -85,6 +90,14 @@ def __init__( self.image_cleaner_type, subarray=subarray, parent=self ) self.check_image = ImageQualityQuery(parent=self) + 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) @@ -94,11 +107,11 @@ def _parameterize_image( tel_id, image, signal_pixels, + geometry, peak_time=None, default=DEFAULT_IMAGE_PARAMETERS, ) -> ImageParametersContainer: """Apply image cleaning and calculate image features - Parameters ---------- tel_id: int @@ -109,15 +122,12 @@ def _parameterize_image( image mask peak_time: np.ndarray peak time image - Returns ------- ImageParametersContainer: cleaning mask, parameters """ - tel = self.subarray.tel[tel_id] - geometry = tel.camera.geometry image_selected = image[signal_pixels] # check if image can be parameterized: @@ -176,9 +186,13 @@ 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(): + if self.use_telescope_frame: + # Use the transformed geometries + geometry = self.telescope_frame_geometries[tel_id] + else: + 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, @@ -191,6 +205,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, ) self.log.debug("params: %s", dl1_camera.parameters.as_dict(recursive=True)) @@ -205,6 +220,7 @@ def _process_telescope_event(self, event): tel_id, image=sim_camera.true_image, signal_pixels=sim_camera.true_image > 0, + 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_hillas.py b/ctapipe/image/tests/test_hillas.py index 51fccb307a7..16fff8c0c27 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -1,14 +1,18 @@ -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 astropy.coordinates import Angle, SkyCoord from astropy import units as u import numpy as np from numpy import isclose, zeros_like 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 +from ctapipe.containers import ( + HillasParametersContainer, + CameraHillasParametersContainer, +) def create_sample_image( @@ -104,6 +108,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) @@ -248,3 +257,84 @@ def test_single_pixel(): assert hillas.length.value == 0 assert hillas.width.value == 0 assert np.isnan(hillas.psi) + + +def test_reconstruction_in_telescope_frame(): + """ + Compare the reconstruction in the telescope + and camera frame. + """ + np.random.seed(42) + + 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 + 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 + + 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: + # 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) + 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, + 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) + + transformed_length = get_transformed_length( + telescope_result, telescope_frame, camera_frame + ) + assert u.isclose(transformed_length, camera_result.length, rtol=0.01) + + transformed_width = get_transformed_width( + telescope_result, telescope_frame, camera_frame + ) + assert u.isclose(transformed_width, camera_result.width, rtol=0.01) diff --git a/ctapipe/image/tests/test_image_processor.py b/ctapipe/image/tests/test_image_processor.py index 86520e06d56..badb32786ca 100644 --- a/ctapipe/image/tests/test_image_processor.py +++ b/ctapipe/image/tests/test_image_processor.py @@ -24,6 +24,36 @@ 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("deg") + 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, + 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("meter") assert isfinite(dl1tel.parameters.timing.slope.value) assert isfinite(dl1tel.parameters.leakage.pixels_width_1) assert isfinite(dl1tel.parameters.concentration.cog) diff --git a/ctapipe/image/tests/test_timing_parameters.py b/ctapipe/image/tests/test_timing_parameters.py index 4eec5891a69..93efd466ea9 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.default_rng(0) 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.default_rng(0) 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.default_rng(0) peak_time = intercept + grad * geom.pix_x.value diff --git a/ctapipe/image/timing.py b/ctapipe/image/timing.py index e6fa99c8942..b7eb81c2f96 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 ( + CameraTimingParametersContainer, + TimingParametersContainer, + CameraHillasParametersContainer, + HillasParametersContainer, +) 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, 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, HillasParametersContainer): + unit = h.fov_lon.unit + pix_x, pix_y, x, y, length, width = all_to_value( + geom.pix_x, geom.pix_y, h.fov_lon, h.fov_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,10 @@ 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) + if unit.is_equivalent(u.m): + return CameraTimingParametersContainer( + slope=beta[0] / unit, intercept=beta[1], deviation=deviation + ) return TimingParametersContainer( slope=beta[0] / unit, intercept=beta[1], deviation=deviation ) diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py index b9413192c98..c45eeafda8d 100644 --- a/ctapipe/io/datawriter.py +++ b/ctapipe/io/datawriter.py @@ -37,8 +37,9 @@ # (meaning readers need to update scripts) # - increase the minor number if new columns or datasets are added # - increase the patch number if there is a small bugfix to the model. -DATA_MODEL_VERSION = "v2.0.0" +DATA_MODEL_VERSION = "v2.1.0" DATA_MODEL_CHANGE_HISTORY = """ +- v2.1.0: hillas and timing parameters are per default saved in telescope frame (degree) as opposed to camera frame (m) - v2.0.0: Match optics and camera tables using indices instead of names - 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 diff --git a/ctapipe/io/dl1eventsource.py b/ctapipe/io/dl1eventsource.py index 1338278446a..38e28bd89e6 100644 --- a/ctapipe/io/dl1eventsource.py +++ b/ctapipe/io/dl1eventsource.py @@ -11,6 +11,7 @@ ArrayEventContainer, DL1CameraContainer, EventIndexContainer, + CameraHillasParametersContainer, HillasParametersContainer, IntensityStatisticsContainer, LeakageContainer, @@ -19,6 +20,7 @@ SimulatedShowerContainer, SimulatedEventContainer, PeakTimeStatisticsContainer, + CameraTimingParametersContainer, TimingParametersContainer, TriggerContainer, ImageParametersContainer, @@ -43,6 +45,7 @@ "v1.1.0", "v1.2.0", "v2.0.0", + "v2.1.0", ] @@ -265,8 +268,16 @@ def _generate_events(self): tel.name: self.reader.read( f"/dl1/event/telescope/parameters/{tel.name}", containers=[ - HillasParametersContainer(), - TimingParametersContainer(), + ( + HillasParametersContainer() + if (self.datamodel_version >= "v2.1.0") + else CameraHillasParametersContainer(prefix="hillas") + ), + ( + TimingParametersContainer() + if (self.datamodel_version >= "v2.1.0") + else CameraTimingParametersContainer(prefix="timing") + ), LeakageContainer(), ConcentrationContainer(), MorphologyContainer(), @@ -282,7 +293,11 @@ def _generate_events(self): tel.name: self.reader.read( f"/simulation/event/telescope/parameters/{tel.name}", containers=[ - HillasParametersContainer(), + ( + HillasParametersContainer() + if (self.datamodel_version >= "v2.1.0") + else CameraHillasParametersContainer(prefix="hillas") + ), LeakageContainer(), ConcentrationContainer(), MorphologyContainer(), @@ -406,7 +421,6 @@ def _generate_events(self): 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( diff --git a/ctapipe/io/tests/test_dl1eventsource.py b/ctapipe/io/tests/test_dl1eventsource.py index 02273f46a64..236bbcc682b 100644 --- a/ctapipe/io/tests/test_dl1eventsource.py +++ b/ctapipe/io/tests/test_dl1eventsource.py @@ -1,16 +1,10 @@ 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 -@pytest.fixture(scope="module") -def dl1_dir(tmp_path_factory): - return tmp_path_factory.mktemp("dl1") - - def test_is_compatible(dl1_file): simtel_path = get_dataset_path("gamma_test_large.simtel.gz") assert not DL1EventSource.is_compatible(simtel_path) @@ -56,15 +50,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 event.simulation.tel[tel].true_parameters.hillas.x != np.nan + reco_lons.append( + event.simulation.tel[tel].true_parameters.hillas.fov_lon.value + ) + 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): @@ -75,18 +82,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 event.dl1.tel[tel].parameters.hillas.x != np.nan + reco_lons.append( + event.simulation.tel[tel].true_parameters.hillas.fov_lon.value + ) + 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 event.dl1.tel[tel].parameters.hillas.x != np.nan + reco_lons.append( + event.simulation.tel[tel].true_parameters.hillas.fov_lon.value + ) + 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): diff --git a/ctapipe/io/tests/test_hdf5.py b/ctapipe/io/tests/test_hdf5.py index bdf69d3ba2e..13d553a2dd9 100644 --- a/ctapipe/io/tests/test_hdf5.py +++ b/ctapipe/io/tests/test_hdf5.py @@ -70,7 +70,7 @@ def test_append_container(tmp_path): def test_read_multiple_containers(tmp_path): path = tmp_path / "test_append.h5" hillas_parameter_container = HillasParametersContainer( - x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m + fov_lon=1 * u.deg, fov_lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) leakage_container = LeakageContainer( @@ -83,7 +83,7 @@ def test_read_multiple_containers(tmp_path): writer.write("params", [hillas_parameter_container, leakage_container]) df = pd.read_hdf(path, key="/dl1/params") - assert "hillas_x" in df.columns + assert "hillas_fov_lon" in df.columns assert "leakage_pixels_width_1" in df.columns # test reading both containers separately @@ -129,7 +129,7 @@ def test_read_without_prefixes(tmp_path): path = tmp_path / "test.h5" hillas_parameter_container = HillasParametersContainer( - x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m + fov_lon=1 * u.deg, fov_lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) leakage_container = LeakageContainer( @@ -143,7 +143,7 @@ def test_read_without_prefixes(tmp_path): writer.write("params", [hillas_parameter_container, leakage_container]) df = pd.read_hdf(path, key="/dl1/params") - assert "x" in df.columns + assert "fov_lon" in df.columns assert "pixels_width_1" in df.columns # call with prefixes=False @@ -189,18 +189,26 @@ def test_read_duplicated_container_types(tmp_path): path = tmp_path / "test.h5" hillas_config_1 = HillasParametersContainer( - x=1 * u.m, y=2 * u.m, length=3 * u.m, width=4 * u.m, prefix="hillas_1" + 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( - x=2 * u.m, y=3 * u.m, length=4 * u.m, width=5 * u.m, prefix="hillas_2" + fov_lon=2 * u.deg, + fov_lat=3 * u.deg, + length=4 * u.deg, + width=5 * u.deg, + prefix="hillas_2", ) with HDF5TableWriter(path, group_name="dl1", add_prefix=True) as writer: writer.write("params", [hillas_config_1, hillas_config_2]) df = pd.read_hdf(path, key="/dl1/params") - assert "hillas_1_x" in df.columns - assert "hillas_2_x" in df.columns + assert "hillas_1_fov_lon" in df.columns + assert "hillas_2_fov_lon" in df.columns with HDF5TableReader(path) as reader: generator = reader.read( @@ -225,7 +233,7 @@ def test_custom_prefix(tmp_path): path = tmp_path / "test.h5" container = HillasParametersContainer( - x=1 * u.m, y=1 * u.m, length=1 * u.m, width=1 * u.m + fov_lon=1 * u.deg, fov_lat=1 * u.deg, length=1 * u.deg, width=1 * u.deg ) container.prefix = "custom" with HDF5TableWriter(path, group_name="dl1", add_prefix=True) as writer: diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py index 362d9bb0597..89f781f777f 100644 --- a/ctapipe/reco/hillas_intersection.py +++ b/ctapipe/reco/hillas_intersection.py @@ -16,18 +16,22 @@ 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, ) -import copy import warnings from ctapipe.core import traits @@ -144,33 +148,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 +205,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 +222,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 +251,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 +265,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 +408,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/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 diff --git a/ctapipe/reco/tests/test_ImPACT.py b/ctapipe/reco/tests/test_ImPACT.py index 1faf086f237..475aa9509d7 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 6e3741dc3e7..e05af2d370c 100644 --- a/ctapipe/reco/tests/test_hillas_intersection.py +++ b/ctapipe/reco/tests/test_hillas_intersection.py @@ -3,7 +3,7 @@ from numpy.testing import assert_allclose import numpy as np from astropy.coordinates import SkyCoord -from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame +from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame, TelescopeFrame from ctapipe.containers import HillasParametersContainer from ctapipe.io import EventSource @@ -82,13 +82,13 @@ def test_intersection_xmax_reco(): hillas_dict = { 1: HillasParametersContainer( - x=-(delta / focal_length) * u.rad, - y=((0 * u.m) / focal_length) * u.rad, + fov_lon=-(delta / focal_length) * u.rad, + fov_lat=((0 * u.m) / focal_length) * u.rad, intensity=1, ), 2: HillasParametersContainer( - x=((0 * u.m) / focal_length) * u.rad, - y=-(delta / focal_length) * u.rad, + fov_lon=((0 * u.m) / focal_length) * u.rad, + fov_lat=-(delta / focal_length) * u.rad, intensity=1, ), } @@ -180,32 +180,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 = HillasParametersContainer( - x=cog_coords_nom_1.fov_lat, - y=cog_coords_nom_1.fov_lon, + fov_lat=cog_coords_nom_1.fov_lat, + fov_lon=cog_coords_nom_1.fov_lon, intensity=100, psi=0 * u.deg, ) hillas_2 = HillasParametersContainer( - x=cog_coords_nom_2.fov_lat, - y=cog_coords_nom_2.fov_lon, + fov_lat=cog_coords_nom_2.fov_lat, + fov_lon=cog_coords_nom_2.fov_lon, intensity=100, psi=45 * u.deg, ) hillas_3 = HillasParametersContainer( - x=cog_coords_nom_3.fov_lat, - y=cog_coords_nom_3.fov_lon, + fov_lat=cog_coords_nom_3.fov_lat, + fov_lon=cog_coords_nom_3.fov_lon, intensity=100, psi=90 * u.deg, ) @@ -257,10 +256,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 +273,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 +320,14 @@ 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() + ) 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 diff --git a/ctapipe/utils/datasets.py b/ctapipe/utils/datasets.py index a4e5af99dc5..42f187f2b82 100644 --- a/ctapipe/utils/datasets.py +++ b/ctapipe/utils/datasets.py @@ -38,7 +38,6 @@ def get_searchpath_dirs(searchpath=os.getenv("CTAPIPE_SVC_PATH"), url=DEFAULT_UR searchpaths = [] else: searchpaths = [Path(p) for p in os.path.expandvars(searchpath).split(":")] - searchpaths.append(get_cache_path(url)) return searchpaths diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index a202715b7d6..8e0de18578b 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -10,7 +10,10 @@ TelescopeDescription, PixelShape, ) -from ctapipe.containers import HillasParametersContainer +from ctapipe.containers import ( + CameraHillasParametersContainer, + HillasParametersContainer, +) import numpy as np from astropy import units as u @@ -49,7 +52,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 ) @@ -141,12 +144,12 @@ 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) # test using hillas params CameraFrame: 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 @@ -161,9 +164,8 @@ def test_array_display(): ) gradient_dict = {1: timing_rot20.slope.value, 2: timing_rot20.slope.value} core_dict = { - tel_id: dl1.parameters.core.psi - for tel_id, dl1 in event.dl1.tel.items() - } + tel_id: dl1.parameters.core.psi for tel_id, dl1 in event.dl1.tel.items() + } ad.set_vector_hillas( hillas_dict=hillas_dict, core_dict=core_dict, @@ -176,16 +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( @@ -200,16 +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( @@ -221,10 +211,7 @@ def test_array_display(): ) # test negative time_gradients - gradient_dict = { - 1: -0.03, - 2: -0.02, - } + gradient_dict = {1: -0.03, 2: -0.02} ad.set_vector_hillas( hillas_dict=hillas_dict, core_dict=core_dict, @@ -233,10 +220,7 @@ def test_array_display(): angle_offset=0 * u.deg, ) # and very small - gradient_dict = { - 1: 0.003, - 2: 0.002, - } + gradient_dict = {1: 0.003, 2: 0.002} ad.set_vector_hillas( hillas_dict=hillas_dict, core_dict=core_dict,