Hillas parameters in telescope frame

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,15 +333,13 @@ class DL1CameraContainer(Container):
         dtype=np.float32,
         ndim=1,
     )
-
     image_mask = Field(
         None,
         "Boolean numpy array where True means the pixel has passed cleaning."
         " Shape: (n_pixel, )",
         dtype=np.bool_,
         ndim=1,
     )
-
     parameters = Field(None, "Image parameters", type=ImageParametersContainer)
 
 

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

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,

ctapipe/image/image_processor.py

@@ -1,15 +1,15 @@
 """
 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
@@ -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(
@@ -64,6 +65,7 @@ def __init__(
         self,
         subarray: SubarrayDescription,
         is_simulation,
+        use_telescope_frame=True,
         config=None,
         parent=None,
         **kwargs,
@@ -93,6 +95,15 @@ def __init__(
         )
         self.check_image = ImageQualityQuery(parent=self)
         self._is_simulation = is_simulation
+        self._use_telescope_frame = use_telescope_frame
+        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)
@@ -102,11 +113,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
@@ -117,15 +128,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:
@@ -184,9 +192,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 = dl1_camera.geometry
             # compute image parameters only if requested to write them
             dl1_camera.image_mask = self.clean(
                 tel_id=tel_id,
@@ -199,6 +211,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))
@@ -212,6 +225,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,
                 )