diff --git a/src/nectarchain/data/container/flatfieldContainer.py b/src/nectarchain/data/container/flatfieldContainer.py new file mode 100644 index 00000000..2ebb8e72 --- /dev/null +++ b/src/nectarchain/data/container/flatfieldContainer.py @@ -0,0 +1,85 @@ +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 ctapipe.core.traits import ComponentNameList, Dict, Integer, List, Tuple + +from .core import NectarCAMContainer + +__all__ = ["FlatFieldContainer"] + + +class FlatFieldContainer(NectarCAMContainer): + """ + Container that holds flat field coefficients and other useful information + + Fields: + run_number (np.uint16): Number of the run + npixels (np.uint16): Number of pixels + pixels_id (np.ndarray): Array of pixel's ID + ucts_timestamp (np.ndarray) : Array of time stamps of each event (UTC) + event_type (np.ndarray): Array of trigger event types (should be all flat field events) + event_id (np.ndarray): Array of the IDs of each event + amp_int_per_pix_per_event (np.ndarray): Array of integrated amplitude of each pulse + t_peak_per_pix_per_event (np.ndarray): Array of samples containing the pulse maximum + FF_coef (np.ndarray): Array of flat field coefficients + bad_pixels (List): List of pixel identified as outliers + """ + + run_number = Field( + type=np.uint16, + description="run number associated to the waveforms", + ) + + npixels = Field( + type=np.uint16, + description="number of effective pixels", + ) + + pixels_id = Field(type=np.ndarray, dtype=np.uint16, ndim=1, description="pixel ids") + + ucts_timestamp = Field( + type=np.ndarray, dtype=np.uint64, ndim=1, description="events ucts timestamp" + ) + + event_type = Field( + type=np.ndarray, dtype=np.uint8, ndim=1, description="trigger event type" + ) + + event_id = Field(type=np.ndarray, dtype=np.uint32, ndim=1, description="event ids") + + amp_int_per_pix_per_event = Field( + type=np.ndarray, + dtype=np.float32, + ndim=3, + description="The amplitude integrated over the window width, per pixel and per event", + ) + + t_peak_per_pix_per_event = Field( + type=np.ndarray, + dtype=np.uint8, + ndim=3, + description="Sample containing the pulse maximum, per pixel and per event", + ) + + FF_coef = Field( + type=np.ndarray, + dtype=np.float32, + ndim=3, + description="The flat field coefficient, per event", + ) + + # masked_wfs = Field( + # type=np.ndarray, dtype=np.uint64, ndim=4, description="Masked array for amplitude integration" + # ) + + bad_pixels = Field( + type=List, + dtype=np.uint16, + ndim=1, + description="Pixels considered as bad in at least one gain channels", + ) diff --git a/src/nectarchain/makers/calibration/flatfieldMakers.py b/src/nectarchain/makers/calibration/flatfieldMakers.py index 01da9958..2a735e70 100644 --- a/src/nectarchain/makers/calibration/flatfieldMakers.py +++ b/src/nectarchain/makers/calibration/flatfieldMakers.py @@ -9,8 +9,30 @@ __all__ = ["FlatfieldNectarCAMCalibrationTool"] -class FlatfieldNectarCAMCalibrationTool(NectarCAMCalibrationTool): - def start(self): - raise NotImplementedError( - "The computation of the flatfield calibration is not yet implemented, feel free to contribute !:)" +#class FlatfieldNectarCAMCalibrationTool(NectarCAMCalibrationTool): +# def start(self): +# raise NotImplementedError( +# "The computation of the flatfield calibration is not yet implemented, feel free to contribute !:)" +# ) + +from nectarchain.makers import EventsLoopNectarCAMCalibrationTool + + +class FlatfieldNectarCAMCalibrationTool(EventsLoopNectarCAMCalibrationTool): + name = "NectarCAM" + + componentsList = ComponentNameList( + NectarCAMComponent, + default_value=["preFlatFieldComponent"], + help="List of Component names to be apply, the order will be respected", + ).tag(config=True) + + def _init_output_path(self): + if self.max_events is None: + filename = f"{self.name}_run{self.run_number}.h5" + else: + filename = f"{self.name}_run{self.run_number}_maxevents{self.max_events}.h5" + self.output_path = pathlib.Path( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/FlatFieldTests/{filename}" ) + diff --git a/src/nectarchain/makers/component/preFlatFieldComponent.py b/src/nectarchain/makers/component/preFlatFieldComponent.py new file mode 100644 index 00000000..20a66249 --- /dev/null +++ b/src/nectarchain/makers/component/preFlatFieldComponent.py @@ -0,0 +1,152 @@ +import logging + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import os +import pathlib + +import matplotlib.pyplot as plt +import numpy as np +from ctapipe.containers import EventType, Field +from ctapipe.coordinates import EngineeringCameraFrame +from ctapipe.core import Component +from ctapipe.core.traits import ComponentNameList, Dict, Float, Integer, List, Tuple +from ctapipe.instrument import CameraGeometry +from ctapipe.io import HDF5TableReader +from ctapipe.visualization import CameraDisplay +from ctapipe_io_nectarcam import constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer + +from nectarchain.data.container import ( + ArrayDataContainer, + NectarCAMContainer, + TriggerMapContainer, +) +from nectarchain.makers import EventsLoopNectarCAMCalibrationTool +from nectarchain.makers.component import ArrayDataComponent, NectarCAMComponent +from nectarchain.utils import ComponentUtils + +__all__ = ["preFlatFieldComponent"] + + +class preFlatFieldComponent(NectarCAMComponent): + window_shift = Integer( + default_value=5, + help="the time in ns before the peak to integrate charge", + ).tag(config=True) + + window_width = Integer( + default_value=12, + help="the duration of the extraction window in ns", + ).tag(config=True) + ## --< final window is 14 samples ! >-- + + g = Float( + default_value=58.0, + help="defaut gain value", + ).tag(config=True) + + hi_lo_ratio = Float( + default_value=13.0, + help="defaut gain value", + ).tag(config=True) + + ## Simple function to substract the pedestal using the 20 first samples of each trace + def substract_pedestal(wfs, window=20): + ped_mean = np.mean(wfs[0][:, :, 0:window], axis=2) + wfs_pedsub = wfs - np.expand_dims(ped_mean, axis=-1) + return wfs_pedsub + + ## Function to make an array that will be used as a mask on the waveform for the calculation of the integrated amplitude of the signal. + def make_masked_array(t_peak, window_shift, window_width): + masked_wfs = [ + np.array([np.zeros(constants.N_SAMPLES)] * constants.N_PIXELS) + ] * constants.N_GAINS + masked_wfs = np.array(masked_wfs) + + t_signal_start = t_peak - window_shift + t_signal_stop = t_peak + window_width - window_shift + + for g in range(0, constants.N_GAINS): + for i in range(0, constants.N_PIXELS): + masked_wfs[g][i][ + t_signal_start[0, g, i] : t_signal_stop[0, g, i] + ] = True + # --< I didn't find a better way to do than using this masked array >-- + + return masked_wfs + + def __init__(self, subarray, config=None, parent=None, *args, **kwargs): + super().__init__( + subarray=subarray, config=config, parent=parent, *args, **kwargs + ) + + self.__ucts_timestamp = [] + self.__event_type = [] + self.__event_id = [] + self.__amp_int_per_pix_per_event = [] + self.__FF_coef = [] + + def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): + wfs = [] + wfs_pedsub = [] + + if event.trigger.event_type.value == EventType.FLATFIELD.value: + # print("event :", (self.__event_id, self.__event_type)) + self.__event_id.append(np.uint32(event.index.event_id)) + self.__event_type.append(event.trigger.event_type.value) + self.__ucts_timestamp.append(event.nectarcam.tel[0].evt.ucts_timestamp) + + wfs.append(event.r0.tel[0].waveform) # not saved + + # substract pedestal using the mean of the 20 first samples + wfs_pedsub = substract_pedestal(wfs, 20) + + # get the masked array for integration window + t_peak = np.argmax(wfs_pedsub, axis=3) + masked_wfs = make_masked_array(t_peak, self.window_shift, self.window_width) + + # get integrated amplitude and mean amplitude over all pixels per event + amp_int_per_pix_per_event = np.sum( + wfs_pedsub[0], axis=2, where=masked_wfs.astype("bool") + ) + self.__amp_int_per_pix_per_event.append(amp_int_per_pix_per_event) + mean_amp_cam_per_event = np.mean(amp_int_per_pix_per_event, axis=-1) + + # get efficiency and flat field coefficient + gain = [self.g, self.g / self.hi_lo_ratio] + + eff = np.divide( + (amp_int_per_pix_per_event[:] / (np.expand_dims(gain[:], axis=-1))), + np.expand_dims((mean_amp_cam_per_event[:] / gain[:]), axis=-1), + ) # efficacité relative + FF_coef = np.ma.array(1.0 / eff, mask=eff == 0) + self.__FF_coef.append(FF_coef) + + def finish(self): + output = FlatFieldContainer( + run_number=FlatFieldContainer.fields["run_number"].type( + self._run_number + ), + npixels=FlatFieldContainer.fields["npixels"].type(self._npixels), + pixels_id=FlatFieldContainer.fields["pixels_id"].dtype.type( + self._pixels_id + ), + ucts_timestamp=FlatFieldContainer.fields["ucts_timestamp"].dtype.type( + self.__ucts_timestamp + ), + event_type=FlatFieldContainer.fields["event_type"].dtype.type( + self.__event_type + ), + event_id=FlatFieldContainer.fields["event_id"].dtype.type( + self.__event_id + ), + amp_int_per_pix_per_event=FlatFieldContainer.fields[ + "amp_int_per_pix_per_event" + ].dtype.type(self.__amp_int_per_pix_per_event), + FF_coef=FlatFieldContainer.fields["FF_coef"].dtype.type(self.__FF_coef), + ) + return output + diff --git a/src/nectarchain/user_scripts/ajardinb/FlatField_test.py b/src/nectarchain/user_scripts/ajardinb/FlatField_test.py new file mode 100644 index 00000000..a121b196 --- /dev/null +++ b/src/nectarchain/user_scripts/ajardinb/FlatField_test.py @@ -0,0 +1,37 @@ +import os +import pathlib + +import matplotlib.pyplot as plt + +# %% +import numpy as np +from ctapipe.containers import Field +from ctapipe.core import Component +from ctapipe.core.traits import ComponentNameList, Integer, List, Dict, Tuple +from ctapipe.io import HDF5TableReader +from ctapipe_io_nectarcam import constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer +from ctapipe.instrument import CameraGeometry +from ctapipe.visualization import CameraDisplay +from ctapipe.coordinates import EngineeringCameraFrame + +from nectarchain.data.container import ( + ArrayDataContainer, + NectarCAMContainer, + TriggerMapContainer, +) +from ctapipe.containers import EventType +from nectarchain.makers import EventsLoopNectarCAMCalibrationTool +from nectarchain.makers.component import ArrayDataComponent, NectarCAMComponent +from nectarchain.utils import ComponentUtils + +from ctapipe.io import EventSource, EventSeeker +from ctapipe_io_nectarcam import NectarCAMEventSource +from astropy import units as u +from astropy.time import Time + +from scipy import signal +from scipy.signal import find_peaks +from scipy.interpolate import Rbf, InterpolatedUnivariateSpline +from scipy.stats import chi2, norm +import scipy