Improvement of the gain calibration + Bugfix of charge extraction (#214)

* change defaut value of n and pp from HHV the free fit

* user script to produce ICRC2025 results

* Bug fix FF SPE combined

* add gain err to photostat method -> propagation of SPE resolution error

* follows photostat error

* logging improbemeltb

* fix name bug in charge makers

* -change pp and n to be taken from config file instead of SPE HHV result in SPE combined algo
-better management of SPE res used for PS method

* bugfix : wrong position of the yeild within a loop -> generator returned was pointing the last slice

* typo

* bugfix : when extracting charge and peak time with ctapipe extractor, even if selected_gain_channel is none, the expected sample axis is wafevorms.shape[-1]. Then if we provide (n_events,n_pix,n_sample) the output is (n_events, n_pix)

* imrove unit test to detect the issue related to charge extraction

* imporove debugging of scripts

* formatting

* fix pedestal makers and component according to last update of the containers I/O methods + unit test fix

* linter

* fix unit test of charge and waveforms makers
improve behavior if event_per_slice > n_events iin data

* fix gain pstat strict arguments

---------

Co-authored-by: Guillaume Grolleron <[email protected]>

Author: guillaumegrolleron

src/nectarchain/data/container/core.py

@@ -269,16 +269,16 @@ def _container_from_hdf5(
         if isinstance(path, str):
             path = Path(path)
         module = importlib.import_module(f"{container_class.__module__}")  # noqa :F841
-        container = eval(f"module.{container_class.__name__}")()
 
         with HDF5TableReader(path) as reader:
             if len(reader._h5file.root.__members__) > 1 and slice_index is None:
                 log.info(
                     f"reading {container_class.__name__} containing"
-                    f"{len(reader._h5file.root.__members__)}"
-                    f"slices, will return a generator"
+                    f" {len(reader._h5file.root.__members__)}"
+                    f" slices, will return a generator"
                 )
                 for data in np.sort(reader._h5file.root.__members__):
+                    container = eval(f"module.{container_class.__name__}")()
                     # container.containers[data] =
                     # eval(f"module.{container_class.__name__}s")()
 
@@ -319,8 +319,9 @@ def _container_from_hdf5(
                             f"table save, unable to load"
                         )
 
-                yield container
+                    yield container
             else:
+                container = eval(f"module.{container_class.__name__}")()
                 if slice_index is None:
                     log.info(
                         f"reading {container_class.__name__} containing"

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

@@ -221,7 +221,7 @@ def _init_output_path(self):
             if "run" in word:
                 HHVrun = int(word.split("run")[-1])
         str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
-            self.extractor_kwargs
+            method=self.method, extractor_kwargs=self.extractor_kwargs
         )
         if self.max_events is None:
             filename = (

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

@@ -122,18 +122,18 @@ def start(
             if len(FF_files) != 1 or len(Ped_files) != 1:
                 self.log.info(
                     f"{len(FF_files)} computed charges FF files found"
-                    f"with max_events >"
+                    f" with max_events >"
                     f"{self.max_events} for run {self.run_number}"
-                    f"with extraction method"
+                    f" with extraction method"
                     f"{self.method} and {str_extractor_kwargs},\n reload charges"
-                    f"from event loop"
+                    f" from event loop"
                 )
                 self.log.info(
                     f"{len(Ped_files)} computed charges FF files found"
-                    f"with max_events >"
+                    f" with max_events >"
                     f"{self.max_events} for run {self.Ped_run_number}"
-                    f"with extraction"
-                    f"method FullWaveformSum,\n reload charges from event loop"
+                    f" with extraction"
+                    f" method FullWaveformSum,\n reload charges from event loop"
                 )
 
             super().start(
@@ -162,7 +162,7 @@ def start(
                     self.components[
                         0
                     ]._FF_chargesContainers = merge_map_ArrayDataContainer(
-                        chargesContainers
+                        next(chargesContainers)
                     )
                 else:
                     self.log.info("merging along slices")
@@ -196,7 +196,7 @@ def start(
                     self.components[
                         0
                     ]._Ped_chargesContainers = merge_map_ArrayDataContainer(
-                        chargesContainers
+                        next(chargesContainers)
                     )
                 else:
                     self.log.info("merging along slices")

src/nectarchain/makers/calibration/pedestal_makers.py

@@ -3,6 +3,7 @@
 import pathlib
 
 import numpy as np
+import tables
 from ctapipe.core.traits import ComponentNameList
 from ctapipe_io_nectarcam.constants import HIGH_GAIN, LOW_GAIN, N_GAINS
 
@@ -55,18 +56,28 @@ def _combine_results(self):
         Method that combines sliced results to reduce memory load
         Can only be called after the file with the sliced results has been saved to disk
         """
+        already_combined = False
+        with tables.open_file(self.output_path, mode="r") as f:
+            keys = list(f.root._v_children)
+            if "data_combined" in keys:
+                log.error(
+                    "Trying to combine results that already contain combined data"
+                )
+                already_combined = True
 
         # re-open results
-        pedestalContainers = next(
-            NectarCAMPedestalContainers.from_hdf5(self.output_path)
-        )
+        if already_combined:
+            pedestalContainers = NectarCAMPedestalContainers.from_hdf5(
+                self.output_path,
+                slice_index="combined",
+            )
+        else:
+            pedestalContainers = NectarCAMPedestalContainers.from_hdf5(self.output_path)
+
         # Loop over sliced results to fill the combined results
-        if "data_combined" in pedestalContainers.containers.keys():
-            log.error("Trying to combine results that already contain combined data")
         self.log.info("Combine sliced results")
-        for i, (_, pedestalContainer) in enumerate(
-            pedestalContainers.containers.items()
-        ):
+        for i, _pedestalContainer in enumerate(pedestalContainers):
+            pedestalContainer = list(_pedestalContainer.containers.values())[0]
             if i == 0:
                 # initialize fields for the combined results based on first slice
                 nsamples = pedestalContainer.nsamples

src/nectarchain/makers/calibration/tests/test_pedestal_tool.py

@@ -1,6 +1,8 @@
+import os
 import tempfile
 
 import numpy as np
+import tables
 from ctapipe.utils import get_dataset_path
 from ctapipe_io_nectarcam.constants import N_SAMPLES
 
@@ -36,6 +38,12 @@ def test_base(self):
             n_pixels = runs["N pixels"][i]
             with tempfile.TemporaryDirectory() as tmpdirname:
                 outfile = tmpdirname + "/pedestal.h5"
+                try:
+                    os.remove(outfile)
+                except FileNotFoundError:
+                    pass
+                except Exception as e:
+                    raise e
 
                 # run tool
                 tool = PedestalNectarCAMCalibrationTool(
@@ -76,18 +84,16 @@ def test_base(self):
                 assert np.allclose(output.pedestal_std_lg, 2.5, atol=2.3)
 
                 # Check output on disk
-                # FIXME: use tables for the moment, update when h5 reader in nectarchain
-                #  is working
-                pedestalContainers = next(
-                    NectarCAMPedestalContainers.from_hdf5(outfile)
-                )
+                pedestalContainers = NectarCAMPedestalContainers.from_hdf5(outfile)
                 j = 0
-                for key, pedestalContainer in pedestalContainers.containers.items():
+                list_slices = []
+                for _pedestalContainer in pedestalContainers:
+                    key = list(_pedestalContainer.containers.keys())[0]
                     if "combined" in key:
                         continue
-                    # Check individual groups
-                    group_name = "data_{}".format(i + 1)
-                    assert group_name in pedestalContainers.containers.keys()
+                    pedestalContainer = _pedestalContainer.containers[key]
+                    list_slices.append(int(key.split("data_")[1]))
+
                     assert pedestalContainer.nsamples == N_SAMPLES
                     assert np.allclose(
                         pedestalContainer.nevents, events_per_slice, atol=7
@@ -110,11 +116,19 @@ def test_base(self):
                         N_SAMPLES,
                     )
                     j += 1
+                assert (np.sort(list_slices) == np.arange(1, n_slices[i] + 1)).all()
+
                 # Check combined results
+                group_name = "data_combined"
+                with tables.open_file(outfile, mode="r") as f:
+                    keys = list(f.root._v_children)
+                    assert group_name in keys
                 pedestalContainers = next(
-                    NectarCAMPedestalContainers.from_hdf5(outfile)
+                    NectarCAMPedestalContainers.from_hdf5(
+                        outfile, slice_index="combined"
+                    )
                 )
-                group_name = "data_combined"
+
                 pedestalContainer = pedestalContainers.containers[group_name]
                 assert pedestalContainer.nsamples == N_SAMPLES
                 assert np.all(pedestalContainer.nevents == max_events[i])

src/nectarchain/makers/charges_makers.py

@@ -56,7 +56,7 @@ def _init_output_path(self):
         if self.max_events is None:
             filename = (
                 f"{self.name}_run{self.run_number}_{self.method}"
-                f"{str_extractor_kwargs}.h5"
+                f"_{str_extractor_kwargs}.h5"
             )
         else:
             filename = (
@@ -87,7 +87,7 @@ def start(
                 files = []
             if len(files) != 1:
                 self.log.info(
-                    f"{len(files)} computed wavforms files found with max_events >="
+                    f"{len(files)} computed waveforms files found with max_events >="
                     f"{self.max_events}  for run {self.run_number}, reload waveforms"
                     f"from event loop"
                 )
@@ -99,7 +99,7 @@ def start(
                 )
             else:
                 self.log.info(
-                    f"{files[0]} is the computed wavforms files found"
+                    f"{files[0]} is the computed waveforms files found"
                     f"with max_events >="
                     f"{self.max_events}  for run {self.run_number}"
                 )

src/nectarchain/makers/component/charges_component.py

@@ -588,23 +588,21 @@ def compute_charges(
             method=method, subarray=subarray, **kwargs
         )
         out = np.array(
-            [
-                CtapipeExtractor.get_image_peak_time(
-                    imageExtractor(
-                        waveforms=waveforms,
-                        tel_id=tel_id,
-                        selected_gain_channel=None,
-                        broken_pixels=broken_pixels,
-                    )
+            CtapipeExtractor.get_image_peak_time(
+                imageExtractor(
+                    waveforms=waveforms,
+                    tel_id=tel_id,
+                    selected_gain_channel=None,
+                    broken_pixels=broken_pixels,
                 )
-            ]
+            )
         )
-
+        # bc waveforms shape is (n_events, n_pix, n_samples),
+        # out is (2, n_events, n_pix)
+        # with the first dimension for charge and the second for peak time
         return ChargesContainer.fields[f"charges_{gain_label}"].dtype.type(
-            out[:, 0, channel, :]
-        ), ChargesContainer.fields[f"peak_{gain_label}"].dtype.type(
-            out[:, 1, channel, :]
-        )
+            out[0, ...]
+        ), ChargesContainer.fields[f"peak_{gain_label}"].dtype.type(out[1, ...])
 
     @staticmethod
     def histo_hg(

src/nectarchain/makers/component/photostatistic_algorithm.py

@@ -94,8 +94,8 @@ def __get_charges_FF_Ped_reshaped(
             ],
             dtype=bool,
         )
-        out["SPE_resolution"] = SPE_result.resolution[mask_SPE].T[0]
-        out["SPE_high_gain"] = SPE_result.high_gain[mask_SPE].T[0]
+        out["SPE_resolution"] = SPE_result.resolution[mask_SPE].T
+        out["SPE_high_gain"] = SPE_result.high_gain[mask_SPE].T
 
         out["pixels_id"] = SPEFFPed_intersection
         out["FFcharge_hg"] = ChargesComponent.select_charges_hg(
@@ -120,7 +120,7 @@ def __check_shape(self) -> None:
         try:
             self.__FFcharge_hg[0] * self.__FFcharge_lg[0] * self.__Pedcharge_hg[
                 0
-            ] * self.__Pedcharge_lg[0] * self.__SPE_resolution * self._pixels_id
+            ] * self.__Pedcharge_lg[0] * self.__SPE_resolution[0] * self._pixels_id
         except Exception as e:
             log.error(e, exc_info=True)
             raise e
@@ -132,19 +132,21 @@ def run(self, pixels_id: np.ndarray = None, **kwargs) -> None:
         mask = np.array(
             [pixel_id in pixels_id for pixel_id in self._pixels_id], dtype=bool
         )
+        gainHG_err = self.gainHG_err * mask
+        gainLG_err = self.gainLG_err * mask
         self._results.high_gain = np.array(
-            (self.gainHG * mask, np.zeros(self.npixels), np.zeros(self.npixels))
+            (self.gainHG * mask, gainHG_err, gainHG_err)
         ).T
         self._results.low_gain = np.array(
-            (self.gainLG * mask, np.zeros(self.npixels), np.zeros(self.npixels))
+            (self.gainLG * mask, gainLG_err, gainLG_err)
         ).T
         self._results.is_valid = mask
 
         figpath = kwargs.get("figpath", False)
         if figpath:
             os.makedirs(figpath, exist_ok=True)
             fig = __class__.plot_correlation(
-                self._results.high_gain.T[0], self.__SPE_high_gain
+                self._results.high_gain.T[0], self.__SPE_high_gain[0]
             )
             fig.savefig(f"{figpath}/plot_correlation_Photo_Stat_SPE.pdf")
             fig.clf()
@@ -286,7 +288,39 @@ def gainHG(self) -> float:
             np.power(self.sigmaChargeHG, 2)
             - np.power(self.sigmaPedHG, 2)
             - np.power(self.BHG * self.meanChargeHG, 2)
-        ) / (self.meanChargeHG * (1 + np.power(self.SPE_resolution, 2)))
+        ) / (self.meanChargeHG * (1 + np.power(self.SPE_resolution[0], 2)))
+
+    @property
+    def gainHG_err(self) -> float:
+        """Calculates and returns the gain for high gain charge data.
+
+        Returns:
+            float: The gain for high gain charge data.
+        """
+        return np.sqrt(
+            np.power(
+                self.gainHG
+                * (-2 * self.SPE_resolution[0])
+                * np.mean(self.SPE_resolution[1:], axis=0),
+                2,
+            )
+        )
+
+    @property
+    def gainLG_err(self) -> float:
+        """Calculates and returns the gain for high gain charge data.
+
+        Returns:
+            float: The gain for high gain charge data.
+        """
+        return np.sqrt(
+            np.power(
+                self.gainLG
+                * (-2 * self.SPE_resolution[0])
+                * np.mean(self.SPE_resolution[1:], axis=0),
+                2,
+            )
+        )
 
     @property
     def sigmaPedLG(self) -> float:
@@ -357,7 +391,7 @@ def gainLG(self) -> float:
             np.power(self.sigmaChargeLG, 2)
             - np.power(self.sigmaPedLG, 2)
             - np.power(self.BLG * self.meanChargeLG, 2)
-        ) / (self.meanChargeLG * (1 + np.power(self.SPE_resolution, 2)))
+        ) / (self.meanChargeLG * (1 + np.power(self.SPE_resolution[0], 2)))
 
     @property
     def results(self):

src/nectarchain/makers/component/spe/parameters_SPECombined_fromHHVFit.yaml

@@ -9,18 +9,18 @@
     min : 0.01,
     max : 5.0
   },
-  pp : { 
-    value : 0.454,
+  pp : { #useless bc use HHVStd result
+    value : 0.42,
     min : .NAN,
     max : .NAN
   }, 
-  resolution : {
-    value : 0.5,
+  resolution : { #useless bc use HHVStd result
+    value : 0.54,
     min : .NAN,
     max : .NAN
   },
-  n : {
-    value : 0.713,
+  n : { #useless bc use HHVStd result
+    value : 0.813,
     min : .NAN,
     max : .NAN
   },

src/nectarchain/makers/component/spe/parameters_SPEHHVStd.yaml

@@ -10,7 +10,7 @@
     max : 5.0
   },
   pp : { 
-    value : 0.454,
+    value : 0.42,
     min : .NAN,
     max : .NAN
   }, 
@@ -20,7 +20,7 @@
     max : 0.7
   },
   n : {
-    value : 0.713,
+    value : 0.813,
     min : .NAN,
     max : .NAN
   },