From 885bfcf1e1d55c7bf66abc59dc3617d06b5aa1a8 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Fri, 4 Oct 2024 16:28:40 +0200 Subject: [PATCH] Switching to only test fitting, not cleaning --- src/ctapipe/calib/camera/pointing.py | 107 ++++++++++++------ .../calib/camera/tests/test_pointing.py | 27 ++++- src/ctapipe/containers.py | 3 +- src/ctapipe/monitoring/interpolation.py | 3 +- 4 files changed, 95 insertions(+), 45 deletions(-) diff --git a/src/ctapipe/calib/camera/pointing.py b/src/ctapipe/calib/camera/pointing.py index 52acea1204a..283ee919a69 100644 --- a/src/ctapipe/calib/camera/pointing.py +++ b/src/ctapipe/calib/camera/pointing.py @@ -146,7 +146,7 @@ def StarImageGenerator( :param dict psf_model_pars: psf model parameters """ camera_geometry = CameraGeometry.from_name(camera_name) - psf = PSFModel.from_name(self.psf_model_type, subarray=self.subarray, parent=self) + psf = PSFModel.from_name(self.psf_model_type) psf.update_model_parameters(psf_model_pars) image = np.zeros(len(camera_geometry)) for r, p, m in zip(radius, phi, magnitude): @@ -353,6 +353,13 @@ def get_expected_star_pixels(self, t, focal_correction=0.0): return res + def get_star_mask(self, t, focal_correction=0.0): + mask = np.full(len(self.camera_geometry), True) + pixels = self.get_expected_star_pixels(t, focal_correction=focal_correction) + mask[pixels] = False + + return mask + class PointingCalculator(TelescopeComponent): """ @@ -418,6 +425,10 @@ class PointingCalculator(TelescopeComponent): focal_length = Float(1.0, help="Camera focal length in meters").tag(config=True) + camera_name = Unicode("LSTCam", help="Name of the camera to be calibrated").tag( + config=True + ) + def __init__( self, subarray, @@ -447,7 +458,7 @@ def __init__( self.stats_aggregator, subarray=self.subarray, parent=self ) - self.set_camera(geometry) + self.set_camera(self.camera_name) def set_camera(self, geometry, focal_lengh): if isinstance(geometry, str): @@ -477,24 +488,6 @@ def ingest_data(self, data_table): self.broken_pixels = np.unique(np.where(data_table["unusable_pixels"])) - for azimuth, altitude, time in zip( - data_table["telescope_pointing_azimuth"], - data_table["telescope_pointing_altitude"], - data_table["time"], - ): - _pointing = SkyCoord( - az=azimuth, - alt=altitude, - frame="altaz", - obstime=time, - location=self.location, - obswl=self.observed_wavelength * u.micron, - relative_humidity=self.meteo_parameters["relative_humidity"], - temperature=self.meteo_parameters["temperature"] * u.deg_C, - pressure=self.meteo_parameters["pressure"] * u.hPa, - ) - self.pointing.append(_pointing) - self.image_size = len( data_table["variance_images"][0].image ) # get the size of images of the camera we are calibrating @@ -635,11 +628,9 @@ def _get_accumulated_images(self, data_table): shower_mask = copy.deepcopy(light_mask) - star_pixels = [ - self.tracer.get_expected_star_pixels(t) for t in data_table["time"] - ] + star_mask = [self.tracer.get_star_mask(t) for t in data_table["time"]] - light_mask[:, star_pixels] = True + light_mask[~star_mask] = True if self.broken_pixels is not None: light_mask[:, self.broken_pixels] = True @@ -672,19 +663,11 @@ def _get_accumulated_images(self, data_table): variance_image_table, col_name="image" ) - self.accumulated_times = np.array( - [x.validity_start for x in variance_statistics] - ) - - accumulated_images = np.array([x.mean for x in variance_statistics]) + self.accumulated_times = variance_statistics["time_start"] - star_pixels = [ - self.tracer.get_expected_star_pixels(t) for t in data_table["time"] - ] + accumulated_images = variance_statistics["mean"] - star_mask = np.ones(self.image_size, dtype=bool) - - star_mask[star_pixels] = False + star_mask = [self.tracer.get_star_mask(t) for t in data_table["time"]] # get NSB values @@ -693,6 +676,56 @@ def _get_accumulated_images(self, data_table): self.clean_images = np.array([x - y for x, y in zip(accumulated_images, nsb)]) + def make_mispointed_data(self, p, pointing_table): + self.tracer = StarTracer.from_lookup( + pointing_table["telescope_pointing_azimuth"], + pointing_table["telescope_pointing_altitude"], + pointing_table["time"], + self.meteo_parameters, + self.observed_wavelength, + self.camera_geometry, + self.focal_length, + self.telescope_location, + ) + + stars = self.tracer.get_star_labels() + + self.accumulated_times = pointing_table["time"] + + self.all_containers = [] + + for t in pointing_table["time"]: + self.all_containers.append([]) + + for star in stars: + x, y = self.tracer.get_position_in_camera(t, star) + r, phi = cart2pol(x, y) + x_mispoint, y_mispoint = self.tracer.get_position_in_camera( + t, star, offset=p + ) + r_mispoint, phi_mispoint = cart2pol(x_mispoint, y_mispoint) + container = StarContainer( + label=star, + magnitude=self.tracer.get_magnitude(star), + expected_x=x, + expected_y=y, + expected_r=r * u.m, + expected_phi=phi * u.rad, + timestamp=t, + reco_x=x_mispoint, + reco_y=y_mispoint, + reco_r=r_mispoint * u.m, + reco_phi=phi_mispoint * u.rad, + reco_dx=0.05 * u.m, + reco_dy=0.05 * u.m, + reco_dr=0.05 * u.m, + reco_dphi=0.15 * u.rad, + ) + + self.all_containers[-1].append(container) + + return self.all_containers + def fit_function(self, p, t): """ Construct the fit function for the pointing correction @@ -715,7 +748,7 @@ def fit_function(self, p, t): return coord_list - def fit(self): + def fit(self, init_mispointing=(0, 0)): """ Performs the ODR fit of stars trajectories and saves the results as self.fit_results @@ -734,7 +767,7 @@ def fit(self): errors.append(star.reco_dx) errors.append(star.reco_dy) - rdata = RealData(x=[t], y=data, sy=self.errors) + rdata = RealData(x=[t], y=data, sy=errors) odr = ODR(rdata, self.fit_function, beta0=init_mispointing) fit_summary = odr.run() fit_results = pd.DataFrame( diff --git a/src/ctapipe/calib/camera/tests/test_pointing.py b/src/ctapipe/calib/camera/tests/test_pointing.py index 1540e9ce3b1..22b1f446f6b 100644 --- a/src/ctapipe/calib/camera/tests/test_pointing.py +++ b/src/ctapipe/calib/camera/tests/test_pointing.py @@ -4,9 +4,11 @@ import astropy.units as u from astropy.coordinates import AltAz, EarthLocation, SkyCoord +from astropy.table import QTable from astropy.time import Time -from ctapipe.calib.camera.pointing import StarTracer +from ctapipe.calib.camera.pointing import PointingCalculator +from ctapipe.instrument.camera.geometry import CameraGeometry # let's prepare a StarTracer to make realistic images # we need a maximum magnitude @@ -35,7 +37,7 @@ # some wavelength -obswl = 0.42 * u.micron +obswl = 0.35 * u.micron # pointing to the crab nebula @@ -44,6 +46,8 @@ alt = [] az = [] +# get the local pointing values + for t in times: contemporary_crab = crab.transform_to( AltAz( @@ -57,8 +61,23 @@ alt.append(contemporary_crab.alt.to_value()) az.append(contemporary_crab.az.to_value()) +# next i make the pointing table for the fake data generator + +pointing_table = QTable( + [alt, az, times], + names=["telescope_pointing_altitude", "telescope_pointing_azimuth", "time"], +) + # the LST geometry -# the LST focal length +geometry = CameraGeometry.from_name(name="LSTCam") + +# now set up the PointingCalculator + +pc = PointingCalculator(geometry) + +# now make fake mispointed data + +pc.make_mispointed_data((1.0, -1.0), pointing_table) -st = StarTracer.from_lookup(max_magnitude, az, alt, times, meteo_parameters, obswl) +pc.fit() diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index a706d50627d..ad4e3199a73 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -425,8 +425,7 @@ class MorphologyContainer(Container): class StatisticsContainer(Container): """Store descriptive statistics of a chunk of images""" - extraction_start = Field(np.float32(nan), "start of the extraction chunk") - extraction_stop = Field(np.float32(nan), "stop of the extraction chunk") + n_events = Field(-1, "number of events used for the extraction of the statistics") mean = Field( None, "mean of a pixel-wise quantity for each channel" diff --git a/src/ctapipe/monitoring/interpolation.py b/src/ctapipe/monitoring/interpolation.py index fda145a12a7..7e216d487b7 100644 --- a/src/ctapipe/monitoring/interpolation.py +++ b/src/ctapipe/monitoring/interpolation.py @@ -9,8 +9,6 @@ from ctapipe.core import Component, traits -from .astropy_helpers import read_table - class ChunkFunction: @@ -175,6 +173,7 @@ def _check_interpolators(self, tel_id): def _read_parameter_table(self, tel_id): # prevent circular import between io and monitoring + from ..io import read_table input_table = read_table( self.h5file,