Add container/component/tool for pedestal estimation (#115)

* create structure for pedestal component/tool

* debugging tool structure

* fix computation of mean, etc

* Added container for pedestals

* Minimal working tool/component

* implemented time selection

* Implemented filtering based on waveforms standard deviation

* Add possibility to filter based on charge distribution

* Restructured to use events_per_slice to reduce memory load

* fully working version based on events_per_slice mechanism

* completed doc strings

* added tests for pedestal container

* First go at unit test for pedestal tool

* Unit tests for pedestal maker

* finished unit test for pedestal tool

* use test file on CC IN2P3 server

---------

Co-authored-by: Luigi Tibaldo <[email protected]>

src/nectarchain/data/container/__init__.py

@@ -9,3 +9,4 @@
 from .eventSource import *
 from .gainContainer import *
 from .waveformsContainer import *
+from .pedestalContainer import *

src/nectarchain/data/container/pedestalContainer.py

@@ -0,0 +1,66 @@
+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 (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_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.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")
+
+    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_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"
+    )

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

src/nectarchain/makers/calibration/pedestalMakers.py

@@ -1,16 +1,148 @@
 import logging
 
+import os
+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"]
 
 
 class PedestalNectarCAMCalibrationTool(NectarCAMCalibrationTool):
-    def start(self):
-        raise NotImplementedError(
-            "The computation of the pedestal calibration is not yet implemented, feel free to contribute !:)"
-        )
+    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)
+
+    def _init_output_path(self):
+        """
+        Initialize output path
+        """
+
+        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 _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):
+        """
+        Redefines 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