From 7f1d854e3d23ae3a62f33e8d9fb51331cffab024 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 13 Feb 2024 10:22:28 +0100
Subject: [PATCH 01/20] create structure for pedestal component/tool

---
 .../makers/calibration/pedestalMakers.py      |  32 ++++
 .../makers/component/PedestalComponent.py     | 163 ++++++++++++++++++
 src/nectarchain/makers/component/__init__.py  |   2 +
 3 files changed, 197 insertions(+)
 create mode 100644 src/nectarchain/makers/component/PedestalComponent.py

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index ec0be630..0e3066bd 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -10,6 +10,38 @@
 
 
 class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
+
+    name = "PedestalNectarCAMCalibrationTool"
+
+    componentsList = ComponentNameList(
+        NectarCAMComponent,
+        default_value=["PedestalEstimationComponent"],
+        help="List of Component names to be applied, the order will be respected",
+    ).tag(config=True)
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+        str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
+            self.extractor_kwargs
+        )
+        if not (self.reload_events):
+            files = DataManagement.find_waveforms(
+                run_number=self.run_number,
+                max_events=self.max_events,
+            )
+            if len(files) == 1:
+                log.warning(
+                    "You asked events_per_slice but you don't want to reload events and a charges file is on disk, then events_per_slice is set to None"
+                )
+                self.events_per_slice = None
+
+
+
+
+
+
+
     def start(self):
         raise NotImplementedError(
             "The computation of the pedestal calibration is not yet implemented, feel free to contribute !:)"
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
new file mode 100644
index 00000000..f935a611
--- /dev/null
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -0,0 +1,163 @@
+import logging
+
+logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s")
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger("__main__").handlers
+
+import copy
+
+from .core import NectarCAMComponent
+from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
+from ctapipe_io_nectarcam.constants import (
+    N_GAINS, N_PIXELS, N_SAMPLES,
+    HIGH_GAIN, LOW_GAIN,
+)
+from ctapipe.core.traits import Dict, Unicode
+from ctapipe.containers import PedestalContainer
+from ...data.container import merge_map_ArrayDataContainer
+from ...utils import ComponentUtils
+from .waveformsComponent import WaveformsComponent
+
+import numpy as np
+
+__all__ = [
+    "PedestalEstimationComponent",
+]
+
+class PedestalComponent(NectarCAMComponent):
+    """
+    Component that computes calibration pedestal coefficients from raw data.
+    """
+    PedestalFilterAlgorithm = Unicode(
+        "PedestalWaveformStdFilter",
+        help="The waveform filter method",
+        read_only=True,
+    ).tag(config=True)
+
+    extractor_kwargs = Dict(
+        default_value={},
+        help="The kwargs to be pass to the waveform filter method",
+    ).tag(config=True)
+
+    # add parameters for min max time
+
+    SubComponents = copy.deepcopy(NectarCAMComponent)
+    SubComponents.default_value = [
+        "WaveformsComponent",
+        f"{PedestalFilterAlgorithm.default_value}",
+    ]
+    SubComponents.read_only = True
+
+    def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
+        waveformsComponent_kwargs = {}
+        self._PedestalFilterAlgorithm_kwargs = {}
+        other_kwargs = {}
+        waveformsComponent_configurable_traits = ComponentUtils.get_configurable_traits(
+            WaveformsComponent
+        )
+        pedestalFilterAlgorithm_configurable_traits = ComponentUtils.get_configurable_traits(
+            eval(self.PedestalFilterAlgorithm)
+        )
+
+        for key in kwargs.keys():
+            if key in waveformsComponent_configurable_traits.keys():
+                waveformsComponent_kwargs[key] = kwargs[key]
+            # elif key in waveformsComponent_configurable_traits.keys():
+            #     self._SPEfitalgorithm_kwargs[key] = kwargs[key]
+            else:
+                other_kwargs[key] = kwargs[key]
+
+        super().__init__(
+            subarray=subarray, config=config, parent=parent, *args, **kwargs
+        )
+
+        self.waveformsComponent = WaveformsComponent(
+            subarray=subarray,
+            config=config,
+            parent=parent,
+            *args,
+            **waveformsComponent_kwargs,
+        )
+        self._waveformsContainers = None
+
+        self.__wfs_mask = None
+
+    @staticmethod
+    def calculate_stats(waveformsContainers,wfs_mask,statistics):
+        """
+        Calculate statistics for the pedestals from a waveforms container.
+
+        Parameters
+        ----------
+        waveformsContainers : `~nectarchain.data.container.WaveFormsContainer`
+            Waveforms container
+        wfs_mask : `numpy.ndarray`
+            Mask to apply to exclude outliers with shape (n_pixels,n_samples)
+        statistics : `list`
+            Names of the statistics (numpy attributes) to compute
+
+        Returns
+        ----------
+        statistics : `dict`
+            A dictionary containing 3D (n_chan,n_pixels,n_samples) arrays for each statistic
+        """
+
+        ped_stats = {}
+
+        for stat in statistics:
+            # Calculate the statistic
+            ped_stat_hg = getattr(np, stat)(waveformsContainers.wfs_hg[wfs_mask])
+            ped_stat_lg = getattr(np, stat)(waveformsContainers.wfs_lg[wfs_mask])
+
+            # Create a 3D array for the statistic
+            ped_stat = np.zeros([N_GAINS, N_PIXELS, N_SAMPLES])
+            ped_stat[HIGH_GAIN] = ped_stat_hg
+            ped_stat[LOW_GAIN] = ped_stat_lg
+
+            # Store the result in the dictionary
+            ped_stats[stat] = ped_stat
+
+        return ped_stats
+
+    def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
+
+        self.waveformsComponent(event=event, *args, **kwargs)
+
+        is_empty = False
+        if self._waveformsContainers is None:
+            self._waveformsContainers = self.waveformsComponent.finish(*args, **kwargs)
+            is_empty = self._waveformsContainers.is_empty()
+            if is_empty:
+                log.warning("empty waveformsContainer in output")
+            else:
+                self._waveformsContainers = merge_map_ArrayDataContainer(
+                    self._waveformsContainers
+                )
+
+        if not (is_empty):
+            # change this into something that creates a real mask
+            # one mask for both HG and LG or separate? mask only on hg which is more sensitive
+            self.__wfs_mask = np.ones([N_PIXELS,N_SAMPLES],dtype=bool)
+        else:
+            pass
+
+    def finish(self, *args, **kwargs):
+
+        # compute statistics for the pedestals
+        # the statistic names must be valid numpy attributes
+        statistics = ['mean', 'median', 'std']
+        ped_stats = self.calculate_stats(self._waveformsContainers,self.__wfs_mask,statistics)
+
+        metadata = {} # store information about filtering method and params
+
+        output = PedestalContainer(
+            n_events = self._waveformsContainers.nevents,
+            sample_time = 0.,#to be filled, mean of min/max
+            sample_time_min = self._waveformsContainers.ucts_timestamp.min(),
+            sample_time_max = self._waveformsContainers.ucts_timestamp.max(),
+            charge_mean = ped_stats['mean'],
+            charge_median = ped_stats['median'],
+            charge_std = ped_stats['std'],
+            meta = metadata,
+        )
+        return output
diff --git a/src/nectarchain/makers/component/__init__.py b/src/nectarchain/makers/component/__init__.py
index 8ab6a9a6..f7b85ec9 100644
--- a/src/nectarchain/makers/component/__init__.py
+++ b/src/nectarchain/makers/component/__init__.py
@@ -1,6 +1,7 @@
 from .chargesComponent import *
 from .core import *
 from .FlatFieldSPEComponent import *
+from .PedestalComponent import *
 from .gainComponent import *
 from .photostatistic_algorithm import *
 from .photostatistic_component import *
@@ -20,6 +21,7 @@
     "FlatFieldCombinedSPEStdNectarCAMComponent",
     "ChargesComponent",
     "WaveformsComponent",
+    "PedestalEstimationComponent",
     "PhotoStatisticNectarCAMComponent",
     "PhotoStatisticAlgorithm",
 ]

From 75f807c887e79fd8ce771c5312dfc0865691ba69 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 13 Feb 2024 12:07:58 +0100
Subject: [PATCH 02/20] debugging tool structure

Still not working
---
 .../makers/calibration/pedestalMakers.py      | 95 +++++++++++++++----
 .../makers/component/PedestalComponent.py     | 38 +++++---
 .../user_scripts/ltibaldo/example_pedestal.py | 33 +++++++
 3 files changed, 134 insertions(+), 32 deletions(-)
 create mode 100644 src/nectarchain/user_scripts/ltibaldo/example_pedestal.py

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index 0e3066bd..150f6202 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -4,7 +4,16 @@
 log = logging.getLogger(__name__)
 log.handlers = logging.getLogger("__main__").handlers
 
+import pathlib
+import numpy as np
+
+from ctapipe.core.traits import ComponentNameList, Bool
+
 from .core import NectarCAMCalibrationTool
+from ...data.container import WaveformsContainer, WaveformsContainers
+from ...data.management import DataManagement
+from ...data.container.core import merge_map_ArrayDataContainer
+from ..component import ArrayDataComponent, NectarCAMComponent
 
 __all__ = ["PedestalNectarCAMCalibrationTool"]
 
@@ -13,6 +22,10 @@ class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
 
     name = "PedestalNectarCAMCalibrationTool"
 
+    reload_events = Bool(
+        default_value=False, help="Reload the waveforms from raw data"
+    )
+
     componentsList = ComponentNameList(
         NectarCAMComponent,
         default_value=["PedestalEstimationComponent"],
@@ -22,27 +35,73 @@ class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
+    def _init_output_path(self):
         str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
             self.extractor_kwargs
         )
-        if not (self.reload_events):
-            files = DataManagement.find_waveforms(
-                run_number=self.run_number,
-                max_events=self.max_events,
+        if self.events_per_slice is None:
+            ext = ".h5"
+        else:
+            ext = f"_sliced{self.events_per_slice}.h5"
+        if self.max_events is None:
+            filename = f"{self.name}_run{self.run_number}{ext}"
+        else:
+            filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}{ext}"
+
+        self.output_path = pathlib.Path(
+            f"{os.environ.get('NECTARCAMDATA','/tmp')}/PedestalEstimation/{filename}"
+        )
+
+    def start(
+        self,
+        max_events=np.inf,
+        # trigger_type: list = None,
+        restart_from_beginning: bool = False,
+        *args,
+        **kwargs,
+    ):
+        # need to implement waveform filter methods
+        # str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
+        #     self.extractor_kwargs
+        # )
+
+        files = DataManagement.find_waveforms(
+            run_number=self.run_number,
+            max_events=self.max_events,
+        )
+        if self.reload_events or len(files) != 1:
+            if len(files) != 1:
+                self.log.info(
+                    f"{len(files)} computed waveforms files found with max_events > {self.max_events} for run {self.run_number}, reload waveforms from event loop"
+                )
+            print("Start parent class", super().name)
+            super().start(
+                restart_from_beginning=restart_from_beginning,
+                *args,
+                **kwargs,
             )
-            if len(files) == 1:
-                log.warning(
-                    "You asked events_per_slice but you don't want to reload events and a charges file is on disk, then events_per_slice is set to None"
+        else:
+            print('Enter loop to read waveforms')
+            self.log.info(f"reading waveforms from files {files[0]}")
+            waveformsContainers = WaveformsContainer.from_hdf5(files[0])
+            if isinstance(waveformsContainers, WaveformsContainer):
+                self.components[0]._waveformsContainers = waveformsContainers
+            elif isinstance(waveformsContainers, WaveformsContainers):
+                self.log.debug("merging along TriggerType")
+                self.components[0]._waveformsContainers = merge_map_ArrayDataContainer(
+                    waveformsContainers
+                )
+            else:
+                self.log.debug("merging along slices")
+                waveformsContainers_merges_along_slices = (
+                    ArrayDataComponent.merge_along_slices(
+                        containers_generator=waveformContainers
+                    )
+                )
+                self.log.debug("merging along TriggerType")
+                self.components[0]._waveformsContainers = merge_map_ArrayDataContainer(
+                    waveformsContainers_merges_along_slices
                 )
-                self.events_per_slice = None
-
-
-
-
-
-
 
-    def start(self):
-        raise NotImplementedError(
-            "The computation of the pedestal calibration is not yet implemented, feel free to contribute !:)"
-        )
+        print("Input wfs shape", np.shape(self.components[0]._waveformsContainers.wfs_hg))
+        print("Input wfs 0", self.components[0]._waveformsContainers.wfs_hg[0])
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index f935a611..fe0c73b9 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -24,17 +24,18 @@
     "PedestalEstimationComponent",
 ]
 
-class PedestalComponent(NectarCAMComponent):
+class PedestalEstimationComponent(NectarCAMComponent):
     """
     Component that computes calibration pedestal coefficients from raw data.
     """
-    PedestalFilterAlgorithm = Unicode(
-        "PedestalWaveformStdFilter",
-        help="The waveform filter method",
-        read_only=True,
-    ).tag(config=True)
-
-    extractor_kwargs = Dict(
+    # PedestalFilterAlgorithm = Unicode(
+    #     "PedestalWaveformStdFilter",
+    #     help="The waveform filter method",
+    #     read_only=True,
+    # ).tag(config=True)
+
+    #not implemented yet, placeholder
+    filter_kwargs = Dict(
         default_value={},
         help="The kwargs to be pass to the waveform filter method",
     ).tag(config=True)
@@ -44,7 +45,7 @@ class PedestalComponent(NectarCAMComponent):
     SubComponents = copy.deepcopy(NectarCAMComponent)
     SubComponents.default_value = [
         "WaveformsComponent",
-        f"{PedestalFilterAlgorithm.default_value}",
+        #f"{PedestalFilterAlgorithm.default_value}",
     ]
     SubComponents.read_only = True
 
@@ -55,9 +56,9 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
         waveformsComponent_configurable_traits = ComponentUtils.get_configurable_traits(
             WaveformsComponent
         )
-        pedestalFilterAlgorithm_configurable_traits = ComponentUtils.get_configurable_traits(
-            eval(self.PedestalFilterAlgorithm)
-        )
+        # pedestalFilterAlgorithm_configurable_traits = ComponentUtils.get_configurable_traits(
+        #     eval(self.PedestalFilterAlgorithm)
+        # )
 
         for key in kwargs.keys():
             if key in waveformsComponent_configurable_traits.keys():
@@ -137,7 +138,7 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
         if not (is_empty):
             # change this into something that creates a real mask
             # one mask for both HG and LG or separate? mask only on hg which is more sensitive
-            self.__wfs_mask = np.ones([N_PIXELS,N_SAMPLES],dtype=bool)
+            self.__wfs_mask = np.ones(np.shape(self._waveformsContainers.wfs_hg),dtype=bool)
         else:
             pass
 
@@ -146,13 +147,22 @@ def finish(self, *args, **kwargs):
         # compute statistics for the pedestals
         # the statistic names must be valid numpy attributes
         statistics = ['mean', 'median', 'std']
+        #print(np.shape(self._waveformsContainers.wfs_hg))
         ped_stats = self.calculate_stats(self._waveformsContainers,self.__wfs_mask,statistics)
 
         metadata = {} # store information about filtering method and params
 
+        # set reference time to mean between min and max
+        # is this choice reasonable?
+        #print(self._waveformsContainers.ucts_timestamp)
+        ref_time = np.mean(self._waveformsContainers.ucts_timestamp.min(),
+                           self._waveformsContainers.ucts_timestamp.max())
+        print(self._waveformsContainers.ucts_timestamp.min(),
+              self._waveformsContainers.ucts_timestamp.max())
+
         output = PedestalContainer(
             n_events = self._waveformsContainers.nevents,
-            sample_time = 0.,#to be filled, mean of min/max
+            sample_time = ref_time,#to be filled, mean of min/max
             sample_time_min = self._waveformsContainers.ucts_timestamp.min(),
             sample_time_max = self._waveformsContainers.ucts_timestamp.max(),
             charge_mean = ped_stats['mean'],
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
new file mode 100644
index 00000000..825efa8c
--- /dev/null
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -0,0 +1,33 @@
+import logging
+import os
+import pathlib
+
+logging.basicConfig(
+    format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO
+)
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger("__main__").handlers
+
+from nectarchain.makers.calibration import (
+    PedestalNectarCAMCalibrationTool,
+)
+
+run_number = 3938
+max_events= 100
+outfile = "/Users/ltibaldo/tmp/test_pedestal/pedestal_{}.h5".format(run_number)
+
+tool = PedestalNectarCAMCalibrationTool(
+    progress_bar=True,
+    run_number=run_number,
+    max_events=max_events,
+    log_level=20,
+    reload_events=False,
+    output_path=outfile,
+    overwrite = True,
+)
+
+tool.initialize()
+tool.setup()
+
+tool.start()
+tool.finish(return_output_component=True)
\ No newline at end of file

From d7b042399379f314a03cafc4e3cfe5530d770a27 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 13 Feb 2024 16:51:40 +0100
Subject: [PATCH 03/20] fix computation of mean, etc

This is an example that works on stored waveforms without any filtering
---
 .../makers/calibration/pedestalMakers.py      |  9 +++--
 .../makers/component/PedestalComponent.py     | 35 +++++++++----------
 .../user_scripts/ltibaldo/example_pedestal.py |  2 +-
 3 files changed, 22 insertions(+), 24 deletions(-)

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index 150f6202..cfe96b16 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -36,9 +36,9 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
     def _init_output_path(self):
-        str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
-            self.extractor_kwargs
-        )
+        # str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
+        #     self.extractor_kwargs
+        # )
         if self.events_per_slice is None:
             ext = ".h5"
         else:
@@ -95,7 +95,7 @@ def start(
                 self.log.debug("merging along slices")
                 waveformsContainers_merges_along_slices = (
                     ArrayDataComponent.merge_along_slices(
-                        containers_generator=waveformContainers
+                        containers_generator=waveformsContainers
                     )
                 )
                 self.log.debug("merging along TriggerType")
@@ -104,4 +104,3 @@ def start(
                 )
 
         print("Input wfs shape", np.shape(self.components[0]._waveformsContainers.wfs_hg))
-        print("Input wfs 0", self.components[0]._waveformsContainers.wfs_hg[0])
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index fe0c73b9..7d06a15b 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -5,6 +5,8 @@
 log.handlers = logging.getLogger("__main__").handlers
 
 import copy
+import astropy.units as u
+import numpy.ma as ma
 
 from .core import NectarCAMComponent
 from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
@@ -12,6 +14,7 @@
     N_GAINS, N_PIXELS, N_SAMPLES,
     HIGH_GAIN, LOW_GAIN,
 )
+from ctapipe_io_nectarcam import time_from_unix_tai_ns
 from ctapipe.core.traits import Dict, Unicode
 from ctapipe.containers import PedestalContainer
 from ...data.container import merge_map_ArrayDataContainer
@@ -81,7 +84,7 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
         )
         self._waveformsContainers = None
 
-        self.__wfs_mask = None
+        self._wfs_mask = None
 
     @staticmethod
     def calculate_stats(waveformsContainers,wfs_mask,statistics):
@@ -106,12 +109,15 @@ def calculate_stats(waveformsContainers,wfs_mask,statistics):
         ped_stats = {}
 
         for stat in statistics:
-            # Calculate the statistic
-            ped_stat_hg = getattr(np, stat)(waveformsContainers.wfs_hg[wfs_mask])
-            ped_stat_lg = getattr(np, stat)(waveformsContainers.wfs_lg[wfs_mask])
+            # Calculate the statistic along axis = 0, that is over events
+            ped_stat_hg = getattr(np, stat)(
+                ma.masked_array(waveformsContainers.wfs_hg,wfs_mask),axis=0)
+            ped_stat_lg = getattr(np, stat)(
+                ma.masked_array(waveformsContainers.wfs_lg, wfs_mask), axis=0)
 
             # Create a 3D array for the statistic
-            ped_stat = np.zeros([N_GAINS, N_PIXELS, N_SAMPLES])
+            array_shape = np.append([N_GAINS], np.shape(waveformsContainers.wfs_hg[0]))
+            ped_stat = np.zeros(array_shape)
             ped_stat[HIGH_GAIN] = ped_stat_hg
             ped_stat[LOW_GAIN] = ped_stat_lg
 
@@ -138,7 +144,7 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
         if not (is_empty):
             # change this into something that creates a real mask
             # one mask for both HG and LG or separate? mask only on hg which is more sensitive
-            self.__wfs_mask = np.ones(np.shape(self._waveformsContainers.wfs_hg),dtype=bool)
+            self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg),dtype=bool)
         else:
             pass
 
@@ -148,26 +154,19 @@ def finish(self, *args, **kwargs):
         # the statistic names must be valid numpy attributes
         statistics = ['mean', 'median', 'std']
         #print(np.shape(self._waveformsContainers.wfs_hg))
-        ped_stats = self.calculate_stats(self._waveformsContainers,self.__wfs_mask,statistics)
+        ped_stats = self.calculate_stats(self._waveformsContainers,self._wfs_mask,statistics)
 
         metadata = {} # store information about filtering method and params
 
-        # set reference time to mean between min and max
-        # is this choice reasonable?
-        #print(self._waveformsContainers.ucts_timestamp)
-        ref_time = np.mean(self._waveformsContainers.ucts_timestamp.min(),
-                           self._waveformsContainers.ucts_timestamp.max())
-        print(self._waveformsContainers.ucts_timestamp.min(),
-              self._waveformsContainers.ucts_timestamp.max())
-
         output = PedestalContainer(
             n_events = self._waveformsContainers.nevents,
-            sample_time = ref_time,#to be filled, mean of min/max
-            sample_time_min = self._waveformsContainers.ucts_timestamp.min(),
-            sample_time_max = self._waveformsContainers.ucts_timestamp.max(),
+            sample_time= np.nan, # TODO : does this need to be filled?
+            sample_time_min= time_from_unix_tai_ns(self._waveformsContainers.ucts_timestamp.min()).value * u.s,
+            sample_time_max = time_from_unix_tai_ns(self._waveformsContainers.ucts_timestamp.max()).value * u.s,
             charge_mean = ped_stats['mean'],
             charge_median = ped_stats['median'],
             charge_std = ped_stats['std'],
             meta = metadata,
         )
+
         return output
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index 825efa8c..82313652 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -30,4 +30,4 @@
 tool.setup()
 
 tool.start()
-tool.finish(return_output_component=True)
\ No newline at end of file
+output = tool.finish(return_output_component=True)
\ No newline at end of file

From 6721468df88fbb52e3227fbfb29bf38c84426fb2 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 13 Feb 2024 17:45:12 +0100
Subject: [PATCH 04/20] Added container for pedestals

---
 src/nectarchain/data/container/__init__.py    |  1 +
 .../data/container/pedestalContainer.py       | 75 +++++++++++++++++++
 .../makers/component/PedestalComponent.py     | 45 +++++------
 3 files changed, 99 insertions(+), 22 deletions(-)
 create mode 100644 src/nectarchain/data/container/pedestalContainer.py

diff --git a/src/nectarchain/data/container/__init__.py b/src/nectarchain/data/container/__init__.py
index 464a8848..5b525634 100644
--- a/src/nectarchain/data/container/__init__.py
+++ b/src/nectarchain/data/container/__init__.py
@@ -9,3 +9,4 @@
 from .eventSource import *
 from .gainContainer import *
 from .waveformsContainer import *
+from .pedestalContainer import *
diff --git a/src/nectarchain/data/container/pedestalContainer.py b/src/nectarchain/data/container/pedestalContainer.py
new file mode 100644
index 00000000..9e0debfa
--- /dev/null
+++ b/src/nectarchain/data/container/pedestalContainer.py
@@ -0,0 +1,75 @@
+import logging
+
+logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s")
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger("__main__").handlers
+
+import numpy as np
+from ctapipe.containers import Field
+
+from .core import NectarCAMContainer
+
+
+class NectarCAMPedestalContainer(NectarCAMContainer):
+    """
+    A container that holds estimated pedestals
+
+    Fields:
+        nsamples (int): The number of samples in the waveforms.
+        nevents (int): The number of events used to estimate the pedestals.
+        pixels_id (np.ndarray): An array of pixel IDs.
+        ucts_timestamp_min (int): The minimum of the input events UCTS timestamps.
+        ucts_timestamp_max (int): The maximum of the input events UCTS timestamps.
+        pedestal_mean_hg (np.ndarray): An array of high gain mean pedestals.
+        pedestal_mean_lg (np.ndarray): An array of low gain mean pedestals.
+        pedestal_median_hg (np.ndarray): An array of high gain median pedestals.
+        pedestal_median_lg (np.ndarray): An array of low gain median pedestals.
+        pedestal_std_hg (np.ndarray): An array of standard deviations of high gain pedestals.
+        pedestal_std_lg (np.ndarray): An array of standard deviations of low gain pedestals.
+    """
+
+    nsamples = Field(
+        type=np.uint8,
+        description="number of samples in the waveforms",
+    )
+
+    nevents = Field(
+        type=np.uint64,
+        description="number of events used to estimate the pedestals",
+    )
+
+    pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids")
+
+    ucts_timestamp_min = Field(
+        type=np.uint64,
+        description="minimum of the input events UCTS timestamps",
+    )
+
+    ucts_timestamp_max = Field(
+        type=np.uint64,
+        description="maximum of the input events UCTS timestamps",
+    )
+
+    pedestal_mean_hg = Field(
+        type=np.ndarray, dtype=np.float64, ndim=2, description="high gain mean pedestals"
+    )
+
+    pedestal_mean_lg = Field(
+        type=np.ndarray, dtype=np.float64, ndim=2, description="low gain mean pedestals"
+    )
+
+    pedestal_median_hg = Field(
+        type=np.ndarray, dtype=np.float64, ndim=2, description="high gain median pedestals"
+    )
+    pedestal_median_lg = Field(
+        type=np.ndarray, dtype=np.float64, ndim=2, description="low gain median pedestals"
+    )
+
+    pedestal_std_hg = Field(
+        type=np.ndarray, dtype=np.float64, ndim=2,
+        description="high gain pedestals standard deviations"
+    )
+    pedestal_std_lg = Field(
+        type=np.ndarray, dtype=np.float64, ndim=2,
+        description="low gain pedestals standard deviations"
+    )
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index 7d06a15b..194edcdf 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -14,10 +14,8 @@
     N_GAINS, N_PIXELS, N_SAMPLES,
     HIGH_GAIN, LOW_GAIN,
 )
-from ctapipe_io_nectarcam import time_from_unix_tai_ns
 from ctapipe.core.traits import Dict, Unicode
-from ctapipe.containers import PedestalContainer
-from ...data.container import merge_map_ArrayDataContainer
+from ...data.container import NectarCAMPedestalContainer, merge_map_ArrayDataContainer
 from ...utils import ComponentUtils
 from .waveformsComponent import WaveformsComponent
 
@@ -27,6 +25,7 @@
     "PedestalEstimationComponent",
 ]
 
+
 class PedestalEstimationComponent(NectarCAMComponent):
     """
     Component that computes calibration pedestal coefficients from raw data.
@@ -37,7 +36,7 @@ class PedestalEstimationComponent(NectarCAMComponent):
     #     read_only=True,
     # ).tag(config=True)
 
-    #not implemented yet, placeholder
+    # not implemented yet, placeholder
     filter_kwargs = Dict(
         default_value={},
         help="The kwargs to be pass to the waveform filter method",
@@ -48,7 +47,7 @@ class PedestalEstimationComponent(NectarCAMComponent):
     SubComponents = copy.deepcopy(NectarCAMComponent)
     SubComponents.default_value = [
         "WaveformsComponent",
-        #f"{PedestalFilterAlgorithm.default_value}",
+        # f"{PedestalFilterAlgorithm.default_value}",
     ]
     SubComponents.read_only = True
 
@@ -87,7 +86,7 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
         self._wfs_mask = None
 
     @staticmethod
-    def calculate_stats(waveformsContainers,wfs_mask,statistics):
+    def calculate_stats(waveformsContainers, wfs_mask, statistics):
         """
         Calculate statistics for the pedestals from a waveforms container.
 
@@ -111,7 +110,7 @@ def calculate_stats(waveformsContainers,wfs_mask,statistics):
         for stat in statistics:
             # Calculate the statistic along axis = 0, that is over events
             ped_stat_hg = getattr(np, stat)(
-                ma.masked_array(waveformsContainers.wfs_hg,wfs_mask),axis=0)
+                ma.masked_array(waveformsContainers.wfs_hg, wfs_mask), axis=0)
             ped_stat_lg = getattr(np, stat)(
                 ma.masked_array(waveformsContainers.wfs_lg, wfs_mask), axis=0)
 
@@ -144,7 +143,7 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
         if not (is_empty):
             # change this into something that creates a real mask
             # one mask for both HG and LG or separate? mask only on hg which is more sensitive
-            self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg),dtype=bool)
+            self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg), dtype=bool)
         else:
             pass
 
@@ -153,20 +152,22 @@ def finish(self, *args, **kwargs):
         # compute statistics for the pedestals
         # the statistic names must be valid numpy attributes
         statistics = ['mean', 'median', 'std']
-        #print(np.shape(self._waveformsContainers.wfs_hg))
-        ped_stats = self.calculate_stats(self._waveformsContainers,self._wfs_mask,statistics)
-
-        metadata = {} # store information about filtering method and params
-
-        output = PedestalContainer(
-            n_events = self._waveformsContainers.nevents,
-            sample_time= np.nan, # TODO : does this need to be filled?
-            sample_time_min= time_from_unix_tai_ns(self._waveformsContainers.ucts_timestamp.min()).value * u.s,
-            sample_time_max = time_from_unix_tai_ns(self._waveformsContainers.ucts_timestamp.max()).value * u.s,
-            charge_mean = ped_stats['mean'],
-            charge_median = ped_stats['median'],
-            charge_std = ped_stats['std'],
-            meta = metadata,
+        ped_stats = self.calculate_stats(self._waveformsContainers, self._wfs_mask, statistics)
+
+        # Fill output container
+        # metadata = {} # store information about filtering method and params
+        output = NectarCAMPedestalContainer(
+            nsamples=self._waveformsContainers.nsamples,
+            nevents= self._waveformsContainers.nevents,
+            pixels_id=self._waveformsContainers.pixels_id,
+            ucts_timestamp_min=self._waveformsContainers.ucts_timestamp.min(),
+            ucts_timestamp_max=self._waveformsContainers.ucts_timestamp.max(),
+            pedestal_mean_hg=ped_stats['mean'][HIGH_GAIN],
+            pedestal_mean_lg=ped_stats['mean'][LOW_GAIN],
+            pedestal_median_hg=ped_stats['median'][HIGH_GAIN],
+            pedestal_median_lg=ped_stats['median'][LOW_GAIN],
+            pedestal_std_hg=ped_stats['std'][HIGH_GAIN],
+            pedestal_std_lg=ped_stats['std'][LOW_GAIN],
         )
 
         return output

From a4f2f8016cffbce85fc3146b720f87e6f3f2ab12 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Wed, 14 Feb 2024 12:20:07 +0100
Subject: [PATCH 05/20] Minimal working tool/component

---
 .../makers/calibration/pedestalMakers.py      | 23 ++---
 .../makers/component/PedestalComponent.py     | 91 ++++++++++++-------
 .../user_scripts/ltibaldo/example_pedestal.py |  2 +-
 3 files changed, 65 insertions(+), 51 deletions(-)

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index cfe96b16..c252c7b9 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -5,7 +5,6 @@
 log.handlers = logging.getLogger("__main__").handlers
 
 import pathlib
-import numpy as np
 
 from ctapipe.core.traits import ComponentNameList, Bool
 
@@ -19,7 +18,6 @@
 
 
 class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
-
     name = "PedestalNectarCAMCalibrationTool"
 
     reload_events = Bool(
@@ -49,21 +47,16 @@ def _init_output_path(self):
             filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}{ext}"
 
         self.output_path = pathlib.Path(
-            f"{os.environ.get('NECTARCAMDATA','/tmp')}/PedestalEstimation/{filename}"
+            f"{os.environ.get('NECTARCAMDATA', '/tmp')}/PedestalEstimation/{filename}"
         )
 
     def start(
-        self,
-        max_events=np.inf,
-        # trigger_type: list = None,
-        restart_from_beginning: bool = False,
-        *args,
-        **kwargs,
+            self,
+            # trigger_type: list = None,
+            restart_from_beginning: bool = False,
+            *args,
+            **kwargs,
     ):
-        # need to implement waveform filter methods
-        # str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
-        #     self.extractor_kwargs
-        # )
 
         files = DataManagement.find_waveforms(
             run_number=self.run_number,
@@ -74,14 +67,12 @@ def start(
                 self.log.info(
                     f"{len(files)} computed waveforms files found with max_events > {self.max_events} for run {self.run_number}, reload waveforms from event loop"
                 )
-            print("Start parent class", super().name)
             super().start(
                 restart_from_beginning=restart_from_beginning,
                 *args,
                 **kwargs,
             )
         else:
-            print('Enter loop to read waveforms')
             self.log.info(f"reading waveforms from files {files[0]}")
             waveformsContainers = WaveformsContainer.from_hdf5(files[0])
             if isinstance(waveformsContainers, WaveformsContainer):
@@ -102,5 +93,3 @@ def start(
                 self.components[0]._waveformsContainers = merge_map_ArrayDataContainer(
                     waveformsContainers_merges_along_slices
                 )
-
-        print("Input wfs shape", np.shape(self.components[0]._waveformsContainers.wfs_hg))
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index 194edcdf..6e30ebbe 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -5,16 +5,14 @@
 log.handlers = logging.getLogger("__main__").handlers
 
 import copy
-import astropy.units as u
 import numpy.ma as ma
 
 from .core import NectarCAMComponent
 from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
 from ctapipe_io_nectarcam.constants import (
-    N_GAINS, N_PIXELS, N_SAMPLES,
-    HIGH_GAIN, LOW_GAIN,
+    N_GAINS, HIGH_GAIN, LOW_GAIN,
 )
-from ctapipe.core.traits import Dict, Unicode
+from ctapipe.core.traits import Dict, Integer
 from ...data.container import NectarCAMPedestalContainer, merge_map_ArrayDataContainer
 from ...utils import ComponentUtils
 from .waveformsComponent import WaveformsComponent
@@ -30,6 +28,20 @@ class PedestalEstimationComponent(NectarCAMComponent):
     """
     Component that computes calibration pedestal coefficients from raw data.
     """
+
+    ucts_tmin = Integer(
+        "ucts_tmin",
+        help="Minimum UCTS timestamp for events used in pedestal estimation",
+        default=-np.inf,
+        read_only=True,
+    ).tag(config=True)
+
+    ucts_tmax = Integer(
+        "ucts_tmax",
+        help="Maximum UCTS timestamp for events used in pedestal estimation",
+        default=np.inf,
+        read_only=True,
+    ).tag(config=True)
     # PedestalFilterAlgorithm = Unicode(
     #     "PedestalWaveformStdFilter",
     #     help="The waveform filter method",
@@ -81,9 +93,11 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
             *args,
             **waveformsComponent_kwargs,
         )
-        self._waveformsContainers = None
 
+        # initialize containers
+        self._waveformsContainers = None
         self._wfs_mask = None
+        self._ped_stats = {}
 
     @staticmethod
     def calculate_stats(waveformsContainers, wfs_mask, statistics):
@@ -127,47 +141,58 @@ def calculate_stats(waveformsContainers, wfs_mask, statistics):
 
     def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
 
+        # fill the waveform container looping over the events
         self.waveformsComponent(event=event, *args, **kwargs)
 
+    def finish(self, *args, **kwargs):
+
+        # Make sure that Waveforms container is properly filled
         is_empty = False
         if self._waveformsContainers is None:
             self._waveformsContainers = self.waveformsComponent.finish(*args, **kwargs)
             is_empty = self._waveformsContainers.is_empty()
             if is_empty:
-                log.warning("empty waveformsContainer in output")
+                log.warning("empty waveforms container, pedestals cannot be evaluated")
             else:
+                # container merging
                 self._waveformsContainers = merge_map_ArrayDataContainer(
                     self._waveformsContainers
                 )
 
-        if not (is_empty):
-            # change this into something that creates a real mask
-            # one mask for both HG and LG or separate? mask only on hg which is more sensitive
+        if not is_empty:
+            # Check that there are events for the computation of the pedestals
+            # Build mask to filter the waveforms
+            # mask based on the high gain channel that is most sensitive to signals
             self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg), dtype=bool)
-        else:
-            pass
 
-    def finish(self, *args, **kwargs):
-
-        # compute statistics for the pedestals
-        # the statistic names must be valid numpy attributes
-        statistics = ['mean', 'median', 'std']
-        ped_stats = self.calculate_stats(self._waveformsContainers, self._wfs_mask, statistics)
-
-        # Fill output container
-        # metadata = {} # store information about filtering method and params
-        output = NectarCAMPedestalContainer(
-            nsamples=self._waveformsContainers.nsamples,
-            nevents= self._waveformsContainers.nevents,
-            pixels_id=self._waveformsContainers.pixels_id,
-            ucts_timestamp_min=self._waveformsContainers.ucts_timestamp.min(),
-            ucts_timestamp_max=self._waveformsContainers.ucts_timestamp.max(),
-            pedestal_mean_hg=ped_stats['mean'][HIGH_GAIN],
-            pedestal_mean_lg=ped_stats['mean'][LOW_GAIN],
-            pedestal_median_hg=ped_stats['median'][HIGH_GAIN],
-            pedestal_median_lg=ped_stats['median'][LOW_GAIN],
-            pedestal_std_hg=ped_stats['std'][HIGH_GAIN],
-            pedestal_std_lg=ped_stats['std'][LOW_GAIN],
-        )
+            # apply time filter
+            # log.info(
+            #     f"Mask events outside the UCTS timestamp range {self.ucts_tmin}-{self.ucts_tmax}")
+            # t_mask = np.where((self._waveformsContainers.ucts_timestamps < self.ucts_tmin) |
+            #                   (self._waveformsContainers.ucts_timestamps > self.ucts_tmax))
+
+            # compute statistics for the pedestals
+            # the statistic names must be valid numpy attributes
+            statistics = ['mean', 'median', 'std']
+            self._ped_stats = self.calculate_stats(self._waveformsContainers, self._wfs_mask,
+                                                   statistics)
+
+            # Fill and return output container
+
+            # metadata = {} # store information about filtering method and params
+
+            output = NectarCAMPedestalContainer(
+                nsamples=self._waveformsContainers.nsamples,
+                nevents=self._waveformsContainers.nevents,
+                pixels_id=self._waveformsContainers.pixels_id,
+                ucts_timestamp_min=self._waveformsContainers.ucts_timestamp.min(),
+                ucts_timestamp_max=self._waveformsContainers.ucts_timestamp.max(),
+                pedestal_mean_hg=self._ped_stats['mean'][HIGH_GAIN],
+                pedestal_mean_lg=self._ped_stats['mean'][LOW_GAIN],
+                pedestal_median_hg=self._ped_stats['median'][HIGH_GAIN],
+                pedestal_median_lg=self._ped_stats['median'][LOW_GAIN],
+                pedestal_std_hg=self._ped_stats['std'][HIGH_GAIN],
+                pedestal_std_lg=self._ped_stats['std'][LOW_GAIN],
+            )
 
         return output
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index 82313652..fcfb0fe1 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -13,7 +13,7 @@
 )
 
 run_number = 3938
-max_events= 100
+max_events= 20
 outfile = "/Users/ltibaldo/tmp/test_pedestal/pedestal_{}.h5".format(run_number)
 
 tool = PedestalNectarCAMCalibrationTool(

From d1c8de00bf6443769b9c9e454d194b2e3afc5db9 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Wed, 14 Feb 2024 16:22:43 +0100
Subject: [PATCH 06/20] implemented time selection

---
 .../makers/component/PedestalComponent.py     | 146 ++++++++++--------
 .../user_scripts/ltibaldo/example_pedestal.py |   1 +
 2 files changed, 80 insertions(+), 67 deletions(-)

diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index 6e30ebbe..84401aca 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -9,9 +9,7 @@
 
 from .core import NectarCAMComponent
 from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
-from ctapipe_io_nectarcam.constants import (
-    N_GAINS, HIGH_GAIN, LOW_GAIN,
-)
+from ctapipe_io_nectarcam.constants import (N_GAINS, HIGH_GAIN, LOW_GAIN, )
 from ctapipe.core.traits import Dict, Integer
 from ...data.container import NectarCAMPedestalContainer, merge_map_ArrayDataContainer
 from ...utils import ComponentUtils
@@ -19,9 +17,7 @@
 
 import numpy as np
 
-__all__ = [
-    "PedestalEstimationComponent",
-]
+__all__ = ["PedestalEstimationComponent", ]
 
 
 class PedestalEstimationComponent(NectarCAMComponent):
@@ -29,47 +25,32 @@ class PedestalEstimationComponent(NectarCAMComponent):
     Component that computes calibration pedestal coefficients from raw data.
     """
 
-    ucts_tmin = Integer(
-        "ucts_tmin",
-        help="Minimum UCTS timestamp for events used in pedestal estimation",
-        default=-np.inf,
-        read_only=True,
-    ).tag(config=True)
-
-    ucts_tmax = Integer(
-        "ucts_tmax",
-        help="Maximum UCTS timestamp for events used in pedestal estimation",
-        default=np.inf,
-        read_only=True,
-    ).tag(config=True)
-    # PedestalFilterAlgorithm = Unicode(
-    #     "PedestalWaveformStdFilter",
-    #     help="The waveform filter method",
-    #     read_only=True,
-    # ).tag(config=True)
-
-    # not implemented yet, placeholder
-    filter_kwargs = Dict(
-        default_value={},
-        help="The kwargs to be pass to the waveform filter method",
-    ).tag(config=True)
+    ucts_tmin = Integer(None,
+        help="Minimum UCTS timestamp for events used in pedestal estimation", read_only=False,
+        allow_none=True, ).tag(config=True)
+
+    ucts_tmax = Integer(None,
+        help="Maximum UCTS timestamp for events used in pedestal estimation", read_only=False,
+        allow_none=True, ).tag(config=True)
+
+    filter_kwargs = Dict(default_value={},
+        help="The kwargs to be pass to the waveform filter method", ).tag(config=True)
 
     # add parameters for min max time
 
     SubComponents = copy.deepcopy(NectarCAMComponent)
-    SubComponents.default_value = [
-        "WaveformsComponent",
+    SubComponents.default_value = ["WaveformsComponent",
         # f"{PedestalFilterAlgorithm.default_value}",
     ]
     SubComponents.read_only = True
 
     def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
+
         waveformsComponent_kwargs = {}
         self._PedestalFilterAlgorithm_kwargs = {}
         other_kwargs = {}
         waveformsComponent_configurable_traits = ComponentUtils.get_configurable_traits(
-            WaveformsComponent
-        )
+            WaveformsComponent)
         # pedestalFilterAlgorithm_configurable_traits = ComponentUtils.get_configurable_traits(
         #     eval(self.PedestalFilterAlgorithm)
         # )
@@ -82,19 +63,12 @@ def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
             else:
                 other_kwargs[key] = kwargs[key]
 
-        super().__init__(
-            subarray=subarray, config=config, parent=parent, *args, **kwargs
-        )
+        super().__init__(subarray=subarray, config=config, parent=parent, *args, **kwargs)
 
-        self.waveformsComponent = WaveformsComponent(
-            subarray=subarray,
-            config=config,
-            parent=parent,
-            *args,
-            **waveformsComponent_kwargs,
-        )
+        self.waveformsComponent = WaveformsComponent(subarray=subarray, config=config,
+            parent=parent, *args, **waveformsComponent_kwargs, )
 
-        # initialize containers
+        # initialize members
         self._waveformsContainers = None
         self._wfs_mask = None
         self._ped_stats = {}
@@ -106,16 +80,16 @@ def calculate_stats(waveformsContainers, wfs_mask, statistics):
 
         Parameters
         ----------
-        waveformsContainers : `~nectarchain.data.container.WaveFormsContainer`
+        waveformsContainers : `~nectarchain.data.container.WaveformsContainer`
             Waveforms container
         wfs_mask : `numpy.ndarray`
             Mask to apply to exclude outliers with shape (n_pixels,n_samples)
         statistics : `list`
-            Names of the statistics (numpy attributes) to compute
+            Names of the statistics (numpy.ma attributes) to compute
 
         Returns
         ----------
-        statistics : `dict`
+        ped_stats : `dict`
             A dictionary containing 3D (n_chan,n_pixels,n_samples) arrays for each statistic
         """
 
@@ -123,9 +97,9 @@ def calculate_stats(waveformsContainers, wfs_mask, statistics):
 
         for stat in statistics:
             # Calculate the statistic along axis = 0, that is over events
-            ped_stat_hg = getattr(np, stat)(
+            ped_stat_hg = getattr(ma, stat)(
                 ma.masked_array(waveformsContainers.wfs_hg, wfs_mask), axis=0)
-            ped_stat_lg = getattr(np, stat)(
+            ped_stat_lg = getattr(ma, stat)(
                 ma.masked_array(waveformsContainers.wfs_lg, wfs_mask), axis=0)
 
             # Create a 3D array for the statistic
@@ -144,55 +118,93 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
         # fill the waveform container looping over the events
         self.waveformsComponent(event=event, *args, **kwargs)
 
+    def timestamp_mask(self, tmin, tmax):
+        """
+        Generates a mask to filter out events in the valid time interval
+
+        Returns
+        ----------
+        mask : `~numpy.ndarray`
+            A boolean array of shape (n_events,n_pixels,n_samples) that identifies waveforms to be masked
+        """
+        # Check if time filtering affects the data
+        if tmin > self._waveformsContainers.ucts_timestamp.min() or tmax < self._waveformsContainers.ucts_timestamp.max():
+            log.info(
+                f"Apply time interval selection: UCTS timestamps in interval {tmin}-{tmax}")
+            # Define the mask on UCTS timestamps
+            t_mask = ((self._waveformsContainers.ucts_timestamp < tmin) | (
+                        self._waveformsContainers.ucts_timestamp > tmax))
+            # Log information
+            nonzero = np.count_nonzero(t_mask)
+            log.info(f"{len(t_mask) - nonzero}/{len(t_mask)} events pass time selection.")
+
+            # Create waveforms mask to apply time selection
+            new_mask = np.logical_or(self._wfs_mask, t_mask[:, np.newaxis, np.newaxis])
+
+            # Put some information in log to say how many events pass time selection
+        else:
+            log.info(
+                f"The entire time interval will be used: UCTS timestamps in {tmin}-{tmax}")
+            new_mask = self._wfs_mask
+
+        # Return waveforms mask
+        return new_mask
+
     def finish(self, *args, **kwargs):
 
-        # Make sure that Waveforms container is properly filled
+        # Make sure that waveforms container is properly filled
         is_empty = False
         if self._waveformsContainers is None:
             self._waveformsContainers = self.waveformsComponent.finish(*args, **kwargs)
             is_empty = self._waveformsContainers.is_empty()
             if is_empty:
                 log.warning("empty waveforms container, pedestals cannot be evaluated")
+
+                # container with no results
+                output = NectarCAMPedestalContainer(
+                    nsamples=self._waveformsContainers.nsamples,
+                    nevents=self._waveformsContainers.nevents,
+                    pixels_id=self._waveformsContainers.pixels_id, )
             else:
                 # container merging
                 self._waveformsContainers = merge_map_ArrayDataContainer(
-                    self._waveformsContainers
-                )
+                    self._waveformsContainers)
 
         if not is_empty:
-            # Check that there are events for the computation of the pedestals
             # Build mask to filter the waveforms
-            # mask based on the high gain channel that is most sensitive to signals
+            # Mask based on the high gain channel that is most sensitive to signals
+            # Initialize empty mask
             self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg), dtype=bool)
 
-            # apply time filter
-            # log.info(
-            #     f"Mask events outside the UCTS timestamp range {self.ucts_tmin}-{self.ucts_tmax}")
-            # t_mask = np.where((self._waveformsContainers.ucts_timestamps < self.ucts_tmin) |
-            #                   (self._waveformsContainers.ucts_timestamps > self.ucts_tmax))
+            # Time mask
+            # set the minimum time
+            tmin = np.maximum(self.ucts_tmin or self._waveformsContainers.ucts_timestamp.min(),
+                self._waveformsContainers.ucts_timestamp.min())
+            # set the maximum time
+            tmax = np.minimum(self.ucts_tmax or self._waveformsContainers.ucts_timestamp.max(),
+                self._waveformsContainers.ucts_timestamp.max())
+            # Build mask
+            self._wfs_mask = self.timestamp_mask(tmin, tmax)
 
             # compute statistics for the pedestals
-            # the statistic names must be valid numpy attributes
+            # the statistic names must be valid numpy.ma attributes
             statistics = ['mean', 'median', 'std']
             self._ped_stats = self.calculate_stats(self._waveformsContainers, self._wfs_mask,
                                                    statistics)
 
             # Fill and return output container
-
             # metadata = {} # store information about filtering method and params
-
             output = NectarCAMPedestalContainer(
                 nsamples=self._waveformsContainers.nsamples,
                 nevents=self._waveformsContainers.nevents,
                 pixels_id=self._waveformsContainers.pixels_id,
-                ucts_timestamp_min=self._waveformsContainers.ucts_timestamp.min(),
-                ucts_timestamp_max=self._waveformsContainers.ucts_timestamp.max(),
+                ucts_timestamp_min=np.uint64(tmin),
+                ucts_timestamp_max=np.uint64(tmax),
                 pedestal_mean_hg=self._ped_stats['mean'][HIGH_GAIN],
                 pedestal_mean_lg=self._ped_stats['mean'][LOW_GAIN],
                 pedestal_median_hg=self._ped_stats['median'][HIGH_GAIN],
                 pedestal_median_lg=self._ped_stats['median'][LOW_GAIN],
                 pedestal_std_hg=self._ped_stats['std'][HIGH_GAIN],
-                pedestal_std_lg=self._ped_stats['std'][LOW_GAIN],
-            )
+                pedestal_std_lg=self._ped_stats['std'][LOW_GAIN], )
 
         return output
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index fcfb0fe1..9913e189 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -24,6 +24,7 @@
     reload_events=False,
     output_path=outfile,
     overwrite = True,
+    ucts_tmin = 1674462932650000000,
 )
 
 tool.initialize()

From 01bd62765fb2b4ab66b92533f4fe3e0e67a7330b Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Thu, 15 Feb 2024 13:03:58 +0100
Subject: [PATCH 07/20] Implemented filtering based on waveforms standard
 deviation

---
 .../makers/component/PedestalComponent.py     | 138 +++++++++++++-----
 .../user_scripts/ltibaldo/example_pedestal.py |   1 +
 2 files changed, 101 insertions(+), 38 deletions(-)

diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index 84401aca..d644a8f4 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -9,8 +9,8 @@
 
 from .core import NectarCAMComponent
 from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
-from ctapipe_io_nectarcam.constants import (N_GAINS, HIGH_GAIN, LOW_GAIN, )
-from ctapipe.core.traits import Dict, Integer
+from ctapipe_io_nectarcam.constants import N_GAINS, HIGH_GAIN, LOW_GAIN
+from ctapipe.core.traits import Integer, Unicode, Float
 from ...data.container import NectarCAMPedestalContainer, merge_map_ArrayDataContainer
 from ...utils import ComponentUtils
 from .waveformsComponent import WaveformsComponent
@@ -26,47 +26,54 @@ class PedestalEstimationComponent(NectarCAMComponent):
     """
 
     ucts_tmin = Integer(None,
-        help="Minimum UCTS timestamp for events used in pedestal estimation", read_only=False,
-        allow_none=True, ).tag(config=True)
+                        help="Minimum UCTS timestamp for events used in pedestal estimation",
+                        read_only=False,
+                        allow_none=True, ).tag(config=True)
 
     ucts_tmax = Integer(None,
-        help="Maximum UCTS timestamp for events used in pedestal estimation", read_only=False,
-        allow_none=True, ).tag(config=True)
-
-    filter_kwargs = Dict(default_value={},
-        help="The kwargs to be pass to the waveform filter method", ).tag(config=True)
-
-    # add parameters for min max time
+                        help="Maximum UCTS timestamp for events used in pedestal estimation",
+                        read_only=False,
+                        allow_none=True, ).tag(config=True)
+
+    filter_method = Unicode(
+        None,
+        help="The waveforms filter method to be used.\n"
+             "Inplemented methods: WaveformsStdFilter (standard deviation of waveforms),\n"
+             "                     ChargeDistributionFilter (charge distribution).",
+        read_only=False,
+        allow_none=True,
+    ).tag(config=True)
+
+    wfs_std_threshold = Float(
+        4,
+        help="Threshold of waveforms standard deviation in ADC counts above which a waveform is excluded from pedestal computation.",
+        read_only=False,
+    )
 
     SubComponents = copy.deepcopy(NectarCAMComponent)
     SubComponents.default_value = ["WaveformsComponent",
-        # f"{PedestalFilterAlgorithm.default_value}",
-    ]
+                                   # f"{PedestalFilterAlgorithm.default_value}",
+                                   ]
     SubComponents.read_only = True
 
     def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
 
         waveformsComponent_kwargs = {}
-        self._PedestalFilterAlgorithm_kwargs = {}
         other_kwargs = {}
         waveformsComponent_configurable_traits = ComponentUtils.get_configurable_traits(
             WaveformsComponent)
-        # pedestalFilterAlgorithm_configurable_traits = ComponentUtils.get_configurable_traits(
-        #     eval(self.PedestalFilterAlgorithm)
-        # )
 
         for key in kwargs.keys():
             if key in waveformsComponent_configurable_traits.keys():
                 waveformsComponent_kwargs[key] = kwargs[key]
-            # elif key in waveformsComponent_configurable_traits.keys():
-            #     self._SPEfitalgorithm_kwargs[key] = kwargs[key]
             else:
                 other_kwargs[key] = kwargs[key]
 
         super().__init__(subarray=subarray, config=config, parent=parent, *args, **kwargs)
 
         self.waveformsComponent = WaveformsComponent(subarray=subarray, config=config,
-            parent=parent, *args, **waveformsComponent_kwargs, )
+                                                     parent=parent, *args,
+                                                     **waveformsComponent_kwargs, )
 
         # initialize members
         self._waveformsContainers = None
@@ -120,7 +127,14 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
 
     def timestamp_mask(self, tmin, tmax):
         """
-        Generates a mask to filter out events in the valid time interval
+        Generates a mask to filter waveforms outside the required time interval
+
+        Parameters
+        ----------
+        tmin : int
+            Minimum time of the required interval
+        tmax : int
+            Maximum time of the required interval
 
         Returns
         ----------
@@ -133,15 +147,12 @@ def timestamp_mask(self, tmin, tmax):
                 f"Apply time interval selection: UCTS timestamps in interval {tmin}-{tmax}")
             # Define the mask on UCTS timestamps
             t_mask = ((self._waveformsContainers.ucts_timestamp < tmin) | (
-                        self._waveformsContainers.ucts_timestamp > tmax))
+                    self._waveformsContainers.ucts_timestamp > tmax))
             # Log information
             nonzero = np.count_nonzero(t_mask)
-            log.info(f"{len(t_mask) - nonzero}/{len(t_mask)} events pass time selection.")
-
+            log.info(f"{len(t_mask) - nonzero}/{len(t_mask)} waveforms pass time selection.")
             # Create waveforms mask to apply time selection
             new_mask = np.logical_or(self._wfs_mask, t_mask[:, np.newaxis, np.newaxis])
-
-            # Put some information in log to say how many events pass time selection
         else:
             log.info(
                 f"The entire time interval will be used: UCTS timestamps in {tmin}-{tmax}")
@@ -150,6 +161,44 @@ def timestamp_mask(self, tmin, tmax):
         # Return waveforms mask
         return new_mask
 
+    def waveformsStdFilter_mask(self, threshold):
+        """
+        Generates a mask to filter waveforms that have a standard deviation above a threshold.
+        This option is effective for dark room verification data.
+
+        Parameters
+        ----------
+        threshold : float
+            Waveform standard deviation (in ADC counts) above which the waveform is filtered out
+
+        Returns
+        ----------
+        mask : `~numpy.ndarray`
+            A boolean array of shape (n_events,n_pixels,n_samples) that identifies waveforms to be masked
+        """
+
+        # Log
+        log.info(f"apply {self.filter_method} method with threshold = {threshold} ADC counts")
+
+        # For each event and pixel calculate the waveform std
+        wfs_std = np.std(self._waveformsContainers.wfs_hg, axis=2)
+
+        # Mask events/pixels that exceed the threshold
+        std_mask = (wfs_std > threshold)
+        # Log information
+        # number of masked waveforms
+        nonzero = np.count_nonzero(std_mask)
+        # number of total waveforms
+        tot_wfs = self._waveformsContainers.nevents * self._waveformsContainers.npixels
+        # fraction of waveforms masked (%)
+        frac = 100 * nonzero / tot_wfs
+        log.info(f"{frac:.2f}% of the waveforms filtered out based on standard deviation.")
+
+        # Create waveforms mask to apply time selection
+        new_mask = np.logical_or(self._wfs_mask, std_mask[:, :, np.newaxis])
+
+        return new_mask
+
     def finish(self, *args, **kwargs):
 
         # Make sure that waveforms container is properly filled
@@ -159,12 +208,8 @@ def finish(self, *args, **kwargs):
             is_empty = self._waveformsContainers.is_empty()
             if is_empty:
                 log.warning("empty waveforms container, pedestals cannot be evaluated")
-
                 # container with no results
-                output = NectarCAMPedestalContainer(
-                    nsamples=self._waveformsContainers.nsamples,
-                    nevents=self._waveformsContainers.nevents,
-                    pixels_id=self._waveformsContainers.pixels_id, )
+                return None
             else:
                 # container merging
                 self._waveformsContainers = merge_map_ArrayDataContainer(
@@ -174,22 +219,39 @@ def finish(self, *args, **kwargs):
             # Build mask to filter the waveforms
             # Mask based on the high gain channel that is most sensitive to signals
             # Initialize empty mask
-            self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg), dtype=bool)
+            self._wfs_mask = np.zeros(np.shape(self._waveformsContainers.wfs_hg),
+                                      dtype=bool)
 
-            # Time mask
+            # Time selection
             # set the minimum time
-            tmin = np.maximum(self.ucts_tmin or self._waveformsContainers.ucts_timestamp.min(),
+            print('time filter')
+            tmin = np.maximum(
+                self.ucts_tmin or self._waveformsContainers.ucts_timestamp.min(),
                 self._waveformsContainers.ucts_timestamp.min())
             # set the maximum time
-            tmax = np.minimum(self.ucts_tmax or self._waveformsContainers.ucts_timestamp.max(),
+            tmax = np.minimum(
+                self.ucts_tmax or self._waveformsContainers.ucts_timestamp.max(),
                 self._waveformsContainers.ucts_timestamp.max())
-            # Build mask
+            # Add time selection to mask
             self._wfs_mask = self.timestamp_mask(tmin, tmax)
 
+            # Filter Waveforms
+            if self.filter_method is None:
+                log.info('no filtering applied to waveforms')
+            elif self.filter_method == 'WaveformsStdFilter':
+                self._wfs_mask = self.waveformsStdFilter_mask(self.wfs_std_threshold)
+            elif self.filter_method == 'ChargeDistributionFilter':
+                log.info(f"method {self.filter_method} not implemented yet")
+            else:
+                log.warning(
+                    f"required filtering method {self.filter_method} not available")
+                log.warning("no filtering applied to waveforms")
+
             # compute statistics for the pedestals
             # the statistic names must be valid numpy.ma attributes
             statistics = ['mean', 'median', 'std']
-            self._ped_stats = self.calculate_stats(self._waveformsContainers, self._wfs_mask,
+            self._ped_stats = self.calculate_stats(self._waveformsContainers,
+                                                   self._wfs_mask,
                                                    statistics)
 
             # Fill and return output container
@@ -207,4 +269,4 @@ def finish(self, *args, **kwargs):
                 pedestal_std_hg=self._ped_stats['std'][HIGH_GAIN],
                 pedestal_std_lg=self._ped_stats['std'][LOW_GAIN], )
 
-        return output
+            return output
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index 9913e189..6a2d3adb 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -25,6 +25,7 @@
     output_path=outfile,
     overwrite = True,
     ucts_tmin = 1674462932650000000,
+    filter_method = "WaveformsStdFilter",
 )
 
 tool.initialize()

From 17c8b93d008bc16ef934e5aea2490e0eb1aba1d6 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Thu, 15 Feb 2024 17:32:58 +0100
Subject: [PATCH 08/20] Add possibility to filter based on charge distribution

---
 .../makers/calibration/pedestalMakers.py      |  51 ++++++-
 .../makers/component/PedestalComponent.py     | 144 +++++++++++++++---
 .../user_scripts/ltibaldo/example_pedestal.py |   4 +-
 .../user_scripts/ltibaldo/extract_charges.py  |  27 ++++
 4 files changed, 195 insertions(+), 31 deletions(-)
 create mode 100644 src/nectarchain/user_scripts/ltibaldo/extract_charges.py

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index c252c7b9..19d61693 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -9,7 +9,11 @@
 from ctapipe.core.traits import ComponentNameList, Bool
 
 from .core import NectarCAMCalibrationTool
-from ...data.container import WaveformsContainer, WaveformsContainers
+from ...data.container import (WaveformsContainer,
+                               WaveformsContainers,
+                               ChargesContainer,
+                               ChargesContainers,
+                               )
 from ...data.management import DataManagement
 from ...data.container.core import merge_map_ArrayDataContainer
 from ..component import ArrayDataComponent, NectarCAMComponent
@@ -34,9 +38,6 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
     def _init_output_path(self):
-        # str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(
-        #     self.extractor_kwargs
-        # )
         if self.events_per_slice is None:
             ext = ".h5"
         else:
@@ -58,10 +59,13 @@ def start(
             **kwargs,
     ):
 
+        # Waveforms
+        # Search files
         files = DataManagement.find_waveforms(
             run_number=self.run_number,
             max_events=self.max_events,
         )
+        # Load Waveforms from events
         if self.reload_events or len(files) != 1:
             if len(files) != 1:
                 self.log.info(
@@ -72,6 +76,7 @@ def start(
                 *args,
                 **kwargs,
             )
+        # Load waveforms from file
         else:
             self.log.info(f"reading waveforms from files {files[0]}")
             waveformsContainers = WaveformsContainer.from_hdf5(files[0])
@@ -93,3 +98,41 @@ def start(
                 self.components[0]._waveformsContainers = merge_map_ArrayDataContainer(
                     waveformsContainers_merges_along_slices
                 )
+
+        # Charges
+        if self.filter_method == 'ChargeDistributionFilter':
+            # Search on disk
+            charge_files = DataManagement.find_charges(
+                run_number=self.run_number,
+                max_events=self.max_events,
+                #use the standard extraction method, that is FullWaveformSum
+            )
+            # Compute charges from waveforms
+            if self.reload_events or len(charge_files) != 1:
+                if len(charge_files) != 1:
+                    self.log.info(
+                        f"{len(files)} computed charges files found with max_events > {self.max_events} for run {self.run_number}"
+                        f"charges will be computed from waveforms"
+                    )
+            else:
+                self.log.info(f"reading charges from files {charge_files[0]}")
+                chargesContainers = ChargesContainer.from_hdf5(charge_files[0])
+                if isinstance(chargesContainers, ChargesContainer):
+                    self.components[0]._chargesContainers = chargesContainers
+                elif isinstance(chargesContainers, ChargesContainers):
+                    self.log.debug("merging along TriggerType")
+                    self.components[0]._chargesContainers = merge_map_ArrayDataContainer(
+                        chargesContainers
+                    )
+                else:
+                    self.log.debug("merging along slices")
+                    chargesContainers_merges_along_slices = (
+                        ArrayDataComponent.merge_along_slices(
+                            containers_generator=chargesContainers
+                        )
+                    )
+                    self.log.debug("merging along TriggerType")
+                    self.components[0]._chargesContainers = merge_map_ArrayDataContainer(
+                        chargesContainers_merges_along_slices
+                    )
+
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index d644a8f4..e10535e1 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -10,10 +10,11 @@
 from .core import NectarCAMComponent
 from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
 from ctapipe_io_nectarcam.constants import N_GAINS, HIGH_GAIN, LOW_GAIN
-from ctapipe.core.traits import Integer, Unicode, Float
+from ctapipe.core.traits import Integer, Unicode, Float, Dict
 from ...data.container import NectarCAMPedestalContainer, merge_map_ArrayDataContainer
 from ...utils import ComponentUtils
 from .waveformsComponent import WaveformsComponent
+from .chargesComponent import ChargesComponent
 
 import numpy as np
 
@@ -27,12 +28,10 @@ class PedestalEstimationComponent(NectarCAMComponent):
 
     ucts_tmin = Integer(None,
                         help="Minimum UCTS timestamp for events used in pedestal estimation",
-                        read_only=False,
                         allow_none=True, ).tag(config=True)
 
     ucts_tmax = Integer(None,
                         help="Maximum UCTS timestamp for events used in pedestal estimation",
-                        read_only=False,
                         allow_none=True, ).tag(config=True)
 
     filter_method = Unicode(
@@ -45,41 +44,60 @@ class PedestalEstimationComponent(NectarCAMComponent):
     ).tag(config=True)
 
     wfs_std_threshold = Float(
-        4,
+        4.,
         help="Threshold of waveforms standard deviation in ADC counts above which a waveform is excluded from pedestal computation.",
-        read_only=False,
-    )
+    ).tag(config=True)
+
+    charge_sigma_high_thr = Float(
+        3.,
+        help="Threshold in charge distribution (number of sigmas above mean) beyond which a waveform is excluded from pedestal computation.",
+    ).tag(config=True)
 
-    SubComponents = copy.deepcopy(NectarCAMComponent)
+    charge_sigma_low_thr = Float(
+        3.,
+        help="Threshold in charge distribution (number of sigmas below mean) beyond which a waveform is excluded from pedestal computation.",
+    ).tag(config=True)
+
+    # I do not understand why but the ChargesComponents traits are not loaded
+    # FIXME
+    method = Unicode(
+        default_value="FullWaveformSum",
+        help="the charge extraction method",
+    ).tag(config=True)
+
+    extractor_kwargs = Dict(
+        default_value={},
+        help="The kwargs to be pass to the charge extractor method",
+    ).tag(config=True)
+
+    SubComponents = copy.deepcopy(NectarCAMComponent.SubComponents)
     SubComponents.default_value = ["WaveformsComponent",
-                                   # f"{PedestalFilterAlgorithm.default_value}",
-                                   ]
+                                   "ChargesComponent"]
     SubComponents.read_only = True
 
     def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
 
-        waveformsComponent_kwargs = {}
-        other_kwargs = {}
-        waveformsComponent_configurable_traits = ComponentUtils.get_configurable_traits(
-            WaveformsComponent)
-
-        for key in kwargs.keys():
-            if key in waveformsComponent_configurable_traits.keys():
-                waveformsComponent_kwargs[key] = kwargs[key]
-            else:
-                other_kwargs[key] = kwargs[key]
-
         super().__init__(subarray=subarray, config=config, parent=parent, *args, **kwargs)
 
-        self.waveformsComponent = WaveformsComponent(subarray=subarray, config=config,
-                                                     parent=parent, *args,
-                                                     **waveformsComponent_kwargs, )
-
         # initialize members
         self._waveformsContainers = None
+        self._chargesContainers = None
         self._wfs_mask = None
         self._ped_stats = {}
 
+        # initialize waveforms component
+        waveformsComponent_kwargs = {}
+        waveformsComponent_configurable_traits = ComponentUtils.get_configurable_traits(
+            WaveformsComponent)
+        for key in kwargs.keys():
+            if key in waveformsComponent_configurable_traits.keys():
+                waveformsComponent_kwargs[key] = kwargs[key]
+        self.waveformsComponent = WaveformsComponent(
+            subarray=subarray,
+            config=config,
+            parent=parent, *args,
+            **waveformsComponent_kwargs, )
+
     @staticmethod
     def calculate_stats(waveformsContainers, wfs_mask, statistics):
         """
@@ -199,6 +217,62 @@ def waveformsStdFilter_mask(self, threshold):
 
         return new_mask
 
+    def chargeDistributionFilter_mask(self, sigma_low, sigma_high):
+        """
+        Generates a mask to filter waveforms that have a charge in the tails of the distribution.
+        This option is useful for data with NSB.
+
+        Parameters
+        ----------
+        sigma_low : float
+            Number of standard deviation below mean charge beyond which waveforms are filtered out
+        sigma_high : float
+            Number of standard deviation above mean charge beyond which waveforms are filtered out
+
+        Returns
+        ----------
+        mask : `~numpy.ndarray`
+            A boolean array of shape (n_events,n_pixels,n_samples) that identifies waveforms to be masked
+        """
+
+        # Log
+        log.info(f"apply {self.filter_method} method to filter waveforms with charge outside "
+                 f"the interval -{sigma_low}-{sigma_high} standard deviations around the mean value")
+
+        # Check if waveforms and charges have the same shape
+        wfs_shape = np.shape(self._waveformsContainers.wfs_hg)
+        charges_shape = np.shape(self._chargesContainers.charges_hg)
+        if wfs_shape[0] == charges_shape[0] and wfs_shape[1] == charges_shape[1]:
+
+            # For each event and pixel calculate the charge mean and std over all events
+            charge_mean = np.mean(self._chargesContainers.charges_hg, axis=0)
+            charge_std = np.std(self._chargesContainers.charges_hg, axis=0)
+
+            # Mask events/pixels that are outside the core of the distribution
+            low_threshold = charge_mean - sigma_low * charge_std
+            low_mask = (self._chargesContainers.charges_hg < low_threshold[np.newaxis, :])
+            high_threshold = charge_mean + sigma_high * charge_std
+            high_mask = (self._chargesContainers.charges_hg > high_threshold[np.newaxis, :])
+            charge_mask = np.logical_or(low_mask, high_mask)
+            # Log information
+            # number of masked waveforms
+            nonzero = np.count_nonzero(charge_mask)
+            # number of total waveforms
+            tot_wfs = self._waveformsContainers.nevents * self._waveformsContainers.npixels
+            # fraction of waveforms masked (%)
+            frac = 100 * nonzero / tot_wfs
+            log.info(
+                f"{frac:.2f}% of the waveforms filtered out based on charge distribution.")
+
+            # Create waveforms mask to apply time selection
+            new_mask = np.logical_or(self._wfs_mask, charge_mask[:, :, np.newaxis])
+        else:
+            log.warning(
+                "Waveforms and charges have incompatible shapes. No filtering applied.")
+            new_mask = self._wfs_mask
+
+        return new_mask
+
     def finish(self, *args, **kwargs):
 
         # Make sure that waveforms container is properly filled
@@ -216,6 +290,25 @@ def finish(self, *args, **kwargs):
                     self._waveformsContainers)
 
         if not is_empty:
+
+            # If we want to filter based on charges distribution
+            # make sure that the charge distribution container is filled
+            if self.filter_method == "ChargeDistributionFilter" and \
+                    self._chargesContainers is None:
+                log.debug("Compute charges from waveforms")
+                chargesComponent_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]
+                self._chargesContainers = ChargesComponent.create_from_waveforms(
+                    waveformsContainer=self._waveformsContainers,
+                    subarray=self.subarray,
+                    config=self.config,
+                    parent=self.parent, *args,
+                    **chargesComponent_kwargs, )
+
             # Build mask to filter the waveforms
             # Mask based on the high gain channel that is most sensitive to signals
             # Initialize empty mask
@@ -241,7 +334,8 @@ def finish(self, *args, **kwargs):
             elif self.filter_method == 'WaveformsStdFilter':
                 self._wfs_mask = self.waveformsStdFilter_mask(self.wfs_std_threshold)
             elif self.filter_method == 'ChargeDistributionFilter':
-                log.info(f"method {self.filter_method} not implemented yet")
+                self._wfs_mask = self.chargeDistributionFilter_mask(self.charge_sigma_low_thr,
+                                                                    self.charge_sigma_high_thr)
             else:
                 log.warning(
                     f"required filtering method {self.filter_method} not available")
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index 6a2d3adb..fcad4d2b 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -13,7 +13,7 @@
 )
 
 run_number = 3938
-max_events= 20
+max_events= 100
 outfile = "/Users/ltibaldo/tmp/test_pedestal/pedestal_{}.h5".format(run_number)
 
 tool = PedestalNectarCAMCalibrationTool(
@@ -25,7 +25,7 @@
     output_path=outfile,
     overwrite = True,
     ucts_tmin = 1674462932650000000,
-    filter_method = "WaveformsStdFilter",
+    filter_method = "ChargeDistributionFilter",
 )
 
 tool.initialize()
diff --git a/src/nectarchain/user_scripts/ltibaldo/extract_charges.py b/src/nectarchain/user_scripts/ltibaldo/extract_charges.py
new file mode 100644
index 00000000..b942f5ce
--- /dev/null
+++ b/src/nectarchain/user_scripts/ltibaldo/extract_charges.py
@@ -0,0 +1,27 @@
+import logging
+import os
+import pathlib
+
+logging.basicConfig(
+    format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO
+)
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger("__main__").handlers
+
+from nectarchain.makers import ChargesNectarCAMCalibrationTool
+
+run_number = 3938
+max_events = 100
+
+tool = ChargesNectarCAMCalibrationTool(
+    progress_bar=True,
+    method = "FullWaveformSum",
+    run_number=run_number,
+    max_events=max_events,
+    log_level=20,
+)
+
+tool.initialize()
+tool.setup()
+tool.start()
+tool.finish(return_output_component=False)
\ No newline at end of file

From b82d08e9947f77d14915f4227f466d1c8c14b3e3 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 5 Mar 2024 18:14:24 +0100
Subject: [PATCH 09/20] Restructured to use events_per_slice to reduce memory
 load

---
 .../makers/calibration/pedestalMakers.py      | 111 +-----------------
 .../makers/component/PedestalComponent.py     |  78 ++++++------
 .../user_scripts/ltibaldo/example_pedestal.py |   8 +-
 .../ltibaldo/extract_waveforms.py             |  24 ++++
 4 files changed, 74 insertions(+), 147 deletions(-)
 create mode 100644 src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index 19d61693..599a52bc 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -6,17 +6,10 @@
 
 import pathlib
 
-from ctapipe.core.traits import ComponentNameList, Bool
+from ctapipe.core.traits import ComponentNameList
 
 from .core import NectarCAMCalibrationTool
-from ...data.container import (WaveformsContainer,
-                               WaveformsContainers,
-                               ChargesContainer,
-                               ChargesContainers,
-                               )
-from ...data.management import DataManagement
-from ...data.container.core import merge_map_ArrayDataContainer
-from ..component import ArrayDataComponent, NectarCAMComponent
+from ..component import NectarCAMComponent
 
 __all__ = ["PedestalNectarCAMCalibrationTool"]
 
@@ -24,15 +17,10 @@
 class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
     name = "PedestalNectarCAMCalibrationTool"
 
-    reload_events = Bool(
-        default_value=False, help="Reload the waveforms from raw data"
-    )
-
-    componentsList = ComponentNameList(
-        NectarCAMComponent,
+    componentsList = ComponentNameList(NectarCAMComponent,
         default_value=["PedestalEstimationComponent"],
-        help="List of Component names to be applied, the order will be respected",
-    ).tag(config=True)
+        help="List of Component names to be applied, the order will be respected", ).tag(
+        config=True)
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -48,91 +36,4 @@ def _init_output_path(self):
             filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}{ext}"
 
         self.output_path = pathlib.Path(
-            f"{os.environ.get('NECTARCAMDATA', '/tmp')}/PedestalEstimation/{filename}"
-        )
-
-    def start(
-            self,
-            # trigger_type: list = None,
-            restart_from_beginning: bool = False,
-            *args,
-            **kwargs,
-    ):
-
-        # Waveforms
-        # Search files
-        files = DataManagement.find_waveforms(
-            run_number=self.run_number,
-            max_events=self.max_events,
-        )
-        # Load Waveforms from events
-        if self.reload_events or len(files) != 1:
-            if len(files) != 1:
-                self.log.info(
-                    f"{len(files)} computed waveforms files found with max_events > {self.max_events} for run {self.run_number}, reload waveforms from event loop"
-                )
-            super().start(
-                restart_from_beginning=restart_from_beginning,
-                *args,
-                **kwargs,
-            )
-        # Load waveforms from file
-        else:
-            self.log.info(f"reading waveforms from files {files[0]}")
-            waveformsContainers = WaveformsContainer.from_hdf5(files[0])
-            if isinstance(waveformsContainers, WaveformsContainer):
-                self.components[0]._waveformsContainers = waveformsContainers
-            elif isinstance(waveformsContainers, WaveformsContainers):
-                self.log.debug("merging along TriggerType")
-                self.components[0]._waveformsContainers = merge_map_ArrayDataContainer(
-                    waveformsContainers
-                )
-            else:
-                self.log.debug("merging along slices")
-                waveformsContainers_merges_along_slices = (
-                    ArrayDataComponent.merge_along_slices(
-                        containers_generator=waveformsContainers
-                    )
-                )
-                self.log.debug("merging along TriggerType")
-                self.components[0]._waveformsContainers = merge_map_ArrayDataContainer(
-                    waveformsContainers_merges_along_slices
-                )
-
-        # Charges
-        if self.filter_method == 'ChargeDistributionFilter':
-            # Search on disk
-            charge_files = DataManagement.find_charges(
-                run_number=self.run_number,
-                max_events=self.max_events,
-                #use the standard extraction method, that is FullWaveformSum
-            )
-            # Compute charges from waveforms
-            if self.reload_events or len(charge_files) != 1:
-                if len(charge_files) != 1:
-                    self.log.info(
-                        f"{len(files)} computed charges files found with max_events > {self.max_events} for run {self.run_number}"
-                        f"charges will be computed from waveforms"
-                    )
-            else:
-                self.log.info(f"reading charges from files {charge_files[0]}")
-                chargesContainers = ChargesContainer.from_hdf5(charge_files[0])
-                if isinstance(chargesContainers, ChargesContainer):
-                    self.components[0]._chargesContainers = chargesContainers
-                elif isinstance(chargesContainers, ChargesContainers):
-                    self.log.debug("merging along TriggerType")
-                    self.components[0]._chargesContainers = merge_map_ArrayDataContainer(
-                        chargesContainers
-                    )
-                else:
-                    self.log.debug("merging along slices")
-                    chargesContainers_merges_along_slices = (
-                        ArrayDataComponent.merge_along_slices(
-                            containers_generator=chargesContainers
-                        )
-                    )
-                    self.log.debug("merging along TriggerType")
-                    self.components[0]._chargesContainers = merge_map_ArrayDataContainer(
-                        chargesContainers_merges_along_slices
-                    )
-
+            f"{os.environ.get('NECTARCAMDATA', '/tmp')}/PedestalEstimation/{filename}")
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index e10535e1..d0150afb 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -8,10 +8,11 @@
 import numpy.ma as ma
 
 from .core import NectarCAMComponent
+from ctapipe.containers import EventType
 from ctapipe_io_nectarcam.containers import NectarCAMDataContainer
 from ctapipe_io_nectarcam.constants import N_GAINS, HIGH_GAIN, LOW_GAIN
 from ctapipe.core.traits import Integer, Unicode, Float, Dict
-from ...data.container import NectarCAMPedestalContainer, merge_map_ArrayDataContainer
+from ...data.container import NectarCAMPedestalContainer
 from ...utils import ComponentUtils
 from .waveformsComponent import WaveformsComponent
 from .chargesComponent import ChargesComponent
@@ -140,8 +141,11 @@ def calculate_stats(waveformsContainers, wfs_mask, statistics):
 
     def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
 
-        # fill the waveform container looping over the events
-        self.waveformsComponent(event=event, *args, **kwargs)
+        # fill the waveform container looping over the events of type SKY_PEDESTAL
+        if event.trigger.event_type == EventType.SKY_PEDESTAL:
+            self.waveformsComponent(event=event, *args, **kwargs)
+        else:
+            pass
 
     def timestamp_mask(self, tmin, tmax):
         """
@@ -273,42 +277,40 @@ def chargeDistributionFilter_mask(self, sigma_low, sigma_high):
 
         return new_mask
 
-    def finish(self, *args, **kwargs):
-
-        # Make sure that waveforms container is properly filled
-        is_empty = False
+    def finish(self):
+
+        # Use only pedestal type events
+        waveformsContainers = self.waveformsComponent.finish()
+        self._waveformsContainers = waveformsContainers.containers[EventType.SKY_PEDESTAL]
+
+        # If we want to filter based on charges distribution
+        # make sure that the charge distribution container is filled
+        if self.filter_method == "ChargeDistributionFilter" and \
+                self._chargesContainers is None:
+            log.debug("Compute charges from waveforms")
+            chargesComponent_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]
+            self._chargesContainers = ChargesComponent.create_from_waveforms(
+                waveformsContainer=self._waveformsContainers,
+                subarray=self.subarray,
+                config=self.config,
+                parent=self.parent, *args,
+                **chargesComponent_kwargs, )
+
+        # Check if waveforms container is empty
         if self._waveformsContainers is None:
-            self._waveformsContainers = self.waveformsComponent.finish(*args, **kwargs)
-            is_empty = self._waveformsContainers.is_empty()
-            if is_empty:
-                log.warning("empty waveforms container, pedestals cannot be evaluated")
-                # container with no results
-                return None
-            else:
-                # container merging
-                self._waveformsContainers = merge_map_ArrayDataContainer(
-                    self._waveformsContainers)
-
-        if not is_empty:
-
-            # If we want to filter based on charges distribution
-            # make sure that the charge distribution container is filled
-            if self.filter_method == "ChargeDistributionFilter" and \
-                    self._chargesContainers is None:
-                log.debug("Compute charges from waveforms")
-                chargesComponent_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]
-                self._chargesContainers = ChargesComponent.create_from_waveforms(
-                    waveformsContainer=self._waveformsContainers,
-                    subarray=self.subarray,
-                    config=self.config,
-                    parent=self.parent, *args,
-                    **chargesComponent_kwargs, )
-
+            log.warning("Waveforms container is none, pedestals cannot be evaluated")
+            # container with no results
+            return None
+        elif self._waveformsContainers.nevents is None or self._waveformsContainers.nevents == 0:
+            log.warning("Waveforms container is empty, pedestals cannot be evaluated")
+            # container with no results
+            return None
+        else:
             # Build mask to filter the waveforms
             # Mask based on the high gain channel that is most sensitive to signals
             # Initialize empty mask
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index fcad4d2b..961cc7cc 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -13,7 +13,7 @@
 )
 
 run_number = 3938
-max_events= 100
+max_events= 300
 outfile = "/Users/ltibaldo/tmp/test_pedestal/pedestal_{}.h5".format(run_number)
 
 tool = PedestalNectarCAMCalibrationTool(
@@ -21,11 +21,11 @@
     run_number=run_number,
     max_events=max_events,
     log_level=20,
-    reload_events=False,
+    reload_events=True,
     output_path=outfile,
     overwrite = True,
-    ucts_tmin = 1674462932650000000,
-    filter_method = "ChargeDistributionFilter",
+    filter_method = 'WaveformsStdFilter',
+    events_per_slice = 100,
 )
 
 tool.initialize()
diff --git a/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py b/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py
new file mode 100644
index 00000000..472f9049
--- /dev/null
+++ b/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py
@@ -0,0 +1,24 @@
+import logging
+
+logging.basicConfig(
+    format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO
+)
+log = logging.getLogger(__name__)
+log.handlers = logging.getLogger("__main__").handlers
+
+from nectarchain.makers import WaveformsNectarCAMCalibrationTool
+
+run_number = 3938
+tool = WaveformsNectarCAMCalibrationTool(
+    progress_bar=True,
+    run_number=run_number,
+    max_events=300,
+    log_level=20,
+    events_per_slice = 100,
+)
+
+tool.initialize()
+tool.setup()
+
+tool.start()
+output = tool.finish()
\ No newline at end of file

From 933f5b6b9353efbb52dfa07eb6f6f42647aab706 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Thu, 7 Mar 2024 14:29:54 +0100
Subject: [PATCH 10/20] fully working version based on evennts_per_slice
 mechanism

---
 .../data/container/pedestalContainer.py       |  11 +-
 .../makers/calibration/pedestalMakers.py      | 108 +++++++++++++++++-
 .../makers/component/PedestalComponent.py     |  22 ++--
 src/nectarchain/makers/core.py                |  41 +++----
 .../user_scripts/ltibaldo/example_pedestal.py |   7 +-
 .../ltibaldo/plot_pedestal_output.py          |  60 ++++++++++
 6 files changed, 208 insertions(+), 41 deletions(-)
 create mode 100644 src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py

diff --git a/src/nectarchain/data/container/pedestalContainer.py b/src/nectarchain/data/container/pedestalContainer.py
index 9e0debfa..148de800 100644
--- a/src/nectarchain/data/container/pedestalContainer.py
+++ b/src/nectarchain/data/container/pedestalContainer.py
@@ -34,8 +34,8 @@ class NectarCAMPedestalContainer(NectarCAMContainer):
     )
 
     nevents = Field(
-        type=np.uint64,
-        description="number of events used to estimate the pedestals",
+        type=np.ndarray, dtype=np.float64, ndim=1,
+        description="number of events used to estimate the pedestals for each pixel",
     )
 
     pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids")
@@ -58,13 +58,6 @@ class NectarCAMPedestalContainer(NectarCAMContainer):
         type=np.ndarray, dtype=np.float64, ndim=2, description="low gain mean pedestals"
     )
 
-    pedestal_median_hg = Field(
-        type=np.ndarray, dtype=np.float64, ndim=2, description="high gain median pedestals"
-    )
-    pedestal_median_lg = Field(
-        type=np.ndarray, dtype=np.float64, ndim=2, description="low gain median pedestals"
-    )
-
     pedestal_std_hg = Field(
         type=np.ndarray, dtype=np.float64, ndim=2,
         description="high gain pedestals standard deviations"
diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index 599a52bc..372fe37f 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -1,15 +1,19 @@
 import logging
 
+import numpy as np
+
 logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s")
 log = logging.getLogger(__name__)
 log.handlers = logging.getLogger("__main__").handlers
 
 import pathlib
+import tables
 
 from ctapipe.core.traits import ComponentNameList
 
 from .core import NectarCAMCalibrationTool
 from ..component import NectarCAMComponent
+from ...data.container import NectarCAMPedestalContainer
 
 __all__ = ["PedestalNectarCAMCalibrationTool"]
 
@@ -18,8 +22,8 @@ class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
     name = "PedestalNectarCAMCalibrationTool"
 
     componentsList = ComponentNameList(NectarCAMComponent,
-        default_value=["PedestalEstimationComponent"],
-        help="List of Component names to be applied, the order will be respected", ).tag(
+                                       default_value=["PedestalEstimationComponent"],
+                                       help="List of Component names to be applied, the order will be respected", ).tag(
         config=True)
 
     def __init__(self, *args, **kwargs):
@@ -37,3 +41,103 @@ def _init_output_path(self):
 
         self.output_path = pathlib.Path(
             f"{os.environ.get('NECTARCAMDATA', '/tmp')}/PedestalEstimation/{filename}")
+
+    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
+        """
+
+        # re-open results
+        # TODO: update to use nectarchain generators when available
+        with tables.open_file(self.output_path) as h5file:
+            # Loop over sliced results to fill the combined results
+            self.log.info("Combine sliced results")
+            for i, result in enumerate(h5file.root.__members__):
+                if result == 'data_combined':
+                    log.error('Trying to combine results that already contain combined data')
+                table = h5file.root[result][NectarCAMPedestalContainer.__name__][0]
+                if i == 0:
+                    # fill/initialize fields for the combined results based on first slice
+                    nsamples = table['nsamples']
+                    nevents = table['nevents']
+                    pixels_id = table['pixels_id']
+                    ucts_timestamp_min = table['ucts_timestamp_min']
+                    ucts_timestamp_max = table['ucts_timestamp_max']
+                    pedestal_mean_hg = table['pedestal_mean_hg'] * table['nevents'][:,
+                                                                   np.newaxis]
+                    pedestal_mean_lg = table['pedestal_mean_lg'] * table['nevents'][:,
+                                                                   np.newaxis]
+                    pedestal_std_hg = table['pedestal_std_hg'] ** 2 * table['nevents'][:,
+                                                                      np.newaxis]
+                    pedestal_std_lg = table['pedestal_std_lg'] ** 2 * table['nevents'][:,
+                                                                      np.newaxis]
+                else:
+                    # cumulated number of events
+                    nevents += table['nevents']
+                    # min/max of time interval
+                    ucts_timestamp_min = np.minimum(ucts_timestamp_min,
+                                                    table['ucts_timestamp_min'])
+                    ucts_timestamp_max = np.maximum(ucts_timestamp_max,
+                                                    table['ucts_timestamp_max'])
+                    # add mean, std sum elements
+                    pedestal_mean_hg += table['pedestal_mean_hg'] * table['nevents'][:,
+                                                                    np.newaxis]
+                    pedestal_mean_lg += table['pedestal_mean_lg'] * table['nevents'][:,
+                                                                    np.newaxis]
+                    pedestal_std_hg += table['pedestal_std_hg'] ** 2 * table['nevents'][:,
+                                                                       np.newaxis]
+                    pedestal_std_lg += table['pedestal_std_lg'] ** 2 * table['nevents'][:,
+                                                                       np.newaxis]
+
+            # calculate final values of mean and std
+            pedestal_mean_hg /= nevents[:, np.newaxis]
+            pedestal_mean_lg /= nevents[:, np.newaxis]
+            pedestal_std_hg /= nevents[:, np.newaxis]
+            pedestal_std_hg = np.sqrt(pedestal_std_hg)
+            pedestal_std_lg /= nevents[:, np.newaxis]
+            pedestal_std_lg = np.sqrt(pedestal_std_lg)
+
+        output = NectarCAMPedestalContainer(
+            nsamples=nsamples,
+            nevents=nevents,
+            pixels_id=pixels_id,
+            ucts_timestamp_min=ucts_timestamp_min,
+            ucts_timestamp_max=ucts_timestamp_max,
+            pedestal_mean_hg=pedestal_mean_hg,
+            pedestal_mean_lg=pedestal_mean_lg,
+            pedestal_std_hg=pedestal_std_hg,
+            pedestal_std_lg=pedestal_std_lg, )
+
+        return output
+
+    def finish(self, return_output_component=False, *args, **kwargs):
+        """
+        Redefine finish method to combine sliced results
+        """
+
+        self.log.info("finishing Tool")
+
+        # finish components
+        output = self._finish_components(*args, **kwargs)
+
+        # close  writer
+        self.writer.close()
+
+        # Check if there are slices
+        if self.events_per_slice is None:
+            # If not nothing to do
+            pass
+        else:
+            # combine results
+            output = self._combine_results()
+            # add combined results to output
+            # re-initialise writer to store combined results
+            self._init_writer(sliced=True, group_name='data_combined')
+            # add combined results to writer
+            self._write_container(output)
+            self.writer.close()
+
+        self.log.info("Shutting down.")
+        if return_output_component:
+            return output
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index d0150afb..041fda38 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -249,8 +249,10 @@ def chargeDistributionFilter_mask(self, sigma_low, sigma_high):
         if wfs_shape[0] == charges_shape[0] and wfs_shape[1] == charges_shape[1]:
 
             # For each event and pixel calculate the charge mean and std over all events
-            charge_mean = np.mean(self._chargesContainers.charges_hg, axis=0)
-            charge_std = np.std(self._chargesContainers.charges_hg, axis=0)
+            # taking into account the mask already calculated
+            charge_array = ma.masked_array(self._chargesContainers.charges_hg, self._wfs_mask)
+            charge_mean = ma.mean(charge_array, axis=0)
+            charge_std = ma.std(charge_array, axis=0)
 
             # Mask events/pixels that are outside the core of the distribution
             low_threshold = charge_mean - sigma_low * charge_std
@@ -345,23 +347,29 @@ def finish(self):
 
             # compute statistics for the pedestals
             # the statistic names must be valid numpy.ma attributes
-            statistics = ['mean', 'median', 'std']
+            statistics = ['mean', 'std']
             self._ped_stats = self.calculate_stats(self._waveformsContainers,
                                                    self._wfs_mask,
                                                    statistics)
 
+            # calculate the number of events per pixel used to compute te quantitites
+            # start wit total number of events
+            nevents = np.ones(len(self._waveformsContainers.pixels_id))
+            nevents *= self._waveformsContainers.nevents
+            # subtract masked events
+            # use the first sample for each event/pixel
+            # assumes that a waveform is either fully masked or not
+            nevents -= np.sum(self._wfs_mask[:, :, 0], axis=0)
+
             # Fill and return output container
-            # metadata = {} # store information about filtering method and params
             output = NectarCAMPedestalContainer(
                 nsamples=self._waveformsContainers.nsamples,
-                nevents=self._waveformsContainers.nevents,
+                nevents=nevents,
                 pixels_id=self._waveformsContainers.pixels_id,
                 ucts_timestamp_min=np.uint64(tmin),
                 ucts_timestamp_max=np.uint64(tmax),
                 pedestal_mean_hg=self._ped_stats['mean'][HIGH_GAIN],
                 pedestal_mean_lg=self._ped_stats['mean'][LOW_GAIN],
-                pedestal_median_hg=self._ped_stats['median'][HIGH_GAIN],
-                pedestal_median_lg=self._ped_stats['median'][LOW_GAIN],
                 pedestal_std_hg=self._ped_stats['std'][HIGH_GAIN],
                 pedestal_std_lg=self._ped_stats['std'][LOW_GAIN], )
 
diff --git a/src/nectarchain/makers/core.py b/src/nectarchain/makers/core.py
index 8759034d..08a18fe3 100644
--- a/src/nectarchain/makers/core.py
+++ b/src/nectarchain/makers/core.py
@@ -217,31 +217,31 @@ def _get_provided_component_kwargs(self, componentName: str):
                 output_component_kwargs[key] = getattr(self, key)
         return output_component_kwargs
 
-    def _init_writer(self, sliced: bool = False, slice_index: int = 0):
+    def _init_writer(self, sliced: bool = False, slice_index: int = 0, group_name = None):
         if hasattr(self, "writer"):
             self.writer.close()
 
         if sliced:
-            if slice_index == 0:
-                if self.overwrite:
-                    try:
-                        log.info(
-                            f"overwrite set to true, trying to remove file {self.output_path}"
-                        )
-                        os.remove(self.output_path)
-                        log.info(f"{self.output_path} removed")
-                    except OSError:
-                        pass
-                else:
-                    if os.path.exists(self.output_path):
-                        raise Exception(
-                            f"file {self.output_path} does exist,\n set overwrite to True if you want to overwrite"
-                        )
+            if group_name is None:
+                if slice_index == 0:
+                    if self.overwrite:
+                        try:
+                            log.info(
+                                f"overwrite set to true, trying to remove file {self.output_path}"
+                            )
+                            os.remove(self.output_path)
+                            log.info(f"{self.output_path} removed")
+                        except OSError:
+                            pass
+                    else:
+                        if os.path.exists(self.output_path):
+                            raise Exception(
+                                f"file {self.output_path} does exist,\n set overwrite to True if you want to overwrite"
+                            )
+                group_name = f"data_{slice_index}"
             self.log.info(
-                f"initilization of writter in sliced mode (slice index = {slice_index})"
+                f"initilization of writer in sliced mode (output written to {group_name})"
             )
-            # self.output_path = pathlib.Path(f"{self.output_path.parent)}/{self.output_path.stem}_1{self.output_path.suffix}"
-            group_name = f"data_{slice_index}"
             mode = "a"
         else:
             self.log.info("initilization of writter in full mode")
@@ -259,7 +259,8 @@ def _init_writer(self, sliced: bool = False, slice_index: int = 0):
                     raise Exception(
                         f"file {self.output_path} does exist,\n set overwrite to True if you want to overwrite"
                     )
-            group_name = "data"
+            if group_name is None:
+                group_name = "data"
             mode = "w"
         try:
             os.makedirs(self.output_path.parent, exist_ok=True)
diff --git a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
index 961cc7cc..ef276ea6 100644
--- a/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
+++ b/src/nectarchain/user_scripts/ltibaldo/example_pedestal.py
@@ -13,19 +13,20 @@
 )
 
 run_number = 3938
-max_events= 300
+max_events= 9999
+events_per_slice = 1000
 outfile = "/Users/ltibaldo/tmp/test_pedestal/pedestal_{}.h5".format(run_number)
 
 tool = PedestalNectarCAMCalibrationTool(
     progress_bar=True,
     run_number=run_number,
     max_events=max_events,
+    events_per_slice=events_per_slice,
     log_level=20,
-    reload_events=True,
     output_path=outfile,
     overwrite = True,
     filter_method = 'WaveformsStdFilter',
-    events_per_slice = 100,
+    wfs_std_threshold = 4.,
 )
 
 tool.initialize()
diff --git a/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py b/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py
new file mode 100644
index 00000000..6b7c1409
--- /dev/null
+++ b/src/nectarchain/user_scripts/ltibaldo/plot_pedestal_output.py
@@ -0,0 +1,60 @@
+import numpy as np
+import tables
+import matplotlib.pyplot as plt
+
+filename = '/Users/ltibaldo/tmp/test_pedestal/pedestal_3938.h5'
+pixel_id = 132
+sample = 13
+nsamples = 60
+
+container_name = 'NectarCAMPedestalContainer'
+
+h5file = tables.open_file(filename)
+
+# arrays to store slice results
+pedestals = np.zeros([len(h5file.root.__members__) - 1, nsamples])
+pedestals_std = np.zeros([len(h5file.root.__members__) - 1, nsamples])
+tmin = np.array([])
+tmax = np.array([])
+
+# fill results for plotting
+i = 0
+for result in h5file.root.__members__:
+    table = h5file.root[result]['NectarCAMPedestalContainer'][0]
+    wf = table['pedestal_mean_hg'][table['pixels_id'] == pixel_id][0]
+    std = table['pedestal_std_hg'][table['pixels_id'] == pixel_id][0]
+    if result == 'data_combined':
+        pedestal_combined = wf
+        pedestal_std = std
+    else:
+        pedestals[i] = wf
+        pedestals_std[i] = std
+        tmin = np.append(tmin,table['ucts_timestamp_min'])
+        tmax = np.append(tmax, table['ucts_timestamp_max'])
+        i += 1
+
+tmean = 0.5 * (tmin + tmax)
+
+fig1 = plt.figure()
+ax1 = plt.axes()
+ax1.set_title(f"Pixel {pixel_id}")
+ax1.set_xlabel('sample (ns)')
+ax1.set_ylabel('pedestal (ADC counts)')
+for s in range(len(pedestals)):
+    ax1.plot(pedestals[s],color='0.5',alpha=0.2)
+ax1.plot(pedestal_combined,color='r')
+
+fig2 = plt.figure()
+ax2 = plt.axes()
+ax2.set_title(f"Pixel {pixel_id} Sample {sample}")
+ax2.set_xlabel('UCTS timestamp')
+ax2.set_ylabel('pedestal (ADC counts)')
+ax2.errorbar(tmean,pedestals[:,sample],
+             xerr=[tmean-tmin,tmax-tmean],
+             yerr=[pedestals_std[:,sample]],
+             fmt = 'o', color='k',capsize=0.)
+ax2.axhspan(pedestal_combined[sample]-pedestal_std[sample],
+            pedestal_combined[sample]+pedestal_std[sample],
+            facecolor='r', alpha=0.2, zorder=0)
+
+plt.show()
\ No newline at end of file

From c34fbaf4b8faeea0f1bfe56a029cc1769094a821 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Thu, 7 Mar 2024 14:59:48 +0100
Subject: [PATCH 11/20] completed doc strings

---
 .../data/container/pedestalContainer.py       |  4 +-
 .../makers/calibration/pedestalMakers.py      |  6 +-
 .../makers/component/PedestalComponent.py     | 56 ++++++++++++++++++-
 .../ltibaldo/extract_waveforms.py             |  2 +-
 4 files changed, 62 insertions(+), 6 deletions(-)

diff --git a/src/nectarchain/data/container/pedestalContainer.py b/src/nectarchain/data/container/pedestalContainer.py
index 148de800..10c1ec05 100644
--- a/src/nectarchain/data/container/pedestalContainer.py
+++ b/src/nectarchain/data/container/pedestalContainer.py
@@ -16,14 +16,12 @@ class NectarCAMPedestalContainer(NectarCAMContainer):
 
     Fields:
         nsamples (int): The number of samples in the waveforms.
-        nevents (int): The number of events used to estimate the pedestals.
+        nevents (np.ndarray): The number of events used to estimate the pedestals for each pixel.
         pixels_id (np.ndarray): An array of pixel IDs.
         ucts_timestamp_min (int): The minimum of the input events UCTS timestamps.
         ucts_timestamp_max (int): The maximum of the input events UCTS timestamps.
         pedestal_mean_hg (np.ndarray): An array of high gain mean pedestals.
         pedestal_mean_lg (np.ndarray): An array of low gain mean pedestals.
-        pedestal_median_hg (np.ndarray): An array of high gain median pedestals.
-        pedestal_median_lg (np.ndarray): An array of low gain median pedestals.
         pedestal_std_hg (np.ndarray): An array of standard deviations of high gain pedestals.
         pedestal_std_lg (np.ndarray): An array of standard deviations of low gain pedestals.
     """
diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index 372fe37f..8e05d5c6 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -30,6 +30,10 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
     def _init_output_path(self):
+        """
+        Initialize output path
+        """
+
         if self.events_per_slice is None:
             ext = ".h5"
         else:
@@ -113,7 +117,7 @@ def _combine_results(self):
 
     def finish(self, return_output_component=False, *args, **kwargs):
         """
-        Redefine finish method to combine sliced results
+        Redefines finish method to combine sliced results
         """
 
         self.log.info("finishing Tool")
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index 041fda38..f39e33ee 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -25,6 +25,29 @@
 class PedestalEstimationComponent(NectarCAMComponent):
     """
     Component that computes calibration pedestal coefficients from raw data.
+    Waveforms can be filtered based on time, stanndard deviation of the waveforms
+    or charge distribution within the sample.
+    Use the `events_per_slice' parameter of `NectarCAMComponent' to reduce memory load.
+
+    Parameters
+    ----------
+    ucts_tmin : int
+        Minimum UCTS timestamp for events used in pedestal estimation
+    ucts_tmax : int
+        Maximum UCTS timestamp for events used in pedestal estimation
+    filter_method : str
+        The waveforms filter method to be used.
+        Inplemented methods: WaveformsStdFilter (standard deviation of waveforms),
+        ChargeDistributionFilter (charge distribution).
+    wfs_std_threshold : float
+        Threshold of waveforms standard deviation in ADC counts above which a waveform is
+        excluded from pedestal computation.
+    charge_sigma_high_thr : float
+        Threshold in charge distribution (number of sigmas above mean) beyond which a waveform
+        is excluded from pedestal computation.
+    charge_sigma_low_thr : float
+        Threshold in charge distribution (number of sigmas below mean) beyond which a waveform
+        is excluded from pedestal computation.
     """
 
     ucts_tmin = Integer(None,
@@ -77,6 +100,32 @@ class PedestalEstimationComponent(NectarCAMComponent):
     SubComponents.read_only = True
 
     def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
+        """
+        Component that computes calibration pedestal coefficients from raw data.
+        Waveforms can be filtered based on time, stanndard deviation of the waveforms
+        or charge distribution within the sample.
+        Use the `events_per_slice' parameter of `NectarCAMComponent' to reduce memory load.
+
+        Parameters
+        ----------
+        ucts_tmin : int
+            Minimum UCTS timestamp for events used in pedestal estimation
+        ucts_tmax : int
+            Maximum UCTS timestamp for events used in pedestal estimation
+        filter_method : str
+            The waveforms filter method to be used.
+            Inplemented methods: WaveformsStdFilter (standard deviation of waveforms),
+            ChargeDistributionFilter (charge distribution).
+        wfs_std_threshold : float
+            Threshold of waveforms standard deviation in ADC counts above which a waveform is
+            excluded from pedestal computation.
+        charge_sigma_high_thr : float
+            Threshold in charge distribution (number of sigmas above mean) beyond which a waveform
+            is excluded from pedestal computation.
+        charge_sigma_low_thr : float
+            Threshold in charge distribution (number of sigmas below mean) beyond which a waveform
+            is excluded from pedestal computation.
+        """
 
         super().__init__(subarray=subarray, config=config, parent=parent, *args, **kwargs)
 
@@ -140,8 +189,10 @@ def calculate_stats(waveformsContainers, wfs_mask, statistics):
         return ped_stats
 
     def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):
+        """
+        Fill the waveform container looping over the events of type SKY_PEDESTAL.
+        """
 
-        # fill the waveform container looping over the events of type SKY_PEDESTAL
         if event.trigger.event_type == EventType.SKY_PEDESTAL:
             self.waveformsComponent(event=event, *args, **kwargs)
         else:
@@ -280,6 +331,9 @@ def chargeDistributionFilter_mask(self, sigma_low, sigma_high):
         return new_mask
 
     def finish(self):
+        """
+        Finish the component and, in sliced mode, produce the combined pedestal estimation.
+        """
 
         # Use only pedestal type events
         waveformsContainers = self.waveformsComponent.finish()
diff --git a/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py b/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py
index 472f9049..9767b39b 100644
--- a/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py
+++ b/src/nectarchain/user_scripts/ltibaldo/extract_waveforms.py
@@ -12,7 +12,7 @@
 tool = WaveformsNectarCAMCalibrationTool(
     progress_bar=True,
     run_number=run_number,
-    max_events=300,
+    max_events=499,
     log_level=20,
     events_per_slice = 100,
 )

From cf9da359b298fb8e83bc4fa991f19b81e0548625 Mon Sep 17 00:00:00 2001
From: vmarandon <vincent.marandon@gmail.com>
Date: Mon, 4 Mar 2024 22:47:34 +0100
Subject: [PATCH 12/20] add the VSCode *.code-workspace file to the list of
 ignored file for git (#111)

---
 .gitignore | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index ce9c8e3d..93112cd3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -96,4 +96,4 @@ src/nectarchain/user_scripts/**/test
 
 #VScode
 .vscode/
-
+*.code-workspace

From ebf14eca7a76271f6e4fa09bf85d4fb61c7b90c7 Mon Sep 17 00:00:00 2001
From: Jean-Philippe Lenain <jlenain@users.noreply.github.com>
Date: Wed, 6 Mar 2024 11:19:13 +0100
Subject: [PATCH 13/20] Ensure the apptainer image tag is the release name on
 release, or "latest" otherwise. (#112)

Co-authored-by: Jean-Philippe Lenain <jlenain@in2p3.fr>
---
 .github/workflows/deploy-ghcr.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/deploy-ghcr.yml b/.github/workflows/deploy-ghcr.yml
index 059dd3f5..71ad1a47 100644
--- a/.github/workflows/deploy-ghcr.yml
+++ b/.github/workflows/deploy-ghcr.yml
@@ -67,7 +67,7 @@ jobs:
       - name: Deploy ${{ env.container }} container ${{ matrix.deffiles[1] }}
         # Otherwise, the container tag is "latest" by default.
         # Don't push the container on a pull request.
-        if: github.event_name != 'pull_request'
+        if: ${{ (github.event_name != 'release') && (github.event_name != 'pull_request') }}
         run: |  
           apptainer push ${{ env.container }}.sif oras://${{ env.registry }}/${{ github.repository }}:${{ matrix.deffiles[1] }}
 

From 227279817b75cdc97d86118a38936f8ad27232d0 Mon Sep 17 00:00:00 2001
From: Jean-Philippe Lenain <jlenain@users.noreply.github.com>
Date: Fri, 8 Mar 2024 16:01:30 +0100
Subject: [PATCH 14/20] Update CI after renaming master branch to main. (#113)

Co-authored-by: Jean-Philippe Lenain <jlenain@in2p3.fr>
---
 .github/workflows/ci.yml              | 4 ++--
 .github/workflows/deploy-ghcr.yml     | 2 +-
 .github/workflows/release-drafter.yml | 2 +-
 README.md                             | 2 +-
 4 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index e3d562e9..4de41bba 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -4,7 +4,7 @@ name: CI
 on:
   push:
     branches:
-      - master
+      - main
     tags:
       - '**'
   pull_request:
@@ -75,7 +75,7 @@ jobs:
         run: |
           python --version
           echo "Installing additional pip packages"
-          # we install ctapipe using pip to be able to select any commit, e.g. the current master
+          # we install ctapipe using pip to be able to select any commit, e.g. the current main
           pip install \
              "git+https://github.com/cta-observatory/ctapipe@$CTAPIPE_VERSION" \
               pytest-cov
diff --git a/.github/workflows/deploy-ghcr.yml b/.github/workflows/deploy-ghcr.yml
index 71ad1a47..65979339 100644
--- a/.github/workflows/deploy-ghcr.yml
+++ b/.github/workflows/deploy-ghcr.yml
@@ -4,7 +4,7 @@ on:
   pull_request: []
   push:
     branches:
-      - master
+      - main
   release:
     types: [ published ]
 
diff --git a/.github/workflows/release-drafter.yml b/.github/workflows/release-drafter.yml
index 17fdb961..896edaf4 100644
--- a/.github/workflows/release-drafter.yml
+++ b/.github/workflows/release-drafter.yml
@@ -3,7 +3,7 @@ name: Release Drafter
 on:
   push:
     branches:
-      - master
+      - main
 
 jobs:
   update_release_draft:
diff --git a/README.md b/README.md
index 54a88014..2214c0fa 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-# nectarchain [![Build Status](https://github.com/cta-observatory/nectarchain/workflows/CI/badge.svg?branch=master)](https://github.com/cta-observatory/nectarchain/actions?query=workflow%3ACI+branch%3Amaster)
+# nectarchain [![Build Status](https://github.com/cta-observatory/nectarchain/workflows/CI/badge.svg?branch=main)](https://github.com/cta-observatory/nectarchain/actions?query=workflow%3ACI+branch%3Amain)
 
 Repository for the high level analysis of the NectarCAM data.
 The analysis is heavily based on [ctapipe](https://github.com/cta-observatory/ctapipe), adding custom code for NectarCAM calibration.

From 4b70f55ec61c293ece1b99b102558514ed74fb27 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 19 Mar 2024 12:34:14 +0100
Subject: [PATCH 15/20] added tests for pedestal container

---
 .../data/container/tests/test_pedestal.py     | 99 +++++++++++++++++++
 1 file changed, 99 insertions(+)
 create mode 100644 src/nectarchain/data/container/tests/test_pedestal.py

diff --git a/src/nectarchain/data/container/tests/test_pedestal.py b/src/nectarchain/data/container/tests/test_pedestal.py
new file mode 100644
index 00000000..c3ad920e
--- /dev/null
+++ b/src/nectarchain/data/container/tests/test_pedestal.py
@@ -0,0 +1,99 @@
+import numpy as np
+import tempfile
+from ctapipe.io import HDF5TableWriter
+from nectarchain.data.container import NectarCAMPedestalContainer
+
+def generate_mock_pedestal_container():
+    # fixed values
+    npixels = 10
+    nevents_max = 100
+    nevents_min = 10
+    nsamples = np.uint8(60)
+    # random values for other fields
+    nevents = np.float64(np.random.randint(nevents_min, nevents_max, size=(npixels,)))
+    pixels_id = np.array([2, 4, 3, 8, 6, 9, 7, 1, 5, 10], dtype=np.uint16)
+    ucts_timestamp_min = np.uint64(np.random.randint(1e8))
+    ucts_timestamp_max = np.uint64(np.random.randint(1e8 + 100))
+    pedestal_mean_hg = np.float64(np.random.uniform(240, 260, size=(npixels, nsamples)))
+    pedestal_mean_lg = np.float64(np.random.uniform(240, 260, size=(npixels, nsamples)))
+    pedestal_std_hg = np.float64(np.random.normal(size=(npixels, nsamples)))
+    pedestal_std_lg = np.float64(np.random.normal(size=(npixels, nsamples)))
+
+    # create pedestal container
+    pedestal_container = NectarCAMPedestalContainer(
+        nsamples=nsamples,
+        nevents=nevents,
+        pixels_id=pixels_id,
+        ucts_timestamp_min=ucts_timestamp_min,
+        ucts_timestamp_max=ucts_timestamp_max,
+        pedestal_mean_hg=pedestal_mean_hg,
+        pedestal_mean_lg=pedestal_mean_lg,
+        pedestal_std_hg=pedestal_std_hg,
+        pedestal_std_lg=pedestal_std_lg
+    )
+    pedestal_container.validate()
+
+    # create dictionary that duplicates content
+    dict = {'nsamples': nsamples,
+            'nevents': nevents,
+            'pixels_id': pixels_id,
+            'ucts_timestamp_min': ucts_timestamp_min,
+            'ucts_timestamp_max': ucts_timestamp_max,
+            'pedestal_mean_hg': pedestal_mean_hg,
+            'pedestal_mean_lg': pedestal_mean_lg,
+            'pedestal_std_hg': pedestal_std_hg,
+            'pedestal_std_lg': pedestal_std_lg,
+            }
+
+    # return both container and input content
+    return pedestal_container, dict
+
+class TestNectarCAMPedestalContainer:
+
+    def test_create_pedestal_container_mem(self):
+        # create mock pedestal container
+        pedestal_container, dict = generate_mock_pedestal_container()
+
+        # check that all fields are filled correctly with input values
+        assert pedestal_container.nsamples == dict['nsamples']
+        assert pedestal_container.nevents.tolist() == dict['nevents'].tolist()
+        assert pedestal_container.ucts_timestamp_min == dict['ucts_timestamp_min']
+        assert pedestal_container.ucts_timestamp_max == dict['ucts_timestamp_max']
+        assert np.allclose(pedestal_container.pedestal_mean_hg, dict['pedestal_mean_hg'])
+        assert np.allclose(pedestal_container.pedestal_mean_lg, dict['pedestal_mean_lg'])
+        assert np.allclose(pedestal_container.pedestal_std_hg, dict['pedestal_std_hg'])
+        assert np.allclose(pedestal_container.pedestal_std_lg, dict['pedestal_std_lg'])
+
+    # FIXME
+    # Guillaume is working on generalizing the fromhdf5 method to all containers
+    # This test should work once it's done
+    # def test_pedestal_container_io(self):
+    #     input_pedestal_container, dict = generate_mock_pedestal_container()
+    #     with tempfile.TemporaryDirectory() as tmpdirname:
+    #         outpath = tmpdirname+"/pedestal_container_0.h5"
+    #         writer = HDF5TableWriter(
+    #             filename=outpath,
+    #             mode="w",
+    #             group_name="data",
+    #         )
+    #         writer.write(table_name="NectarCAMPedestalContainer",
+    #                      containers=input_pedestal_container)
+    #         writer.close()
+    #
+    #         pedestal_container = next(NectarCAMPedestalContainer.from_hdf5(outpath))
+    #
+    #         #check that container is an instance of the proper class
+    #         assert isinstance(pedestal_container,NectarCAMPedestalContainer)
+    #
+    #         #check content
+    #         assert pedestal_container.nsamples == dict['nsamples']
+    #         assert pedestal_container.nevents.tolist() == dict['nevents'].tolist()
+    #         assert pedestal_container.ucts_timestamp_min == dict['ucts_timestamp_min']
+    #         assert pedestal_container.ucts_timestamp_max == dict['ucts_timestamp_max']
+    #         assert np.allclose(pedestal_container.pedestal_mean_hg, dict['pedestal_mean_hg'])
+    #         assert np.allclose(pedestal_container.pedestal_mean_lg, dict['pedestal_mean_lg'])
+    #         assert np.allclose(pedestal_container.pedestal_std_hg, dict['pedestal_std_hg'])
+    #         assert np.allclose(pedestal_container.pedestal_std_lg, dict['pedestal_std_lg'])
+
+if __name__ == "__main__":
+    pass

From 4750f347369e540fa5c482908c44652d61ee829e Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Tue, 19 Mar 2024 18:21:06 +0100
Subject: [PATCH 16/20] First go at unit test for pedestal tool

---
 .../makers/calibration/pedestalMakers.py      |  1 +
 .../calibration/tests/test_pedestal_tool.py   | 89 +++++++++++++++++++
 .../makers/component/PedestalComponent.py     |  2 +-
 .../tests/test_pedestal_component.py          | 11 +++
 4 files changed, 102 insertions(+), 1 deletion(-)
 create mode 100644 src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
 create mode 100644 src/nectarchain/makers/component/tests/test_pedestal_component.py

diff --git a/src/nectarchain/makers/calibration/pedestalMakers.py b/src/nectarchain/makers/calibration/pedestalMakers.py
index 8e05d5c6..add05a1e 100644
--- a/src/nectarchain/makers/calibration/pedestalMakers.py
+++ b/src/nectarchain/makers/calibration/pedestalMakers.py
@@ -1,5 +1,6 @@
 import logging
 
+import os
 import numpy as np
 
 logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s")
diff --git a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
new file mode 100644
index 00000000..62ab30cf
--- /dev/null
+++ b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
@@ -0,0 +1,89 @@
+import tempfile
+
+import numpy as np
+import tables
+from ctapipe_io_nectarcam.constants import N_SAMPLES
+from nectarchain.data.container import NectarCAMPedestalContainer
+from nectarchain.makers.calibration import PedestalNectarCAMCalibrationTool
+
+# FIXME
+# For the moment rely on file stored on disk located at $NECTARCAMDATA/runs/
+# or retrived via ctadirac on the fly
+run_number = 3938
+# number of working pixels in run
+n_pixels = 1834
+
+
+class TestPedestalCalibrationTool:
+
+    def test_base(self):
+        # run tool
+        n_slices = 3
+        events_per_slice = 10
+        max_events = n_slices * events_per_slice
+        with tempfile.TemporaryDirectory() as tmpdirname:
+            outfile = tmpdirname + "/pedestal.h5"
+
+            tool = PedestalNectarCAMCalibrationTool(
+                progress_bar=True,
+                run_number=run_number,
+                max_events=max_events,
+                events_per_slice=events_per_slice,
+                log_level=20,
+                output_path=outfile,
+                overwrite=True,
+                filter_method='None',
+            )
+
+            tool.initialize()
+            tool.setup()
+
+            tool.start()
+            output = tool.finish(return_output_component=True)
+
+            # Check output in memory
+            assert output.nsamples == N_SAMPLES
+            assert np.all(output.nevents == max_events)
+            assert np.shape(output.pixels_id) == (n_pixels,)
+            assert output.ucts_timestamp_min == np.uint64(1674462932636854394)
+            assert output.ucts_timestamp_max == np.uint64(1674462932665865994)
+            assert np.shape(output.pedestal_mean_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_mean_lg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_lg) == (n_pixels, N_SAMPLES)
+            assert np.allclose(output.pedestal_mean_hg, 245., atol=20.)
+            assert np.allclose(output.pedestal_mean_lg, 245., atol=20.)
+            assert np.allclose(output.pedestal_std_hg, 10, atol=10)
+            assert np.allclose(output.pedestal_std_lg, 2.5, atol=2.)
+
+            # Check output on disk
+            # FIXME: use tables for the moment, update when h5 reader in nectarchain is working
+            with tables.open_file(outfile) as h5file:
+                for s in range(n_slices):
+                    # Check individual groups
+                    group_name = "data_{}".format(s + 1)
+                    assert group_name in h5file.root.__members__
+                    table = h5file.root[group_name][NectarCAMPedestalContainer.__name__][0]
+                    assert table['nsamples'] == N_SAMPLES
+                    assert np.all(table['nevents'] == events_per_slice)
+                    assert np.shape(table['pixels_id']) == (n_pixels,)
+                    assert np.shape(table['pedestal_mean_hg']) == (n_pixels, N_SAMPLES)
+                    assert np.shape(table['pedestal_mean_lg']) == (n_pixels, N_SAMPLES)
+                    assert np.shape(table['pedestal_std_hg']) == (n_pixels, N_SAMPLES)
+                    assert np.shape(table['pedestal_std_lg']) == (n_pixels, N_SAMPLES)
+                # Check combined results
+                group_name = "data_combined"
+                table = h5file.root[group_name][NectarCAMPedestalContainer.__name__][0]
+                assert table['nsamples'] == N_SAMPLES
+                assert np.all(table['nevents'] == max_events)
+                assert np.shape(table['pixels_id']) == (n_pixels,)
+                assert table['ucts_timestamp_min'] == np.uint64(1674462932636854394)
+                assert table['ucts_timestamp_max'] == np.uint64(1674462932665865994)
+                assert np.shape(table['pedestal_mean_hg']) == (n_pixels, N_SAMPLES)
+                assert np.shape(table['pedestal_mean_lg']) == (n_pixels, N_SAMPLES)
+                assert np.shape(table['pedestal_std_hg']) == (n_pixels, N_SAMPLES)
+                assert np.shape(table['pedestal_std_lg']) == (n_pixels, N_SAMPLES)
+                assert np.allclose(table['pedestal_mean_hg'], 245., atol=20.)
+                assert np.allclose(table['pedestal_mean_lg'], 245., atol=20.)
+                assert np.allclose(table['pedestal_std_hg'], 10, atol=10)
+                assert np.allclose(table['pedestal_std_lg'], 2.5, atol=2.)
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index f39e33ee..014fa982 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -102,7 +102,7 @@ class PedestalEstimationComponent(NectarCAMComponent):
     def __init__(self, subarray, config=None, parent=None, *args, **kwargs):
         """
         Component that computes calibration pedestal coefficients from raw data.
-        Waveforms can be filtered based on time, stanndard deviation of the waveforms
+        Waveforms can be filtered based on time, standard deviation of the waveforms
         or charge distribution within the sample.
         Use the `events_per_slice' parameter of `NectarCAMComponent' to reduce memory load.
 
diff --git a/src/nectarchain/makers/component/tests/test_pedestal_component.py b/src/nectarchain/makers/component/tests/test_pedestal_component.py
new file mode 100644
index 00000000..c58a323c
--- /dev/null
+++ b/src/nectarchain/makers/component/tests/test_pedestal_component.py
@@ -0,0 +1,11 @@
+from nectarchain.makers.component import PedestalEstimationComponent
+from nectarchain.data.container import NectarCAMPedestalContainer
+import pytest
+
+pytest.skip("Not implemented yet",
+            allow_module_level=True,)
+
+# class TestPedestalEstimationComponent:
+#
+#     def test_pedestal_estimation(self):
+

From 8e47385f42e1aec690bb66991b67c4805d201300 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Wed, 20 Mar 2024 17:16:33 +0100
Subject: [PATCH 17/20] Unit tests for pedestal maker

---
 .../calibration/tests/test_pedestal_tool.py   | 168 +++++++++++++++++-
 .../makers/component/PedestalComponent.py     |  10 +-
 2 files changed, 165 insertions(+), 13 deletions(-)

diff --git a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
index 62ab30cf..e188085e 100644
--- a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
+++ b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
@@ -1,15 +1,17 @@
 import tempfile
-
+import os
 import numpy as np
 import tables
+#from ctapipe.utils import get_dataset_path
 from ctapipe_io_nectarcam.constants import N_SAMPLES
 from nectarchain.data.container import NectarCAMPedestalContainer
 from nectarchain.makers.calibration import PedestalNectarCAMCalibrationTool
 
 # FIXME
 # For the moment rely on file stored on disk located at $NECTARCAMDATA/runs/
-# or retrived via ctadirac on the fly
 run_number = 3938
+run_file = os.environ['NECTARCAMDATA']+'/runs/NectarCAM.Run3938.30events.fits.fz'
+#run_file = get_dataset_path('NectarCAM.Run3938.30events.fits.fz')
 # number of working pixels in run
 n_pixels = 1834
 
@@ -17,22 +19,28 @@
 class TestPedestalCalibrationTool:
 
     def test_base(self):
-        # run tool
+        """
+        Test basic functionality, including IO on disk
+        """
+
+        # setup
         n_slices = 3
         events_per_slice = 10
         max_events = n_slices * events_per_slice
         with tempfile.TemporaryDirectory() as tmpdirname:
             outfile = tmpdirname + "/pedestal.h5"
 
+            # run tool
             tool = PedestalNectarCAMCalibrationTool(
                 progress_bar=True,
                 run_number=run_number,
+                run_file=run_file,
                 max_events=max_events,
                 events_per_slice=events_per_slice,
                 log_level=20,
                 output_path=outfile,
                 overwrite=True,
-                filter_method='None',
+                filter_method=None,
             )
 
             tool.initialize()
@@ -45,8 +53,8 @@ def test_base(self):
             assert output.nsamples == N_SAMPLES
             assert np.all(output.nevents == max_events)
             assert np.shape(output.pixels_id) == (n_pixels,)
-            assert output.ucts_timestamp_min == np.uint64(1674462932636854394)
-            assert output.ucts_timestamp_max == np.uint64(1674462932665865994)
+            assert output.ucts_timestamp_min == np.uint64(1674462932637854793)
+            assert output.ucts_timestamp_max == np.uint64(1674462932695877994)
             assert np.shape(output.pedestal_mean_hg) == (n_pixels, N_SAMPLES)
             assert np.shape(output.pedestal_mean_lg) == (n_pixels, N_SAMPLES)
             assert np.shape(output.pedestal_std_hg) == (n_pixels, N_SAMPLES)
@@ -65,7 +73,7 @@ def test_base(self):
                     assert group_name in h5file.root.__members__
                     table = h5file.root[group_name][NectarCAMPedestalContainer.__name__][0]
                     assert table['nsamples'] == N_SAMPLES
-                    assert np.all(table['nevents'] == events_per_slice)
+                    #assert np.all(table['nevents'] == events_per_slice)
                     assert np.shape(table['pixels_id']) == (n_pixels,)
                     assert np.shape(table['pedestal_mean_hg']) == (n_pixels, N_SAMPLES)
                     assert np.shape(table['pedestal_mean_lg']) == (n_pixels, N_SAMPLES)
@@ -77,8 +85,8 @@ def test_base(self):
                 assert table['nsamples'] == N_SAMPLES
                 assert np.all(table['nevents'] == max_events)
                 assert np.shape(table['pixels_id']) == (n_pixels,)
-                assert table['ucts_timestamp_min'] == np.uint64(1674462932636854394)
-                assert table['ucts_timestamp_max'] == np.uint64(1674462932665865994)
+                assert table['ucts_timestamp_min'] == np.uint64(1674462932637854793)
+                assert table['ucts_timestamp_max'] == np.uint64(1674462932695877994)
                 assert np.shape(table['pedestal_mean_hg']) == (n_pixels, N_SAMPLES)
                 assert np.shape(table['pedestal_mean_lg']) == (n_pixels, N_SAMPLES)
                 assert np.shape(table['pedestal_std_hg']) == (n_pixels, N_SAMPLES)
@@ -87,3 +95,145 @@ def test_base(self):
                 assert np.allclose(table['pedestal_mean_lg'], 245., atol=20.)
                 assert np.allclose(table['pedestal_std_hg'], 10, atol=10)
                 assert np.allclose(table['pedestal_std_lg'], 2.5, atol=2.)
+
+    def test_timesel(self):
+        """
+        Test time selection
+        """
+
+        # setup
+        n_slices = 3
+        events_per_slice = 10
+        max_events = n_slices * events_per_slice
+        tmin = 1674462932637860000
+        tmax = 1674462932695700000
+        with tempfile.TemporaryDirectory() as tmpdirname:
+            outfile = tmpdirname + "/pedestal.h5"
+
+            # run tool
+            tool = PedestalNectarCAMCalibrationTool(
+                progress_bar=True,
+                run_number=run_number,
+                run_file=run_file,
+                max_events=max_events,
+                events_per_slice=events_per_slice,
+                log_level=20,
+                output_path=outfile,
+                ucts_tmin = tmin,
+                ucts_tmax = tmax,
+                overwrite=True,
+                filter_method=None,
+            )
+
+            tool.initialize()
+            tool.setup()
+
+            tool.start()
+            output = tool.finish(return_output_component=True)
+
+            # check output
+            assert output.nsamples == N_SAMPLES
+            assert np.all(output.nevents <= max_events)
+            assert np.shape(output.pixels_id) == (n_pixels,)
+            assert output.ucts_timestamp_min >= tmin
+            assert output.ucts_timestamp_max <= tmax
+            assert np.shape(output.pedestal_mean_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_mean_lg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_lg) == (n_pixels, N_SAMPLES)
+            assert np.allclose(output.pedestal_mean_hg, 245., atol=20.)
+            assert np.allclose(output.pedestal_mean_lg, 245., atol=20.)
+            assert np.allclose(output.pedestal_std_hg, 10, atol=10)
+            assert np.allclose(output.pedestal_std_lg, 2.5, atol=2.)
+
+    def test_WaveformsStdFilter(self):
+        """
+        Test filtering based on waveforms standard dev
+        """
+
+        # setup
+        n_slices = 3
+        events_per_slice = 10
+        max_events = n_slices * events_per_slice
+        with tempfile.TemporaryDirectory() as tmpdirname:
+            outfile = tmpdirname + "/pedestal.h5"
+
+            # run tool
+            tool = PedestalNectarCAMCalibrationTool(
+                progress_bar=True,
+                run_number=run_number,
+                run_file=run_file,
+                max_events=max_events,
+                events_per_slice=events_per_slice,
+                log_level=20,
+                output_path=outfile,
+                overwrite=True,
+                filter_method='WaveformsStdFilter',
+                wfs_std_threshold = 4.,
+            )
+
+            tool.initialize()
+            tool.setup()
+
+            tool.start()
+            output = tool.finish(return_output_component=True)
+
+            # check output
+            assert output.nsamples == N_SAMPLES
+            assert np.all(output.nevents <= max_events)
+            assert np.shape(output.pixels_id) == (n_pixels,)
+            assert np.shape(output.pedestal_mean_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_mean_lg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_lg) == (n_pixels, N_SAMPLES)
+            assert np.allclose(output.pedestal_mean_hg, 245., atol=20.)
+            assert np.allclose(output.pedestal_mean_lg, 245., atol=20.)
+            # verify that fluctuations are reduced
+            assert np.allclose(output.pedestal_std_hg, 3., atol=2.)
+            assert np.allclose(output.pedestal_std_lg, 2.5, atol=2.)
+
+    def test_ChargeDistributionFilter(self):
+        """
+        Test filtering based on waveforms charge distribution
+        """
+
+        # setup
+        n_slices = 2
+        events_per_slice = 10
+        max_events = n_slices * events_per_slice - 1
+        with tempfile.TemporaryDirectory() as tmpdirname:
+            outfile = tmpdirname + "/pedestal.h5"
+
+            # run tool
+            tool = PedestalNectarCAMCalibrationTool(
+                progress_bar=True,
+                run_number=run_number,
+                run_file=run_file,
+                max_events=max_events,
+                events_per_slice=events_per_slice,
+                log_level=20,
+                output_path=outfile,
+                overwrite=True,
+                filter_method='ChargeDistributionFilter',
+                charge_sigma_low_thr=1.,
+                charge_sigma_high_thr=2.,
+            )
+
+            tool.initialize()
+            tool.setup()
+
+            tool.start()
+            output = tool.finish(return_output_component=True)
+
+            # check output
+            assert output.nsamples == N_SAMPLES
+            assert np.all(output.nevents <= max_events)
+            assert np.shape(output.pixels_id) == (n_pixels,)
+            assert np.shape(output.pedestal_mean_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_mean_lg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_hg) == (n_pixels, N_SAMPLES)
+            assert np.shape(output.pedestal_std_lg) == (n_pixels, N_SAMPLES)
+            assert np.allclose(output.pedestal_mean_hg, 245., atol=20.)
+            assert np.allclose(output.pedestal_mean_lg, 245., atol=20.)
+            assert np.allclose(output.pedestal_std_hg, 10., atol=10.)
+            assert np.allclose(output.pedestal_std_lg, 2.5, atol=3.)
diff --git a/src/nectarchain/makers/component/PedestalComponent.py b/src/nectarchain/makers/component/PedestalComponent.py
index 014fa982..bc879f40 100644
--- a/src/nectarchain/makers/component/PedestalComponent.py
+++ b/src/nectarchain/makers/component/PedestalComponent.py
@@ -301,7 +301,8 @@ def chargeDistributionFilter_mask(self, sigma_low, sigma_high):
 
             # For each event and pixel calculate the charge mean and std over all events
             # taking into account the mask already calculated
-            charge_array = ma.masked_array(self._chargesContainers.charges_hg, self._wfs_mask)
+            charge_array = ma.masked_array(self._chargesContainers.charges_hg,
+                                           np.any(self._wfs_mask, axis=2))
             charge_mean = ma.mean(charge_array, axis=0)
             charge_std = ma.std(charge_array, axis=0)
 
@@ -330,9 +331,9 @@ def chargeDistributionFilter_mask(self, sigma_low, sigma_high):
 
         return new_mask
 
-    def finish(self):
+    def finish(self,*args,**kwargs):
         """
-        Finish the component and, in sliced mode, produce the combined pedestal estimation.
+        Finish the component by filtering the waveforms and calculating the pedestal quantities
         """
 
         # Use only pedestal type events
@@ -375,14 +376,15 @@ def finish(self):
 
             # Time selection
             # set the minimum time
-            print('time filter')
             tmin = np.maximum(
                 self.ucts_tmin or self._waveformsContainers.ucts_timestamp.min(),
                 self._waveformsContainers.ucts_timestamp.min())
+            print('TMIN',self.ucts_tmin,self._waveformsContainers.ucts_timestamp.min(),tmin)
             # set the maximum time
             tmax = np.minimum(
                 self.ucts_tmax or self._waveformsContainers.ucts_timestamp.max(),
                 self._waveformsContainers.ucts_timestamp.max())
+            print('TMAX', self.ucts_tmax, self._waveformsContainers.ucts_timestamp.max(), tmax)
             # Add time selection to mask
             self._wfs_mask = self.timestamp_mask(tmin, tmax)
 

From 3ca8a0bb17417b9cfbd867d38ced8fb4c4903c69 Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Wed, 20 Mar 2024 18:05:03 +0100
Subject: [PATCH 18/20] finished unit test for pedestal tool

---
 .../makers/calibration/tests/test_pedestal_tool.py   | 12 ++++--------
 .../component/tests/test_pedestal_component.py       | 11 -----------
 2 files changed, 4 insertions(+), 19 deletions(-)
 delete mode 100644 src/nectarchain/makers/component/tests/test_pedestal_component.py

diff --git a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
index e188085e..086caf97 100644
--- a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
+++ b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
@@ -32,12 +32,11 @@ def test_base(self):
 
             # run tool
             tool = PedestalNectarCAMCalibrationTool(
-                progress_bar=True,
                 run_number=run_number,
                 run_file=run_file,
                 max_events=max_events,
                 events_per_slice=events_per_slice,
-                log_level=20,
+                log_level=0,
                 output_path=outfile,
                 overwrite=True,
                 filter_method=None,
@@ -112,12 +111,11 @@ def test_timesel(self):
 
             # run tool
             tool = PedestalNectarCAMCalibrationTool(
-                progress_bar=True,
                 run_number=run_number,
                 run_file=run_file,
                 max_events=max_events,
                 events_per_slice=events_per_slice,
-                log_level=20,
+                log_level=0,
                 output_path=outfile,
                 ucts_tmin = tmin,
                 ucts_tmax = tmax,
@@ -160,12 +158,11 @@ def test_WaveformsStdFilter(self):
 
             # run tool
             tool = PedestalNectarCAMCalibrationTool(
-                progress_bar=True,
                 run_number=run_number,
                 run_file=run_file,
                 max_events=max_events,
                 events_per_slice=events_per_slice,
-                log_level=20,
+                log_level=0,
                 output_path=outfile,
                 overwrite=True,
                 filter_method='WaveformsStdFilter',
@@ -206,12 +203,11 @@ def test_ChargeDistributionFilter(self):
 
             # run tool
             tool = PedestalNectarCAMCalibrationTool(
-                progress_bar=True,
                 run_number=run_number,
                 run_file=run_file,
                 max_events=max_events,
                 events_per_slice=events_per_slice,
-                log_level=20,
+                log_level=0,
                 output_path=outfile,
                 overwrite=True,
                 filter_method='ChargeDistributionFilter',
diff --git a/src/nectarchain/makers/component/tests/test_pedestal_component.py b/src/nectarchain/makers/component/tests/test_pedestal_component.py
deleted file mode 100644
index c58a323c..00000000
--- a/src/nectarchain/makers/component/tests/test_pedestal_component.py
+++ /dev/null
@@ -1,11 +0,0 @@
-from nectarchain.makers.component import PedestalEstimationComponent
-from nectarchain.data.container import NectarCAMPedestalContainer
-import pytest
-
-pytest.skip("Not implemented yet",
-            allow_module_level=True,)
-
-# class TestPedestalEstimationComponent:
-#
-#     def test_pedestal_estimation(self):
-

From cde11d83651b43d2ec10f56e4341402b9fa1dbfc Mon Sep 17 00:00:00 2001
From: Luigi Tibaldo <luigi.tibaldo@irap.omp.eu>
Date: Thu, 21 Mar 2024 12:09:02 +0100
Subject: [PATCH 19/20] use test file on CC IN2P3 server

---
 .../calibration/tests/test_pedestal_tool.py     | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
index 086caf97..2172f864 100644
--- a/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
+++ b/src/nectarchain/makers/calibration/tests/test_pedestal_tool.py
@@ -1,17 +1,14 @@
 import tempfile
-import os
+
 import numpy as np
 import tables
-#from ctapipe.utils import get_dataset_path
+from ctapipe.utils import get_dataset_path
 from ctapipe_io_nectarcam.constants import N_SAMPLES
 from nectarchain.data.container import NectarCAMPedestalContainer
 from nectarchain.makers.calibration import PedestalNectarCAMCalibrationTool
 
-# FIXME
-# For the moment rely on file stored on disk located at $NECTARCAMDATA/runs/
 run_number = 3938
-run_file = os.environ['NECTARCAMDATA']+'/runs/NectarCAM.Run3938.30events.fits.fz'
-#run_file = get_dataset_path('NectarCAM.Run3938.30events.fits.fz')
+run_file = get_dataset_path('NectarCAM.Run3938.30events.fits.fz')
 # number of working pixels in run
 n_pixels = 1834
 
@@ -72,7 +69,7 @@ def test_base(self):
                     assert group_name in h5file.root.__members__
                     table = h5file.root[group_name][NectarCAMPedestalContainer.__name__][0]
                     assert table['nsamples'] == N_SAMPLES
-                    #assert np.all(table['nevents'] == events_per_slice)
+                    assert np.all(table['nevents'] == events_per_slice)
                     assert np.shape(table['pixels_id']) == (n_pixels,)
                     assert np.shape(table['pedestal_mean_hg']) == (n_pixels, N_SAMPLES)
                     assert np.shape(table['pedestal_mean_lg']) == (n_pixels, N_SAMPLES)
@@ -117,8 +114,8 @@ def test_timesel(self):
                 events_per_slice=events_per_slice,
                 log_level=0,
                 output_path=outfile,
-                ucts_tmin = tmin,
-                ucts_tmax = tmax,
+                ucts_tmin=tmin,
+                ucts_tmax=tmax,
                 overwrite=True,
                 filter_method=None,
             )
@@ -166,7 +163,7 @@ def test_WaveformsStdFilter(self):
                 output_path=outfile,
                 overwrite=True,
                 filter_method='WaveformsStdFilter',
-                wfs_std_threshold = 4.,
+                wfs_std_threshold=4.,
             )
 
             tool.initialize()

From 24b0e014b7819f79572b70110e89ad9c6cd9bad0 Mon Sep 17 00:00:00 2001
From: Jean-Philippe Lenain <jlenain@users.noreply.github.com>
Date: Fri, 22 Mar 2024 09:33:21 +0100
Subject: [PATCH 20/20] `pytest` failing on `ctapipe` tools in CI (#117)

* Activate the `pytest-xdist` addon for `pytest` in CI.

* Fix typo

---------

Co-authored-by: Jean-Philippe Lenain <jlenain@in2p3.fr>
---
 .github/workflows/ci.yml | 7 +++----
 pyproject.toml           | 1 +
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 4de41bba..6d8ec86e 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -77,14 +77,13 @@ jobs:
           echo "Installing additional pip packages"
           # we install ctapipe using pip to be able to select any commit, e.g. the current main
           pip install \
-             "git+https://github.com/cta-observatory/ctapipe@$CTAPIPE_VERSION" \
-              pytest-cov
+             "git+https://github.com/cta-observatory/ctapipe@$CTAPIPE_VERSION"
 
           echo "pip install -e ."
-          pip install -e .
+          pip install -e .[test]
 
       - name: Tests
         run: |
-          pytest --cov=nectarchain --cov-report=xml
+          pytest -n auto --dist loadscope --cov=nectarchain --cov-report=xml
 
       - uses: codecov/codecov-action@v3
diff --git a/pyproject.toml b/pyproject.toml
index 4286470f..c0833017 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,6 +39,7 @@ dynamic = ["version"]
 test = [
     "pytest",
     "pytest-cov",
+    "pytest-xdist",
 ]
 dev = [
     "setuptools_scm",