diff --git a/lstchain/scripts/lstchain_convert_drs4_pedestal_to_evb.py b/lstchain/scripts/lstchain_convert_drs4_pedestal_to_evb.py index 0e0b38d187..c0df690946 100644 --- a/lstchain/scripts/lstchain_convert_drs4_pedestal_to_evb.py +++ b/lstchain/scripts/lstchain_convert_drs4_pedestal_to_evb.py @@ -59,8 +59,8 @@ def main(): pedestal = read_pedestal(args.drs4_baseline_file, args.tel_id) - if pedestal.dtype != np.uint16: - print(f'Input pedestal must have dtype uint16, got {pedestal.dtype}.', file=sys.stderr) + if not (pedestal.dtype == np.uint16 or pedestal.dtype == np.int16): + print(f'Input pedestal must have dtype (u)int16, got {pedestal.dtype}.', file=sys.stderr) print('Make sure the pedestal file was generated by lstchain > 0.8.2', file=sys.stderr) sys.exit(1) diff --git a/lstchain/tools/lstchain_create_drs4_pedestal_file.py b/lstchain/tools/lstchain_create_drs4_pedestal_file.py index d332acded2..96d64d15ab 100644 --- a/lstchain/tools/lstchain_create_drs4_pedestal_file.py +++ b/lstchain/tools/lstchain_create_drs4_pedestal_file.py @@ -14,7 +14,7 @@ flag, ) from ctapipe.io.hdf5tableio import HDF5TableWriter -from ctapipe_io_lst import LSTEventSource +from ctapipe_io_lst import LSTEventSource, EVBPreprocessing, TriggerBits from ctapipe_io_lst.calibration import get_spike_A_positions from ctapipe_io_lst.constants import ( N_CAPACITORS_PIXEL, @@ -31,13 +31,13 @@ class DRS4CalibrationContainer(Container): baseline_mean = Field( None, "Mean baseline of each capacitor, shape (N_GAINS, N_PIXELS, N_CAPACITORS_PIXEL)", - dtype=np.uint16, + dtype=np.int16, ndim=3, ) baseline_std = Field( None, "Std Dev. of the baseline calculation, shape (N_GAINS, N_PIXELS, N_CAPACITORS_PIXEL)", - dtype=np.uint16, + dtype=np.int16, ndim=3, ) baseline_counts = Field( @@ -51,14 +51,16 @@ class DRS4CalibrationContainer(Container): None, "Mean spike height for each pixel, shape (N_GAINS, N_PIXELS, 3)", ndim=3, - dtype=np.uint16, + dtype=np.int16, ) -def convert_to_uint16(array): - '''Convert an array to uint16, rounding and clipping before to avoid under/overflow''' - array = np.clip(np.round(array), 0, np.iinfo(np.uint16).max) - return array.astype(np.uint16) +def convert_to_int16(array): + '''Convert an array to int16, rounding and clipping before to avoid under/overflow''' + dtype = np.int16 + info = np.iinfo(dtype) + array = np.clip(np.round(array), info.min, info.max) + return array.astype(dtype) @numba.njit(cache=True, inline='always') @@ -257,6 +259,65 @@ def mean_spike_height(self): return spike_heights + def check_baseline_values(self, baseline_mean): + """ + Check the values of the drs4 baseline for issues. + + In case of no EVBv5 data or no pre-calibration by EVBv6, + we expect no negative and no small values. + + In case of EVBv6 applied baseline corrections, values should + be close to 0, but negative values are ok. + """ + + check_large_values = False + check_negative = True + check_small = True + + if self.source.cta_r1: + preprocessing = self.source.evb_preprocessing[TriggerBits.MONO] + if EVBPreprocessing.BASELINE_SUBTRACTION in preprocessing: + check_large_values = True + check_negative = False + check_small = False + + + if check_negative: + negative = baseline_mean < 0 + n_negative = np.count_nonzero(negative) + if n_negative > 0: + gain, pixel, capacitor = np.nonzero(negative) + self.log.critical(f'{n_negative} baseline values are smaller than 0') + self.log.info("Gain | Pixel | Capacitor | Baseline ") + for g, p, c in zip(gain, pixel, capacitor): + self.log.info( + f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}" + ) + + if check_small: + small = baseline_mean < 25 + n_small = np.count_nonzero(small) + if n_small > 0: + gain, pixel, capacitor = np.nonzero(small) + self.log.warning(f'{n_small} baseline values are smaller than 25') + self.log.info("Gain | Pixel | Capacitor | Baseline ") + for g, p, c in zip(gain, pixel, capacitor): + self.log.info( + f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}" + ) + + if check_large_values: + large = np.abs(baseline_mean) > 25 + n_large = np.count_nonzero(large) + if n_large > 0: + gain, pixel, capacitor = np.nonzero(large) + self.log.warning(f'{n_large} baseline values have an abs value >= 25') + self.log.info("Gain | Pixel | Capacitor | Baseline ") + for g, p, c in zip(gain, pixel, capacitor): + self.log.info( + f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}" + ) + def finish(self): tel_id = self.source.tel_id self.log.info('Writing output to %s', self.output_path) @@ -267,28 +328,13 @@ def finish(self): baseline_std = self.baseline_stats.std.reshape(shape) baseline_counts = self.baseline_stats.counts.reshape(shape).astype(np.uint16) - n_negative = np.count_nonzero(baseline_mean < 0) - if n_negative > 0: - gain, pixel, capacitor = np.nonzero(baseline_mean < 0) - self.log.critical(f'{n_negative} baseline values are smaller than 0') - self.log.info("Gain | Pixel | Capacitor | Baseline ") - for g, p, c in zip(gain, pixel, capacitor): - self.log.info(f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}") - - n_small = np.count_nonzero(baseline_mean < 25) - if n_small > 0: - gain, pixel, capacitor = np.nonzero(baseline_mean < 25) - self.log.warning(f'{n_small} baseline values are smaller than 25') - self.log.info("Gain | Pixel | Capacitor | Baseline ") - for g, p, c in zip(gain, pixel, capacitor): - self.log.info(f"{g:4d} | {p:4d} | {c:9d} | {baseline_mean[g][p][c]:6.1f}") - + self.check_baseline_values(baseline_mean) # Convert baseline mean and spike heights to uint16, handle missing # values and values smaller 0, larger maxint - baseline_mean = convert_to_uint16(baseline_mean) - baseline_std = convert_to_uint16(baseline_std) - spike_height = convert_to_uint16(self.mean_spike_height()) + baseline_mean = convert_to_int16(baseline_mean) + baseline_std = convert_to_int16(baseline_std) + spike_height = convert_to_int16(self.mean_spike_height()) with HDF5TableWriter(self.output_path) as writer: Provenance().add_output_file(str(self.output_path)) diff --git a/lstchain/tools/tests/test_create_drs4_pedestal_file.py b/lstchain/tools/tests/test_create_drs4_pedestal_file.py index a7e0477c08..e83cdaeca5 100644 --- a/lstchain/tools/tests/test_create_drs4_pedestal_file.py +++ b/lstchain/tools/tests/test_create_drs4_pedestal_file.py @@ -35,12 +35,12 @@ def test_create_drs4_pedestal_file(tmp_path): baseline_mean = drs4_data['baseline_mean'] baseline_counts = drs4_data['baseline_std'] - assert baseline_mean.dtype == np.uint16 + assert baseline_mean.dtype == np.int16 assert baseline_mean.shape == (N_GAINS, N_PIXELS, N_CAPACITORS_PIXEL) assert np.isclose(np.average(baseline_mean[baseline_counts > 0], weights=baseline_counts[baseline_counts > 0]), 400, rtol=0.05) spike_height = drs4_data['spike_height'] - assert spike_height.dtype == np.uint16 + assert spike_height.dtype == np.int16 mean_spike_height = np.nanmean(spike_height, axis=(0, 1)) # these are the expected spike heights, but due to the low statistics,