Skip to content

Commit

Permalink
implement of comments and cleaning
Browse files Browse the repository at this point in the history
  • Loading branch information
ArmelleJB committed Jan 14, 2025
1 parent fdbbeca commit 457e057
Showing 1 changed file with 96 additions and 65 deletions.
161 changes: 96 additions & 65 deletions src/nectarchain/makers/component/preFlatFieldComponent.py
Original file line number Diff line number Diff line change
@@ -1,37 +1,41 @@
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.containers import EventType
from ctapipe.core.traits import Float, Integer
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
from nectarchain.data.container import FlatFieldContainer
from nectarchain.makers.component import NectarCAMComponent

logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s")
log = logging.getLogger(__name__)
log.handlers = logging.getLogger("__main__").handlers

__all__ = ["preFlatFieldComponent"]


class preFlatFieldComponent(NectarCAMComponent):
"""
Component that computes flat field coefficients from raw data.
Parameters
----------
window_shift: int
time in ns before the peak to integrate charge (default value = 5)
window_width: int
duration of the extraction window in ns (default value = 12)
g: float
default gain value (default value = 58)
hi_lo_ratio: float
default high gain to low gain ratio (default value = 13)
"""

window_shift = Integer(
default_value=5,
help="the time in ns before the peak to integrate charge",
Expand All @@ -41,43 +45,18 @@ class preFlatFieldComponent(NectarCAMComponent):
default_value=12,
help="the duration of the extraction window in ns",
).tag(config=True)
## --< final window is 14 samples ! >--
# --< final window is 14 samples ! >--

g = Float(
default_value=58.0,
help="defaut gain value",
help="default gain value",
).tag(config=True)

hi_lo_ratio = Float(
default_value=13.0,
help="defaut gain value",
help="default high gain to low gain ratio",
).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
Expand All @@ -101,35 +80,90 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs):

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)
# subtract pedestal using the mean of the 20 first samples
wfs_pedsub = self.subtract_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)
masked_wfs = self.make_masked_array(
t_peak, self.window_shift, self.window_width
)
# --< I didn't find a better way to do than using this masked array >--

# 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)
# 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]

amp_int_per_pix_per_event_pe = amp_int_per_pix_per_event[:] / (
np.expand_dims(gain[:], axis=-1)
)
mean_amp_cam_per_event_pe = np.mean(amp_int_per_pix_per_event_pe, axis=-1)

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
amp_int_per_pix_per_event_pe,
np.expand_dims(mean_amp_cam_per_event_pe, axis=-1),
)

FF_coef = np.ma.array(1.0 / eff, mask=eff == 0)
self.__FF_coef.append(FF_coef)

@staticmethod
def subtract_pedestal(wfs, window=20):
"""
Subtract the pedestal defined as the average of the first samples of each trace
Args:
wfs: raw wavefroms
window: number of samples n to calculate the pedestal (default value is 20)
Returns:
wfs_pedsub: wavefroms subtracted from the pedestal
"""

ped_mean = np.mean(wfs[0][:, :, 0:window], axis=2)
wfs_pedsub = wfs - np.expand_dims(ped_mean, axis=-1)

return wfs_pedsub

@staticmethod
def make_masked_array(t_peak, window_shift, window_width):
"""
Define an array that will be used as a mask on the waveforms for the calculation
of the integrated amplitude of the signal
Args:
t_peak: time corresponding the the highest peak of the trace
window_shift: time in ns before the peak to integrate charge
window_width: duration of the extraction window in ns
Returns:
masked_wfs: a mask array
"""

masked_wfs = np.zeros(
shape=(constants.N_GAINS, constants.N_PIXELS, constants.N_SAMPLES),
dtype=bool,
)
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

return masked_wfs

def finish(self):
output = FlatFieldContainer(
run_number=FlatFieldContainer.fields["run_number"].type(
self._run_number
),
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
Expand All @@ -140,13 +174,10 @@ def finish(self):
event_type=FlatFieldContainer.fields["event_type"].dtype.type(
self.__event_type
),
event_id=FlatFieldContainer.fields["event_id"].dtype.type(
self.__event_id
),
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

0 comments on commit 457e057

Please sign in to comment.