From 12d373850d8b12e4550d21aefe807ba6174a1929 Mon Sep 17 00:00:00 2001 From: nieves Date: Fri, 15 Dec 2023 15:54:06 +0100 Subject: [PATCH] removed camera frame --- ctapipe/containers.py | 36 +-- ctapipe/image/ellipsoid.py | 91 +++---- ctapipe/image/tests/test_ellipsoid.py | 333 +++++++------------------- 3 files changed, 116 insertions(+), 344 deletions(-) diff --git a/ctapipe/containers.py b/ctapipe/containers.py index b2ed441ff99..ae8d97101b0 100644 --- a/ctapipe/containers.py +++ b/ctapipe/containers.py @@ -34,7 +34,6 @@ "MorphologyContainer", "BaseHillasParametersContainer", "CameraHillasParametersContainer", - "CameraImageFitParametersContainer", "ImageFitParametersContainer", "CameraTimingParametersContainer", "ParticleClassificationContainer", @@ -282,37 +281,6 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer): psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg) -class CameraImageFitParametersContainer(CameraHillasParametersContainer): - """ - Hillas Parameters after fitting. The cog position - is given in meter from the camera center. - """ - - skewness_uncertainty = Field(nan, "measure of skewness uncertainty") - likelihood = Field(nan, "measure of likelihood") - goodness_of_fit = Field( - nan, "measure of goodness of fit, mean likelihood subtracted to the likelihood" - ) - n_pixels = Field(nan, "number of pixels used in the fit") - free_parameters = Field(nan, "number of free parameters") - is_valid = Field(False, "True if the fit is valid") - is_accurate = Field( - False, "returns True if the fit is accurate. If False, the fit is not reliable." - ) - - x_uncertainty = Field(nan * u.m, "centroid x uncertainty", unit=u.m) - y_uncertainty = Field(nan * u.m, "centroid y uncertainty", unit=u.m) - r_uncertainty = Field(nan * u.m, "centroid r uncertainty", unit=u.m) - phi_uncertainty = Field( - nan * u.deg, "polar coordinate of centroid uncertainty", unit=u.deg - ) - psi_uncertainty = Field( - nan * u.deg, "Uncertainty in rotation angle of ellipse", unit=u.deg - ) - amplitude = Field(nan, "Amplitude of the fitted model") - amplitude_uncertainty = Field(nan, "error in amplitude from the fit") - - class HillasParametersContainer(BaseHillasParametersContainer): """ Hillas Parameters in a spherical system centered on the pointing position @@ -352,8 +320,8 @@ class ImageFitParametersContainer(HillasParametersContainer): goodness_of_fit = Field( nan, "measure of goodness of fit, mean likelihood subtracted to the likelihood" ) - n_pixels = Field(nan, "number of pixels used in the fit") - free_parameters = Field(nan, "number of free parameters") + n_pixels = Field(-1, "number of pixels used in the fit") + free_parameters = Field(-1, "number of free parameters") is_valid = Field(False, "True if the fit is valid") is_accurate = Field( False, "returns True if the fit is accurate. If False, the fit is not reliable." diff --git a/ctapipe/image/ellipsoid.py b/ctapipe/image/ellipsoid.py index 2e048f411b7..7f3eb1444aa 100644 --- a/ctapipe/image/ellipsoid.py +++ b/ctapipe/image/ellipsoid.py @@ -19,16 +19,14 @@ ) from ctapipe.image.toymodel import SkewedGaussian -from ..containers import CameraImageFitParametersContainer, ImageFitParametersContainer +from ..containers import ImageFitParametersContainer PED_TABLE = { "LSTCam": 2.8, "NectarCam": 2.3, "FlashCam": 2.3, - "SSTCam": 0.5, - "CHEC": 0.5, - "DUMMY": 0, - "testcam": 0, + "SST-Camera": 0.5, + "testcam": 0.5, } SPE_WIDTH = 0.5 @@ -46,7 +44,6 @@ class PDFType(Enum): def create_initial_guess(geometry, hillas, size): """ - This function computes the seeds of the fit with the Hillas parameters Parameters @@ -66,12 +63,8 @@ def create_initial_guess(geometry, hillas, size): unit = geometry.pix_x.unit initial_guess = {} - if unit.is_equivalent(u.m): - initial_guess["cog_x"] = hillas.x.to_value(unit) - initial_guess["cog_y"] = hillas.y.to_value(unit) - else: - initial_guess["cog_x"] = hillas.fov_lon.to_value(unit) - initial_guess["cog_y"] = hillas.fov_lat.to_value(unit) + initial_guess["cog_x"] = hillas.fov_lon.to_value(unit) + initial_guess["cog_y"] = hillas.fov_lat.to_value(unit) initial_guess["length"] = hillas.length.to_value(unit) initial_guess["width"] = hillas.width.to_value(unit) @@ -79,7 +72,7 @@ def create_initial_guess(geometry, hillas, size): initial_guess["skewness"] = hillas.skewness - if (hillas.width.to_value(unit) == 0) or (hillas.length.to_value(unit) == 0): + if (hillas.width.value == 0) or (hillas.length.value == 0): raise ImageFitParameterizationError("Hillas width and/or length is zero") initial_guess["amplitude"] = size @@ -88,7 +81,6 @@ def create_initial_guess(geometry, hillas, size): def boundaries(geometry, image, dilated_mask, hillas, pdf): """ - Computes the boundaries of the fit. Parameters @@ -113,18 +105,13 @@ def boundaries(geometry, image, dilated_mask, hillas, pdf): unit = geometry.pix_x.unit camera_radius = geometry.radius.to_value(unit) - # Using dilated image - row_image = image.copy() - row_image[~dilated_mask] = 0.0 - row_image[row_image < 0] = 0.0 - max_x = np.max(x[dilated_mask]) min_x = np.min(x[dilated_mask]) max_y = np.max(y[dilated_mask]) min_y = np.min(y[dilated_mask]) ang_unit = hillas.psi.unit - psi_unc = 10 * u.deg + psi_unc = 10 * unit psi_min, psi_max = max( hillas.psi.value - psi_unc.to_value(ang_unit), -np.pi / 2 ), min(hillas.psi.value + psi_unc.to_value(ang_unit), np.pi / 2) @@ -141,10 +128,7 @@ def boundaries(geometry, image, dilated_mask, hillas, pdf): (max_x - min_x) ** 2 + (max_y - min_y) ** 2 ) # maximum distance of the shower in the longitudinal direction - if unit.is_equivalent(u.m): - width_unc = 0.05 * u.m - else: - width_unc = 0.1 * u.deg + width_unc = 0.1 * unit length_min, length_max = hillas.length.to_value(unit), long_dis width_min, width_max = hillas.width.to_value(unit), hillas.width.to_value( unit @@ -156,9 +140,15 @@ def boundaries(geometry, image, dilated_mask, hillas, pdf): ) # Guess from Hillas unit tests if pdf == PDFType.gaussian: - amplitude = np.sum(row_image) / (2 * np.pi * width_min * length_min) + amplitude = np.sum(image[(dilated_mask) & (image > 0)]) / ( + 2 * np.pi * width_min * length_min + ) else: - amplitude = np.sum(row_image) / scale / (2 * np.pi * width_min) + amplitude = ( + np.sum(image[(dilated_mask) & (image > 0)]) + / scale + / (2 * np.pi * width_min) + ) bounds = [ (cogx_min, cogx_max), @@ -182,11 +172,10 @@ def image_fit_parameters( image, n_row, cleaned_mask, - pdf=PDFType("skewed"), + pdf=PDFType.skewed, bounds=None, ): """ - Computes image parameters for a given shower image. Implementation based on https://arxiv.org/pdf/1211.0254.pdf @@ -326,6 +315,7 @@ def likelihood(cog_x, cog_y, psi, length, width, skewness, amplitude): fit_rcog = np.linalg.norm([pars[0], pars[1]]) fit_phi = np.arctan2(pars[1], pars[0]) + # The uncertainty in r and phi is calculated by propagating errors of the x and y coordinates b = pars[1] ** 2 + pars[0] ** 2 A = (-pars[1] / b) ** 2 B = (pars[0] / b) ** 2 @@ -337,46 +327,25 @@ def likelihood(cog_x, cog_y, psi, length, width, skewness, amplitude): delta_x = pix_x.value - pars[0] delta_y = pix_y.value - pars[1] + # calculate higher order moments along shower axes longitudinal = delta_x * np.cos(pars[2]) + delta_y * np.sin(pars[2]) - m4_long = np.average(longitudinal**4, weights=dilated_image) + m3_long = np.average(longitudinal**3, weights=image) + skewness_ = m3_long / pars[3] ** 3 + + m4_long = np.average(longitudinal**4, weights=image) kurtosis_long = m4_long / pars[3] ** 4 - skewness_long = pars[5] - skewness_uncertainty = errors[5] + if pdf == PDFType.gaussian: + skewness_long = skewness_ + skewness_uncertainty = np.nan + else: + skewness_long = pars[5] + skewness_uncertainty = errors[5] + amplitude = pars[6] amplitude_uncertainty = errors[6] - if unit.is_equivalent(u.m): - return CameraImageFitParametersContainer( - x=u.Quantity(pars[0], unit), - x_uncertainty=u.Quantity(errors[0], unit), - y=u.Quantity(pars[1], unit), - y_uncertainty=u.Quantity(errors[1], unit), - r=u.Quantity(fit_rcog, unit), - r_uncertainty=u.Quantity(fit_rcog_err, unit), - phi=Angle(fit_phi, unit=u.rad), - phi_uncertainty=Angle(fit_phi_err, unit=u.rad), - intensity=size, - amplitude=amplitude, - amplitude_uncertainty=amplitude_uncertainty, - length=u.Quantity(pars[3], unit), - length_uncertainty=u.Quantity(errors[3], unit), - width=u.Quantity(pars[4], unit), - width_uncertainty=u.Quantity(errors[4], unit), - psi=Angle(pars[2], unit=u.rad), - psi_uncertainty=Angle(errors[2], unit=u.rad), - skewness=skewness_long, - skewness_uncertainty=skewness_uncertainty, - kurtosis=kurtosis_long, - likelihood=likelihood, - goodness_of_fit=goodness_of_fit, - n_pixels=np.count_nonzero(cleaned_image), - free_parameters=m.nfit, - is_valid=m.valid, - is_accurate=m.accurate, - ) - return ImageFitParametersContainer( fov_lon=u.Quantity(pars[0], unit), fov_lon_uncertainty=u.Quantity(errors[0], unit), diff --git a/ctapipe/image/tests/test_ellipsoid.py b/ctapipe/image/tests/test_ellipsoid.py index a16a324e244..25b1c8806a8 100644 --- a/ctapipe/image/tests/test_ellipsoid.py +++ b/ctapipe/image/tests/test_ellipsoid.py @@ -2,15 +2,12 @@ import numpy.ma as ma import pytest from astropy import units as u -from astropy.coordinates import Angle, SkyCoord +from astropy.coordinates import Angle from numpy import isclose from numpy.testing import assert_allclose from pytest import approx -from ctapipe.containers import ( - CameraImageFitParametersContainer, - ImageFitParametersContainer, -) +from ctapipe.containers import ImageFitParametersContainer from ctapipe.coordinates import TelescopeFrame from ctapipe.image import hillas_parameters, tailcuts_clean, toymodel from ctapipe.image.concentration import concentration_parameters @@ -55,13 +52,15 @@ def create_sample_image( def test_boundaries(prod5_lst): # Test default functin for finding the boundaries of the fit geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) image, clean_mask = create_sample_image(geometry=geom) cleaned_image = image.copy() cleaned_image[~clean_mask] = 0.0 - hillas = hillas_parameters(geom, cleaned_image) - bounds = boundaries(geom, image, clean_mask, hillas, pdf=PDFType("gaussian")) + hillas = hillas_parameters(geom_nom, cleaned_image) + bounds = boundaries(geom_nom, image, clean_mask, hillas, pdf=PDFType("gaussian")) for i in range(len(bounds)): assert bounds[i][1] > bounds[i][0] # upper limit > lower limit @@ -86,6 +85,8 @@ def compare_fit_params(fit1, fit2): def test_fit_selected(prod5_lst): geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) image, clean_mask = create_sample_image(geometry=geom) image_zeros = image.copy() @@ -94,13 +95,13 @@ def test_fit_selected(prod5_lst): image_selected = ma.masked_array(image.copy(), mask=~clean_mask) results = image_fit_parameters( - geom, + geom_nom, image_zeros, n_row=2, cleaned_mask=clean_mask, ) results_selected = image_fit_parameters( - geom, + geom_nom, image_selected, n_row=2, cleaned_mask=clean_mask, @@ -110,10 +111,12 @@ def test_fit_selected(prod5_lst): def test_dilation(prod5_lst): geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) image, clean_mask = create_sample_image(geometry=geom) results = image_fit_parameters( - geom, + geom_nom, image, n_row=0, cleaned_mask=clean_mask, @@ -122,7 +125,7 @@ def test_dilation(prod5_lst): assert_allclose(results.intensity, np.sum(image[clean_mask]), rtol=1e-4, atol=1e-4) results = image_fit_parameters( - geom, + geom_nom, image, n_row=2, cleaned_mask=clean_mask, @@ -133,11 +136,13 @@ def test_dilation(prod5_lst): def test_imagefit_failure(prod5_lst): geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) blank_image = np.zeros(geom.n_pixels) with pytest.raises(ImageFitParameterizationError): image_fit_parameters( - geom, + geom_nom, blank_image, n_row=2, cleaned_mask=(blank_image == 1), @@ -148,14 +153,6 @@ def test_fit_container(prod5_lst): geom = prod5_lst.camera.geometry image, clean_mask = create_sample_image(psi="0d", geometry=geom) - params = image_fit_parameters( - geom, - image, - n_row=2, - cleaned_mask=clean_mask, - ) - assert isinstance(params, CameraImageFitParametersContainer) - geom_telescope_frame = geom.transform_to(TelescopeFrame()) params = image_fit_parameters( geom_telescope_frame, @@ -168,43 +165,45 @@ def test_fit_container(prod5_lst): def test_truncated(prod5_lst): geom = prod5_lst.camera.geometry - widths = u.Quantity([0.03, 0.05], u.m) - lengths = u.Quantity([0.3, 0.2], u.m) + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) + widths = u.Quantity([0.05, 0.06], u.deg) + lengths = u.Quantity([0.4, 0.5], u.deg) intensity = 5000 - xs = u.Quantity([0.8, 0.7, -0.7, -0.8], u.m) - ys = u.Quantity([-0.7, -0.8, 0.8, 0.7], u.m) + xs = u.Quantity([0.8, 0.9, -0.9, -0.8], u.deg) + ys = u.Quantity([-0.9, -0.8, 0.8, 0.9], u.deg) for x, y, width, length in zip(xs, ys, widths, lengths): image, clean_mask = create_sample_image( - geometry=geom, x=x, y=y, length=length, width=width, intensity=intensity + geometry=geom_nom, x=x, y=y, length=length, width=width, intensity=intensity ) cleaned_image = image.copy() cleaned_image[~clean_mask] = 0.0 # Gaussian result = image_fit_parameters( - geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("gaussian") + geom_nom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("gaussian") ) - hillas = hillas_parameters(geom, cleaned_image) + hillas = hillas_parameters(geom_nom, cleaned_image) assert result.length.value > hillas.length.value - conc_fit = concentration_parameters(geom, cleaned_image, result).core - conc_hillas = concentration_parameters(geom, cleaned_image, hillas).core + conc_fit = concentration_parameters(geom_nom, cleaned_image, result).core + conc_hillas = concentration_parameters(geom_nom, cleaned_image, hillas).core assert conc_fit > conc_hillas assert conc_fit > 0.4 # Skewed result = image_fit_parameters( - geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("skewed") + geom_nom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("skewed") ) assert result.length.value > hillas.length.value - conc_fit = concentration_parameters(geom, cleaned_image, result).core + conc_fit = concentration_parameters(geom_nom, cleaned_image, result).core assert conc_fit > conc_hillas assert conc_fit > 0.4 @@ -212,168 +211,38 @@ def test_truncated(prod5_lst): def test_percentage(prod5_lst): geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) # Gaussian image, clean_mask = create_sample_image(psi="0d", geometry=geom) fit = image_fit_parameters( - geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("gaussian") + geom_nom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("gaussian") ) cleaned_image = image.copy() cleaned_image[~clean_mask] = 0.0 - conc = concentration_parameters(geom, image, fit) + conc = concentration_parameters(geom_nom, image, fit) signal_inside_ellipse = conc.core if fit.is_valid and fit.is_accurate: assert signal_inside_ellipse > 0.3 -def test_with_toy_mst_tel(prod5_mst_flashcam): - rng = np.random.default_rng(42) - - geom = prod5_mst_flashcam.camera.geometry - - width = 0.03 * u.m - length = 0.15 * u.m - intensity = 500 - - xs = u.Quantity([0.2, 0.2, -0.2, -0.2], u.m) - ys = u.Quantity([0.2, -0.2, 0.2, -0.2], u.m) - psis = Angle([-60, -45, 0, 45, 60], unit="deg") - - for x, y in zip(xs, ys): - for psi in psis: - - # make a toy shower model - model_gaussian = toymodel.Gaussian( - x=x, y=y, width=width, length=length, psi=psi - ) - model_skewed = toymodel.SkewedGaussian( - x=x, y=y, width=width, length=length, psi=psi, skewness=0.5 - ) - - image, signal, noise = model_gaussian.generate_image( - geom, intensity=intensity, nsb_level_pe=0, rng=rng - ) - - clean_mask = np.array(signal) > 0 - result = image_fit_parameters( - geom, - signal, - n_row=0, - cleaned_mask=clean_mask, - pdf=PDFType("gaussian"), - ) - - if result.is_valid or result.is_accurate: - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) - - assert u.isclose(result.width, width, rtol=0.1) - assert u.isclose(result.length, length, rtol=0.1) - assert (result.psi.to_value(u.deg) == approx(psi.deg, abs=2)) or abs( - result.psi.to_value(u.deg) - psi.deg - ) == approx(180.0, abs=2) - - image, signal, noise = model_skewed.generate_image( - geom, intensity=intensity, nsb_level_pe=0, rng=rng - ) - - clean_mask = np.array(signal) > 0 - result = image_fit_parameters( - geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("skewed") - ) - - if result.is_valid or result.is_accurate: - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) - - assert u.isclose(result.width, width, rtol=0.1) - assert u.isclose(result.length, length, rtol=0.1) - assert (result.psi.to_value(u.deg) == approx(psi.deg, abs=2)) or abs( - result.psi.to_value(u.deg) - psi.deg - ) == approx(180.0, abs=2) - - -def test_with_toy_lst(prod5_lst): - rng = np.random.default_rng(42) - - geom = prod5_lst.camera.geometry - - width = 0.03 * u.m - length = 0.15 * u.m - intensity = 500 - - xs = u.Quantity([0.2, 0.2, -0.2, -0.2], u.m) - ys = u.Quantity([0.2, -0.2, 0.2, -0.2], u.m) - psis = Angle([-60, -45, 0, 45, 60], unit="deg") - - for x, y in zip(xs, ys): - for psi in psis: - - # make a toy shower model - model_gaussian = toymodel.Gaussian( - x=x, y=y, width=width, length=length, psi=psi - ) - model_skewed = toymodel.SkewedGaussian( - x=x, y=y, width=width, length=length, psi=psi, skewness=0.5 - ) - - image, signal, noise = model_gaussian.generate_image( - geom, intensity=intensity, nsb_level_pe=0, rng=rng - ) - - clean_mask = np.array(signal) > 0 - result = image_fit_parameters( - geom, - signal, - n_row=0, - cleaned_mask=clean_mask, - pdf=PDFType("gaussian"), - ) - - if result.is_valid or result.is_accurate: - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) - - assert u.isclose(result.width, width, rtol=0.1) - assert u.isclose(result.length, length, rtol=0.1) - assert (result.psi.to_value(u.deg) == approx(psi.deg, abs=2)) or abs( - result.psi.to_value(u.deg) - psi.deg - ) == approx(180.0, abs=2) - - image, signal, noise = model_skewed.generate_image( - geom, intensity=intensity, nsb_level_pe=0, rng=rng - ) - - clean_mask = np.array(signal) > 0 - result = image_fit_parameters( - geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("skewed") - ) - - if result.is_valid or result.is_accurate: - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) - - assert u.isclose(result.width, width, rtol=0.1) - assert u.isclose(result.length, length, rtol=0.1) - assert (result.psi.to_value(u.deg) == approx(psi.deg, abs=2)) or abs( - result.psi.to_value(u.deg) - psi.deg - ) == approx(180.0, abs=2) - - def test_skewness(prod5_lst): rng = np.random.default_rng(42) geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) - widths = u.Quantity([0.02, 0.03, 0.04], u.m) - lengths = u.Quantity([0.15, 0.2, 0.3], u.m) + widths = u.Quantity([0.04, 0.05, 0.06], u.deg) + lengths = u.Quantity([0.2, 0.3, 0.4], u.deg) intensities = np.array([500, 1500, 2000]) - xs = u.Quantity([0.2, 0.2, -0.2, -0.2], u.m) - ys = u.Quantity([0.2, -0.2, 0.2, -0.2], u.m) + xs = u.Quantity([0.9, 0.9, -0.9, -0.9], u.deg) + ys = u.Quantity([0.9, -0.9, 0.9, -0.9], u.deg) psis = Angle([-60, -45, 0, 45, 60, 90], unit="deg") skews = np.array([0, 0.5, -0.5, 0.8, -0.8, 0.9]) @@ -387,24 +256,28 @@ def test_skewness(prod5_lst): ) image, signal, noise = model_skewed.generate_image( - geom, intensity=intensity, nsb_level_pe=0, rng=rng + geom_nom, intensity=intensity, nsb_level_pe=0, rng=rng ) clean_mask = np.array(signal) > 0 result = image_fit_parameters( - geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("skewed") + geom_nom, + signal, + n_row=0, + cleaned_mask=clean_mask, + pdf=PDFType("skewed"), ) if result.is_valid or result.is_accurate: - assert u.isclose(result.x, x, rtol=0.1) - assert u.isclose(result.y, y, rtol=0.1) + assert u.isclose(result.fov_lon, x, rtol=0.1) + assert u.isclose(result.fov_lat, y, rtol=0.1) assert u.isclose(result.width, width, rtol=0.1) assert u.isclose(result.length, length, rtol=0.1) - psi_same = result.psi.to_value(u.deg) == approx(psi.deg, abs=1) + psi_same = result.psi.to_value(u.deg) == approx(psi.deg, abs=2) psi_opposite = abs(result.psi.to_value(u.deg) - psi.deg) == approx( - 180.0, abs=1 + 180.0, abs=2 ) assert psi_same or psi_opposite @@ -419,27 +292,29 @@ def test_skewness(prod5_lst): def test_gaussian_skewness(prod5_lst): rng = np.random.default_rng(42) geom = prod5_lst.camera.geometry + telescope_frame = TelescopeFrame() + geom_nom = geom.transform_to(telescope_frame) model_gaussian = toymodel.Gaussian( - x=0 * u.m, - y=0 * u.m, - width=0.02 * u.m, - length=0.1 * u.m, + x=1.0 * u.deg, + y=1.0 * u.deg, + width=0.05 * u.deg, + length=0.3 * u.deg, psi=20 * u.deg, ) image, signal, noise = model_gaussian.generate_image( - geom, intensity=1500, nsb_level_pe=0, rng=rng + geom_nom, intensity=1500, nsb_level_pe=0, rng=rng ) clean_mask = np.array(signal) > 0 result = image_fit_parameters( - geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("gaussian") + geom_nom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("gaussian") ) result_skew = image_fit_parameters( - geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("skewed") + geom_nom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("skewed") ) - assert result.skewness == 0 + assert result.skewness == approx(0, abs=0.1) assert result_skew.skewness == approx(0, abs=0.1) @@ -451,10 +326,10 @@ def test_single_pixel(): geom = CameraGeometry( name="testcam", pix_id=np.arange(9), - pix_x=x.ravel() * u.cm, - pix_y=y.ravel() * u.cm, + pix_x=x.ravel() * u.deg, + pix_y=y.ravel() * u.deg, pix_type="rectangular", - pix_area=np.full(9, 1.0) * u.cm**2, + pix_area=np.full(9, 1.0) * u.deg**2, ) image = np.zeros((3, 3)) @@ -483,45 +358,17 @@ def test_reconstruction_in_telescope_frame(prod5_lst): geom = prod5_lst.camera.geometry telescope_frame = TelescopeFrame() - camera_frame = geom.frame geom_nom = geom.transform_to(telescope_frame) - width = 0.03 * u.m - length = 0.15 * u.m + width = 0.04 * u.deg + length = 0.40 * u.deg intensity = 5000 - 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) + xs = u.Quantity([0.9, 0.9, -0.9, -0.9], u.deg) + ys = u.Quantity([0.9, -0.9, 0.9, -0.9], u.deg) psis = Angle([-90, -45, 0, 45, 90], unit="deg") skews = 0.0, 0.2, 0.5 - 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: for skew in skews: @@ -530,7 +377,7 @@ def get_transformed_width(telescope_hillas, telescope_frame, camera_frame): x=x, y=y, width=width, length=length, psi=psi, skewness=skew ) image, signal, noise = model.generate_image( - geom, intensity=intensity, nsb_level_pe=5 + geom_nom, intensity=intensity, nsb_level_pe=5 ) telescope_result = image_fit_parameters( @@ -541,35 +388,23 @@ def get_transformed_width(telescope_hillas, telescope_frame, camera_frame): pdf=PDFType("skewed"), ) - camera_result = image_fit_parameters( - geom, - signal, - n_row=0, - cleaned_mask=(np.array(signal) > 0), - pdf=PDFType("skewed"), - ) - - 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) - if telescope_result.is_valid or telescope_result.is_accurate: - assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01) - assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01) + assert u.isclose(telescope_result.fov_lon, x, rtol=0.02) + assert u.isclose(telescope_result.fov_lat, y, rtol=0.02) + assert u.isclose(telescope_result.length, length, rtol=0.1) + assert u.isclose(telescope_result.width, width, rtol=0.1) - transformed_length = get_transformed_length( - telescope_result, telescope_frame, camera_frame - ) - assert u.isclose( - transformed_length, camera_result.length, rtol=0.01 + psi_same = telescope_result.psi.to_value(u.deg) == approx( + psi.deg, abs=1 ) + psi_opposite = abs( + telescope_result.psi.to_value(u.deg) - psi.deg + ) == approx(180.0, abs=1) + assert psi_same or psi_opposite - transformed_width = get_transformed_width( - telescope_result, telescope_frame, camera_frame - ) - assert u.isclose(transformed_width, camera_result.width, rtol=0.01) + if psi_same: + assert telescope_result.skewness == approx(skew, abs=0.3) + else: + assert telescope_result.skewness == approx(-skew, abs=0.3) + + assert signal.sum() == telescope_result.intensity