Add HiLo correction (#303)

* Add linearity range for HiLo correction (https://arxiv.org/abs/2311.11631)

* Initial commit of `HiLoComponent`

* Initial commit of `HiLoNectarCAMCalibrationTool`

* update `__init__` files

* add example file

* remove `import logging` since we can use `self.logging` for any `Tool` or `Component`

* update logger format

* add unit tests

* We assume the `$NECTARCAMDATA` environment variable is already set up in one's configuration

* Update example script so that it can be run stand-alone. Made sure the `extractor_kwargs` for the charge are the same for the `gain_tool` and `hilo_tool`.

* Add a description to the tool

---------

Co-authored-by: Jean-Philippe Lenain <jlenain@in2p3.fr>

Author: pcorcam

src/nectarchain/makers/calibration/__init__.py

@@ -5,6 +5,7 @@
     FlatFieldSPEHHVStdNectarCAMCalibrationTool,
     FlatFieldSPENominalNectarCAMCalibrationTool,
     FlatFieldSPENominalStdNectarCAMCalibrationTool,
+    HiLoNectarCAMCalibrationTool,
     PhotoStatisticNectarCAMCalibrationTool,
 )
 from .pedestal_makers import PedestalNectarCAMCalibrationTool
@@ -18,4 +19,5 @@
     "FlatFieldSPENominalStdNectarCAMCalibrationTool",
     "PedestalNectarCAMCalibrationTool",
     "PhotoStatisticNectarCAMCalibrationTool",
+    "HiLoNectarCAMCalibrationTool",
 ]

src/nectarchain/makers/calibration/gain/__init__.py

@@ -5,6 +5,7 @@
     FlatFieldSPENominalNectarCAMCalibrationTool,
     FlatFieldSPENominalStdNectarCAMCalibrationTool,
 )
+from .hilo_makers import HiLoNectarCAMCalibrationTool
 from .photostat_makers import PhotoStatisticNectarCAMCalibrationTool
 
 # from .white_target_spe_makers import *
@@ -16,4 +17,5 @@
     "FlatFieldSPEHHVStdNectarCAMCalibrationTool",
     "FlatFieldSPECombinedStdNectarCAMCalibrationTool",
     "PhotoStatisticNectarCAMCalibrationTool",
+    "HiLoNectarCAMCalibrationTool",
 ]

src/nectarchain/makers/calibration/gain/hilo_makers.py

@@ -0,0 +1,106 @@
+import numpy as np
+from ctapipe.core.traits import ComponentNameList
+
+from nectarchain.makers.component import NectarCAMComponent
+
+from ....data.container import (
+    ChargesContainer,
+    ChargesContainers,
+    merge_map_ArrayDataContainer,
+)
+from ....data.management import DataManagement
+from ...component import ArrayDataComponent
+from ...extractor.utils import CtapipeExtractor
+from .core import GainNectarCAMCalibrationTool
+
+__all__ = ["HiLoNectarCAMCalibrationTool"]
+
+
+class HiLoNectarCAMCalibrationTool(GainNectarCAMCalibrationTool):
+    name = "HiLoNectarCAMCalibrationTool"
+    description = (
+        "Calibrate the gain of the low-gain channel. "
+        "Do so by computing the HiLo ratio in the regime where the gain is linear."
+    )
+
+    componentsList = ComponentNameList(
+        NectarCAMComponent,
+        default_value=["HiLoComponent"],
+        help="List of Component names to be apply, the order will be respected",
+    ).tag(config=True)
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def _init_output_path(self):
+        if self.gain_file is not None:
+            self.output_path = self.gain_file.with_name(
+                f"{self.gain_file.stem}_hilo_corrected{self.gain_file.suffix}"
+            )
+        else:
+            # The HiLoComponent will raise an error if no gain_file is provided
+            super()._init_output_path()
+
+    def start(self, n_events=np.inf, *args, **kwargs):
+        str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
+            method=self.method,
+            extractor_kwargs=self.extractor_kwargs,
+        )
+        try:
+            files = DataManagement.find_charges(
+                run_number=self.run_number,
+                method=self.method,
+                str_extractor_kwargs=str_extractor_kwargs,
+                max_events=self.max_events,
+            )
+        except Exception as e:
+            self.log.warning(e)
+            files = []
+        if self.reload_events or len(files) != 1:
+            if len(files) != 1:
+                self.log.info(
+                    f"{len(files)} computed charges files found with max_events >"
+                    f"{self.max_events} for run {self.run_number} with extraction"
+                    f"method {self.method} and {str_extractor_kwargs},\n reload"
+                    f"charges from event loop"
+                )
+            super().start(
+                n_events=n_events,
+                restart_from_begining=False,
+                *args,
+                **kwargs,
+            )
+        else:
+            self.log.info(f"reading computed charge from files {files[0]}")
+            chargesContainers = ChargesContainers.from_hdf5(files[0])
+            if isinstance(chargesContainers, ChargesContainer):
+                self.components[0]._chargesContainers = chargesContainers
+            else:
+                n_slices = 0
+                try:
+                    while True:
+                        next(chargesContainers)
+                        n_slices += 1
+                except StopIteration:
+                    pass
+                chargesContainers = ChargesContainers.from_hdf5(files[0])
+                if n_slices == 1:
+                    self.log.info("merging along TriggerType")
+                    self.components[
+                        0
+                    ]._chargesContainers = merge_map_ArrayDataContainer(
+                        next(chargesContainers)
+                    )
+                else:
+                    self.log.info("merging along slices")
+                    chargesContainers_merged_along_slices = (
+                        ArrayDataComponent.merge_along_slices(
+                            containers_generator=chargesContainers
+                        )
+                    )
+                    self.log.info("merging along TriggerType")
+                    self.components[
+                        0
+                    ]._chargesContainers = merge_map_ArrayDataContainer(
+                        chargesContainers_merged_along_slices
+                    )

src/nectarchain/makers/calibration/gain/tests/test_hilo_makers.py

@@ -0,0 +1,17 @@
+from pathlib import Path
+
+import pytest
+
+from nectarchain.makers.calibration.gain import HiLoNectarCAMCalibrationTool
+
+
+class TestHiLoNectarCAMCalibrationTool:
+    @pytest.fixture
+    def instance(self):
+        tool = HiLoNectarCAMCalibrationTool()
+        return tool
+
+    def test_init_output_path(self, instance):
+        instance.gain_file = Path("/tmp/gain.h5")
+        instance._init_output_path()
+        assert instance.output_path == Path("/tmp/gain_hilo_corrected.h5")

src/nectarchain/makers/component/__init__.py

@@ -9,6 +9,7 @@
     FlatFieldSingleNominalSPEStdNectarCAMComponent,
 )
 from .gain_component import GainNectarCAMComponent
+from .hilo_component import HiLoComponent
 from .pedestal_component import PedestalEstimationComponent
 from .photostatistic_algorithm import PhotoStatisticAlgorithm
 from .photostatistic_component import PhotoStatisticNectarCAMComponent
@@ -42,4 +43,5 @@
     "PhotoStatisticAlgorithm",
     "GainNectarCAMComponent",
     "FlatFieldComponent",
+    "HiLoComponent",
 ]

src/nectarchain/makers/component/hilo_component.py

@@ -0,0 +1,129 @@
+import copy
+
+import numpy as np
+from ctapipe.containers import EventType
+from ctapipe.core.traits import Path
+from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
+
+from ...data.container import (
+    PhotostatContainer,
+    SPEfitContainer,
+    merge_map_ArrayDataContainer,
+)
+from ...makers.component import ChargesComponent, GainNectarCAMComponent
+from ...utils import ComponentUtils, ContainerUtils
+from ...utils.constants import GAIN_DEFAULT, GAIN_LINEAR_RANGE, HILO_DEFAULT
+
+GAIN_CONTAINER_CLASSES = [PhotostatContainer, SPEfitContainer]
+
+
+class HiLoComponent(GainNectarCAMComponent):
+    """
+    Component that computes the HiLo ratio.
+    """
+
+    gain_file = Path(
+        default_value=None,
+        help="Path to h5 file with gain calibration coefficients",
+        allow_none=True,
+    ).tag(config=True)
+
+    SubComponents = copy.deepcopy(GainNectarCAMComponent.SubComponents)
+    SubComponents.default_value = [
+        "ChargesComponent",
+    ]
+    SubComponents.read_only = True
+
+    def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
+        chargesComponent_kwargs = {}
+        other_kwargs = {}
+        chargesComponent_configurable_traits = ComponentUtils.get_configurable_traits(
+            ChargesComponent
+        )
+        for key in kwargs.keys():
+            if key in chargesComponent_configurable_traits.keys():
+                chargesComponent_kwargs[key] = kwargs[key]
+            else:
+                other_kwargs[key] = kwargs[key]
+
+        super().__init__(
+            subarray=subarray, config=config, parent=parent, *args, **other_kwargs
+        )
+
+        self.chargesComponent = ChargesComponent(
+            subarray=subarray,
+            config=config,
+            parent=parent,
+            *args,
+            **chargesComponent_kwargs,
+        )
+        self._chargesContainer = None
+        self.log.debug(f"{kwargs.keys()}")
+
+        self._init_gain_container()
+
+    def _init_gain_container(self):
+        self.__gain_container = None
+
+        if self.gain_file is None:
+            raise ValueError("Need to provide a gain_file to compute HiLo ratio")
+
+        try:
+            self.__gain_container = ContainerUtils.get_container_from_hdf5(
+                self.gain_file,
+                GAIN_CONTAINER_CLASSES,
+            )
+            ContainerUtils.add_missing_pixels_to_container(
+                self.__gain_container, pad_value=GAIN_DEFAULT
+            )
+        except Exception as e:
+            raise RuntimeError(
+                f"Failed to initialize gain container from {self.gain_file}"
+            ) from e
+
+    def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
+        # For now only flat-field events, to be updated for e.g. white-target
+        if event.trigger.event_type == EventType.FLATFIELD:
+            self.chargesComponent(event=event, *args, **kwargs)
+
+    def finish(self, *args, **kwargs):
+        if self._chargesContainer is None:
+            chargesContainers = self.chargesComponent.finish(*args, **kwargs)
+            self._chargesContainer = merge_map_ArrayDataContainer(chargesContainers)
+
+        self._compute_low_gain()
+
+        return self.__gain_container
+
+    def _compute_low_gain(self):
+        ContainerUtils.add_missing_pixels_to_container(self._chargesContainer)
+        charges_hg = self._chargesContainer["charges_hg"]
+        charges_lg = self._chargesContainer["charges_lg"]
+        high_gain = self.__gain_container["high_gain"][:, 0]
+
+        # Mask the linear regime between low gain and high gain
+        charges_hg_pe = charges_hg / high_gain
+        mask_linearity = np.logical_and(
+            charges_hg_pe > np.min(GAIN_LINEAR_RANGE),
+            charges_hg_pe < np.max(GAIN_LINEAR_RANGE),
+        )
+
+        # Add failsafe to not divide by 0
+        mask = np.logical_and(mask_linearity, charges_lg > 0)
+
+        # Compute HiLo ratio (per pixel) for each event
+        hilo_ratio_all_events = np.divide(
+            charges_hg,
+            charges_lg,
+            where=mask,
+            out=np.full_like(charges_hg, np.nan, dtype=float),
+        )
+
+        # Compute HiLo ratio (per pixel) averaged over all events
+        hilo_ratio = np.nanmean(hilo_ratio_all_events, axis=0)
+
+        # Set default values if all events were masked
+        hilo_ratio = np.where(np.isnan(hilo_ratio), HILO_DEFAULT, hilo_ratio)
+
+        # Fill gain container
+        self.__gain_container["low_gain"][:, 0] = high_gain / hilo_ratio