Merge pull request #2715 from cta-observatory/stats_tool_shape_bugfix

Processing gain-selected dl1 data with the pixel stats tool

Author: maxnoe

docs/changes/2715.bugfix.rst

@@ -0,0 +1 @@
+Fix processing gain-selected dl1 data with the ``PixelStatisticsCalculatorTool``.

src/ctapipe/monitoring/outlier.py

@@ -65,6 +65,12 @@ class RangeOutlierDetector(OutlierDetector):
     ).tag(config=True)
 
     def __call__(self, column):
+        # Validate that the shape of the statistic values has three dimensions
+        if column.ndim != 3:
+            raise ValueError(
+                f"Invalid shape of the column '{column.name}': '{column.shape}'. "
+                "Expected the statistic values of shape (n_entries, n_channels, n_pixels)."
+            )
         # Remove outliers is statistical values out a given range
         outliers = np.logical_or(
             column < self.validity_range[0],
@@ -93,6 +99,12 @@ class MedianOutlierDetector(OutlierDetector):
     ).tag(config=True)
 
     def __call__(self, column):
+        # Validate that the shape of the statistic values has three dimensions
+        if column.ndim != 3:
+            raise ValueError(
+                f"Invalid shape of the column '{column.name}': '{column.shape}'. "
+                "Expected the statistic values of shape (n_entries, n_channels, n_pixels)."
+            )
         # Camera median
         camera_median = np.ma.median(column, axis=2)
         # Detect outliers based on the deviation of the median distribution
@@ -124,6 +136,12 @@ class StdOutlierDetector(OutlierDetector):
     ).tag(config=True)
 
     def __call__(self, column):
+        # Validate that the shape of the statistic values has three dimensions
+        if column.ndim != 3:
+            raise ValueError(
+                f"Invalid shape of the column '{column.name}': '{column.shape}'. "
+                "Expected the statistic values of shape (n_entries, n_channels, n_pixels)."
+            )
         # Camera median
         camera_median = np.ma.median(column, axis=2)
         # Camera std

src/ctapipe/monitoring/tests/test_outlier.py

@@ -3,10 +3,12 @@
 """
 
 import numpy as np
+import pytest
 from astropy.table import Table
 
 from ctapipe.monitoring.outlier import (
     MedianOutlierDetector,
+    OutlierDetector,
     RangeOutlierDetector,
     StdOutlierDetector,
 )
@@ -102,3 +104,28 @@ def test_std_detection(example_subarray):
     # Check if outliers where detected correctly
     np.testing.assert_array_equal(ff_outliers, ff_expected_outliers)
     np.testing.assert_array_equal(ped_outliers, ped_expected_outliers)
+
+
+def test_check_for_column_shape(example_subarray):
+    """test the check for the shape of the column"""
+
+    # Create dummy data for testing
+    rng = np.random.default_rng(0)
+    # Distribution mimics the median values of gain-selected charge images of flat-field events
+    median = rng.normal(77.0, 0.6, size=(50, 1855))
+    # Create astropy table
+    table = Table([median], names=("median",))
+    # Initialize the outlier detector based on the deviation from the camera median
+    outlier_detectors = [
+        "RangeOutlierDetector",
+        "MedianOutlierDetector",
+        "StdOutlierDetector",
+    ]
+    for outlier_detector_name in outlier_detectors:
+        outlier_detector = OutlierDetector.from_name(
+            outlier_detector_name,
+            subarray=example_subarray,
+        )
+        # Check if ValueError is raised when the provided column shape is invalid
+        with pytest.raises(ValueError, match="Invalid shape of the column"):
+            outlier_detector(table["median"])

src/ctapipe/tools/calculate_pixel_stats.py

@@ -82,6 +82,8 @@ class PixelStatisticsCalculatorTool(Tool):
         TableLoader,
     ] + classes_with_traits(PixelStatisticsCalculator)
 
+    DL1_COLUMN_NAMES = ["image", "peak_time"]
+
     def setup(self):
         # Read the input data with the 'TableLoader'
         self.input_data = self.enter_context(
@@ -137,6 +139,10 @@ def start(self):
                     f"Column '{self.input_column_name}' not found "
                     f"in the input data for telescope 'tel_id={tel_id}'."
                 )
+            # Check if the dl1 data is gain selected and add an extra dimension for n_channels
+            for col_name in self.DL1_COLUMN_NAMES:
+                if col_name in dl1_table.colnames and dl1_table[col_name].ndim == 2:
+                    dl1_table[col_name] = dl1_table[col_name][:, np.newaxis]
             # Perform the first pass of the statistics calculation
             aggregated_stats = self.stats_calculator.first_pass(
                 table=dl1_table,

src/ctapipe/tools/tests/test_calculate_pixel_stats.py

@@ -29,6 +29,13 @@ def test_calculate_pixel_stats_tool(tmp_path, dl1_image_file):
                 "stats_aggregator_type": [
                     ("id", tel_id, "PlainAggregator"),
                 ],
+                "outlier_detector_list": [
+                    {
+                        "apply_to": "mean",
+                        "name": "MedianOutlierDetector",
+                        "config": {"median_range_factors": [-2.0, 2.0]},
+                    }
+                ],
             },
             "PlainAggregator": {
                 "chunk_size": 1,
@@ -50,13 +57,13 @@ def test_calculate_pixel_stats_tool(tmp_path, dl1_image_file):
     )
     # Check that the output file has been created
     assert monitoring_file.exists()
-    # Check that the output file is not empty
+    # Check if the shape of the aggregated statistic values has three dimension
     assert (
         read_table(
             monitoring_file,
             path=f"/dl1/monitoring/telescope/statistics/tel_{tel_id:03d}",
-        )["mean"]
-        is not None
+        )["mean"].ndim
+        == 3
     )
     # Read subarray description from the created monitoring file
     subarray = SubarrayDescription.from_hdf(monitoring_file)