diff --git a/docs/changes/3003.feature.rst b/docs/changes/3003.feature.rst new file mode 100644 index 00000000000..5100dc14872 --- /dev/null +++ b/docs/changes/3003.feature.rst @@ -0,0 +1,5 @@ +Added calculation and storage of ``true_disp_norm`` and ``true_disp_sign`` per-telescope +when running ``ctapipe-process``. These values are computed in ``ImageProcessor`` +(telescope frame only) and written to ``simulation/event/telescope/parameters``. +A new ``TrueDispContainer`` is added to ``SimulatedCameraContainer``, and +``HDF5EventSource`` reads these columns back with NaN defaults for older files. diff --git a/docs/conf.py b/docs/conf.py index 6e6d3f8cfd1..329bff3cdb8 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -131,6 +131,8 @@ def add_reference_type(prefix, objs): "traitlets.traitlets.ObserveHandler", "traitlets.traitlets.T", "traitlets.traitlets.G", + "traitlets.traitlets.K", + "traitlets.traitlets.V", "Sentinel", "ObserveHandler", "dict[K, V]", diff --git a/docs/user-guide/data_format/index.rst b/docs/user-guide/data_format/index.rst index ba0c56a6679..f8dea13b606 100644 --- a/docs/user-guide/data_format/index.rst +++ b/docs/user-guide/data_format/index.rst @@ -128,8 +128,10 @@ Simulation Data Format - simulated camera images - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.SimulatedCameraContainer` * - ``/simulation/event/telescope/parameters/tel_{TEL_ID:03d}`` - - Parameters derived form the simulated camera images - - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ImageParametersContainer` + - Parameters derived from the simulated camera images, including true Hillas parameters + and true disp values (:py:class:`~ctapipe.containers.TrueDispContainer`) for training + and evaluating angular reconstruction algorithms + - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ImageParametersContainer`, :py:class:`~ctapipe.containers.TrueDispContainer` * - ``/simulation/service/shower_distribution`` - simulated shower distribution histograms - :py:class:`~ctapipe.containers.SimulatedShowerDistribution` diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 1d49146ce13..e8c8153b33f 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -50,6 +50,7 @@ "ReconstructedEnergyContainer", "ReconstructedGeometryContainer", "DispContainer", + "TrueDispContainer", "SimulatedCameraContainer", "SimulatedShowerContainer", "SimulatedShowerDistribution", @@ -764,6 +765,26 @@ class TelescopeImpactParameterContainer(Container): # Simulation containers +class TrueDispContainer(Container): + """ + Container for the true (simulated) disp values, used for debugging angular + reconstruction and for training the DispReconstructor. + """ + + default_prefix = "true_disp" + + norm = Field( + nan * u.deg, + "True norm of the disp parameter (distance from shower COG to true source position)", + unit=u.deg, + ) + sign = Field( + np.float32(nan), + "True sign of the disp parameter (+1, -1, or 0 in degenerate case where " + "true source position coincides exactly with the shower CoG)", + ) + + class SimulatedShowerContainer(Container): default_prefix = "true" energy = Field(nan * u.TeV, "Simulated Energy", unit=u.TeV) @@ -811,6 +832,12 @@ class SimulatedCameraContainer(Container): type=ImageParametersContainer, ) + true_disp = Field( + default_factory=TrueDispContainer, + description="True disp parameters (norm and sign) derived from the true shower direction", + type=TrueDispContainer, + ) + impact = Field( default_factory=TelescopeImpactParameterContainer, description="true impact parameter", diff --git a/src/ctapipe/image/image_processor.py b/src/ctapipe/image/image_processor.py index 0483bac7049..fe2b0166ff1 100644 --- a/src/ctapipe/image/image_processor.py +++ b/src/ctapipe/image/image_processor.py @@ -2,11 +2,13 @@ High level image processing (ImageProcessor Component) """ +import warnings from copy import deepcopy +import astropy.units as u import numpy as np -from ctapipe.coordinates import TelescopeFrame +from ctapipe.coordinates import MissingFrameAttributeWarning, TelescopeFrame from ..containers import ( ArrayEventContainer, @@ -16,6 +18,7 @@ IntensityStatisticsContainer, PeakTimeStatisticsContainer, TimingParametersContainer, + TrueDispContainer, ) from ..core import QualityQuery, TelescopeComponent from ..core.traits import Bool, BoolTelescopeParameter, ComponentName, List @@ -213,6 +216,72 @@ def _parameterize_image( # parameterization return default + def _calculate_true_disp(self, event, sim_camera): + """ + Calculate the true disp parameters (norm and sign) for a simulated telescope event. + + This method should only be called when ``use_telescope_frame=True``, as it + requires the hillas parameters in the telescope frame (fov_lon, fov_lat, psi). + + Parameters + ---------- + event : ArrayEventContainer + The event container with simulation and monitoring data. + sim_camera : SimulatedCameraContainer + The simulated camera container with true_parameters already set. + + Returns + ------- + TrueDispContainer + Container with true_disp_norm and true_disp_sign values, or default + NaN-filled container if the calculation is not possible (missing shower + info or pointing data). + """ + from astropy.coordinates import AltAz + + shower = event.simulation.shower + if shower is None: + return TrueDispContainer() + + pointing_alt = event.monitoring.pointing.array_altitude + pointing_az = event.monitoring.pointing.array_azimuth + + if np.isnan(pointing_alt.value) or np.isnan(pointing_az.value): + return TrueDispContainer() + + # Clamp altitude to [-90, 90] deg to avoid floating-point rounding errors + # that can push values slightly outside the valid range for astropy AltAz + pointing_alt = np.clip(pointing_alt.to_value(u.deg), -90.0, 90.0) * u.deg + shower_alt = np.clip(shower.alt.to_value(u.deg), -90.0, 90.0) * u.deg + + hillas = sim_camera.true_parameters.hillas + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", MissingFrameAttributeWarning) + horizontal_coord = AltAz(alt=shower_alt, az=shower.az) + pointing = AltAz(alt=pointing_alt, az=pointing_az) + tel_frame = TelescopeFrame(telescope_pointing=pointing) + tel_coord = horizontal_coord.transform_to(tel_frame) + + fov_lon = tel_coord.fov_lon.to(u.deg) + fov_lat = tel_coord.fov_lat.to(u.deg) + + psi = hillas.psi.to_value(u.rad) + cog_lon = hillas.fov_lon.to(u.deg) + cog_lat = hillas.fov_lat.to(u.deg) + + delta_lon = fov_lon - cog_lon + delta_lat = fov_lat - cog_lat + + # Project displacement vector onto the shower major axis (defined by psi) + # to determine on which end of the shower axis the true source lies + true_disp_projected = np.cos(psi) * delta_lon + np.sin(psi) * delta_lat + true_sign = np.float32(np.sign(true_disp_projected.value)) + # Norm is the full Euclidean distance from COG to true source position + true_norm = np.sqrt(delta_lon**2 + delta_lat**2) + + return TrueDispContainer(norm=true_norm, sign=true_sign) + def _process_telescope_event(self, event): """ Loop over telescopes and process the calibrated images into parameters @@ -255,6 +324,9 @@ def _process_telescope_event(self, event): if not container.prefix.startswith("true_"): container.prefix = f"true_{container.prefix}" + if self.use_telescope_frame: + sim_camera.true_disp = self._calculate_true_disp(event, sim_camera) + self.log.debug( "sim params: %s", event.simulation.tel[tel_id].true_parameters.as_dict( diff --git a/src/ctapipe/image/tests/test_image_processor.py b/src/ctapipe/image/tests/test_image_processor.py index 7cb74bc5130..5f6d3db75f0 100644 --- a/src/ctapipe/image/tests/test_image_processor.py +++ b/src/ctapipe/image/tests/test_image_processor.py @@ -7,14 +7,99 @@ import astropy.units as u import numpy as np import pytest +from astropy.coordinates import EarthLocation from numpy import isfinite from ctapipe.calib import CameraCalibrator from ctapipe.containers import ( + ArrayEventContainer, CameraHillasParametersContainer, CameraTimingParametersContainer, + DL1CameraContainer, + SimulatedCameraContainer, + SimulatedEventContainer, + SimulatedShowerContainer, ) +from ctapipe.coordinates import CameraFrame from ctapipe.image import ImageCleaner, ImageProcessor +from ctapipe.instrument import CameraGeometry, SubarrayDescription, TelescopeDescription +from ctapipe.instrument.camera import CameraDescription, CameraReadout +from ctapipe.instrument.optics import OpticsDescription, ReflectorShape, SizeType + + +@pytest.fixture(scope="module") +def synthetic_subarray(): + """ + Create a fully synthetic subarray (no network access needed). + Uses a rectangular camera with a 1m focal length. + """ + focal_length = 1.0 * u.m + geom = CameraGeometry.make_rectangular(10, 10) + geom.frame = CameraFrame(focal_length=focal_length) + n_pix = geom.n_pixels + + readout = CameraReadout( + name="TestCam", + sampling_rate=1 * u.GHz, + reference_pulse_shape=np.ones((1, 3)), + reference_pulse_sample_width=1 * u.ns, + n_channels=1, + n_pixels=n_pix, + n_samples=50, + ) + cam_desc = CameraDescription("TestCam", geom, readout) + optics = OpticsDescription( + name="testoptics", + size_type=SizeType.MST, + n_mirrors=1, + n_mirror_tiles=1, + mirror_area=1.0 * u.m**2, + equivalent_focal_length=focal_length, + effective_focal_length=focal_length, + reflector_shape=ReflectorShape.PARABOLIC, + ) + tel_desc = TelescopeDescription(name="TestTel", optics=optics, camera=cam_desc) + reference_location = EarthLocation( + lon=-17 * u.deg, lat=28 * u.deg, height=2200 * u.m + ) + return SubarrayDescription( + name="test", + tel_positions={1: np.array([0.0, 0.0, 0.0]) * u.m}, + tel_descriptions={1: tel_desc}, + reference_location=reference_location, + ) + + +@pytest.fixture +def synthetic_event_with_simulation(synthetic_subarray): + """ + Create a synthetic ArrayEventContainer with simulation data for tel_id=1. + The image has a simple signal in the central pixels so that Hillas parameters + can be computed. + """ + n_pix = synthetic_subarray.tel[1].camera.geometry.n_pixels + image = np.zeros(n_pix) + # Put signal in the central pixels (rows 3-7, cols 3-7 of the 10x10 grid) + # Signal falls off in one direction to give a clear major axis + for i in range(5): + image[30 + i * 10 + 3 : 30 + i * 10 + 7] = 100.0 * np.exp(-((i - 2) ** 2) / 2) + + true_image = (image * 0.5).astype(np.int32) + + event = ArrayEventContainer() + event.monitoring.pointing.array_altitude = 70.0 * u.deg + event.monitoring.pointing.array_azimuth = 0.0 * u.deg + + event.simulation = SimulatedEventContainer() + event.simulation.shower = SimulatedShowerContainer( + alt=70.5 * u.deg, + az=0.5 * u.deg, + ) + + event.dl1.tel[1] = DL1CameraContainer(image=image, peak_time=None) + event.monitoring.tel[1] # initialize monitoring + event.simulation.tel[1] = SimulatedCameraContainer(true_image=true_image) + return event @pytest.mark.parametrize("cleaner", ImageCleaner.non_abstract_subclasses().values()) @@ -94,3 +179,48 @@ def test_image_processor_camera_frame(cleaner, example_event, example_subarray): assert isinstance(dl1.parameters.timing, CameraTimingParametersContainer) assert np.isnan(dl1.parameters.hillas.length.value) assert dl1.parameters.hillas.length.unit == u.m + + +def test_true_disp_calculation(synthetic_subarray, synthetic_event_with_simulation): + """Test that true_disp_norm and true_disp_sign are calculated for simulation events.""" + process_images = ImageProcessor(subarray=synthetic_subarray) + process_images(synthetic_event_with_simulation) + + sim_camera = synthetic_event_with_simulation.simulation.tel[1] + true_disp = sim_camera.true_disp + + # Hillas parameters should have been computed for the non-trivial image + if isfinite(sim_camera.true_parameters.hillas.fov_lon.value): + assert isfinite(true_disp.norm.value) + assert true_disp.norm >= 0 * u.deg + assert true_disp.norm.unit == u.deg + assert true_disp.sign in {-1.0, 0.0, 1.0} + + +def test_true_disp_requires_telescope_frame( + synthetic_subarray, synthetic_event_with_simulation +): + """Test that true_disp remains NaN when use_telescope_frame=False.""" + process_images = ImageProcessor( + subarray=synthetic_subarray, use_telescope_frame=False + ) + process_images(synthetic_event_with_simulation) + + sim_camera = synthetic_event_with_simulation.simulation.tel[1] + # When not in telescope frame, true_disp should remain at defaults (NaN) + assert np.isnan(sim_camera.true_disp.norm.value) + assert np.isnan(sim_camera.true_disp.sign) + + +def test_true_disp_no_simulation(synthetic_subarray, synthetic_event_with_simulation): + """Test that _calculate_true_disp handles missing simulation shower gracefully.""" + event = deepcopy(synthetic_event_with_simulation) + event.simulation.shower = None + + process_images = ImageProcessor(subarray=synthetic_subarray) + process_images(event) + + sim_camera = event.simulation.tel[1] + # Without shower, true_disp should remain at defaults (NaN) + assert np.isnan(sim_camera.true_disp.norm.value) + assert np.isnan(sim_camera.true_disp.sign) diff --git a/src/ctapipe/io/datawriter.py b/src/ctapipe/io/datawriter.py index 920cafa4177..13066277649 100644 --- a/src/ctapipe/io/datawriter.py +++ b/src/ctapipe/io/datawriter.py @@ -662,6 +662,7 @@ def _write_dl1_telescope_events(self, event: ArrayEventContainer): true_parameters.concentration, true_parameters.morphology, true_parameters.intensity_statistics, + event.simulation.tel[tel_id].true_disp, ], ) diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index 9674c87db70..f5a76eb1cbb 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -44,6 +44,7 @@ TelEventIndexContainer, TimingParametersContainer, TriggerContainer, + TrueDispContainer, ) from ..core import Container, Field, Provenance from ..core.traits import UseEnum @@ -663,6 +664,7 @@ def _init_dl1_true_parameter_readers(self): ConcentrationContainer, MorphologyContainer, IntensityStatisticsContainer, + TrueDispContainer, ], prefixes=[ f"true_{hillas_prefix}", @@ -670,6 +672,7 @@ def _init_dl1_true_parameter_readers(self): "true_concentration", "true_morphology", "true_intensity", + "true_disp", ], ) for table in self.file_.root.dl1.event.telescope.parameters @@ -1002,6 +1005,7 @@ def _fill_parameters_for_tel( morphology=simulated_params[3], intensity_statistics=simulated_params[4], ) + simulated.true_disp = simulated_params[5] def _fill_muons_for_tel(self, data, tel_id, key, muon_readers): if not self.has_muon_parameters: