Skip to content

Commit

Permalink
Switching to only test fitting, not cleaning
Browse files Browse the repository at this point in the history
  • Loading branch information
Christoph Toennis committed Oct 4, 2024
1 parent 666e1d1 commit 885bfcf
Show file tree
Hide file tree
Showing 4 changed files with 95 additions and 45 deletions.
107 changes: 70 additions & 37 deletions src/ctapipe/calib/camera/pointing.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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):
"""
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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(
Expand Down
27 changes: 23 additions & 4 deletions src/ctapipe/calib/camera/tests/test_pointing.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -35,7 +37,7 @@

# some wavelength

obswl = 0.42 * u.micron
obswl = 0.35 * u.micron

# pointing to the crab nebula

Expand All @@ -44,6 +46,8 @@
alt = []
az = []

# get the local pointing values

for t in times:
contemporary_crab = crab.transform_to(
AltAz(
Expand All @@ -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()
3 changes: 1 addition & 2 deletions src/ctapipe/containers.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
3 changes: 1 addition & 2 deletions src/ctapipe/monitoring/interpolation.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,6 @@

from ctapipe.core import Component, traits

from .astropy_helpers import read_table


class ChunkFunction:

Expand Down Expand Up @@ -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,
Expand Down

0 comments on commit 885bfcf

Please sign in to comment.