From b1c862a81043efdbacdd686e20be1ebfcb44e63a Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 19 Mar 2024 18:23:37 +0100 Subject: [PATCH 01/42] Initial commit --- src/ctapipe/containers.py | 7 +- src/ctapipe/image/extractor.py | 151 ++++++++++++++++------ src/ctapipe/image/reducer.py | 2 +- src/ctapipe/image/tests/test_extractor.py | 82 ++++++------ src/ctapipe/image/tests/test_reducer.py | 14 +- src/ctapipe/image/toymodel.py | 4 +- src/ctapipe/io/simteleventsource.py | 8 +- 7 files changed, 172 insertions(+), 96 deletions(-) diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 067cd8b6c6e..19045a971e1 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -1,6 +1,7 @@ """ Container structures for data that should be read or written to disk """ + import enum from functools import partial @@ -559,7 +560,7 @@ class R0CameraContainer(Container): """ waveform = Field( - None, ("numpy array containing ADC samples" "(n_channels, n_pixels, n_samples)") + None, ("numpy array containing ADC samples: (n_channels, n_pixels, n_samples)") ) @@ -586,7 +587,7 @@ class R1CameraContainer(Container): None, ( "numpy array containing a set of images, one per ADC sample" - "Shape: (n_pixels, n_samples)" + "Shape: (n_channels, n_pixels, n_samples)" ), ) @@ -660,7 +661,7 @@ class DL0CameraContainer(Container): ( "numpy array containing data volume reduced " "p.e. samples" - "(n_pixels, n_samples). Note this may be a masked array, " + "(n_channels, n_pixels, n_samples). Note this may be a masked array, " "if pixels or time slices are zero-suppressed" ), ) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index be920e0a0d4..caa8772cd31 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -84,7 +84,7 @@ def extract_around_peak( ---------- waveforms : ndarray Waveforms stored in a numpy array. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix, n_samples) peak_index : ndarray or int Peak index for each pixel. width : ndarray or int @@ -105,7 +105,10 @@ def extract_around_peak( ------- charge : ndarray Extracted charge. - Shape: (n_pix) + Shape: (n_channels, n_pix) + peak_time: ndarray + Extracted peak time. + Shape: (n_channels, n_pix) """ n_samples = waveforms.size @@ -162,7 +165,7 @@ def extract_sliding_window(waveforms, width, sampling_rate_ghz, sum_, peak_time) ---------- waveforms : ndarray Waveforms stored in a numpy array. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix, n_samples) width : ndarray or int Window size of integration window for each pixel. sampling_rate_ghz : float @@ -179,7 +182,10 @@ def extract_sliding_window(waveforms, width, sampling_rate_ghz, sum_, peak_time) ------- charge : ndarray Extracted charge. - Shape: (n_pix) + Shape: (n_channels, n_pix) + peak_time: ndarray + Extracted peak time. + Shape: (n_channels, n_pix) """ @@ -215,7 +221,7 @@ def neighbor_average_maximum( ---------- waveforms : ndarray Waveforms stored in a numpy array. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix, n_samples) neighbors_indices : ndarray indices of a scipy csr sparse matrix of neighbors, i.e. ``ctapipe.instrument.CameraGeometry.neighbor_matrix_sparse.indices``. @@ -233,29 +239,30 @@ def neighbor_average_maximum( ------- average_wf : ndarray Average of neighbor waveforms for each pixel. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix) """ - n_pixels = waveforms.shape[0] + n_channels, n_pixels, _ = waveforms.shape indptr = neighbors_indptr indices = neighbors_indices # initialize to waveforms weighted with local_weight # so the value of the pixel itself is already taken into account - peak_pos = np.empty(n_pixels, dtype=np.int64) + peak_pos = np.empty((n_channels, n_pixels), dtype=np.int64) for pixel in prange(n_pixels): - average = waveforms[pixel] * local_weight + average = waveforms[:, pixel] * local_weight neighbors = indices[indptr[pixel] : indptr[pixel + 1]] for neighbor in neighbors: if broken_pixels[neighbor]: continue - average += waveforms[neighbor] + for ichannel in range(n_channels): + average[ichannel] += waveforms[ichannel][neighbor] - peak_pos[pixel] = np.argmax(average) + peak_pos[:, pixel] = np.argmax(average, axis=-1) return peak_pos @@ -269,7 +276,7 @@ def subtract_baseline(waveforms, baseline_start, baseline_end): ---------- waveforms : ndarray Waveforms stored in a numpy array. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix, n_samples) baseline_start : int Sample where the baseline window starts baseline_end : int @@ -279,6 +286,7 @@ def subtract_baseline(waveforms, baseline_start, baseline_end): ------- baseline_corrected : ndarray Waveform with the baseline subtracted + Shape: (n_channels, n_pix, n_samples) """ baseline_corrected = ( waveforms @@ -397,7 +405,7 @@ def __call__( ---------- waveforms : ndarray Waveforms stored in a numpy array of shape - (n_pix, n_samples). + (n_channels, n_pix, n_samples). tel_id : int The telescope id. Used to obtain to correct traitlet configuration and instrument properties @@ -424,7 +432,9 @@ def __call__( charge, peak_time = extract_around_peak( waveforms, 0, waveforms.shape[-1], 0, self.sampling_rate_ghz[tel_id] ) - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1CameraContainer( + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True + ) class FixedWindowSum(ImageExtractor): @@ -451,7 +461,7 @@ class FixedWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): """ - Calculate the correction for the extracted change such that the value + Calculate the correction for the extracted charge such that the value returned would equal 1 for a noise-less unit pulse. This method is decorated with @lru_cache to ensure it is only @@ -488,8 +498,17 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + if selected_gain_channel is not None and waveforms.shape[-3] == 1: + charge *= self._calculate_correction(tel_id=tel_id)[ + selected_gain_channel + ] + else: + charge *= self._calculate_correction(tel_id=tel_id)[ + :, np.newaxis, np.newaxis + ] + return DL1CameraContainer( + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True + ) class GlobalPeakWindowSum(ImageExtractor): @@ -530,7 +549,7 @@ class GlobalPeakWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): """ - Calculate the correction for the extracted change such that the value + Calculate the correction for the extracted charge such that the value returned would equal 1 for a noise-less unit pulse. This method is decorated with @lru_cache to ensure it is only @@ -561,26 +580,35 @@ def __call__( ) -> DL1CameraContainer: if self.pixel_fraction.tel[tel_id] == 1.0: # average over pixels then argmax over samples - peak_index = waveforms[~broken_pixels].mean(axis=-2).argmax() + peak_index = waveforms[:, ~broken_pixels].mean(axis=-2).argmax(axis=-1) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - brightest = np.argsort(waveforms.max(axis=-1))[~broken_pixels][ + brightest = np.argsort(waveforms.max(axis=-1))[:, ~broken_pixels][ ..., -n_pixels: ] # average over brightest pixels then argmax over samples - peak_index = waveforms[brightest].mean(axis=-2).argmax() + peak_index = waveforms[:, brightest].mean(axis=-2).argmax(axis=-1) charge, peak_time = extract_around_peak( waveforms, - peak_index, + peak_index[:, np.newaxis], self.window_width.tel[tel_id], self.window_shift.tel[tel_id], self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + if selected_gain_channel is not None and waveforms.shape[-3] == 1: + charge *= self._calculate_correction(tel_id=tel_id)[ + selected_gain_channel + ] + else: + charge *= self._calculate_correction(tel_id=tel_id)[ + :, np.newaxis, np.newaxis + ] + return DL1CameraContainer( + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True + ) class LocalPeakWindowSum(ImageExtractor): @@ -606,7 +634,7 @@ class LocalPeakWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): """ - Calculate the correction for the extracted change such that the value + Calculate the correction for the extracted charge such that the value returned would equal 1 for a noise-less unit pulse. This method is decorated with @lru_cache to ensure it is only @@ -644,8 +672,17 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + if selected_gain_channel is not None and waveforms.shape[-3] == 1: + charge *= self._calculate_correction(tel_id=tel_id)[ + selected_gain_channel + ] + else: + charge *= self._calculate_correction(tel_id=tel_id)[ + :, np.newaxis, np.newaxis + ] + return DL1CameraContainer( + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True + ) class SlidingWindowMaxSum(ImageExtractor): @@ -720,8 +757,17 @@ def __call__( waveforms, self.window_width.tel[tel_id], self.sampling_rate_ghz[tel_id] ) if self.apply_integration_correction.tel[tel_id]: - charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + if selected_gain_channel is not None and waveforms.shape[-3] == 1: + charge *= self._calculate_correction(tel_id=tel_id)[ + selected_gain_channel + ] + else: + charge *= self._calculate_correction(tel_id=tel_id)[ + :, np.newaxis, np.newaxis + ] + return DL1CameraContainer( + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True + ) class NeighborPeakWindowSum(ImageExtractor): @@ -753,7 +799,7 @@ class NeighborPeakWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): """ - Calculate the correction for the extracted change such that the value + Calculate the correction for the extracted charge such that the value returned would equal 1 for a noise-less unit pulse. This method is decorated with @lru_cache to ensure it is only @@ -798,8 +844,17 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - charge *= self._calculate_correction(tel_id=tel_id)[selected_gain_channel] - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + if selected_gain_channel is not None and waveforms.shape[-3] == 1: + charge *= self._calculate_correction(tel_id=tel_id)[ + selected_gain_channel + ] + else: + charge *= self._calculate_correction(tel_id=tel_id)[ + :, np.newaxis, np.newaxis + ] + return DL1CameraContainer( + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True + ) class BaselineSubtractedNeighborPeakWindowSum(NeighborPeakWindowSum): @@ -836,7 +891,7 @@ class TwoPassWindowSum(ImageExtractor): the range of the sliding is the one allowing extension from 3 to 5; add 1 sample on each side and integrate charge in the 5-sample window; time is obtained as a charge-weighted average of the sample numbers; - No information from neighboouring pixels is used. + No information from neighbouring pixels is used. #. Preliminary image cleaning via simple tailcut with minimum number of core neighbours set at 1, #. Only the brightest cluster of pixels is kept. @@ -1265,6 +1320,13 @@ def _apply_second_pass( def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels ) -> DL1CameraContainer: + if waveforms.shape[-3] != 1: + raise AttributeError( + "The data needs to be gain selected to use the TwoPassWindowSum." + ) + else: + waveforms = waveforms[0, :, :] + charge1, pulse_time1, correction1 = self._apply_first_pass(waveforms, tel_id) # FIXME: properly make sure that output is 32Bit instead of downcasting here @@ -1422,8 +1484,8 @@ def __filtfilt_fast(signal, filt): filtfilt has some speed issues (https://github.com/scipy/scipy/issues/17080) """ forward = convolve1d(signal, filt, axis=-1, mode="nearest") - backward = convolve1d(forward[:, ::-1], filt, axis=-1, mode="nearest") - return backward[:, ::-1] + backward = convolve1d(forward[..., ::-1], filt, axis=-1, mode="nearest") + return backward[..., ::-1] def deconvolve( @@ -1441,7 +1503,7 @@ def deconvolve( ---------- waveforms : ndarray Waveforms stored in a numpy array. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix, n_samples) baselines : ndarray or float Baseline estimates for each pixel that are subtracted from the waveforms before deconvolution. @@ -1459,8 +1521,8 @@ def deconvolve( Shape: (n_pix, upsampling * n_samples) """ deconvolved_waveforms = np.atleast_2d(waveforms) - np.atleast_2d(baselines).T - deconvolved_waveforms[:, 1:] -= pole_zero * deconvolved_waveforms[:, :-1] - deconvolved_waveforms[:, 0] = 0 + deconvolved_waveforms[..., 1:] -= pole_zero * deconvolved_waveforms[..., :-1] + deconvolved_waveforms[..., 0] = 0 if upsampling > 1: filt = np.ones(upsampling) @@ -1491,7 +1553,7 @@ def adaptive_centroid(waveforms, peak_index, rel_descend_limit, centroids): ---------- waveforms : ndarray Waveforms stored in a numpy array. - Shape: (n_pix, n_samples) + Shape: (n_channels, n_pix, n_samples) peak_index : ndarray or int Peak index for each pixel. rel_descend_limit : ndarray or float @@ -1644,6 +1706,13 @@ def clip(x, lo=0.0, hi=np.inf): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels ) -> DL1CameraContainer: + if waveforms.shape[-3] != 1: + raise AttributeError( + "The FlashCam concept foresees the use of preamplifiers with non-linear " + "gain characteristics, instead of splitting the signal into two separate " + f"channels but the data consists of {waveforms.shape[-3]} channels." + ) + upsampling = self.upsampling.tel[tel_id] integration_window_width = self.window_width.tel[tel_id] integration_window_shift = self.window_shift.tel[tel_id] @@ -1696,4 +1765,6 @@ def __call__( if shift != 0: peak_time -= shift - return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) + return DL1CameraContainer( + image=charge[0, :], peak_time=peak_time[0, :], is_valid=True + ) diff --git a/src/ctapipe/image/reducer.py b/src/ctapipe/image/reducer.py index c79f1602216..17bb5e77e46 100644 --- a/src/ctapipe/image/reducer.py +++ b/src/ctapipe/image/reducer.py @@ -187,8 +187,8 @@ def select_pixels(self, waveforms, tel_id=None, selected_gain_channel=None): dl1: DL1CameraContainer = extractor( waveforms, tel_id=tel_id, - selected_gain_channel=selected_gain_channel, broken_pixels=broken_pixels, + selected_gain_channel=selected_gain_channel, ) # 1) Step: TailcutCleaning at first diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 23c67abd5d5..28289169414 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -72,7 +72,7 @@ def subarray_1_LST(prod3_lst, reference_location): subarray = SubarrayDescription( "One LST", tel_positions={1: np.zeros(3) * u.m}, - tel_descriptions={1: prod3_lst}, + tel_descriptions={1: deepcopy(prod3_lst)}, reference_location=reference_location, ) return subarray @@ -83,7 +83,7 @@ def subarray_mst_fc(prod5_mst_flashcam, reference_location): subarray = SubarrayDescription( "One MST with FlashCam", tel_positions={1: np.zeros(3) * u.m}, - tel_descriptions={1: prod5_mst_flashcam}, + tel_descriptions={1: deepcopy(prod5_mst_flashcam)}, reference_location=reference_location, ) return subarray @@ -158,7 +158,7 @@ def toymodel_mst_fc(subarray_mst_fc: object) -> object: def test_extract_around_peak(toymodel): waveforms, _, _, _, _, _ = toymodel - n_pixels, n_samples = waveforms.shape + _, n_pixels, n_samples = waveforms.shape rand = np.random.RandomState(1) peak_index = rand.uniform(0, n_samples, n_pixels).astype(np.int64) charge, peak_time = extract_around_peak(waveforms, peak_index, 7, 3, 1) @@ -180,7 +180,7 @@ def test_extract_around_peak(toymodel): def test_extract_sliding_window(toymodel): waveforms, _, _, _, _, _ = toymodel - n_pixels, n_samples = waveforms.shape + n_samples = waveforms.shape[-1] charge, peak_time = extract_sliding_window(waveforms, 7, 1) assert (charge >= 0).all() assert (peak_time >= 0).all() and (peak_time <= n_samples).all() @@ -199,7 +199,7 @@ def test_extract_sliding_window(toymodel): def test_extract_around_peak_charge_expected(toymodel): - waveforms = np.ones((2048, 96)) + waveforms = np.ones((1, 2048, 96)) n_samples = waveforms.shape[-1] sampling_rate_ghz = 1 @@ -255,18 +255,18 @@ def test_extract_around_peak_charge_expected(toymodel): def test_neighbor_average_peakpos(toymodel): waveforms, subarray, tel_id, _, _, _ = toymodel neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) peak_pos = neighbor_average_maximum( waveforms, neighbors_indices=neighbors.indices, neighbors_indptr=neighbors.indptr, local_weight=0, broken_pixels=broken_pixels, - ) + )[0] pixel = 0 _, nei_pixel = np.where(neighbors[pixel].A) - expected_average = waveforms[nei_pixel].sum(0) / len(nei_pixel) + expected_average = waveforms[0][nei_pixel].sum(0) / len(nei_pixel) expected_peak_pos = np.argmax(expected_average, axis=-1) assert (peak_pos[pixel] == expected_peak_pos).all() @@ -277,12 +277,12 @@ def test_neighbor_average_peakpos(toymodel): neighbors_indptr=neighbors.indptr, local_weight=local_weight, broken_pixels=broken_pixels, - ) + )[0] pixel = 1 _, nei_pixel = np.where(neighbors[pixel].A) nei_pixel = np.concatenate([nei_pixel, [pixel] * local_weight]) - expected_average = waveforms[nei_pixel].sum(0) / len(nei_pixel) + expected_average = waveforms[0][nei_pixel].sum(0) / len(nei_pixel) expected_peak_pos = np.argmax(expected_average, axis=-1) assert (peak_pos[pixel] == expected_peak_pos).all() @@ -298,11 +298,11 @@ def test_extract_peak_time_within_range(): def test_baseline_subtractor(toymodel): waveforms, _, _, _, _, _ = toymodel - n_pixels, _ = waveforms.shape + n_pixels = waveforms.shape[-2] rand = np.random.RandomState(1) offset = np.arange(n_pixels)[:, np.newaxis] waveforms = rand.normal(0, 0.1, waveforms.shape) + offset - assert_allclose(waveforms[3].mean(), 3, rtol=1e-2) + assert_allclose(waveforms[0][3].mean(), 3, rtol=1e-2) baseline_subtracted = subtract_baseline(waveforms, 0, 10) assert_allclose(baseline_subtracted.mean(), 0, atol=1e-3) @@ -366,7 +366,7 @@ def test_extractors(Extractor, toymodel): true_time, ) = toymodel extractor = Extractor(subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) @@ -388,7 +388,7 @@ def test_integration_correction_off(Extractor, toymodel): true_time, ) = toymodel extractor = Extractor(subarray=subarray, apply_integration_correction=False) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid @@ -410,7 +410,7 @@ def test_fixed_window_sum(toymodel): true_time, ) = toymodel extractor = FixedWindowSum(subarray=subarray, peak_index=47) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) @@ -427,7 +427,7 @@ def test_sliding_window_max_sum(toymodel): true_time, ) = toymodel extractor = SlidingWindowMaxSum(subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) print(true_charge, dl1.image, true_time, dl1.peak_time) assert_allclose(dl1.image, true_charge, rtol=0.1) @@ -446,7 +446,7 @@ def test_neighbor_peak_window_sum_local_weight(toymodel): ) = toymodel extractor = NeighborPeakWindowSum(subarray=subarray, local_weight=4) assert extractor.local_weight.tel[tel_id] == 4 - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) @@ -523,7 +523,7 @@ def test_Two_pass_window_sum_no_noise(subarray_1_LST): # Test only the 1st pass extractor.disable_second_pass = True - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1_pass1 = extractor(waveforms, 1, selected_gain_channel, broken_pixels) assert_allclose( dl1_pass1.image[signal_pixels & integration_window_inside], @@ -577,7 +577,7 @@ def test_waveform_extractor_factory(toymodel): ) = toymodel extractor = ImageExtractor.from_name("LocalPeakWindowSum", subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) @@ -603,7 +603,7 @@ def test_waveform_extractor_factory_args(subarray): def test_extractor_tel_param(toymodel): waveforms, subarray, _, _, _, _ = toymodel - _, n_samples = waveforms.shape + n_samples = waveforms.shape[-1] config = Config( { @@ -615,7 +615,7 @@ def test_extractor_tel_param(toymodel): ) waveforms, subarray, _, _, _, _ = toymodel - n_pixels, n_samples = waveforms.shape + _, n_pixels, n_samples = waveforms.shape extractor = ImageExtractor.from_name( "FixedWindowSum", subarray=subarray, config=config ) @@ -634,9 +634,9 @@ def test_dtype(Extractor, subarray): n_pixels = subarray.tel[tel_id].camera.geometry.n_pixels selected_gain_channel = np.zeros(n_pixels, dtype=int) - waveforms = np.ones((n_pixels, 50), dtype="float64") + waveforms = np.ones((1, n_pixels, 50), dtype="float64") extractor = Extractor(subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.image.dtype == np.float32 assert dl1.peak_time.dtype == np.float32 @@ -655,13 +655,13 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray): bright_pixels[np.random.choice(n_pixels, size=int(0.1 * n_pixels))] = True # signal in dim pixels is in slice 10, signal in bright pixels is in slice 30 - waveforms = np.zeros((n_pixels, 50), dtype="float64") - waveforms[~bright_pixels, 9] = 3 - waveforms[~bright_pixels, 10] = 5 - waveforms[~bright_pixels, 11] = 2 - waveforms[bright_pixels, 29] = 5 - waveforms[bright_pixels, 30] = 10 - waveforms[bright_pixels, 31] = 3 + waveforms = np.zeros((1, n_pixels, 50), dtype="float64") + waveforms[0][~bright_pixels, 9] = 3 + waveforms[0][~bright_pixels, 10] = 5 + waveforms[0][~bright_pixels, 11] = 2 + waveforms[0][bright_pixels, 29] = 5 + waveforms[0][bright_pixels, 30] = 10 + waveforms[0][bright_pixels, 31] = 3 extractor = GlobalPeakWindowSum( subarray=subarray, @@ -671,7 +671,7 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray): apply_integration_correction=False, ) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert np.allclose(dl1.image[bright_pixels], 18) @@ -692,7 +692,7 @@ def test_adaptive_centroid(toymodel_mst_fc): ) = toymodel_mst_fc neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) trig_time = np.argmax(waveforms, axis=-1) peak_time = adaptive_centroid( @@ -703,7 +703,7 @@ def test_adaptive_centroid(toymodel_mst_fc): assert (peak_time == trig_time).all() - waveforms = waveforms[np.min(waveforms, axis=-1) > 0.0] + waveforms = waveforms[np.min(waveforms, axis=-1) > 0.0].reshape((0, 0, 0)) peak_pos = neighbor_average_maximum( waveforms, neighbors_indices=neighbors.indices, @@ -732,11 +732,11 @@ def test_deconvolve(toymodel_mst_fc): deconvolved_waveforms_0 = deconvolve(waveforms, 0, 0, 0.0) - assert (deconvolved_waveforms_0[:, 1:] == waveforms[:, 1:]).all() + assert (deconvolved_waveforms_0[..., 1:] == waveforms[..., 1:]).all() deconvolved_waveforms_1 = deconvolve(waveforms, 0, 0, 1.0) - assert (deconvolved_waveforms_1[:, 1:] == np.diff(waveforms, axis=-1)).all() + assert (deconvolved_waveforms_1[..., 1:] == np.diff(waveforms, axis=-1)).all() def test_upsampling(toymodel_mst_fc): @@ -796,7 +796,7 @@ def test_FC_time(toymodel_mst_fc_time): ) = toymodel_mst_fc_time extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid @@ -804,7 +804,7 @@ def test_FC_time(toymodel_mst_fc_time): extractor = FlashCamExtractor( subarray=subarray, upsampling=1, leading_edge_timing=True ) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid @@ -812,13 +812,13 @@ def test_FC_time(toymodel_mst_fc_time): extractor = FlashCamExtractor( subarray=subarray, upsampling=1, leading_edge_timing=False ) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=False) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid @@ -835,7 +835,7 @@ def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): true_time, ) = toymodel_mst_fc extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) - broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) + broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert dl1.is_valid @@ -853,7 +853,7 @@ def is_flashcam(tel_id): true_charge = event.simulation.tel[tel_id].true_image waveforms = event.r1.tel[tel_id].waveform - selected_gain_channel = np.zeros(waveforms.shape[0], dtype=np.int64) + selected_gain_channel = np.zeros(waveforms.shape[-2], dtype=np.int64) broken_pixels = event.mon.tel[ tel_id ].pixel_status.hardware_failing_pixels[0] diff --git a/src/ctapipe/image/tests/test_reducer.py b/src/ctapipe/image/tests/test_reducer.py index ca6e8fcaba6..d70709403b2 100644 --- a/src/ctapipe/image/tests/test_reducer.py +++ b/src/ctapipe/image/tests/test_reducer.py @@ -31,7 +31,7 @@ def subarray_lst(prod3_lst, reference_location): def test_null_data_volume_reducer(subarray_lst): subarray, _, _, _, _ = subarray_lst rng = np.random.default_rng(0) - waveforms = rng.uniform(0, 1, (2048, 96)) + waveforms = rng.uniform(0, 1, (1, 2048, 96)) reducer = NullDataVolumeReducer(subarray=subarray) reduced_waveforms_mask = reducer(waveforms) reduced_waveforms = waveforms.copy() @@ -43,22 +43,24 @@ def test_tailcuts_data_volume_reducer(subarray_lst): subarray, tel_id, selected_gain_channel, n_pixels, n_samples = subarray_lst # create signal - waveforms_signal = np.zeros((n_pixels, n_samples), dtype=np.float64) + waveforms_signal = np.zeros((1, n_pixels, n_samples), dtype=np.float64) # Should be selected as core-pixel from Step 1) tailcuts_clean - waveforms_signal[9] = 100 + waveforms_signal[0][9] = 100 # 10 and 8 as boundary-pixel from Step 1) tailcuts_clean # 6 and 5 as iteration-pixel in Step 2) - waveforms_signal[[10, 8, 6, 5]] = 50 + waveforms_signal[0][[10, 8, 6, 5]] = 50 # pixels from dilate at the end in Step 3) - waveforms_signal[[0, 1, 4, 7, 11, 13, 121, 122, 136, 137, 257, 258, 267, 272]] = 25 + waveforms_signal[0][ + [0, 1, 4, 7, 11, 13, 121, 122, 136, 137, 257, 258, 267, 272] + ] = 25 expected_waveforms = waveforms_signal.copy() # add some random pixels, which should not be selected - waveforms_signal[[50, 51, 135, 138, 54, 170, 210, 400]] = 25 + waveforms_signal[0][[50, 51, 135, 138, 54, 170, 210, 400]] = 25 # Reduction parameters reduction_param = Config( diff --git a/src/ctapipe/image/toymodel.py b/src/ctapipe/image/toymodel.py index 7c1575e75c4..da27c33e9aa 100644 --- a/src/ctapipe/image/toymodel.py +++ b/src/ctapipe/image/toymodel.py @@ -135,7 +135,7 @@ def get_waveform(self, charge, time, n_samples): ------- waveform : ndarray Toy model waveform - Shape (n_pixels, n_samples) + Shape (n_channels, n_pixels, n_samples) """ n_pixels = charge.size @@ -156,7 +156,7 @@ def get_waveform(self, charge, time, n_samples): ).sum(-1) * self.ref_width_ns # Waveform units: p.e. ) - return sampled + return sampled[np.newaxis, :] @classmethod def from_camera_readout(cls, readout, gain_channel=0): diff --git a/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py index 4cda9ae6b59..a1b89ecbe36 100644 --- a/src/ctapipe/io/simteleventsource.py +++ b/src/ctapipe/io/simteleventsource.py @@ -272,7 +272,7 @@ def apply_simtel_r1_calibration( ------- r1_waveforms : ndarray Calibrated waveforms - Shape: (n_pixels, n_samples) + Shape: (n_channels, n_pixels, n_samples) selected_gain_channel : ndarray The gain channel selected for each pixel Shape: (n_pixels) @@ -286,10 +286,12 @@ def apply_simtel_r1_calibration( if n_channels == 1: selected_gain_channel = np.zeros(n_pixels, dtype=np.int8) - r1_waveforms = r1_waveforms[0] + r1_waveforms = r1_waveforms else: selected_gain_channel = gain_selector(r0_waveforms) - r1_waveforms = r1_waveforms[selected_gain_channel, np.arange(n_pixels)] + r1_waveforms = r1_waveforms[ + np.newaxis, selected_gain_channel, np.arange(n_pixels) + ] return r1_waveforms, selected_gain_channel From fd7fc3a070b7118b60f5efff1a2d5a6ab4eaaf0e Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 20 Mar 2024 19:00:57 +0100 Subject: [PATCH 02/42] Adapting CameraCalibrator - Fixing some unit tests --- src/ctapipe/calib/camera/calibrator.py | 14 ++++----- .../calib/camera/tests/test_calibrator.py | 18 ++++++----- src/ctapipe/image/extractor.py | 31 +++++-------------- src/ctapipe/image/tests/test_extractor.py | 2 +- src/ctapipe/image/tests/test_reducer.py | 7 +++-- .../io/tests/test_simteleventsource.py | 17 +++++----- 6 files changed, 38 insertions(+), 51 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index c7864ca2e23..23328458f96 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -188,7 +188,7 @@ def _calibrate_dl0(self, event, tel_id): ) dl0_waveform = r1.waveform.copy() - dl0_waveform[~signal_pixels] = 0 + dl0_waveform[:, ~signal_pixels] = 0 dl0_pixel_status = r1.pixel_status.copy() # set dvr pixel bit in pixel_status for pixels kept by DVR @@ -211,7 +211,7 @@ def _calibrate_dl1(self, event, tel_id): if self._check_dl0_empty(waveforms): return - n_pixels, n_samples = waveforms.shape + _, n_pixels, n_samples = waveforms.shape selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel broken_pixels = _get_invalid_pixels( @@ -227,8 +227,8 @@ def _calibrate_dl1(self, event, tel_id): # subtract any remaining pedestal before extraction if dl1_calib.pedestal_offset is not None: # this copies intentionally, we don't want to modify the dl0 data - # waveforms have shape (n_pixel, n_samples), pedestals (n_pixels, ) - waveforms = waveforms - dl1_calib.pedestal_offset[:, np.newaxis] + # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels, ) + waveforms = waveforms - dl1_calib.pedestal_offset[np.newaxis, :, np.newaxis] if n_samples == 1: # To handle ASTRI and dst @@ -238,7 +238,7 @@ def _calibrate_dl1(self, event, tel_id): # - Don't do anything if dl1 container already filled # - Update on SST review decision dl1 = DL1CameraContainer( - image=waveforms[..., 0].astype(np.float32), + image=waveforms[0, ..., 0].astype(np.float32), peak_time=np.zeros(n_pixels, dtype=np.float32), is_valid=True, ) @@ -309,7 +309,7 @@ def shift_waveforms(waveforms, time_shift_samples): Parameters ---------- - waveforms: ndarray of shape (n_pixels, n_samples) + waveforms: ndarray of shape (n_channels, n_pixels, n_samples) The waveforms to shift time_shift_samples: ndarray of shape (n_pixels, ) The shift to apply in units of samples. @@ -318,7 +318,7 @@ def shift_waveforms(waveforms, time_shift_samples): Returns ------- - shifted_waveforms: ndarray of shape (n_pixels, n_samples) + shifted_waveforms: ndarray of shape (n_channels, n_pixels, n_samples) The shifted waveforms remaining_shift: ndarray of shape (n_pixels, ) The remaining shift after applying the integer shift to the waveforms. diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index 1e22d4c28c9..27bebb4f02a 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -109,7 +109,7 @@ def test_check_r1_empty(example_event, example_subarray): image_extractor=FullWaveformSum(subarray=example_subarray), ) event = ArrayEventContainer() - event.dl0.tel[tel_id].waveform = np.full((2048, 128), 2) + event.dl0.tel[tel_id].waveform = np.full((1, 2048, 128), 2) calibrator(event) assert (event.dl0.tel[tel_id].waveform == 2).all() assert (event.dl1.tel[tel_id].image == 2 * 128).all() @@ -148,22 +148,24 @@ def test_dl1_charge_calib(example_subarray): # Randomize times and create pulses time_offset = random.uniform(-10, +10, n_pixels) - y = norm.pdf(x, mid + time_offset[:, np.newaxis], pulse_sigma).astype("float32") + y = norm.pdf(x, mid + time_offset[np.newaxis, :, np.newaxis], pulse_sigma).astype( + "float32" + ) camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma)[np.newaxis, :] camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate # Define absolute calibration coefficients absolute = random.uniform(100, 1000, n_pixels).astype("float32") - y *= absolute[:, np.newaxis] + y *= absolute[np.newaxis, :, np.newaxis] # Define relative coefficients relative = random.normal(1, 0.01, n_pixels) - y /= relative[:, np.newaxis] + y /= relative[np.newaxis, :, np.newaxis] # Define pedestal pedestal = random.uniform(-4, 4, n_pixels) - y += pedestal[:, np.newaxis] + y += pedestal[np.newaxis, :, np.newaxis] event = ArrayEventContainer() tel_id = list(subarray.tel.keys())[0] @@ -176,7 +178,7 @@ def test_dl1_charge_calib(example_subarray): subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) ) calibrator(event) - np.testing.assert_allclose(event.dl1.tel[tel_id].image, y.sum(1), rtol=1e-4) + np.testing.assert_allclose(event.dl1.tel[tel_id].image, y.sum(2)[0, :], rtol=1e-4) event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal event.calibration.tel[tel_id].dl1.absolute_factor = absolute @@ -276,8 +278,8 @@ def test_invalid_pixels(example_event, example_subarray): event.mon.tel[tel_id].pixel_status.flatfield_failing_pixels[:, 0] = True event.r1.tel[tel_id].waveform.fill(0.0) - event.r1.tel[tel_id].waveform[1:, 20] = 1.0 - event.r1.tel[tel_id].waveform[0, 10] = 9999 + event.r1.tel[tel_id].waveform[:, 1:, 20] = 1.0 + event.r1.tel[tel_id].waveform[:, 0, 10] = 9999 calibrator = CameraCalibrator( subarray=example_subarray, diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index caa8772cd31..998d9cd4852 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -503,9 +503,7 @@ def __call__( selected_gain_channel ] else: - charge *= self._calculate_correction(tel_id=tel_id)[ - :, np.newaxis, np.newaxis - ] + charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -603,9 +601,7 @@ def __call__( selected_gain_channel ] else: - charge *= self._calculate_correction(tel_id=tel_id)[ - :, np.newaxis, np.newaxis - ] + charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -677,9 +673,7 @@ def __call__( selected_gain_channel ] else: - charge *= self._calculate_correction(tel_id=tel_id)[ - :, np.newaxis, np.newaxis - ] + charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -762,9 +756,7 @@ def __call__( selected_gain_channel ] else: - charge *= self._calculate_correction(tel_id=tel_id)[ - :, np.newaxis, np.newaxis - ] + charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -849,9 +841,7 @@ def __call__( selected_gain_channel ] else: - charge *= self._calculate_correction(tel_id=tel_id)[ - :, np.newaxis, np.newaxis - ] + charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -1518,7 +1508,7 @@ def deconvolve( ------- deconvolved_waveforms : ndarray Deconvolved and upsampled waveforms stored in a numpy array. - Shape: (n_pix, upsampling * n_samples) + Shape: (n_channels, n_pix, upsampling * n_samples) """ deconvolved_waveforms = np.atleast_2d(waveforms) - np.atleast_2d(baselines).T deconvolved_waveforms[..., 1:] -= pole_zero * deconvolved_waveforms[..., :-1] @@ -1706,13 +1696,6 @@ def clip(x, lo=0.0, hi=np.inf): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels ) -> DL1CameraContainer: - if waveforms.shape[-3] != 1: - raise AttributeError( - "The FlashCam concept foresees the use of preamplifiers with non-linear " - "gain characteristics, instead of splitting the signal into two separate " - f"channels but the data consists of {waveforms.shape[-3]} channels." - ) - upsampling = self.upsampling.tel[tel_id] integration_window_width = self.window_width.tel[tel_id] integration_window_shift = self.window_shift.tel[tel_id] @@ -1766,5 +1749,5 @@ def __call__( peak_time -= shift return DL1CameraContainer( - image=charge[0, :], peak_time=peak_time[0, :], is_valid=True + image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 28289169414..603be1ec423 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -198,7 +198,7 @@ def test_extract_sliding_window(toymodel): assert charge.dtype == np.float32 -def test_extract_around_peak_charge_expected(toymodel): +def test_extract_around_peak_charge_expected(): waveforms = np.ones((1, 2048, 96)) n_samples = waveforms.shape[-1] sampling_rate_ghz = 1 diff --git a/src/ctapipe/image/tests/test_reducer.py b/src/ctapipe/image/tests/test_reducer.py index d70709403b2..49cd06776be 100644 --- a/src/ctapipe/image/tests/test_reducer.py +++ b/src/ctapipe/image/tests/test_reducer.py @@ -35,7 +35,7 @@ def test_null_data_volume_reducer(subarray_lst): reducer = NullDataVolumeReducer(subarray=subarray) reduced_waveforms_mask = reducer(waveforms) reduced_waveforms = waveforms.copy() - reduced_waveforms[~reduced_waveforms_mask] = 0 + reduced_waveforms[:, ~reduced_waveforms_mask] = 0 assert_array_equal(waveforms, reduced_waveforms) @@ -43,7 +43,8 @@ def test_tailcuts_data_volume_reducer(subarray_lst): subarray, tel_id, selected_gain_channel, n_pixels, n_samples = subarray_lst # create signal - waveforms_signal = np.zeros((1, n_pixels, n_samples), dtype=np.float64) + n_channels = 1 + waveforms_signal = np.zeros((n_channels, n_pixels, n_samples), dtype=np.float64) # Should be selected as core-pixel from Step 1) tailcuts_clean waveforms_signal[0][9] = 100 @@ -87,7 +88,7 @@ def test_tailcuts_data_volume_reducer(subarray_lst): reduced_waveforms_mask = reducer( waveforms_signal, tel_id=tel_id, selected_gain_channel=selected_gain_channel ) - reduced_waveforms[~reduced_waveforms_mask] = 0 + reduced_waveforms[:, ~reduced_waveforms_mask] = 0 assert (reduced_waveforms != 0).sum() == (1 + 4 + 14) * n_samples assert_array_equal(expected_waveforms, reduced_waveforms) diff --git a/src/ctapipe/io/tests/test_simteleventsource.py b/src/ctapipe/io/tests/test_simteleventsource.py index 43c146c837c..c614bab2145 100644 --- a/src/ctapipe/io/tests/test_simteleventsource.py +++ b/src/ctapipe/io/tests/test_simteleventsource.py @@ -1,4 +1,5 @@ """ tests of SimTelEventSource """ + # pylint: disable=import-outside-toplevel import copy from itertools import zip_longest @@ -285,12 +286,12 @@ def test_apply_simtel_r1_calibration_1_channel(): ) assert (selected_gain_channel == 0).all() - assert r1_waveforms.ndim == 2 - assert r1_waveforms.shape == (n_pixels, n_samples) + assert r1_waveforms.ndim == 3 + assert r1_waveforms.shape == (1, n_pixels, n_samples) ped = pedestal - assert r1_waveforms[0, 0] == (r0_waveforms[0, 0, 0] - ped[0, 0]) * dc_to_pe[0, 0] - assert r1_waveforms[1, 0] == (r0_waveforms[0, 1, 0] - ped[0, 1]) * dc_to_pe[0, 1] + assert r1_waveforms[0, 0, 0] == (r0_waveforms[0, 0, 0] - ped[0, 0]) * dc_to_pe[0, 0] + assert r1_waveforms[0, 1, 0] == (r0_waveforms[0, 1, 0] - ped[0, 1]) * dc_to_pe[0, 1] def test_apply_simtel_r1_calibration_2_channel(): @@ -317,12 +318,12 @@ def test_apply_simtel_r1_calibration_2_channel(): assert selected_gain_channel[0] == 1 assert (selected_gain_channel[np.arange(1, 2048)] == 0).all() - assert r1_waveforms.ndim == 2 - assert r1_waveforms.shape == (n_pixels, n_samples) + assert r1_waveforms.ndim == 3 + assert r1_waveforms.shape == (1, n_pixels, n_samples) ped = pedestal - assert r1_waveforms[0, 0] == (r0_waveforms[1, 0, 0] - ped[1, 0]) * dc_to_pe[1, 0] - assert r1_waveforms[1, 0] == (r0_waveforms[0, 1, 0] - ped[0, 1]) * dc_to_pe[0, 1] + assert r1_waveforms[0, 0, 0] == (r0_waveforms[1, 0, 0] - ped[1, 0]) * dc_to_pe[1, 0] + assert r1_waveforms[0, 1, 0] == (r0_waveforms[0, 1, 0] - ped[0, 1]) * dc_to_pe[0, 1] def test_focal_length_choice(): From 2fb7f06dd87e8c019cbfa9c34dee4047a4e7f216 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 22 Mar 2024 10:26:27 +0100 Subject: [PATCH 03/42] Enhance Toymodel - Adapt toymodel to the current needs - Add possibility to create 2 gain toymodel --- src/ctapipe/image/extractor.py | 2 +- src/ctapipe/image/tests/test_extractor.py | 85 ++++++++++++++++------- src/ctapipe/image/toymodel.py | 62 ++++++++++++----- 3 files changed, 104 insertions(+), 45 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 998d9cd4852..2a46736d268 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -90,7 +90,7 @@ def extract_around_peak( width : ndarray or int Window size of integration window for each pixel. shift : ndarray or int - Window size of integration window for each pixel. + Shift of the integration window from the peak_index. sampling_rate_ghz : float Sampling rate of the camera, in units of GHz Astropy units should have to_value('GHz') applied before being passed diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 603be1ec423..a85d4eed326 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -37,6 +37,8 @@ extractors.remove(FixedWindowSum) extractors.remove(FlashCamExtractor) +camera_toymodels = ["toymodel", "toymodel_lst", "toymodel_mst_fc"] + @pytest.fixture(scope="module") def subarray(prod5_sst, reference_location): @@ -92,8 +94,8 @@ def subarray_mst_fc(prod5_mst_flashcam, reference_location): def get_test_toymodel(subarray, minCharge=100, maxCharge=1000): tel_id = list(subarray.tel.keys())[0] n_pixels = subarray.tel[tel_id].camera.geometry.n_pixels - n_samples = 96 readout = subarray.tel[tel_id].camera.readout + n_samples = readout.n_samples random = np.random.RandomState(1) charge = random.uniform(minCharge, maxCharge, n_pixels) @@ -109,10 +111,20 @@ def get_test_toymodel(subarray, minCharge=100, maxCharge=1000): @pytest.fixture(scope="module") -def toymodel(subarray): +def toymodel(subarray: object) -> object: return get_test_toymodel(subarray) +@pytest.fixture(scope="module") +def toymodel_lst(subarray_1_LST: object) -> object: + return get_test_toymodel(subarray_1_LST) + + +@pytest.fixture(scope="module") +def toymodel_mst_fc(subarray_mst_fc: object) -> object: + return get_test_toymodel(subarray_mst_fc) + + def get_test_toymodel_gradient(subarray, minCharge=100, maxCharge=1000): tel_id = list(subarray.tel.keys())[0] n_pixels = subarray.tel[tel_id].camera.geometry.n_pixels @@ -129,7 +141,7 @@ def get_test_toymodel_gradient(subarray, minCharge=100, maxCharge=1000): mask = geometry.pix_id == pix_id dilated_mask = mask.copy() - for row in range(4): + for _ in range(4): dilated_mask = dilate(geometry, dilated_mask) x_d = subarray.tel[tel_id].camera.geometry.pix_x.value @@ -151,11 +163,6 @@ def toymodel_mst_fc_time(subarray_mst_fc: object) -> object: return get_test_toymodel_gradient(subarray_mst_fc) -@pytest.fixture(scope="module") -def toymodel_mst_fc(subarray_mst_fc: object) -> object: - return get_test_toymodel(subarray_mst_fc) - - def test_extract_around_peak(toymodel): waveforms, _, _, _, _, _ = toymodel _, n_pixels, n_samples = waveforms.shape @@ -252,8 +259,9 @@ def test_extract_around_peak_charge_expected(): assert_equal(charge, n_samples) -def test_neighbor_average_peakpos(toymodel): - waveforms, subarray, tel_id, _, _, _ = toymodel +@pytest.mark.parametrize("toymodels", camera_toymodels) +def test_neighbor_average_peakpos(toymodels, request): + waveforms, subarray, tel_id, _, _, _ = request.getfixturevalue(toymodels) neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) peak_pos = neighbor_average_maximum( @@ -262,13 +270,14 @@ def test_neighbor_average_peakpos(toymodel): neighbors_indptr=neighbors.indptr, local_weight=0, broken_pixels=broken_pixels, - )[0] + ) pixel = 0 _, nei_pixel = np.where(neighbors[pixel].A) - expected_average = waveforms[0][nei_pixel].sum(0) / len(nei_pixel) + expected_average = waveforms[:, nei_pixel].sum(1) / len(nei_pixel) expected_peak_pos = np.argmax(expected_average, axis=-1) - assert (peak_pos[pixel] == expected_peak_pos).all() + for ichannel in range(waveforms.shape[-3]): + assert (peak_pos[ichannel][pixel] == expected_peak_pos).all() local_weight = 4 peak_pos = neighbor_average_maximum( @@ -277,14 +286,15 @@ def test_neighbor_average_peakpos(toymodel): neighbors_indptr=neighbors.indptr, local_weight=local_weight, broken_pixels=broken_pixels, - )[0] + ) pixel = 1 _, nei_pixel = np.where(neighbors[pixel].A) nei_pixel = np.concatenate([nei_pixel, [pixel] * local_weight]) - expected_average = waveforms[0][nei_pixel].sum(0) / len(nei_pixel) + expected_average = waveforms[:, nei_pixel].sum(1) / len(nei_pixel) expected_peak_pos = np.argmax(expected_average, axis=-1) - assert (peak_pos[pixel] == expected_peak_pos).all() + for ichannel in range(waveforms.shape[-3]): + assert (peak_pos[ichannel][pixel] == expected_peak_pos).all() def test_extract_peak_time_within_range(): @@ -355,8 +365,9 @@ def test_integration_correction_outofbounds(subarray): np.testing.assert_allclose(full_integral, window_integral * correction) +@pytest.mark.parametrize("toymodels", camera_toymodels) @pytest.mark.parametrize("Extractor", extractors) -def test_extractors(Extractor, toymodel): +def test_extractors(Extractor, toymodels, request): ( waveforms, subarray, @@ -364,17 +375,30 @@ def test_extractors(Extractor, toymodel): selected_gain_channel, true_charge, true_time, - ) = toymodel + ) = request.getfixturevalue(toymodels) extractor = Extractor(subarray=subarray) broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + + if Extractor is TwoPassWindowSum and waveforms.shape[-3] != 1: + with pytest.raises(AttributeError): + extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) + return + dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) - assert_allclose(dl1.image, true_charge, rtol=0.1) - assert_allclose(dl1.peak_time, true_time, rtol=0.1) assert dl1.is_valid + breakpoint() + if dl1.image.ndim == 1: + assert_allclose(dl1.image, true_charge, rtol=0.2) + assert_allclose(dl1.peak_time, true_time, rtol=0.1) + else: + for ichannel in range(dl1.image.shape[-2]): + assert_allclose(dl1.image[ichannel], true_charge, rtol=0.2) + assert_allclose(dl1.peak_time[ichannel], true_time, rtol=0.1) +@pytest.mark.parametrize("toymodels", camera_toymodels) @pytest.mark.parametrize("Extractor", extractors) -def test_integration_correction_off(Extractor, toymodel): +def test_integration_correction_off(Extractor, toymodels, request): # full waveform extractor does not have an integration correction if Extractor is FullWaveformSum: return @@ -386,18 +410,27 @@ def test_integration_correction_off(Extractor, toymodel): selected_gain_channel, true_charge, true_time, - ) = toymodel + ) = request.getfixturevalue(toymodels) extractor = Extractor(subarray=subarray, apply_integration_correction=False) broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) - dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) + if Extractor is TwoPassWindowSum and waveforms.shape[-3] != 1: + with pytest.raises(AttributeError): + extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) + return + + dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid # peak time should stay the same - assert_allclose(dl1.peak_time, true_time, rtol=0.1) - # charge should be too small without correction - assert np.all(dl1.image <= true_charge) + if dl1.image.ndim == 1: + assert np.all(dl1.image <= true_charge) + assert_allclose(dl1.peak_time, true_time, rtol=0.1) + else: + for ichannel in range(dl1.image.shape[-2]): + assert np.all(dl1.image[ichannel] <= true_charge) + assert_allclose(dl1.peak_time[ichannel], true_time, rtol=0.1) def test_fixed_window_sum(toymodel): diff --git a/src/ctapipe/image/toymodel.py b/src/ctapipe/image/toymodel.py index da27c33e9aa..ba9cc167600 100644 --- a/src/ctapipe/image/toymodel.py +++ b/src/ctapipe/image/toymodel.py @@ -18,6 +18,7 @@ >>> print(image.shape) (400,) """ + from abc import ABCMeta, abstractmethod import astropy.units as u @@ -26,6 +27,7 @@ from scipy.ndimage import convolve1d from scipy.stats import multivariate_normal, norm, skewnorm +from ctapipe.calib.camera.gainselection import GainChannel from ctapipe.image.hillas import camera_to_shower_coordinates from ctapipe.utils import linalg @@ -39,6 +41,7 @@ TOYMODEL_RNG = default_rng(0) +VALID_GAIN_CHANNEL = {"HIGH", "LOW", "ALL"} def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_intercept): @@ -87,8 +90,8 @@ def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_int class WaveformModel: @u.quantity_input(reference_pulse_sample_width=u.ns, sample_width=u.ns) def __init__(self, reference_pulse, reference_pulse_sample_width, sample_width): - """Generate a toy model waveform using the reference pulse shape of a - camera. Useful for testing image extraction algorithms. + """Generate a toy model waveform for each gain channel using the reference + pulse shape of a camera. Useful for testing image extraction algorithms. Does not include the electronic noise and the Excess Noise Factor of the photosensor, and therefore should not be used to make charge @@ -99,23 +102,28 @@ def __init__(self, reference_pulse, reference_pulse_sample_width, sample_width): reference_pulse_sample_width : u.Quantity[time] Sample width of the reference pulse shape reference_pulse : ndarray - Reference pulse shape + Reference pulse shape for each channel sample_width : u.Quantity[time] Sample width of the waveform """ + self.n_channels = reference_pulse.shape[-2] self.upsampling = 10 reference_pulse_sample_width = reference_pulse_sample_width.to_value(u.ns) sample_width_ns = sample_width.to_value(u.ns) - ref_max_sample = reference_pulse.size * reference_pulse_sample_width + ref_max_sample = reference_pulse[0].size * reference_pulse_sample_width reference_pulse_x = np.arange(0, ref_max_sample, reference_pulse_sample_width) self.ref_width_ns = sample_width_ns / self.upsampling self.ref_interp_x = np.arange(0, reference_pulse_x.max(), self.ref_width_ns) - self.ref_interp_y = np.interp( - self.ref_interp_x, reference_pulse_x, reference_pulse - ) - self.ref_interp_y /= self.ref_interp_y.sum() * self.ref_width_ns - self.origin = self.ref_interp_y.argmax() - self.ref_interp_y.size // 2 + self.ref_interp_y = np.zeros((self.n_channels, self.ref_interp_x.size)) + for channel in range(self.n_channels): + self.ref_interp_y[channel] = np.interp( + self.ref_interp_x, reference_pulse_x, reference_pulse[channel] + ) + self.ref_interp_y = ( + self.ref_interp_y.T / (self.ref_interp_y.sum(-1) * self.ref_width_ns) + ).T + self.origin = self.ref_interp_y.argmax(-1) - self.ref_interp_y[0].size // 2 def get_waveform(self, charge, time, n_samples): """Obtain the waveform toy model. @@ -147,34 +155,52 @@ def get_waveform(self, charge, time, n_samples): sample[outofrange] = 0 charge[outofrange] = 0 readout[np.arange(n_pixels), sample] = charge - convolved = convolve1d( - readout, self.ref_interp_y, mode="constant", origin=self.origin - ) + convolved = np.zeros((self.n_channels, n_pixels, n_upsampled_samples)) + for channel in range(self.n_channels): + convolved[channel] = convolve1d( + readout, + self.ref_interp_y[channel], + mode="constant", + origin=self.origin[channel], + ) sampled = ( convolved.reshape( - (n_pixels, convolved.shape[-1] // self.upsampling, self.upsampling) + ( + self.n_channels, + n_pixels, + convolved.shape[-1] // self.upsampling, + self.upsampling, + ) ).sum(-1) * self.ref_width_ns # Waveform units: p.e. ) - return sampled[np.newaxis, :] + return sampled @classmethod - def from_camera_readout(cls, readout, gain_channel=0): + def from_camera_readout(cls, readout, gain_channel="ALL"): """Create class from a `ctapipe.instrument.CameraReadout`. Parameters ---------- readout : `ctapipe.instrument.CameraReadout` - gain_channel : int - The reference pulse gain channel to use + gain_channel : str + The reference pulse gain channel to use. + Choose between `HIGH`, `LOW` and `ALL`. Returns ------- WaveformModel """ + if gain_channel not in VALID_GAIN_CHANNEL: + raise ValueError(f"gain_channel must be one of {VALID_GAIN_CHANNEL}") + + ref_pulse_shape = readout.reference_pulse_shape + if gain_channel != "ALL": + ref_pulse_shape = ref_pulse_shape[np.newaxis, GainChannel[gain_channel]] + return cls( - readout.reference_pulse_shape[gain_channel], + ref_pulse_shape, readout.reference_pulse_sample_width, (1 / readout.sampling_rate).to(u.ns), ) From f547cfb5ff442b9f616e132fd7c8ad289e9a5d14 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 26 Mar 2024 13:02:36 +0100 Subject: [PATCH 04/42] Fix tests --- src/ctapipe/containers.py | 12 +-- src/ctapipe/image/extractor.py | 68 +++++------- src/ctapipe/image/tests/test_extractor.py | 101 ++++++++++++------ .../tests/test_sliding_window_correction.py | 62 ----------- src/ctapipe/image/tests/test_toy.py | 20 ++-- 5 files changed, 116 insertions(+), 147 deletions(-) delete mode 100644 src/ctapipe/image/tests/test_sliding_window_correction.py diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 19045a971e1..11a17107a71 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -484,16 +484,16 @@ class DL1CameraContainer(Container): image = Field( None, - "Numpy array of camera image, after waveform extraction." "Shape: (n_pixel)", - dtype=np.float32, - ndim=1, + "Numpy array of camera image, after waveform extraction." + "Shape: (n_pixel) if n_channels is 1 or data is gain selected" + "else: (n_channels, n_pixel)", ) peak_time = Field( None, "Numpy array containing position of the peak of the pulse as determined by " - "the extractor. Shape: (n_pixel, )", - dtype=np.float32, - ndim=1, + "the extractor." + "Shape: (n_pixel) if n_channels is 1 or data is gain selected" + "else: (n_channels, n_pixel)", ) image_mask = Field( None, diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 2a46736d268..c58ce61928d 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -251,18 +251,17 @@ def neighbor_average_maximum( # so the value of the pixel itself is already taken into account peak_pos = np.empty((n_channels, n_pixels), dtype=np.int64) - for pixel in prange(n_pixels): - average = waveforms[:, pixel] * local_weight - neighbors = indices[indptr[pixel] : indptr[pixel + 1]] + for ichannel in prange(n_channels): + for pixel in prange(n_pixels): + average = waveforms[ichannel, pixel] * local_weight + neighbors = indices[indptr[pixel] : indptr[pixel + 1]] - for neighbor in neighbors: - if broken_pixels[neighbor]: - continue + for neighbor in neighbors: + if broken_pixels[neighbor]: + continue + average += waveforms[ichannel][neighbor] - for ichannel in range(n_channels): - average[ichannel] += waveforms[ichannel][neighbor] - - peak_pos[:, pixel] = np.argmax(average, axis=-1) + peak_pos[ichannel, pixel] = np.argmax(average) return peak_pos @@ -498,12 +497,11 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - if selected_gain_channel is not None and waveforms.shape[-3] == 1: - charge *= self._calculate_correction(tel_id=tel_id)[ - selected_gain_channel - ] + correction = self._calculate_correction(tel_id=tel_id) + if selected_gain_channel is None: + charge *= correction[:, np.newaxis] else: - charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] + charge *= correction[selected_gain_channel] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -596,12 +594,11 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - if selected_gain_channel is not None and waveforms.shape[-3] == 1: - charge *= self._calculate_correction(tel_id=tel_id)[ - selected_gain_channel - ] + correction = self._calculate_correction(tel_id=tel_id) + if selected_gain_channel is None: + charge *= correction[:, np.newaxis] else: - charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] + charge *= correction[selected_gain_channel] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -668,12 +665,11 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - if selected_gain_channel is not None and waveforms.shape[-3] == 1: - charge *= self._calculate_correction(tel_id=tel_id)[ - selected_gain_channel - ] + correction = self._calculate_correction(tel_id=tel_id) + if selected_gain_channel is None: + charge *= correction[:, np.newaxis] else: - charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] + charge *= correction[selected_gain_channel] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -751,12 +747,11 @@ def __call__( waveforms, self.window_width.tel[tel_id], self.sampling_rate_ghz[tel_id] ) if self.apply_integration_correction.tel[tel_id]: - if selected_gain_channel is not None and waveforms.shape[-3] == 1: - charge *= self._calculate_correction(tel_id=tel_id)[ - selected_gain_channel - ] + correction = self._calculate_correction(tel_id=tel_id) + if selected_gain_channel is None: + charge *= correction[:, np.newaxis] else: - charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] + charge *= correction[selected_gain_channel] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -836,12 +831,11 @@ def __call__( self.sampling_rate_ghz[tel_id], ) if self.apply_integration_correction.tel[tel_id]: - if selected_gain_channel is not None and waveforms.shape[-3] == 1: - charge *= self._calculate_correction(tel_id=tel_id)[ - selected_gain_channel - ] + correction = self._calculate_correction(tel_id=tel_id) + if selected_gain_channel is None: + charge *= correction[:, np.newaxis] else: - charge *= self._calculate_correction(tel_id=tel_id)[:, np.newaxis] + charge *= correction[selected_gain_channel] return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -1630,10 +1624,6 @@ class FlashCamExtractor(ImageExtractor): "(peak_index - shift)", ).tag(config=True) - apply_integration_correction = BoolTelescopeParameter( - default_value=True, help="Apply the integration window correction" - ).tag(config=True) - local_weight = IntTelescopeParameter( default_value=0, help="Weight of the local pixel (0: peak from neighbors only, " diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index a85d4eed326..165d5949312 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -37,7 +37,7 @@ extractors.remove(FixedWindowSum) extractors.remove(FlashCamExtractor) -camera_toymodels = ["toymodel", "toymodel_lst", "toymodel_mst_fc"] +camera_toymodels = ["toymodel", "toymodel_sst", "toymodel_lst", "toymodel_mst_fc"] @pytest.fixture(scope="module") @@ -66,6 +66,21 @@ def subarray(prod5_sst, reference_location): reference_pulse_sample_width, u.ns ) readout.sampling_rate = u.Quantity(1 / sample_width, u.GHz) + readout.n_samples = 96 + return subarray + + +@pytest.fixture(scope="module") +def subarray_2_sst(prod5_sst, reference_location): + subarray = SubarrayDescription( + "test array", + tel_positions={1: np.zeros(3) * u.m, 2: np.zeros(3) * u.m}, + tel_descriptions={ + 1: deepcopy(prod5_sst), + 2: deepcopy(prod5_sst), + }, + reference_location=reference_location, + ) return subarray @@ -81,7 +96,7 @@ def subarray_1_LST(prod3_lst, reference_location): @pytest.fixture(scope="module") -def subarray_mst_fc(prod5_mst_flashcam, reference_location): +def subarray_1_mst_fc(prod5_mst_flashcam, reference_location): subarray = SubarrayDescription( "One MST with FlashCam", tel_positions={1: np.zeros(3) * u.m}, @@ -106,6 +121,8 @@ def get_test_toymodel(subarray, minCharge=100, maxCharge=1000): waveform = waveform_model.get_waveform(charge, time, n_samples) selected_gain_channel = np.zeros(charge.size, dtype=np.int64) + if waveform.shape[-3] != 1: + selected_gain_channel = None return waveform, subarray, tel_id, selected_gain_channel, charge, time @@ -115,14 +132,19 @@ def toymodel(subarray: object) -> object: return get_test_toymodel(subarray) +@pytest.fixture(scope="module") +def toymodel_sst(subarray_2_sst: object) -> object: + return get_test_toymodel(subarray_2_sst) + + @pytest.fixture(scope="module") def toymodel_lst(subarray_1_LST: object) -> object: return get_test_toymodel(subarray_1_LST) @pytest.fixture(scope="module") -def toymodel_mst_fc(subarray_mst_fc: object) -> object: - return get_test_toymodel(subarray_mst_fc) +def toymodel_mst_fc(subarray_1_mst_fc: object) -> object: + return get_test_toymodel(subarray_1_mst_fc) def get_test_toymodel_gradient(subarray, minCharge=100, maxCharge=1000): @@ -159,8 +181,8 @@ def get_test_toymodel_gradient(subarray, minCharge=100, maxCharge=1000): @pytest.fixture(scope="module") -def toymodel_mst_fc_time(subarray_mst_fc: object) -> object: - return get_test_toymodel_gradient(subarray_mst_fc) +def toymodel_mst_fc_time(subarray_1_mst_fc: object) -> object: + return get_test_toymodel_gradient(subarray_1_mst_fc) def test_extract_around_peak(toymodel): @@ -386,14 +408,13 @@ def test_extractors(Extractor, toymodels, request): dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid - breakpoint() if dl1.image.ndim == 1: assert_allclose(dl1.image, true_charge, rtol=0.2) - assert_allclose(dl1.peak_time, true_time, rtol=0.1) + assert_allclose(dl1.peak_time, true_time, rtol=0.2) else: for ichannel in range(dl1.image.shape[-2]): assert_allclose(dl1.image[ichannel], true_charge, rtol=0.2) - assert_allclose(dl1.peak_time[ichannel], true_time, rtol=0.1) + assert_allclose(dl1.peak_time[ichannel], true_time, rtol=0.2) @pytest.mark.parametrize("toymodels", camera_toymodels) @@ -442,15 +463,23 @@ def test_fixed_window_sum(toymodel): true_charge, true_time, ) = toymodel + extractor = FixedWindowSum(subarray=subarray, peak_index=47) broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) + assert dl1.is_valid assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) - assert dl1.is_valid + + waveforms = np.append(waveforms, waveforms, axis=-3) + dl1 = extractor(waveforms, tel_id, None, broken_pixels) + for ichannel in range(dl1.image.shape[-2]): + assert_allclose(dl1.image[ichannel], true_charge, rtol=0.1) + assert_allclose(dl1.peak_time[ichannel], true_time, rtol=0.1) -def test_sliding_window_max_sum(toymodel): +@pytest.mark.parametrize("toymodels", camera_toymodels) +def test_sliding_window_max_sum(toymodels, request): ( waveforms, subarray, @@ -458,14 +487,18 @@ def test_sliding_window_max_sum(toymodel): selected_gain_channel, true_charge, true_time, - ) = toymodel + ) = request.getfixturevalue(toymodels) extractor = SlidingWindowMaxSum(subarray=subarray) broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) - print(true_charge, dl1.image, true_time, dl1.peak_time) - assert_allclose(dl1.image, true_charge, rtol=0.1) - assert_allclose(dl1.peak_time, true_time, rtol=0.1) assert dl1.is_valid + if dl1.image.ndim == 1: + assert_allclose(dl1.image, true_charge, rtol=0.1) + assert_allclose(dl1.peak_time, true_time, rtol=0.1) + else: + for ichannel in range(dl1.image.shape[-2]): + assert_allclose(dl1.image[ichannel], true_charge, rtol=0.1) + assert_allclose(dl1.peak_time[ichannel], true_time, rtol=0.1) def test_neighbor_peak_window_sum_local_weight(toymodel): @@ -535,7 +568,7 @@ def test_Two_pass_window_sum_no_noise(subarray_1_LST): # Define the model for the waveforms to fill with the information from # the simulated image - waveform_model = WaveformModel.from_camera_readout(readout) + waveform_model = WaveformModel.from_camera_readout(readout, gain_channel="HIGH") waveforms = waveform_model.get_waveform(true_charge, true_time, n_samples) selected_gain_channel = np.zeros(true_charge.size, dtype=np.int64) @@ -575,11 +608,11 @@ def test_Two_pass_window_sum_no_noise(subarray_1_LST): # Check that we have gained signal charge by using the 2nd pass # This also checks that the test image has triggered the 2nd pass - # (i.e. it is not so bad to have <3 pixels in the preliminary cleaned image) + # (i.e. it is not so bad to have < 3 pixels in the preliminary cleaned image) reco_charge1 = np.sum(dl1_pass1.image[signal_pixels & integration_window_inside]) reco_charge2 = np.sum(dl1_pass2.image[signal_pixels & integration_window_inside]) # since there is no noise in this test, 1st pass will find the peak and 2nd - # can at most to the same + # can at most do the same assert (reco_charge2 / reco_charge1) < 1 assert dl1_pass1.is_valid @@ -648,7 +681,7 @@ def test_extractor_tel_param(toymodel): ) waveforms, subarray, _, _, _, _ = toymodel - _, n_pixels, n_samples = waveforms.shape + _, _, n_samples = waveforms.shape extractor = ImageExtractor.from_name( "FixedWindowSum", subarray=subarray, config=config ) @@ -719,9 +752,9 @@ def test_adaptive_centroid(toymodel_mst_fc): waveforms, subarray, tel_id, - selected_gain_channel, - true_charge, - true_time, + _, + _, + _, ) = toymodel_mst_fc neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse @@ -756,11 +789,11 @@ def test_adaptive_centroid(toymodel_mst_fc): def test_deconvolve(toymodel_mst_fc): ( waveforms, - subarray, - tel_id, - selected_gain_channel, - true_charge, - true_time, + _, + _, + _, + _, + _, ) = toymodel_mst_fc deconvolved_waveforms_0 = deconvolve(waveforms, 0, 0, 0.0) @@ -775,11 +808,11 @@ def test_deconvolve(toymodel_mst_fc): def test_upsampling(toymodel_mst_fc): ( waveforms, - subarray, - tel_id, - selected_gain_channel, - true_charge, - true_time, + _, + _, + _, + _, + _, ) = toymodel_mst_fc upsampling_even = 4 upsampling_odd = 3 @@ -823,7 +856,7 @@ def test_FC_time(toymodel_mst_fc_time): subarray, tel_id, selected_gain_channel, - true_charge, + _, true_time, mask, ) = toymodel_mst_fc_time @@ -865,7 +898,7 @@ def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): tel_id, selected_gain_channel, true_charge, - true_time, + _, ) = toymodel_mst_fc extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) diff --git a/src/ctapipe/image/tests/test_sliding_window_correction.py b/src/ctapipe/image/tests/test_sliding_window_correction.py deleted file mode 100644 index 73625511c95..00000000000 --- a/src/ctapipe/image/tests/test_sliding_window_correction.py +++ /dev/null @@ -1,62 +0,0 @@ -""" -Test of sliding window extractor for LST camera pulse shape with -the correction for the integration window completeness -""" - -import astropy.units as u -import numpy as np -from numpy.testing import assert_allclose -from traitlets.config.loader import Config - -from ctapipe.containers import DL1CameraContainer -from ctapipe.image.extractor import ImageExtractor, SlidingWindowMaxSum -from ctapipe.image.toymodel import WaveformModel -from ctapipe.instrument import SubarrayDescription - - -def test_sw_pulse_lst(prod5_lst, reference_location): - """ - Test function of sliding window extractor for LST camera pulse shape with - the correction for the integration window completeness - """ - - # prepare array with 1 LST - subarray = SubarrayDescription( - "LST1", - tel_positions={1: np.zeros(3) * u.m}, - tel_descriptions={1: prod5_lst}, - reference_location=reference_location, - ) - - tel_id = list(subarray.tel.keys())[0] - - n_pixels = subarray.tel[tel_id].camera.geometry.n_pixels - n_samples = 40 - readout = subarray.tel[tel_id].camera.readout - - random = np.random.RandomState(1) - min_charge = 100 - max_charge = 1000 - charge_true = random.uniform(min_charge, max_charge, n_pixels) - time_true = random.uniform( - n_samples // 2 - 1, n_samples // 2 + 1, n_pixels - ) / readout.sampling_rate.to_value(u.GHz) - - waveform_model = WaveformModel.from_camera_readout(readout) - waveform = waveform_model.get_waveform(charge_true, time_true, n_samples) - selected_gain_channel = np.zeros(charge_true.size, dtype=np.int8) - - # define extractor - config = Config({"SlidingWindowMaxSum": {"window_width": 8}}) - extractor = SlidingWindowMaxSum(subarray=subarray) - extractor = ImageExtractor.from_name( - "SlidingWindowMaxSum", subarray=subarray, config=config - ) - - broken_pixels = np.zeros(n_pixels, dtype=bool) - dl1: DL1CameraContainer = extractor( - waveform, tel_id, selected_gain_channel, broken_pixels - ) - print(dl1.image / charge_true) - assert_allclose(dl1.image, charge_true, rtol=0.02) - assert dl1.is_valid diff --git a/src/ctapipe/image/tests/test_toy.py b/src/ctapipe/image/tests/test_toy.py index 8229163fc64..f5b6fa96c52 100644 --- a/src/ctapipe/image/tests/test_toy.py +++ b/src/ctapipe/image/tests/test_toy.py @@ -224,7 +224,7 @@ def test_waveform_model(frame, prod5_sst): ref_x_norm = np.linspace(0, ref_duration, n_ref_samples) ref_y_norm = norm.pdf(ref_x_norm, ref_duration / 2, pulse_sigma) - readout.reference_pulse_shape = ref_y_norm[np.newaxis, :] + readout.reference_pulse_shape = np.array([ref_y_norm, ref_y_norm]) readout.reference_pulse_sample_width = u.Quantity( ref_x_norm[1] - ref_x_norm[0], u.ns ) @@ -254,20 +254,28 @@ def test_waveform_model(frame, prod5_sst): waveform_model = WaveformModel.from_camera_readout(readout) waveform = waveform_model.get_waveform(charge, time, 96) - np.testing.assert_allclose(waveform.sum(axis=1), charge, rtol=1e-3) + assert waveform.shape[-3] == 2 + np.testing.assert_allclose(waveform.sum(axis=-1)[0], charge, rtol=1e-3) np.testing.assert_allclose( - waveform.argmax(axis=1) / readout.sampling_rate.to_value(u.GHz), time, rtol=1e-1 + waveform.argmax(axis=-1)[0] / readout.sampling_rate.to_value(u.GHz), + time, + rtol=1e-1, ) time_2 = time + 1 time_2[charge == 0] = 0 waveform_2 = waveform_model.get_waveform(charge, time_2, 96) - np.testing.assert_allclose(waveform_2.sum(axis=1), charge, rtol=1e-3) + np.testing.assert_allclose(waveform_2.sum(axis=-1)[0], charge, rtol=1e-3) np.testing.assert_allclose( - waveform_2.argmax(axis=1) / readout.sampling_rate.to_value(u.GHz), + waveform_2.argmax(axis=-1)[0] / readout.sampling_rate.to_value(u.GHz), time_2, rtol=1e-1, ) assert ( - waveform_2.argmax(axis=1)[charge != 0] > waveform.argmax(axis=1)[charge != 0] + waveform_2.argmax(axis=-1)[0, charge != 0] + > waveform.argmax(axis=-1)[0, charge != 0] ).all() + + waveform_model = WaveformModel.from_camera_readout(readout, gain_channel="HIGH") + waveform = waveform_model.get_waveform(charge, time, 96) + assert waveform.shape[-3] == 1 From 6268cbb2866f5d18da1dbcf44fdf7336020f2949 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 27 Mar 2024 17:03:28 +0100 Subject: [PATCH 05/42] Adapting shift_waveforms test - now testing that it also works for n_channels --- .../calib/camera/tests/test_calibrator.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index 27bebb4f02a..c43458c0971 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -240,20 +240,20 @@ def test_dl1_charge_calib(example_subarray): def test_shift_waveforms(): from ctapipe.calib.camera.calibrator import shift_waveforms - # 5 pixels, 40 samples - waveforms = np.zeros((5, 40)) - waveforms[:, 10] = 1 + # 2 channels, 5 pixels, 40 samples + waveforms = np.zeros((2, 5, 40)) + waveforms[:, :, 10] = 1 shifts = np.array([1.4, 2.1, -1.8, 3.1, -4.4]) shifted_waveforms, remaining_shift = shift_waveforms(waveforms, shifts) assert np.allclose(remaining_shift, [0.4, 0.1, 0.2, 0.1, -0.4]) - assert shifted_waveforms[0, 9] == 1 - assert shifted_waveforms[1, 8] == 1 - assert shifted_waveforms[2, 12] == 1 - assert shifted_waveforms[3, 7] == 1 - assert shifted_waveforms[4, 14] == 1 + assert shifted_waveforms[:, 0, 9].all() == 1 + assert shifted_waveforms[:, 1, 8].all() == 1 + assert shifted_waveforms[:, 2, 12].all() == 1 + assert shifted_waveforms[:, 3, 7].all() == 1 + assert shifted_waveforms[:, 4, 14].all() == 1 def test_invalid_pixels(example_event, example_subarray): From f4a2916f4e6cafb39d3d81521ffa9d30bcb4d52a Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 27 Mar 2024 17:10:01 +0100 Subject: [PATCH 06/42] Fix reducer --- src/ctapipe/image/reducer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/image/reducer.py b/src/ctapipe/image/reducer.py index 17bb5e77e46..c79f1602216 100644 --- a/src/ctapipe/image/reducer.py +++ b/src/ctapipe/image/reducer.py @@ -187,8 +187,8 @@ def select_pixels(self, waveforms, tel_id=None, selected_gain_channel=None): dl1: DL1CameraContainer = extractor( waveforms, tel_id=tel_id, - broken_pixels=broken_pixels, selected_gain_channel=selected_gain_channel, + broken_pixels=broken_pixels, ) # 1) Step: TailcutCleaning at first From d57dedd863b704dc87717c815fead22af7ef6ede Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 28 Mar 2024 14:01:07 +0100 Subject: [PATCH 07/42] Change GainSelector handle new waveforms shape - With the current API of the GainSelector its not possible anymore - to distinguish between gain selected waveforms and waveforms - that only have one gain channel --- src/ctapipe/calib/camera/gainselection.py | 11 ++++++----- src/ctapipe/calib/camera/tests/test_gainselection.py | 9 --------- src/ctapipe/io/simteleventsource.py | 12 +++--------- 3 files changed, 9 insertions(+), 23 deletions(-) diff --git a/src/ctapipe/calib/camera/gainselection.py b/src/ctapipe/calib/camera/gainselection.py index 6ba3e550076..39dad985ebc 100644 --- a/src/ctapipe/calib/camera/gainselection.py +++ b/src/ctapipe/calib/camera/gainselection.py @@ -1,6 +1,7 @@ """ Algorithms to select correct gain channel """ + from abc import abstractmethod from enum import IntEnum @@ -47,16 +48,16 @@ def __call__(self, waveforms): Shape: n_pix Dtype: int8 """ - if waveforms.ndim == 2: # Return None if already gain selected - return None - elif waveforms.ndim == 3: + if waveforms.ndim == 3: n_channels, n_pixels, _ = waveforms.shape - if n_channels == 1: # Must be first channel if only one channel + if n_channels == 1: return np.zeros(n_pixels, dtype=np.int8) else: return self.select_channel(waveforms) else: - raise ValueError(f"Cannot handle waveform array of shape: {waveforms.ndim}") + raise ValueError( + f"Cannot handle waveform array of shape: {waveforms.shape}" + ) @abstractmethod def select_channel(self, waveforms): diff --git a/src/ctapipe/calib/camera/tests/test_gainselection.py b/src/ctapipe/calib/camera/tests/test_gainselection.py index 254fc0b31b9..e06d8d3027d 100644 --- a/src/ctapipe/calib/camera/tests/test_gainselection.py +++ b/src/ctapipe/calib/camera/tests/test_gainselection.py @@ -23,15 +23,6 @@ def test_gain_selector(): np.testing.assert_equal(selected_gain_channel, 0) -def test_pre_selected(): - shape = (2048, 128) - waveforms = np.zeros(shape) - - gain_selector = DummyGainSelector() - selected_gain_channel = gain_selector(waveforms) - assert selected_gain_channel is None - - def test_single_channel(): shape = (1, 2048, 128) waveforms = np.zeros(shape) diff --git a/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py index a1b89ecbe36..5853f2b9277 100644 --- a/src/ctapipe/io/simteleventsource.py +++ b/src/ctapipe/io/simteleventsource.py @@ -277,21 +277,15 @@ def apply_simtel_r1_calibration( The gain channel selected for each pixel Shape: (n_pixels) """ - n_channels, n_pixels, _ = r0_waveforms.shape + n_pixels = r0_waveforms.shape[-2] ped = pedestal[..., np.newaxis] DC_to_PHE = dc_to_pe[..., np.newaxis] gain = DC_to_PHE * calib_scale r1_waveforms = (r0_waveforms - ped) * gain + calib_shift - if n_channels == 1: - selected_gain_channel = np.zeros(n_pixels, dtype=np.int8) - r1_waveforms = r1_waveforms - else: - selected_gain_channel = gain_selector(r0_waveforms) - r1_waveforms = r1_waveforms[ - np.newaxis, selected_gain_channel, np.arange(n_pixels) - ] + selected_gain_channel = gain_selector(r0_waveforms) + r1_waveforms = r1_waveforms[np.newaxis, selected_gain_channel, np.arange(n_pixels)] return r1_waveforms, selected_gain_channel From d7d450c689d5d9234eca4ef077406345f8eb78ff Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 28 Mar 2024 16:37:50 +0100 Subject: [PATCH 08/42] Adapt some comment --- src/ctapipe/image/tests/test_extractor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 165d5949312..b29d80e8eba 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -444,7 +444,8 @@ def test_integration_correction_off(Extractor, toymodels, request): assert dl1.is_valid # peak time should stay the same - # charge should be too small without correction + # charge should be too small without correction for the used reference pulse + # shapes (not in general). if dl1.image.ndim == 1: assert np.all(dl1.image <= true_charge) assert_allclose(dl1.peak_time, true_time, rtol=0.1) From 1411bdd4501b4f0c293de65419a4fd9956fdecc4 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 2 Apr 2024 13:31:49 +0200 Subject: [PATCH 09/42] Minor fixes --- src/ctapipe/calib/camera/calibrator.py | 2 +- src/ctapipe/calib/camera/gainselection.py | 3 +-- src/ctapipe/calib/camera/tests/test_calibrator.py | 10 +++++----- src/ctapipe/image/extractor.py | 3 +-- src/ctapipe/image/tests/test_extractor.py | 8 ++++---- src/ctapipe/image/toymodel.py | 2 +- 6 files changed, 13 insertions(+), 15 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 23328458f96..58bc0b899a7 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -227,7 +227,7 @@ def _calibrate_dl1(self, event, tel_id): # subtract any remaining pedestal before extraction if dl1_calib.pedestal_offset is not None: # this copies intentionally, we don't want to modify the dl0 data - # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels, ) + # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels,) waveforms = waveforms - dl1_calib.pedestal_offset[np.newaxis, :, np.newaxis] if n_samples == 1: diff --git a/src/ctapipe/calib/camera/gainselection.py b/src/ctapipe/calib/camera/gainselection.py index 39dad985ebc..b32a0cdced4 100644 --- a/src/ctapipe/calib/camera/gainselection.py +++ b/src/ctapipe/calib/camera/gainselection.py @@ -52,8 +52,7 @@ def __call__(self, waveforms): n_channels, n_pixels, _ = waveforms.shape if n_channels == 1: return np.zeros(n_pixels, dtype=np.int8) - else: - return self.select_channel(waveforms) + return self.select_channel(waveforms) else: raise ValueError( f"Cannot handle waveform array of shape: {waveforms.shape}" diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index c43458c0971..6a59530d792 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -249,11 +249,11 @@ def test_shift_waveforms(): assert np.allclose(remaining_shift, [0.4, 0.1, 0.2, 0.1, -0.4]) - assert shifted_waveforms[:, 0, 9].all() == 1 - assert shifted_waveforms[:, 1, 8].all() == 1 - assert shifted_waveforms[:, 2, 12].all() == 1 - assert shifted_waveforms[:, 3, 7].all() == 1 - assert shifted_waveforms[:, 4, 14].all() == 1 + assert (shifted_waveforms[:, 0, 9] == 1).all() + assert (shifted_waveforms[:, 1, 8] == 1).all() + assert (shifted_waveforms[:, 2, 12] == 1).all() + assert (shifted_waveforms[:, 3, 7] == 1).all() + assert (shifted_waveforms[:, 4, 14] == 1).all() def test_invalid_pixels(example_event, example_subarray): diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index c58ce61928d..51524dd7574 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -1308,8 +1308,7 @@ def __call__( raise AttributeError( "The data needs to be gain selected to use the TwoPassWindowSum." ) - else: - waveforms = waveforms[0, :, :] + waveforms = waveforms[0, :, :] charge1, pulse_time1, correction1 = self._apply_first_pass(waveforms, tel_id) diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index b29d80e8eba..a15c7c1587d 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -128,22 +128,22 @@ def get_test_toymodel(subarray, minCharge=100, maxCharge=1000): @pytest.fixture(scope="module") -def toymodel(subarray: object) -> object: +def toymodel(subarray): return get_test_toymodel(subarray) @pytest.fixture(scope="module") -def toymodel_sst(subarray_2_sst: object) -> object: +def toymodel_sst(subarray_2_sst): return get_test_toymodel(subarray_2_sst) @pytest.fixture(scope="module") -def toymodel_lst(subarray_1_LST: object) -> object: +def toymodel_lst(subarray_1_LST): return get_test_toymodel(subarray_1_LST) @pytest.fixture(scope="module") -def toymodel_mst_fc(subarray_1_mst_fc: object) -> object: +def toymodel_mst_fc(subarray_1_mst_fc): return get_test_toymodel(subarray_1_mst_fc) diff --git a/src/ctapipe/image/toymodel.py b/src/ctapipe/image/toymodel.py index ba9cc167600..ae482ae9b7d 100644 --- a/src/ctapipe/image/toymodel.py +++ b/src/ctapipe/image/toymodel.py @@ -185,7 +185,7 @@ def from_camera_readout(cls, readout, gain_channel="ALL"): readout : `ctapipe.instrument.CameraReadout` gain_channel : str The reference pulse gain channel to use. - Choose between `HIGH`, `LOW` and `ALL`. + Choose between 'HIGH', 'LOW' and 'ALL'. Returns ------- From 9e7b2b5cb73debf03d2d8dd9be48333d9221cbc9 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 14:17:26 +0200 Subject: [PATCH 10/42] Make CameraCalibrator work for not gain selected data --- src/ctapipe/calib/camera/calibrator.py | 16 ++++++++++------ .../calib/camera/tests/test_calibrator.py | 4 ++-- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 58bc0b899a7..993725bf441 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -228,7 +228,9 @@ def _calibrate_dl1(self, event, tel_id): if dl1_calib.pedestal_offset is not None: # this copies intentionally, we don't want to modify the dl0 data # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels,) - waveforms = waveforms - dl1_calib.pedestal_offset[np.newaxis, :, np.newaxis] + waveforms = ( + waveforms - np.atleast_2d(dl1_calib.pedestal_offset)[..., np.newaxis] + ) if n_samples == 1: # To handle ASTRI and dst @@ -238,7 +240,7 @@ def _calibrate_dl1(self, event, tel_id): # - Don't do anything if dl1 container already filled # - Update on SST review decision dl1 = DL1CameraContainer( - image=waveforms[0, ..., 0].astype(np.float32), + image=np.squeeze(waveforms).astype(np.float32), peak_time=np.zeros(n_pixels, dtype=np.float32), is_valid=True, ) @@ -265,10 +267,12 @@ def _calibrate_dl1(self, event, tel_id): # correct non-integer remainder of the shift if given if self.apply_peak_time_shift.tel[tel_id] and time_shift is not None: - dl1.peak_time -= remaining_shift + dl1.peak_time = (dl1.peak_time.T - remaining_shift).T # Calibrate extracted charge - dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor + dl1.image = ( + dl1.image.T * (dl1_calib.relative_factor / dl1_calib.absolute_factor) + ).T # handle invalid pixels if self.invalid_pixel_handler is not None: @@ -323,8 +327,8 @@ def shift_waveforms(waveforms, time_shift_samples): remaining_shift: ndarray of shape (n_pixels, ) The remaining shift after applying the integer shift to the waveforms. """ - mean_shift = time_shift_samples.mean() - integer_shift = np.round(time_shift_samples - mean_shift).astype("int16") + mean_shift = time_shift_samples.mean(-1) + integer_shift = np.round((time_shift_samples.T - mean_shift).T).astype("int16") remaining_shift = time_shift_samples - integer_shift shifted_waveforms = _shift_waveforms_by_integer(waveforms, integer_shift) return shifted_waveforms, remaining_shift diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index 6a59530d792..ce1c4b13107 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -240,8 +240,8 @@ def test_dl1_charge_calib(example_subarray): def test_shift_waveforms(): from ctapipe.calib.camera.calibrator import shift_waveforms - # 2 channels, 5 pixels, 40 samples - waveforms = np.zeros((2, 5, 40)) + # 1 channel, 5 pixels, 40 samples + waveforms = np.zeros((1, 5, 40)) waveforms[:, :, 10] = 1 shifts = np.array([1.4, 2.1, -1.8, 3.1, -4.4]) From 151e9c0c9cba0446161c47b5b2b67d591d0a91d7 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 17:07:41 +0200 Subject: [PATCH 11/42] Change datamodel version --- src/ctapipe/io/datawriter.py | 4 +++- src/ctapipe/io/hdf5eventsource.py | 4 +--- src/ctapipe/io/tests/test_hdf5eventsource.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/io/datawriter.py b/src/ctapipe/io/datawriter.py index eb353d43620..d0f93905f66 100644 --- a/src/ctapipe/io/datawriter.py +++ b/src/ctapipe/io/datawriter.py @@ -48,8 +48,10 @@ def _get_tel_index(event, tel_id): # (meaning readers need to update scripts) # - increase the minor number if new columns or datasets are added # - increase the patch number if there is a small bugfix to the model. -DATA_MODEL_VERSION = "v5.1.0" +DATA_MODEL_VERSION = "v6.0.0" DATA_MODEL_CHANGE_HISTORY = """ +- v6.0.0: - Change R1- and DL0-waveform shape from (n_pixels, n_samples) to be always + (n_channels, n_pixels, n_samples). - v5.1.0: - Remove redundant 'is_valid' column in ``DispContainer``. - Rename content of ``DispContainer`` from 'norm' to 'parameter' and use the same default prefix ('disp') for all containers filled by ``DispReconstructor``. diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index a451bde4116..4b352b5953f 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -71,9 +71,7 @@ COMPATIBLE_DATA_MODEL_VERSIONS = [ - "v4.0.0", - "v5.0.0", - "v5.1.0", + "v6.0.0", ] diff --git a/src/ctapipe/io/tests/test_hdf5eventsource.py b/src/ctapipe/io/tests/test_hdf5eventsource.py index bc8b1aa33d9..eef8d2b0fa2 100644 --- a/src/ctapipe/io/tests/test_hdf5eventsource.py +++ b/src/ctapipe/io/tests/test_hdf5eventsource.py @@ -30,7 +30,7 @@ def test_is_compatible(compatible_file, request): def test_metadata(dl1_file): with HDF5EventSource(input_url=dl1_file) as source: assert source.is_simulation - assert source.datamodel_version == (5, 1, 0) + assert source.datamodel_version == (6, 0, 0) assert set(source.datalevels) == { DataLevel.DL1_IMAGES, DataLevel.DL1_PARAMETERS, From e93b35cc3ba223b6a8f4ced0814daa2da1dc13ea Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 17:17:39 +0200 Subject: [PATCH 12/42] Test shift_waveforms for both use cases - 1 gain channel and 2 gain channel --- .../calib/camera/tests/test_calibrator.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index ce1c4b13107..f768ba69316 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -249,6 +249,23 @@ def test_shift_waveforms(): assert np.allclose(remaining_shift, [0.4, 0.1, 0.2, 0.1, -0.4]) + assert shifted_waveforms[0, 0, 9] == 1 + assert shifted_waveforms[0, 1, 8] == 1 + assert shifted_waveforms[0, 2, 12] == 1 + assert shifted_waveforms[0, 3, 7] == 1 + assert shifted_waveforms[0, 4, 14] == 1 + + # 2 channel, 5 pixels, 40 samples + waveforms = np.zeros((2, 5, 40)) + waveforms[:, :, 10] = 1 + shifts = np.array([[1.4, 2.1, -1.8, 3.1, -4.4], [1.4, 2.1, -1.8, 3.1, -4.4]]) + + shifted_waveforms, remaining_shift = shift_waveforms(waveforms, shifts) + + assert np.allclose( + remaining_shift, [[0.4, 0.1, 0.2, 0.1, -0.4], [0.4, 0.1, 0.2, 0.1, -0.4]] + ) + assert (shifted_waveforms[:, 0, 9] == 1).all() assert (shifted_waveforms[:, 1, 8] == 1).all() assert (shifted_waveforms[:, 2, 12] == 1).all() From a69308876f30a70d84f429e5b6ce68e55fb0115e Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 17:26:27 +0200 Subject: [PATCH 13/42] Test wrong ndim in GainSelector --- src/ctapipe/calib/camera/gainselection.py | 12 ++++++------ src/ctapipe/calib/camera/tests/test_gainselection.py | 11 +++++++++++ 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/calib/camera/gainselection.py b/src/ctapipe/calib/camera/gainselection.py index b32a0cdced4..c1954f3e3dc 100644 --- a/src/ctapipe/calib/camera/gainselection.py +++ b/src/ctapipe/calib/camera/gainselection.py @@ -48,16 +48,16 @@ def __call__(self, waveforms): Shape: n_pix Dtype: int8 """ - if waveforms.ndim == 3: - n_channels, n_pixels, _ = waveforms.shape - if n_channels == 1: - return np.zeros(n_pixels, dtype=np.int8) - return self.select_channel(waveforms) - else: + if waveforms.ndim != 3: raise ValueError( f"Cannot handle waveform array of shape: {waveforms.shape}" ) + n_channels, n_pixels, _ = waveforms.shape + if n_channels == 1: + return np.zeros(n_pixels, dtype=np.int8) + return self.select_channel(waveforms) + @abstractmethod def select_channel(self, waveforms): """ diff --git a/src/ctapipe/calib/camera/tests/test_gainselection.py b/src/ctapipe/calib/camera/tests/test_gainselection.py index e06d8d3027d..6966791673d 100644 --- a/src/ctapipe/calib/camera/tests/test_gainselection.py +++ b/src/ctapipe/calib/camera/tests/test_gainselection.py @@ -1,4 +1,5 @@ import numpy as np +import pytest from ctapipe.calib.camera.gainselection import ( GainChannel, @@ -23,6 +24,16 @@ def test_gain_selector(): np.testing.assert_equal(selected_gain_channel, 0) +def test_ndim(): + shape = (5, 2, 2048, 128) + waveforms = np.indices(shape)[1] + waveforms[1] *= 2 + + gain_selector = DummyGainSelector() + with pytest.raises(ValueError): + gain_selector(waveforms) + + def test_single_channel(): shape = (1, 2048, 128) waveforms = np.zeros(shape) From 7c5a4d0272d303a6567834743c38b2e0d11ed8e9 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 17:32:52 +0200 Subject: [PATCH 14/42] Add toymodel test for wrong gain channel --- src/ctapipe/image/tests/test_toy.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ctapipe/image/tests/test_toy.py b/src/ctapipe/image/tests/test_toy.py index f5b6fa96c52..e7c5bbe7168 100644 --- a/src/ctapipe/image/tests/test_toy.py +++ b/src/ctapipe/image/tests/test_toy.py @@ -279,3 +279,6 @@ def test_waveform_model(frame, prod5_sst): waveform_model = WaveformModel.from_camera_readout(readout, gain_channel="HIGH") waveform = waveform_model.get_waveform(charge, time, 96) assert waveform.shape[-3] == 1 + + with pytest.raises(ValueError): + WaveformModel.from_camera_readout(readout, gain_channel=0) From 7a3cf785199b74c704ecc2624780d59d58d5603c Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 19:16:48 +0200 Subject: [PATCH 15/42] Update test for neighbor_average_maximum --- src/ctapipe/image/tests/test_extractor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index a15c7c1587d..512dadaa391 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -299,7 +299,7 @@ def test_neighbor_average_peakpos(toymodels, request): expected_average = waveforms[:, nei_pixel].sum(1) / len(nei_pixel) expected_peak_pos = np.argmax(expected_average, axis=-1) for ichannel in range(waveforms.shape[-3]): - assert (peak_pos[ichannel][pixel] == expected_peak_pos).all() + assert peak_pos[ichannel][pixel] == expected_peak_pos[ichannel] local_weight = 4 peak_pos = neighbor_average_maximum( @@ -316,7 +316,7 @@ def test_neighbor_average_peakpos(toymodels, request): expected_average = waveforms[:, nei_pixel].sum(1) / len(nei_pixel) expected_peak_pos = np.argmax(expected_average, axis=-1) for ichannel in range(waveforms.shape[-3]): - assert (peak_pos[ichannel][pixel] == expected_peak_pos).all() + assert peak_pos[ichannel][pixel] == expected_peak_pos[ichannel] def test_extract_peak_time_within_range(): From 279e6bfda27edf524a97c61c4a269c644880b487 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 3 Apr 2024 20:16:07 +0200 Subject: [PATCH 16/42] Fix GlobalPeakWindowSum for 2 gains --- src/ctapipe/image/extractor.py | 2 +- src/ctapipe/image/tests/test_extractor.py | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 51524dd7574..377706ea340 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -584,7 +584,7 @@ def __call__( ] # average over brightest pixels then argmax over samples - peak_index = waveforms[:, brightest].mean(axis=-2).argmax(axis=-1) + peak_index = waveforms[:, brightest].mean(axis=(-2, -3)).argmax(axis=-1) charge, peak_time = extract_around_peak( waveforms, diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 512dadaa391..be50c08cb19 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -723,12 +723,12 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray): # signal in dim pixels is in slice 10, signal in bright pixels is in slice 30 waveforms = np.zeros((1, n_pixels, 50), dtype="float64") - waveforms[0][~bright_pixels, 9] = 3 - waveforms[0][~bright_pixels, 10] = 5 - waveforms[0][~bright_pixels, 11] = 2 - waveforms[0][bright_pixels, 29] = 5 - waveforms[0][bright_pixels, 30] = 10 - waveforms[0][bright_pixels, 31] = 3 + waveforms[:, ~bright_pixels, 9] = 3 + waveforms[:, ~bright_pixels, 10] = 5 + waveforms[:, ~bright_pixels, 11] = 2 + waveforms[:, bright_pixels, 29] = 5 + waveforms[:, bright_pixels, 30] = 10 + waveforms[:, bright_pixels, 31] = 3 extractor = GlobalPeakWindowSum( subarray=subarray, From 9a65df2bbe5524b89ca5c11c003f68af911fdfa3 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 4 Apr 2024 11:20:42 +0200 Subject: [PATCH 17/42] Add Changelog --- docs/changes/2529.datamodel.rst | 6 ++++++ src/ctapipe/image/tests/test_extractor.py | 10 +++++----- 2 files changed, 11 insertions(+), 5 deletions(-) create mode 100644 docs/changes/2529.datamodel.rst diff --git a/docs/changes/2529.datamodel.rst b/docs/changes/2529.datamodel.rst new file mode 100644 index 00000000000..e08c54a2c61 --- /dev/null +++ b/docs/changes/2529.datamodel.rst @@ -0,0 +1,6 @@ +Change R1- and DL0-waveforms datamodel shape from (n_pixels, n_samples) +to be always (n_channels, n_pixels, n_samples). This breaks reading in single +gain data of older versions (< v6.0.0) with the ``HDF5EventSource``. + +Re-introduce also the possibility of running ``ImageExtractor``s on data +consisting of multiple gain channels. diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index be50c08cb19..dd12ac4a065 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -85,7 +85,7 @@ def subarray_2_sst(prod5_sst, reference_location): @pytest.fixture(scope="module") -def subarray_1_LST(prod3_lst, reference_location): +def subarray_1_lst(prod3_lst, reference_location): subarray = SubarrayDescription( "One LST", tel_positions={1: np.zeros(3) * u.m}, @@ -138,8 +138,8 @@ def toymodel_sst(subarray_2_sst): @pytest.fixture(scope="module") -def toymodel_lst(subarray_1_LST): - return get_test_toymodel(subarray_1_LST) +def toymodel_lst(subarray_1_lst): + return get_test_toymodel(subarray_1_lst) @pytest.fixture(scope="module") @@ -520,10 +520,10 @@ def test_neighbor_peak_window_sum_local_weight(toymodel): assert dl1.is_valid -def test_Two_pass_window_sum_no_noise(subarray_1_LST): +def test_Two_pass_window_sum_no_noise(subarray_1_lst): rng = np.random.default_rng(0) - subarray = subarray_1_LST + subarray = subarray_1_lst camera = subarray.tel[1].camera geometry = camera.geometry From 7633667ec97a0eda896c9c6b62f734619e14c670 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 4 Apr 2024 13:49:20 +0200 Subject: [PATCH 18/42] Fix docs - close remaining open file - fix changelog --- docs/changes/2529.datamodel.rst | 2 +- examples/core/table_writer_reader.py | 10 +++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/docs/changes/2529.datamodel.rst b/docs/changes/2529.datamodel.rst index e08c54a2c61..06fb49fb081 100644 --- a/docs/changes/2529.datamodel.rst +++ b/docs/changes/2529.datamodel.rst @@ -2,5 +2,5 @@ Change R1- and DL0-waveforms datamodel shape from (n_pixels, n_samples) to be always (n_channels, n_pixels, n_samples). This breaks reading in single gain data of older versions (< v6.0.0) with the ``HDF5EventSource``. -Re-introduce also the possibility of running ``ImageExtractor``s on data +Re-introduce also the possibility of running ``ImageExtractor``\s on data consisting of multiple gain channels. diff --git a/examples/core/table_writer_reader.py b/examples/core/table_writer_reader.py index 7a4ce4d1bce..7c8e54cbc2f 100644 --- a/examples/core/table_writer_reader.py +++ b/examples/core/table_writer_reader.py @@ -15,7 +15,6 @@ """ - ###################################################################### # Caveats to think about: \* vector columns in Containers *can* be # written, but some lilbraries like Pandas can not read those (so you must @@ -228,6 +227,7 @@ def create_stream(n_event): # note that here we can still access the metadata # + table.attrs ###################################################################### @@ -236,6 +236,14 @@ def create_stream(n_event): h5.close() +###################################################################### +# close the file afterwards +# + + +h5.close() + + ###################################################################### # Reading one-row-at-a-time: # ~~~~~~~~~~~~~~~~~~~~~~~~~~ From 074f40507868900e4d43fe5318c0b9c4d7500dae Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 4 Apr 2024 14:07:59 +0200 Subject: [PATCH 19/42] Minor beauty changes --- src/ctapipe/image/extractor.py | 5 ++++- src/ctapipe/image/tests/test_extractor.py | 27 +++-------------------- 2 files changed, 7 insertions(+), 25 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 377706ea340..3be48dc5f27 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -584,8 +584,11 @@ def __call__( ] # average over brightest pixels then argmax over samples - peak_index = waveforms[:, brightest].mean(axis=(-2, -3)).argmax(axis=-1) + peak_index = ( + waveforms[:, brightest][:, 0, ...].mean(axis=-2).argmax(axis=-1) + ) + breakpoint() charge, peak_time = extract_around_peak( waveforms, peak_index[:, np.newaxis], diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index dd12ac4a065..3f7f08ab303 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -749,14 +749,7 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray): def test_adaptive_centroid(toymodel_mst_fc): - ( - waveforms, - subarray, - tel_id, - _, - _, - _, - ) = toymodel_mst_fc + waveforms, subarray, tel_id, _, _, _ = toymodel_mst_fc neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) @@ -788,14 +781,7 @@ def test_adaptive_centroid(toymodel_mst_fc): def test_deconvolve(toymodel_mst_fc): - ( - waveforms, - _, - _, - _, - _, - _, - ) = toymodel_mst_fc + waveforms, _, _, _, _, _ = toymodel_mst_fc deconvolved_waveforms_0 = deconvolve(waveforms, 0, 0, 0.0) @@ -807,14 +793,7 @@ def test_deconvolve(toymodel_mst_fc): def test_upsampling(toymodel_mst_fc): - ( - waveforms, - _, - _, - _, - _, - _, - ) = toymodel_mst_fc + waveforms, _, _, _, _, _ = toymodel_mst_fc upsampling_even = 4 upsampling_odd = 3 filt_even = np.ones(upsampling_even) From e50eb7282cd2ddf4e87525aff5959bbdde50f468 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Thu, 4 Apr 2024 14:56:09 +0200 Subject: [PATCH 20/42] Delete breakpoint ...-.- --- src/ctapipe/image/extractor.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 3be48dc5f27..bc306122f61 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -588,7 +588,6 @@ def __call__( waveforms[:, brightest][:, 0, ...].mean(axis=-2).argmax(axis=-1) ) - breakpoint() charge, peak_time = extract_around_peak( waveforms, peak_index[:, np.newaxis], From 89bf1ddebf41a4e730a101aef7ce08c4af72722e Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 5 Apr 2024 16:33:18 +0200 Subject: [PATCH 21/42] Keep compatibility for older datamodel versions --- src/ctapipe/io/hdf5eventsource.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index 4b352b5953f..afbe2153d2c 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -71,6 +71,9 @@ COMPATIBLE_DATA_MODEL_VERSIONS = [ + "v4.0.0", + "v5.0.0", + "v5.1.0", "v6.0.0", ] @@ -650,6 +653,11 @@ def _generator(self): if DataLevel.R1 in self.datalevels: data.r1.tel[tel_id] = next(waveform_readers[key]) + # TODO: Delete if we drop support for datamodel versions < v6.0.0 + r1_waveform = data.r1.tel[tel_id].waveform + if r1_waveform.ndim == 2: + data.r1.tel[tel_id].waveform = r1_waveform[np.newaxis, ...] + if self.has_simulated_dl1: simulated = data.simulation.tel[tel_id] From 0bad55cf4f3198c1105a124c0cc69af7bb0d3d03 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Fri, 5 Apr 2024 17:45:28 +0200 Subject: [PATCH 22/42] Create helper function for applying integration correction --- src/ctapipe/image/extractor.py | 34 ++++++++++++++-------------------- 1 file changed, 14 insertions(+), 20 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index bc306122f61..ce8b335d568 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -358,6 +358,15 @@ def integration_correction( return correction +def _apply_correction(charge, correction, selected_gain_channel): + """ + Helper function for applying the integration correction for certain `ImageExtractor`s. + """ + if selected_gain_channel is None: + return charge * correction[:, np.newaxis] + return charge * correction[selected_gain_channel] + + class ImageExtractor(TelescopeComponent): def __init__(self, subarray, config=None, parent=None, **kwargs): """ @@ -498,10 +507,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - if selected_gain_channel is None: - charge *= correction[:, np.newaxis] - else: - charge *= correction[selected_gain_channel] + charge = _apply_correction(charge, correction, selected_gain_channel) return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -597,10 +603,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - if selected_gain_channel is None: - charge *= correction[:, np.newaxis] - else: - charge *= correction[selected_gain_channel] + charge = _apply_correction(charge, correction, selected_gain_channel) return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -668,10 +671,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - if selected_gain_channel is None: - charge *= correction[:, np.newaxis] - else: - charge *= correction[selected_gain_channel] + charge = _apply_correction(charge, correction, selected_gain_channel) return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -750,10 +750,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - if selected_gain_channel is None: - charge *= correction[:, np.newaxis] - else: - charge *= correction[selected_gain_channel] + charge = _apply_correction(charge, correction, selected_gain_channel) return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) @@ -834,10 +831,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - if selected_gain_channel is None: - charge *= correction[:, np.newaxis] - else: - charge *= correction[selected_gain_channel] + charge = _apply_correction(charge, correction, selected_gain_channel) return DL1CameraContainer( image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True ) From a76ea8f31a2c06aec5b7af8b37f7ccb8a04aea1b Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 8 Apr 2024 11:33:33 +0200 Subject: [PATCH 23/42] Keep dtype after correction --- src/ctapipe/image/extractor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index ce8b335d568..3d8d1e09616 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -363,8 +363,8 @@ def _apply_correction(charge, correction, selected_gain_channel): Helper function for applying the integration correction for certain `ImageExtractor`s. """ if selected_gain_channel is None: - return charge * correction[:, np.newaxis] - return charge * correction[selected_gain_channel] + return (charge * correction[:, np.newaxis]).astype(charge.dtype) + return (charge * correction[selected_gain_channel]).astype(charge.dtype) class ImageExtractor(TelescopeComponent): From 3b82775ff44a1415e9b3d1f6034c5d43ca9c523a Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 8 Apr 2024 14:18:34 +0200 Subject: [PATCH 24/42] Handle broken pixels with 2 gain channels - Adapted extractors and calibrator to handel broken pixels for multiple gain channels - InvalidPixelHandel cannot handle multiple gain channels for now --- src/ctapipe/calib/camera/calibrator.py | 5 ++++- src/ctapipe/image/extractor.py | 17 ++++++++++++----- src/ctapipe/image/invalid_pixels.py | 3 ++- src/ctapipe/image/tests/test_extractor.py | 20 ++++++++++++++++---- 4 files changed, 34 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 993725bf441..f9e5957f97c 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -31,7 +31,10 @@ def _get_invalid_pixels(n_pixels, pixel_status, selected_gain_channel): ) for mask in masks: if mask is not None: - broken_pixels |= mask[selected_gain_channel, index] + if selected_gain_channel is not None: + broken_pixels |= mask[selected_gain_channel, index] + else: + broken_pixels |= mask return broken_pixels diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 3d8d1e09616..ed029cce674 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -257,7 +257,7 @@ def neighbor_average_maximum( neighbors = indices[indptr[pixel] : indptr[pixel + 1]] for neighbor in neighbors: - if broken_pixels[neighbor]: + if np.atleast_2d(broken_pixels)[ichannel, neighbor]: continue average += waveforms[ichannel][neighbor] @@ -580,14 +580,21 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels ) -> DL1CameraContainer: + # Make broken_pixels atleast 2D to allow multiple gains + broken_pix_2d = np.atleast_2d(broken_pixels) if self.pixel_fraction.tel[tel_id] == 1.0: # average over pixels then argmax over samples - peak_index = waveforms[:, ~broken_pixels].mean(axis=-2).argmax(axis=-1) + peak_index = ( + waveforms[~broken_pix_2d] + .reshape(waveforms.shape) + .mean(axis=-2) + .argmax(axis=-1) + ) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - brightest = np.argsort(waveforms.max(axis=-1))[:, ~broken_pixels][ - ..., -n_pixels: - ] + brightest = np.argsort(waveforms.max(axis=-1))[~broken_pix_2d].reshape( + waveforms.shape[-3], waveforms.shape[-2] + )[..., -n_pixels:] # average over brightest pixels then argmax over samples peak_index = ( diff --git a/src/ctapipe/image/invalid_pixels.py b/src/ctapipe/image/invalid_pixels.py index 10d22cb50a0..7bb71ff9aeb 100644 --- a/src/ctapipe/image/invalid_pixels.py +++ b/src/ctapipe/image/invalid_pixels.py @@ -1,6 +1,7 @@ """ Methods to interpolate broken pixels """ + from abc import ABCMeta, abstractmethod import numpy as np @@ -51,7 +52,7 @@ class NeighborAverage(InvalidPixelHandler): def __call__(self, tel_id, image, peak_time, pixel_mask): """Interpolate pixels in dl1 images and peak_times - Pixels to be interpolated are replaced by the average value their neighbors. + Pixels to be interpolated are replaced by the average value of their neighbors. Pixels where no valid neighbors are available are filled with zeros. diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 3f7f08ab303..2c4e389183a 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -284,8 +284,12 @@ def test_extract_around_peak_charge_expected(): @pytest.mark.parametrize("toymodels", camera_toymodels) def test_neighbor_average_peakpos(toymodels, request): waveforms, subarray, tel_id, _, _, _ = request.getfixturevalue(toymodels) + n_channels, n_pixels, _ = waveforms.shape neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + if n_channels == 1: + broken_pixels = np.zeros(n_pixels, dtype=bool) + else: + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) peak_pos = neighbor_average_maximum( waveforms, neighbors_indices=neighbors.indices, @@ -399,7 +403,11 @@ def test_extractors(Extractor, toymodels, request): true_time, ) = request.getfixturevalue(toymodels) extractor = Extractor(subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + if n_channels == 1: + broken_pixels = np.zeros(n_pixels, dtype=bool) + else: + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) if Extractor is TwoPassWindowSum and waveforms.shape[-3] != 1: with pytest.raises(AttributeError): @@ -433,7 +441,11 @@ def test_integration_correction_off(Extractor, toymodels, request): true_time, ) = request.getfixturevalue(toymodels) extractor = Extractor(subarray=subarray, apply_integration_correction=False) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + if n_channels == 1: + broken_pixels = np.zeros(n_pixels, dtype=bool) + else: + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) if Extractor is TwoPassWindowSum and waveforms.shape[-3] != 1: with pytest.raises(AttributeError): @@ -738,7 +750,7 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray): apply_integration_correction=False, ) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + broken_pixels = np.zeros((waveforms.shape[-2]), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert np.allclose(dl1.image[bright_pixels], 18) From dfff5667a8865621e63b7f7ed38e8cdde35738ef Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 8 Apr 2024 14:37:18 +0200 Subject: [PATCH 25/42] Minor Fix --- src/ctapipe/image/extractor.py | 2 +- src/ctapipe/image/tests/test_extractor.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index ed029cce674..8f01c8fd125 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -580,7 +580,7 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels ) -> DL1CameraContainer: - # Make broken_pixels atleast 2D to allow multiple gains + # Make broken_pixels at least 2D to allow multiple gains broken_pix_2d = np.atleast_2d(broken_pixels) if self.pixel_fraction.tel[tel_id] == 1.0: # average over pixels then argmax over samples diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 2c4e389183a..7349b94c287 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -76,8 +76,8 @@ def subarray_2_sst(prod5_sst, reference_location): "test array", tel_positions={1: np.zeros(3) * u.m, 2: np.zeros(3) * u.m}, tel_descriptions={ - 1: deepcopy(prod5_sst), - 2: deepcopy(prod5_sst), + 1: prod5_sst, + 2: prod5_sst, }, reference_location=reference_location, ) @@ -89,7 +89,7 @@ def subarray_1_lst(prod3_lst, reference_location): subarray = SubarrayDescription( "One LST", tel_positions={1: np.zeros(3) * u.m}, - tel_descriptions={1: deepcopy(prod3_lst)}, + tel_descriptions={1: prod3_lst}, reference_location=reference_location, ) return subarray @@ -100,7 +100,7 @@ def subarray_1_mst_fc(prod5_mst_flashcam, reference_location): subarray = SubarrayDescription( "One MST with FlashCam", tel_positions={1: np.zeros(3) * u.m}, - tel_descriptions={1: deepcopy(prod5_mst_flashcam)}, + tel_descriptions={1: prod5_mst_flashcam}, reference_location=reference_location, ) return subarray From fa561266771a9bd08861e1df88b54761bddad961 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 8 Apr 2024 16:04:31 +0200 Subject: [PATCH 26/42] Change broken_pixels to have always shape (n_channels, n_pixels) --- src/ctapipe/calib/camera/calibrator.py | 7 +-- src/ctapipe/image/extractor.py | 12 ++--- src/ctapipe/image/tests/test_extractor.py | 53 +++++++++++------------ 3 files changed, 37 insertions(+), 35 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index f9e5957f97c..e3881c2a1fe 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -21,8 +21,8 @@ __all__ = ["CameraCalibrator"] -def _get_invalid_pixels(n_pixels, pixel_status, selected_gain_channel): - broken_pixels = np.zeros(n_pixels, dtype=bool) +def _get_invalid_pixels(n_channels, n_pixels, pixel_status, selected_gain_channel): + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) index = np.arange(n_pixels) masks = ( pixel_status.hardware_failing_pixels, @@ -214,10 +214,11 @@ def _calibrate_dl1(self, event, tel_id): if self._check_dl0_empty(waveforms): return - _, n_pixels, n_samples = waveforms.shape + n_channels, n_pixels, n_samples = waveforms.shape selected_gain_channel = event.dl0.tel[tel_id].selected_gain_channel broken_pixels = _get_invalid_pixels( + n_channels, n_pixels, event.mon.tel[tel_id].pixel_status, selected_gain_channel, diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 8f01c8fd125..6a6cdec1c58 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -234,6 +234,7 @@ def neighbor_average_maximum( broken_pixels : ndarray Mask of broken pixels. Broken pixels are ignored in the sum over the neighbors. + Shape: (n_channels, n_pix) Returns ------- @@ -257,7 +258,7 @@ def neighbor_average_maximum( neighbors = indices[indptr[pixel] : indptr[pixel + 1]] for neighbor in neighbors: - if np.atleast_2d(broken_pixels)[ichannel, neighbor]: + if broken_pixels[ichannel, neighbor]: continue average += waveforms[ichannel][neighbor] @@ -421,6 +422,9 @@ def __call__( The channel selected in the gain selection, per pixel. Required in some cases to calculate the correct correction for the charge extraction. + broken_pixels : ndarray + Mask of broken pixels used for certain `Image Extractor` types. + Shape: (n_channels, n_pix) Returns ------- @@ -580,19 +584,17 @@ def _calculate_correction(self, tel_id): def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels ) -> DL1CameraContainer: - # Make broken_pixels at least 2D to allow multiple gains - broken_pix_2d = np.atleast_2d(broken_pixels) if self.pixel_fraction.tel[tel_id] == 1.0: # average over pixels then argmax over samples peak_index = ( - waveforms[~broken_pix_2d] + waveforms[~broken_pixels] .reshape(waveforms.shape) .mean(axis=-2) .argmax(axis=-1) ) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - brightest = np.argsort(waveforms.max(axis=-1))[~broken_pix_2d].reshape( + brightest = np.argsort(waveforms.max(axis=-1))[~broken_pixels].reshape( waveforms.shape[-3], waveforms.shape[-2] )[..., -n_pixels:] diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 7349b94c287..a49d2b1324b 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -286,10 +286,7 @@ def test_neighbor_average_peakpos(toymodels, request): waveforms, subarray, tel_id, _, _, _ = request.getfixturevalue(toymodels) n_channels, n_pixels, _ = waveforms.shape neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse - if n_channels == 1: - broken_pixels = np.zeros(n_pixels, dtype=bool) - else: - broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) peak_pos = neighbor_average_maximum( waveforms, neighbors_indices=neighbors.indices, @@ -404,10 +401,7 @@ def test_extractors(Extractor, toymodels, request): ) = request.getfixturevalue(toymodels) extractor = Extractor(subarray=subarray) n_channels, n_pixels, _ = waveforms.shape - if n_channels == 1: - broken_pixels = np.zeros(n_pixels, dtype=bool) - else: - broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) if Extractor is TwoPassWindowSum and waveforms.shape[-3] != 1: with pytest.raises(AttributeError): @@ -442,10 +436,7 @@ def test_integration_correction_off(Extractor, toymodels, request): ) = request.getfixturevalue(toymodels) extractor = Extractor(subarray=subarray, apply_integration_correction=False) n_channels, n_pixels, _ = waveforms.shape - if n_channels == 1: - broken_pixels = np.zeros(n_pixels, dtype=bool) - else: - broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) if Extractor is TwoPassWindowSum and waveforms.shape[-3] != 1: with pytest.raises(AttributeError): @@ -478,7 +469,8 @@ def test_fixed_window_sum(toymodel): ) = toymodel extractor = FixedWindowSum(subarray=subarray, peak_index=47) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid assert_allclose(dl1.image, true_charge, rtol=0.1) @@ -502,7 +494,8 @@ def test_sliding_window_max_sum(toymodels, request): true_time, ) = request.getfixturevalue(toymodels) extractor = SlidingWindowMaxSum(subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid if dl1.image.ndim == 1: @@ -525,7 +518,8 @@ def test_neighbor_peak_window_sum_local_weight(toymodel): ) = toymodel extractor = NeighborPeakWindowSum(subarray=subarray, local_weight=4) assert extractor.local_weight.tel[tel_id] == 4 - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) @@ -602,7 +596,8 @@ def test_Two_pass_window_sum_no_noise(subarray_1_lst): # Test only the 1st pass extractor.disable_second_pass = True - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1_pass1 = extractor(waveforms, 1, selected_gain_channel, broken_pixels) assert_allclose( dl1_pass1.image[signal_pixels & integration_window_inside], @@ -656,7 +651,8 @@ def test_waveform_extractor_factory(toymodel): ) = toymodel extractor = ImageExtractor.from_name("LocalPeakWindowSum", subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert_allclose(dl1.peak_time, true_time, rtol=0.1) @@ -715,7 +711,8 @@ def test_dtype(Extractor, subarray): waveforms = np.ones((1, n_pixels, 50), dtype="float64") extractor = Extractor(subarray=subarray) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.image.dtype == np.float32 assert dl1.peak_time.dtype == np.float32 @@ -750,7 +747,8 @@ def test_global_peak_window_sum_with_pixel_fraction(subarray): apply_integration_correction=False, ) - broken_pixels = np.zeros((waveforms.shape[-2]), dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert np.allclose(dl1.image[bright_pixels], 18) @@ -764,7 +762,8 @@ def test_adaptive_centroid(toymodel_mst_fc): waveforms, subarray, tel_id, _, _, _ = toymodel_mst_fc neighbors = subarray.tel[tel_id].camera.geometry.neighbor_matrix_sparse - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) trig_time = np.argmax(waveforms, axis=-1) peak_time = adaptive_centroid( @@ -853,8 +852,10 @@ def test_FC_time(toymodel_mst_fc_time): mask, ) = toymodel_mst_fc_time + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) + extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid @@ -862,7 +863,6 @@ def test_FC_time(toymodel_mst_fc_time): extractor = FlashCamExtractor( subarray=subarray, upsampling=1, leading_edge_timing=True ) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid @@ -870,13 +870,11 @@ def test_FC_time(toymodel_mst_fc_time): extractor = FlashCamExtractor( subarray=subarray, upsampling=1, leading_edge_timing=False ) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=False) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.peak_time[mask], true_time[mask], rtol=0.1) assert dl1.is_valid @@ -893,7 +891,8 @@ def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): _, ) = toymodel_mst_fc extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) - broken_pixels = np.zeros(waveforms.shape[-2], dtype=bool) + n_channels, n_pixels, _ = waveforms.shape + broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert_allclose(dl1.image, true_charge, rtol=0.1) assert dl1.is_valid @@ -914,13 +913,13 @@ def is_flashcam(tel_id): selected_gain_channel = np.zeros(waveforms.shape[-2], dtype=np.int64) broken_pixels = event.mon.tel[ tel_id - ].pixel_status.hardware_failing_pixels[0] + ].pixel_status.hardware_failing_pixels dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) assert dl1.is_valid bright_pixels = ( - (true_charge > 30) & (true_charge < 3000) & (~broken_pixels) + (true_charge > 30) & (true_charge < 3000) & (~broken_pixels[0]) ) assert_allclose( dl1.image[bright_pixels], true_charge[bright_pixels], rtol=0.35 From 82a956a9099eae21e55dd1dd444ff9cb5b79f1cf Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Mon, 8 Apr 2024 17:38:08 +0200 Subject: [PATCH 27/42] Adapt InvalidPixelHandler to work with 2 gain channels --- src/ctapipe/calib/camera/flatfield.py | 4 +- src/ctapipe/calib/camera/pedestals.py | 4 +- src/ctapipe/image/invalid_pixels.py | 58 ++++++++++--------- src/ctapipe/image/reducer.py | 4 +- .../image/tests/test_invalid_pixels.py | 24 ++++---- 5 files changed, 53 insertions(+), 41 deletions(-) diff --git a/src/ctapipe/calib/camera/flatfield.py b/src/ctapipe/calib/camera/flatfield.py index ecd805c2d3c..9e9a008e0c0 100644 --- a/src/ctapipe/calib/camera/flatfield.py +++ b/src/ctapipe/calib/camera/flatfield.py @@ -182,9 +182,11 @@ def _extract_charge(self, event) -> DL1CameraContainer: """ waveforms = event.r1.tel[self.tel_id].waveform + n_channels, n_pixels, _ = waveforms.shape selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel broken_pixels = _get_invalid_pixels( - n_pixels=waveforms.shape[-2], + n_channels=n_channels, + n_pixels=n_pixels, pixel_status=event.mon.tel[self.tel_id].pixel_status, selected_gain_channel=selected_gain_channel, ) diff --git a/src/ctapipe/calib/camera/pedestals.py b/src/ctapipe/calib/camera/pedestals.py index 094ce7645cb..383719d9654 100644 --- a/src/ctapipe/calib/camera/pedestals.py +++ b/src/ctapipe/calib/camera/pedestals.py @@ -211,9 +211,11 @@ def _extract_charge(self, event) -> DL1CameraContainer: DL1CameraContainer """ waveforms = event.r1.tel[self.tel_id].waveform + n_channels, n_pixels, _ = waveforms.shape selected_gain_channel = event.r1.tel[self.tel_id].selected_gain_channel broken_pixels = _get_invalid_pixels( - n_pixels=waveforms.shape[-2], + n_channels=n_channels, + n_pixels=n_pixels, pixel_status=event.mon.tel[self.tel_id].pixel_status, selected_gain_channel=selected_gain_channel, ) diff --git a/src/ctapipe/image/invalid_pixels.py b/src/ctapipe/image/invalid_pixels.py index 7bb71ff9aeb..6d427e890ab 100644 --- a/src/ctapipe/image/invalid_pixels.py +++ b/src/ctapipe/image/invalid_pixels.py @@ -38,6 +38,7 @@ def __call__( Array of pixel peak_time values pixel_mask : np.ndarray Boolean mask of the pixels to be interpolated + Shape: (n_channels, n_pixels) Returns ------- @@ -66,6 +67,7 @@ def __call__(self, tel_id, image, peak_time, pixel_mask): Array of pixel peak_time values pixel_mask : np.ndarray Boolean mask of the pixels to be interpolated + Shape: (n_channels, n_pixels) Returns ------- @@ -76,31 +78,35 @@ def __call__(self, tel_id, image, peak_time, pixel_mask): """ geometry = self.subarray.tel[tel_id].camera.geometry - n_interpolated = np.count_nonzero(pixel_mask) - if n_interpolated == 0: + n_interpolated = np.count_nonzero(pixel_mask, axis=-1, keepdims=True) + if (n_interpolated == 0).all(): return image, peak_time - # exclude to-be-interpolated pixels from neighbors - neighbors = geometry.neighbor_matrix[pixel_mask] & ~pixel_mask - - index, neighbor = np.nonzero(neighbors) - image_sum = np.zeros(n_interpolated, dtype=image.dtype) - count = np.zeros(n_interpolated, dtype=int) - peak_time_sum = np.zeros(n_interpolated, dtype=peak_time.dtype) - - # calculate average of image and peak_time - np.add.at(count, index, 1) - np.add.at(image_sum, index, image[neighbor]) - np.add.at(peak_time_sum, index, peak_time[neighbor]) - - valid = count > 0 - np.divide(image_sum, count, out=image_sum, where=valid) - np.divide(peak_time_sum, count, out=peak_time_sum, where=valid) - - peak_time = peak_time.copy() - peak_time[pixel_mask] = peak_time_sum - - image = image.copy() - image[pixel_mask] = image_sum - - return image, peak_time + image = np.atleast_2d(image) + peak_time = np.atleast_2d(peak_time) + for ichannel in range(image.shape[-2]): + if n_interpolated[ichannel] == 0: + continue + # exclude to-be-interpolated pixels from neighbors + neighbors = ( + geometry.neighbor_matrix[pixel_mask[ichannel]] & ~pixel_mask[ichannel] + ) + + index, neighbor = np.nonzero(neighbors) + image_sum = np.zeros(n_interpolated[ichannel], dtype=image.dtype) + count = np.zeros(n_interpolated[ichannel], dtype=int) + peak_time_sum = np.zeros(n_interpolated[ichannel], dtype=peak_time.dtype) + + # calculate average of image and peak_time + np.add.at(count, index, 1) + np.add.at(image_sum, index, image[ichannel, neighbor]) + np.add.at(peak_time_sum, index, peak_time[ichannel, neighbor]) + + valid = count > 0 + np.divide(image_sum, count, out=image_sum, where=valid) + np.divide(peak_time_sum, count, out=peak_time_sum, where=valid) + + peak_time[ichannel, pixel_mask[ichannel]] = peak_time_sum + image[ichannel, pixel_mask[ichannel]] = image_sum + + return np.squeeze(image), np.squeeze(peak_time) diff --git a/src/ctapipe/image/reducer.py b/src/ctapipe/image/reducer.py index c79f1602216..535aca63355 100644 --- a/src/ctapipe/image/reducer.py +++ b/src/ctapipe/image/reducer.py @@ -183,7 +183,9 @@ def select_pixels(self, waveforms, tel_id=None, selected_gain_channel=None): # Pulse-integrate waveforms extractor = self.image_extractors[self.image_extractor_type.tel[tel_id]] # do not treat broken pixels in data volume reduction - broken_pixels = np.zeros(camera_geom.n_pixels, dtype=bool) + broken_pixels = np.zeros( + (waveforms.shape[-3], camera_geom.n_pixels), dtype=bool + ) dl1: DL1CameraContainer = extractor( waveforms, tel_id=tel_id, diff --git a/src/ctapipe/image/tests/test_invalid_pixels.py b/src/ctapipe/image/tests/test_invalid_pixels.py index 280c84f886f..43a6c1312cf 100644 --- a/src/ctapipe/image/tests/test_invalid_pixels.py +++ b/src/ctapipe/image/tests/test_invalid_pixels.py @@ -12,16 +12,16 @@ def test_neighbor_average(prod5_gamma_simtel_path): neighbor_average = NeighborAverage(subarray) - broken_pixels = np.zeros(geometry.n_pixels, dtype=bool) + broken_pixels = np.zeros((2, geometry.n_pixels), dtype=bool) # one border pixel, one non-border pixel border = geometry.get_border_pixel_mask(1) - broken_pixels[np.nonzero(border)[0][0]] = True - broken_pixels[np.nonzero(~border)[0][0]] = True + broken_pixels[:, np.nonzero(border)[0][0]] = True + broken_pixels[:, np.nonzero(~border)[0][0]] = True - image = np.ones(geometry.n_pixels) + image = np.ones((2, geometry.n_pixels)) image[broken_pixels] = 9999 - peak_time = np.full(geometry.n_pixels, 20.0) + peak_time = np.full((2, geometry.n_pixels), 20.0) peak_time[broken_pixels] = -1 interpolated_image, interpolated_peaktime = neighbor_average( @@ -42,20 +42,20 @@ def test_negative_image(prod5_gamma_simtel_path): neighbor_average = NeighborAverage(subarray) - broken_pixels = np.zeros(geometry.n_pixels, dtype=bool) + broken_pixels = np.zeros((1, geometry.n_pixels), dtype=bool) # one border pixel, one non-border pixel border = geometry.get_border_pixel_mask(1) - broken_pixels[np.nonzero(border)[0][0]] = True - broken_pixels[np.nonzero(~border)[0][0]] = True + broken_pixels[:, np.nonzero(border)[0][0]] = True + broken_pixels[:, np.nonzero(~border)[0][0]] = True image = np.full(geometry.n_pixels, -10.0) - image[broken_pixels] = 9999 + image[broken_pixels[0]] = 9999 peak_time = np.full(geometry.n_pixels, 20.0) - peak_time[broken_pixels] = -1 + peak_time[broken_pixels[0]] = -1 interpolated_image, interpolated_peaktime = neighbor_average( tel_id=1, image=image, peak_time=peak_time, pixel_mask=broken_pixels ) - assert np.allclose(interpolated_image[broken_pixels], -10) - assert np.allclose(interpolated_peaktime[broken_pixels], 20.0) + assert np.allclose(interpolated_image[broken_pixels[0]], -10) + assert np.allclose(interpolated_peaktime[broken_pixels[0]], 20.0) From 887eb6db7b1530c4415180c8dfd76f6efc23fb78 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 9 Apr 2024 15:17:28 +0200 Subject: [PATCH 28/42] Fix reference --- src/ctapipe/image/extractor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 6a6cdec1c58..56a739013f6 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -423,7 +423,7 @@ def __call__( some cases to calculate the correct correction for the charge extraction. broken_pixels : ndarray - Mask of broken pixels used for certain `Image Extractor` types. + Mask of broken pixels used for certain `ImageExtractor` types. Shape: (n_channels, n_pix) Returns From 9ebbdd589fa032d1bf5d961382f02a1c48291c46 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 9 Apr 2024 14:25:16 +0200 Subject: [PATCH 29/42] Change charge and peak_time calibration broadcasting - I actually cant remember why i used transpose to broadcast the calib coeffs to the charges and peak times. I think it was due to a usecase which is actually not needed. - I've added a test especially for 2 gain data to check that the calibrator also works this way --- src/ctapipe/calib/camera/calibrator.py | 25 ++-- .../calib/camera/tests/test_calibrator.py | 117 +++++++++++++++++- 2 files changed, 126 insertions(+), 16 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index e3881c2a1fe..8648066a7a0 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -225,16 +225,15 @@ def _calibrate_dl1(self, event, tel_id): ) dl1_calib = event.calibration.tel[tel_id].dl1 - time_shift = event.calibration.tel[tel_id].dl1.time_shift + ped_offset = dl1_calib.pedestal_offset + time_shift = dl1_calib.time_shift readout = self.subarray.tel[tel_id].camera.readout # subtract any remaining pedestal before extraction - if dl1_calib.pedestal_offset is not None: + if ped_offset is not None: # this copies intentionally, we don't want to modify the dl0 data - # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels,) - waveforms = ( - waveforms - np.atleast_2d(dl1_calib.pedestal_offset)[..., np.newaxis] - ) + # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels) + waveforms = waveforms - np.atleast_2d(ped_offset)[..., np.newaxis] if n_samples == 1: # To handle ASTRI and dst @@ -271,12 +270,10 @@ def _calibrate_dl1(self, event, tel_id): # correct non-integer remainder of the shift if given if self.apply_peak_time_shift.tel[tel_id] and time_shift is not None: - dl1.peak_time = (dl1.peak_time.T - remaining_shift).T + dl1.peak_time -= remaining_shift # Calibrate extracted charge - dl1.image = ( - dl1.image.T * (dl1_calib.relative_factor / dl1_calib.absolute_factor) - ).T + dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor # handle invalid pixels if self.invalid_pixel_handler is not None: @@ -319,7 +316,7 @@ def shift_waveforms(waveforms, time_shift_samples): ---------- waveforms: ndarray of shape (n_channels, n_pixels, n_samples) The waveforms to shift - time_shift_samples: ndarray of shape (n_pixels, ) + time_shift_samples: ndarray The shift to apply in units of samples. Waveforms are shifted to the left by the smallest integer that minimizes inter-pixel differences. @@ -328,11 +325,11 @@ def shift_waveforms(waveforms, time_shift_samples): ------- shifted_waveforms: ndarray of shape (n_channels, n_pixels, n_samples) The shifted waveforms - remaining_shift: ndarray of shape (n_pixels, ) + remaining_shift: ndarray The remaining shift after applying the integer shift to the waveforms. """ - mean_shift = time_shift_samples.mean(-1) - integer_shift = np.round((time_shift_samples.T - mean_shift).T).astype("int16") + mean_shift = time_shift_samples.mean(axis=-1, keepdims=True) + integer_shift = np.round(time_shift_samples - mean_shift).astype("int16") remaining_shift = time_shift_samples - integer_shift shifted_waveforms = _shift_waveforms_by_integer(waveforms, integer_shift) return shifted_waveforms, remaining_shift diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index f768ba69316..c28aa8783a0 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -130,7 +130,7 @@ def test_check_dl0_empty(example_event, example_subarray): assert (event.dl1.tel[tel_id].image == 2).all() -def test_dl1_charge_calib(example_subarray): +def test_dl1_charge_calib_1_gain(example_subarray): # copy because we mutate the camera, should not affect other tests subarray = deepcopy(example_subarray) camera = subarray.tel[1].camera @@ -178,7 +178,7 @@ def test_dl1_charge_calib(example_subarray): subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) ) calibrator(event) - np.testing.assert_allclose(event.dl1.tel[tel_id].image, y.sum(2)[0, :], rtol=1e-4) + np.testing.assert_allclose(event.dl1.tel[tel_id].image, y.sum(-1)[0], rtol=1e-4) event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal event.calibration.tel[tel_id].dl1.absolute_factor = absolute @@ -237,6 +237,119 @@ def test_dl1_charge_calib(example_subarray): assert not np.allclose(event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1) +def test_dl1_charge_calib_2_gain(example_subarray): + # copy because we mutate the camera, should not affect other tests + subarray = deepcopy(example_subarray) + camera = subarray.tel[1].camera + # test with a sampling_rate different than 1 to + # test if we handle time vs. slices correctly + sampling_rate = 2 + camera.readout.sampling_rate = sampling_rate * u.GHz + + n_pixels = camera.geometry.n_pixels + n_samples = 96 + mid = n_samples // 2 + pulse_sigma = 6 + random = np.random.default_rng(1) + x = np.arange(n_samples) + + # Randomize times and create pulses + time_offset_2_gain = random.uniform(-10, +10, (2, n_pixels)) + y_2_gain = np.ones((2, n_pixels, n_samples)) + y_2_gain[:] = norm.pdf( + x, mid + time_offset_2_gain[..., np.newaxis], pulse_sigma + ).astype("float32") + + camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma) + camera.readout.reference_pulse_shape = np.repeat( + camera.readout.reference_pulse_shape[np.newaxis, :], 2, axis=0 + ) + camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate + + # Define absolute calibration coefficients + absolute_2_gain = random.uniform(100, 1000, (2, n_pixels)).astype("float32") + y_2_gain *= absolute_2_gain[..., np.newaxis] + + # Define relative coefficients + relative_2_gain = random.normal(1, 0.01, (2, n_pixels)) + y_2_gain /= relative_2_gain[..., np.newaxis] + + # Define pedestal + pedestal_2_gain = random.uniform(-4, 4, (2, n_pixels)) + y_2_gain += pedestal_2_gain[..., np.newaxis] + + event = ArrayEventContainer() + tel_id = list(subarray.tel.keys())[0] + event.dl0.tel[tel_id].waveform = y_2_gain + event.dl0.tel[tel_id].selected_gain_channel = None + event.r1.tel[tel_id].selected_gain_channel = None + + # Test default + calibrator = CameraCalibrator( + subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) + ) + calibrator(event) + np.testing.assert_allclose( + event.dl1.tel[tel_id].image[0], y_2_gain.sum(2)[0, :], rtol=1e-4 + ) + + event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal_2_gain + event.calibration.tel[tel_id].dl1.absolute_factor = absolute_2_gain + event.calibration.tel[tel_id].dl1.relative_factor = relative_2_gain + + # Test without timing corrections + calibrator(event) + dl1 = event.dl1.tel[tel_id] + np.testing.assert_allclose(dl1.image[0], 1, rtol=1e-5) + expected_peak_time = (mid + time_offset_2_gain) / sampling_rate + np.testing.assert_allclose(dl1.peak_time, expected_peak_time, rtol=1e-5) + + # test with timing corrections + event.calibration.tel[tel_id].dl1.time_shift = time_offset_2_gain / sampling_rate + calibrator(event) + + # more rtol since shifting might lead to reduced integral + np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-5) + np.testing.assert_allclose( + event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + ) + + # test not applying time shifts + # now we should be back to the result without setting time shift + calibrator.apply_peak_time_shift = False + calibrator.apply_waveform_time_shift = False + calibrator(event) + + np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-4) + np.testing.assert_allclose( + event.dl1.tel[tel_id].peak_time, expected_peak_time, atol=1 + ) + + # We now use GlobalPeakWindowSum to see the effect of missing charge + # due to not correcting time offsets. + calibrator = CameraCalibrator( + subarray=subarray, + image_extractor=GlobalPeakWindowSum(subarray=subarray), + apply_waveform_time_shift=True, + ) + calibrator(event) + # test with timing corrections, should work + # higher rtol because we cannot shift perfectly + np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=0.01) + np.testing.assert_allclose( + event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + ) + + # test deactivating timing corrections + calibrator.apply_waveform_time_shift = False + calibrator(event) + + # make sure we chose an example where the time shifts matter + # charges should be quite off due to summing around global shift + assert not np.allclose(event.dl1.tel[tel_id].image, 1, rtol=0.1) + assert not np.allclose(event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1) + + def test_shift_waveforms(): from ctapipe.calib.camera.calibrator import shift_waveforms From 732c5cf0b71876c79dc5a6994d6c28f6ea1910f6 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 9 Apr 2024 14:31:49 +0200 Subject: [PATCH 30/42] Adapting changelog --- docs/changes/2529.datamodel.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/changes/2529.datamodel.rst b/docs/changes/2529.datamodel.rst index 06fb49fb081..eb8231c348a 100644 --- a/docs/changes/2529.datamodel.rst +++ b/docs/changes/2529.datamodel.rst @@ -1,6 +1,6 @@ Change R1- and DL0-waveforms datamodel shape from (n_pixels, n_samples) -to be always (n_channels, n_pixels, n_samples). This breaks reading in single -gain data of older versions (< v6.0.0) with the ``HDF5EventSource``. +to be always (n_channels, n_pixels, n_samples). ``HDF5EventSource`` was adjusted +accordingly to support also older datamodel versions. Re-introduce also the possibility of running ``ImageExtractor``\s on data consisting of multiple gain channels. From a85cc583d72c981cf3fe69d37eb436f0402b898a Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 9 Apr 2024 15:06:46 +0200 Subject: [PATCH 31/42] Add test for 2 gains in InvaliPixelHandler --- .../image/tests/test_invalid_pixels.py | 22 ++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/tests/test_invalid_pixels.py b/src/ctapipe/image/tests/test_invalid_pixels.py index 43a6c1312cf..0aaed1dd7c8 100644 --- a/src/ctapipe/image/tests/test_invalid_pixels.py +++ b/src/ctapipe/image/tests/test_invalid_pixels.py @@ -11,11 +11,31 @@ def test_neighbor_average(prod5_gamma_simtel_path): geometry = subarray.tel[1].camera.geometry neighbor_average = NeighborAverage(subarray) + border = geometry.get_border_pixel_mask(1) + + ### 1 GAIN ### + broken_pixels = np.zeros((1, geometry.n_pixels), dtype=bool) + + # one border pixel, one non-border pixel + broken_pixels[:, np.nonzero(border)[0][0]] = True + broken_pixels[:, np.nonzero(~border)[0][0]] = True + + image = np.ones((geometry.n_pixels)) + image[broken_pixels[0]] = 9999 + peak_time = np.full(geometry.n_pixels, 20.0) + peak_time[broken_pixels[0]] = -1 + interpolated_image, interpolated_peaktime = neighbor_average( + tel_id=1, image=image, peak_time=peak_time, pixel_mask=broken_pixels + ) + + assert np.allclose(interpolated_image[broken_pixels[0]], 1.0) + assert np.allclose(interpolated_peaktime[broken_pixels[0]], 20.0) + + ### 2 GAIN ### broken_pixels = np.zeros((2, geometry.n_pixels), dtype=bool) # one border pixel, one non-border pixel - border = geometry.get_border_pixel_mask(1) broken_pixels[:, np.nonzero(border)[0][0]] = True broken_pixels[:, np.nonzero(~border)[0][0]] = True From bd59741684ab0f0cc48bdb9223960e2f556b1816 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 9 Apr 2024 17:25:52 +0200 Subject: [PATCH 32/42] Add select_gain option to SimTelEventSource --- src/ctapipe/io/simteleventsource.py | 31 +++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py index 5853f2b9277..2969396a1ed 100644 --- a/src/ctapipe/io/simteleventsource.py +++ b/src/ctapipe/io/simteleventsource.py @@ -284,8 +284,13 @@ def apply_simtel_r1_calibration( r1_waveforms = (r0_waveforms - ped) * gain + calib_shift - selected_gain_channel = gain_selector(r0_waveforms) - r1_waveforms = r1_waveforms[np.newaxis, selected_gain_channel, np.arange(n_pixels)] + if gain_selector is not None: + selected_gain_channel = gain_selector(r0_waveforms) + r1_waveforms = r1_waveforms[ + np.newaxis, selected_gain_channel, np.arange(n_pixels) + ] + else: + selected_gain_channel = None return r1_waveforms, selected_gain_channel @@ -487,6 +492,15 @@ class SimTelEventSource(EventSource): help="Use the given obs_id instead of the run number from sim_telarray", ).tag(config=True) + select_gain = Bool( + default_value=None, + allow_none=True, + help=( + "Whether to perform gain selection. The default (None) means only" + " select gain of physics events, not of calibration events." + ), + ).tag(config=True) + def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs): """ EventSource for simtelarray files using the pyeventio library. @@ -626,7 +640,7 @@ def prepare_subarray_info(self, telescope_descriptions, header): # TODO: switch to warning or even an exception once # we start relying on this. - self.log.info( + self.log.debug( "Could not determine telescope from sim_telarray metadata," " guessing using builtin lookup-table: %d: %s", tel_id, @@ -835,11 +849,20 @@ def _generate_events(self): mon.calibration.pedestal_per_sample = pedestal mon.pixel_status = self._fill_mon_pixels_status(tel_id) + select_gain = self.select_gain is True or ( + self.select_gain is None + and trigger.event_type is EventType.SUBARRAY + ) + if select_gain: + gain_selector = self.gain_selector + else: + gain_selector = None + r1_waveform, selected_gain_channel = apply_simtel_r1_calibration( adc_samples, pedestal, dc_to_pe, - self.gain_selector, + gain_selector, self.calib_scale, self.calib_shift, ) From 25c3bce211094917fc825c7078a6022961d01eb3 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 9 Apr 2024 17:26:34 +0200 Subject: [PATCH 33/42] Calibration coefficients always have shape (n_channels, n_pixels) Co-Authored-By: Jonas Hackfeld --- src/ctapipe/calib/camera/calibrator.py | 18 +++- src/ctapipe/calib/camera/pedestals.py | 6 +- .../calib/camera/tests/test_calibrator.py | 47 +++++++--- .../calib/camera/tests/test_flatfield.py | 1 - .../calib/camera/tests/test_pedestals.py | 1 - src/ctapipe/containers.py | 21 ++--- src/ctapipe/image/extractor.py | 86 ++++++++++++------- src/ctapipe/image/invalid_pixels.py | 26 +++--- src/ctapipe/io/simteleventsource.py | 4 +- 9 files changed, 136 insertions(+), 74 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 8648066a7a0..447189e5d36 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -23,6 +23,7 @@ def _get_invalid_pixels(n_channels, n_pixels, pixel_status, selected_gain_channel): broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) + index = np.arange(n_pixels) masks = ( pixel_status.hardware_failing_pixels, @@ -223,6 +224,7 @@ def _calibrate_dl1(self, event, tel_id): event.mon.tel[tel_id].pixel_status, selected_gain_channel, ) + pixel_index = np.arange(n_pixels) dl1_calib = event.calibration.tel[tel_id].dl1 ped_offset = dl1_calib.pedestal_offset @@ -250,6 +252,9 @@ def _calibrate_dl1(self, event, tel_id): else: # shift waveforms if time_shift calibration is available if time_shift is not None: + if selected_gain_channel is not None: + time_shift = time_shift[selected_gain_channel, pixel_index] + if self.apply_waveform_time_shift.tel[tel_id]: sampling_rate = readout.sampling_rate.to_value(u.GHz) time_shift_samples = time_shift * sampling_rate @@ -273,7 +278,18 @@ def _calibrate_dl1(self, event, tel_id): dl1.peak_time -= remaining_shift # Calibrate extracted charge - dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor + if ( + dl1_calib.relative_factor is not None + and dl1_calib.absolute_factor is not None + ): + if selected_gain_channel is None: + dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor + else: + corr = ( + dl1_calib.relative_factor[selected_gain_channel, pixel_index] + / dl1_calib.absolute_factor[selected_gain_channel, pixel_index] + ) + dl1.image *= corr # handle invalid pixels if self.invalid_pixel_handler is not None: diff --git a/src/ctapipe/calib/camera/pedestals.py b/src/ctapipe/calib/camera/pedestals.py index 383719d9654..ba8ee336838 100644 --- a/src/ctapipe/calib/camera/pedestals.py +++ b/src/ctapipe/calib/camera/pedestals.py @@ -219,7 +219,6 @@ def _extract_charge(self, event) -> DL1CameraContainer: pixel_status=event.mon.tel[self.tel_id].pixel_status, selected_gain_channel=selected_gain_channel, ) - # Extract charge and time if self.extractor: return self.extractor( @@ -249,10 +248,7 @@ def calculate_pedestals(self, event): # real data trigger_time = event.trigger.time - if event.meta["origin"] != "hessio": - pixel_mask = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels - else: # patches for MC data - pixel_mask = np.zeros(waveform.shape[1], dtype=bool) + pixel_mask = event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels if self.n_events_seen == 0: self.time_start = trigger_time diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index c28aa8783a0..2579fbf5f10 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -139,6 +139,7 @@ def test_dl1_charge_calib_1_gain(example_subarray): sampling_rate = 2 camera.readout.sampling_rate = sampling_rate * u.GHz + n_channels = camera.readout.n_channels n_pixels = camera.geometry.n_pixels n_samples = 96 mid = n_samples // 2 @@ -147,31 +148,29 @@ def test_dl1_charge_calib_1_gain(example_subarray): x = np.arange(n_samples) # Randomize times and create pulses - time_offset = random.uniform(-10, +10, n_pixels) - y = norm.pdf(x, mid + time_offset[np.newaxis, :, np.newaxis], pulse_sigma).astype( - "float32" - ) + time_offset = random.uniform(-10, +10, (n_channels, n_pixels)) + y = norm.pdf(x, mid + time_offset[:, :, np.newaxis], pulse_sigma).astype("float32") camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma)[np.newaxis, :] camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate # Define absolute calibration coefficients - absolute = random.uniform(100, 1000, n_pixels).astype("float32") - y *= absolute[np.newaxis, :, np.newaxis] + absolute = random.uniform(100, 1000, (n_channels, n_pixels)).astype("float32") + y *= absolute[..., np.newaxis] # Define relative coefficients - relative = random.normal(1, 0.01, n_pixels) - y /= relative[np.newaxis, :, np.newaxis] + relative = random.normal(1, 0.01, (n_channels, n_pixels)) + y /= relative[..., np.newaxis] # Define pedestal - pedestal = random.uniform(-4, 4, n_pixels) - y += pedestal[np.newaxis, :, np.newaxis] + pedestal = random.uniform(-4, 4, (n_channels, n_pixels)) + y += pedestal[..., np.newaxis] event = ArrayEventContainer() tel_id = list(subarray.tel.keys())[0] event.dl0.tel[tel_id].waveform = y - event.dl0.tel[tel_id].selected_gain_channel = np.zeros(len(y), dtype=int) - event.r1.tel[tel_id].selected_gain_channel = np.zeros(len(y), dtype=int) + event.dl0.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=int) + event.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=int) # Test default calibrator = CameraCalibrator( @@ -188,7 +187,9 @@ def test_dl1_charge_calib_1_gain(example_subarray): calibrator(event) dl1 = event.dl1.tel[tel_id] np.testing.assert_allclose(dl1.image, 1, rtol=1e-5) - expected_peak_time = (mid + time_offset) / sampling_rate + + # we use only the first gain + expected_peak_time = (mid + time_offset[0]) / sampling_rate np.testing.assert_allclose(dl1.peak_time, expected_peak_time, rtol=1e-5) # test with timing corrections @@ -428,3 +429,23 @@ def test_invalid_pixels(example_event, example_subarray): calibrator(event) assert event.dl1.tel[tel_id].image[0] == 9999 assert event.dl1.tel[tel_id].peak_time[0] == 10.0 / sampling_rate + + +def test_no_gain_selection(prod5_gamma_simtel_path): + from ctapipe.io import SimTelEventSource + + with SimTelEventSource(prod5_gamma_simtel_path, select_gain=False) as source: + event = next(iter(source)) + + tel_id = next(iter(event.r1.tel)) + readout = source.subarray.tel[tel_id].camera.readout + + calibrator = CameraCalibrator(subarray=source.subarray) + calibrator(event) + + image = event.dl1.tel[tel_id].image + peak_time = event.dl1.tel[tel_id].peak_time + assert image.ndim == 2 + assert peak_time.ndim == 2 + assert image.shape == (readout.n_channels, readout.n_pixels) + assert peak_time.shape == (readout.n_channels, readout.n_pixels) diff --git a/src/ctapipe/calib/camera/tests/test_flatfield.py b/src/ctapipe/calib/camera/tests/test_flatfield.py index 3b5de6139a1..7c68f613af9 100644 --- a/src/ctapipe/calib/camera/tests/test_flatfield.py +++ b/src/ctapipe/calib/camera/tests/test_flatfield.py @@ -53,7 +53,6 @@ def test_flasherflatfieldcalculator(prod5_sst, reference_location): (n_gain, n_pixels), dtype=bool ) data.r1.tel[tel_id].waveform = np.zeros((n_gain, n_pixels, 40)) - data.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8) # flat-field signal put == delta function of height ff_level at sample 20 data.r1.tel[tel_id].waveform[:, :, 20] = ff_level diff --git a/src/ctapipe/calib/camera/tests/test_pedestals.py b/src/ctapipe/calib/camera/tests/test_pedestals.py index 7f17c9f1703..af37677d79f 100644 --- a/src/ctapipe/calib/camera/tests/test_pedestals.py +++ b/src/ctapipe/calib/camera/tests/test_pedestals.py @@ -49,7 +49,6 @@ def test_pedestal_integrator(prod5_sst, reference_location): (n_gain, n_pixels), dtype=bool ) data.r1.tel[tel_id].waveform = np.full((2, n_pixels, 40), ped_level) - data.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=np.uint8) while ped_calculator.n_events_seen < n_events: if ped_calculator.calculate_pedestals(data): diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py index 11a17107a71..9a5b39f4da8 100644 --- a/src/ctapipe/containers.py +++ b/src/ctapipe/containers.py @@ -534,23 +534,24 @@ class DL1CameraCalibrationContainer(Container): None, "Residual mean pedestal of the waveforms for each pixel." " This value is subtracted from the waveforms of each pixel before" - " the pulse extraction.", + " the pulse extraction. Shape: (n_channels, n_pixels)", ) absolute_factor = Field( - 1, - "Multiplicative coefficients for the absolute calibration of extracted charge into " - "physical units (e.g. photoelectrons or photons) for each pixel", + None, + "Multiplicative coefficients for the absolute calibration of extracted charge" + " into physical units (e.g. photoelectrons or photons) for each pixel." + " Shape: (n_channels, n_pixels)", ) relative_factor = Field( - 1, - "Multiplicative Coefficients for the relative correction between pixels to achieve a " - "uniform charge response (post absolute calibration) from a " - "uniform illumination.", + None, + "Multiplicative Coefficients for the relative correction between pixels to" + " achieve a uniform charge response (post absolute calibration) from a" + " uniform illumination. Shape: (n_channels, n_pixels)", ) time_shift = Field( None, - "Additive coefficients for the timing correction before charge extraction " - "for each pixel", + "Additive coefficients for the timing correction before charge extraction" + " for each pixel. Shape: (n_channels, n_pixels)", ) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 56a739013f6..8f2c4f07e9a 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -444,9 +444,13 @@ def __call__( charge, peak_time = extract_around_peak( waveforms, 0, waveforms.shape[-1], 0, self.sampling_rate_ghz[tel_id] ) - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) class FixedWindowSum(ImageExtractor): @@ -512,9 +516,13 @@ def __call__( if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) charge = _apply_correction(charge, correction, selected_gain_channel) - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) class GlobalPeakWindowSum(ImageExtractor): @@ -586,16 +594,15 @@ def __call__( ) -> DL1CameraContainer: if self.pixel_fraction.tel[tel_id] == 1.0: # average over pixels then argmax over samples - peak_index = ( - waveforms[~broken_pixels] - .reshape(waveforms.shape) - .mean(axis=-2) - .argmax(axis=-1) - ) + peak_index = waveforms.mean( + axis=-2, where=~broken_pixels[..., np.newaxis] + ).argmax(axis=-1) else: n_pixels = int(self.pixel_fraction.tel[tel_id] * waveforms.shape[-2]) - brightest = np.argsort(waveforms.max(axis=-1))[~broken_pixels].reshape( - waveforms.shape[-3], waveforms.shape[-2] + brightest = np.argsort( + waveforms.max( + axis=-1, where=~broken_pixels[..., np.newaxis], initial=-np.inf + ) )[..., -n_pixels:] # average over brightest pixels then argmax over samples @@ -613,9 +620,13 @@ def __call__( if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) charge = _apply_correction(charge, correction, selected_gain_channel) - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) class LocalPeakWindowSum(ImageExtractor): @@ -681,9 +692,13 @@ def __call__( if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) charge = _apply_correction(charge, correction, selected_gain_channel) - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) class SlidingWindowMaxSum(ImageExtractor): @@ -757,12 +772,17 @@ def __call__( charge, peak_time = extract_sliding_window( waveforms, self.window_width.tel[tel_id], self.sampling_rate_ghz[tel_id] ) + if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) charge = _apply_correction(charge, correction, selected_gain_channel) - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) class NeighborPeakWindowSum(ImageExtractor): @@ -838,12 +858,17 @@ def __call__( self.window_shift.tel[tel_id], self.sampling_rate_ghz[tel_id], ) + if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) charge = _apply_correction(charge, correction, selected_gain_channel) - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) class BaselineSubtractedNeighborPeakWindowSum(NeighborPeakWindowSum): @@ -1742,6 +1767,9 @@ def __call__( if shift != 0: peak_time -= shift - return DL1CameraContainer( - image=np.squeeze(charge), peak_time=np.squeeze(peak_time), is_valid=True - ) + # reduce dimensions for gain selected data to (n_pixels, ) + if selected_gain_channel is not None: + charge = charge[0] + peak_time = peak_time[0] + + return DL1CameraContainer(image=charge, peak_time=peak_time, is_valid=True) diff --git a/src/ctapipe/image/invalid_pixels.py b/src/ctapipe/image/invalid_pixels.py index 6d427e890ab..5b483f7895e 100644 --- a/src/ctapipe/image/invalid_pixels.py +++ b/src/ctapipe/image/invalid_pixels.py @@ -82,31 +82,35 @@ def __call__(self, tel_id, image, peak_time, pixel_mask): if (n_interpolated == 0).all(): return image, peak_time + input_1d = image.ndim == 1 image = np.atleast_2d(image) peak_time = np.atleast_2d(peak_time) - for ichannel in range(image.shape[-2]): - if n_interpolated[ichannel] == 0: + n_channels = image.shape[0] + for channel in range(n_channels): + if n_interpolated[channel] == 0: continue # exclude to-be-interpolated pixels from neighbors neighbors = ( - geometry.neighbor_matrix[pixel_mask[ichannel]] & ~pixel_mask[ichannel] + geometry.neighbor_matrix[pixel_mask[channel]] & ~pixel_mask[channel] ) index, neighbor = np.nonzero(neighbors) - image_sum = np.zeros(n_interpolated[ichannel], dtype=image.dtype) - count = np.zeros(n_interpolated[ichannel], dtype=int) - peak_time_sum = np.zeros(n_interpolated[ichannel], dtype=peak_time.dtype) + image_sum = np.zeros(n_interpolated[channel], dtype=image.dtype) + count = np.zeros(n_interpolated[channel], dtype=int) + peak_time_sum = np.zeros(n_interpolated[channel], dtype=peak_time.dtype) # calculate average of image and peak_time np.add.at(count, index, 1) - np.add.at(image_sum, index, image[ichannel, neighbor]) - np.add.at(peak_time_sum, index, peak_time[ichannel, neighbor]) + np.add.at(image_sum, index, image[channel, neighbor]) + np.add.at(peak_time_sum, index, peak_time[channel, neighbor]) valid = count > 0 np.divide(image_sum, count, out=image_sum, where=valid) np.divide(peak_time_sum, count, out=peak_time_sum, where=valid) - peak_time[ichannel, pixel_mask[ichannel]] = peak_time_sum - image[ichannel, pixel_mask[ichannel]] = image_sum + peak_time[channel, pixel_mask[channel]] = peak_time_sum + image[channel, pixel_mask[channel]] = image_sum - return np.squeeze(image), np.squeeze(peak_time) + if input_1d: + return image[0], peak_time[0] + return image, peak_time diff --git a/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py index 2969396a1ed..e4ff831dacb 100644 --- a/src/ctapipe/io/simteleventsource.py +++ b/src/ctapipe/io/simteleventsource.py @@ -880,10 +880,8 @@ def _generate_events(self): # get time_shift from laser calibration time_calib = array_event["laser_calibrations"][tel_id]["tm_calib"] - pix_index = np.arange(n_pixels) - dl1_calib = data.calibration.tel[tel_id].dl1 - dl1_calib.time_shift = time_calib[selected_gain_channel, pix_index] + dl1_calib.time_shift = time_calib yield data From 9fe0018fbe6d603b38617f453ed94a9f8ab10a43 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 10 Apr 2024 14:42:53 +0200 Subject: [PATCH 34/42] Test all tels in first event --- .../calib/camera/tests/test_calibrator.py | 25 +++++++++++-------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index 2579fbf5f10..f9fc4f04f56 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -437,15 +437,20 @@ def test_no_gain_selection(prod5_gamma_simtel_path): with SimTelEventSource(prod5_gamma_simtel_path, select_gain=False) as source: event = next(iter(source)) - tel_id = next(iter(event.r1.tel)) - readout = source.subarray.tel[tel_id].camera.readout + tested_n_channels = set() - calibrator = CameraCalibrator(subarray=source.subarray) - calibrator(event) + for tel_id in event.r1.tel: + readout = source.subarray.tel[tel_id].camera.readout + tested_n_channels.add(readout.n_channels) + + calibrator = CameraCalibrator(subarray=source.subarray) + calibrator(event) + + image = event.dl1.tel[tel_id].image + peak_time = event.dl1.tel[tel_id].peak_time + assert image.ndim == 2 + assert peak_time.ndim == 2 + assert image.shape == (readout.n_channels, readout.n_pixels) + assert peak_time.shape == (readout.n_channels, readout.n_pixels) - image = event.dl1.tel[tel_id].image - peak_time = event.dl1.tel[tel_id].peak_time - assert image.ndim == 2 - assert peak_time.ndim == 2 - assert image.shape == (readout.n_channels, readout.n_pixels) - assert peak_time.shape == (readout.n_channels, readout.n_pixels) + assert tested_n_channels == {1, 2} From 0a2e3098b3b66099c0348a7d94fa0df521384cce Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 10 Apr 2024 14:51:28 +0200 Subject: [PATCH 35/42] Remove unneeded differentiation between simtel and "real" data in flatfield --- src/ctapipe/calib/camera/flatfield.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/src/ctapipe/calib/camera/flatfield.py b/src/ctapipe/calib/camera/flatfield.py index 9e9a008e0c0..c3f74678aae 100644 --- a/src/ctapipe/calib/camera/flatfield.py +++ b/src/ctapipe/calib/camera/flatfield.py @@ -217,20 +217,15 @@ def calculate_relative_gain(self, event): if self.n_events_seen == self.sample_size: self.n_events_seen = 0 - # real data trigger_time = event.trigger.time - if event.meta["origin"] != "hessio": - hardware_or_pedestal_mask = np.logical_or( - event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels, - event.mon.tel[self.tel_id].pixel_status.pedestal_failing_pixels, - ) - pixel_mask = np.logical_or( - hardware_or_pedestal_mask, - event.mon.tel[self.tel_id].pixel_status.flatfield_failing_pixels, - ) - - else: # patches for MC data - pixel_mask = np.zeros(waveform.shape[1], dtype=bool) + hardware_or_pedestal_mask = np.logical_or( + event.mon.tel[self.tel_id].pixel_status.hardware_failing_pixels, + event.mon.tel[self.tel_id].pixel_status.pedestal_failing_pixels, + ) + pixel_mask = np.logical_or( + hardware_or_pedestal_mask, + event.mon.tel[self.tel_id].pixel_status.flatfield_failing_pixels, + ) if self.n_events_seen == 0: self.time_start = trigger_time From 0402d1bee977de15e4147be1c99662532035bef4 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 11 Apr 2024 17:27:05 +0200 Subject: [PATCH 36/42] Improve readability --- src/ctapipe/calib/camera/calibrator.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 447189e5d36..3678e61fb33 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -227,15 +227,14 @@ def _calibrate_dl1(self, event, tel_id): pixel_index = np.arange(n_pixels) dl1_calib = event.calibration.tel[tel_id].dl1 - ped_offset = dl1_calib.pedestal_offset - time_shift = dl1_calib.time_shift readout = self.subarray.tel[tel_id].camera.readout # subtract any remaining pedestal before extraction - if ped_offset is not None: + if dl1_calib.pedestal_offset is not None: # this copies intentionally, we don't want to modify the dl0 data # waveforms have shape (n_channels, n_pixel, n_samples), pedestals (n_pixels) - waveforms = waveforms - np.atleast_2d(ped_offset)[..., np.newaxis] + waveforms = waveforms.copy() + waveforms -= np.atleast_2d(dl1_calib.pedestal_offset)[..., np.newaxis] if n_samples == 1: # To handle ASTRI and dst @@ -251,6 +250,8 @@ def _calibrate_dl1(self, event, tel_id): ) else: # shift waveforms if time_shift calibration is available + time_shift = dl1_calib.time_shift + remaining_shift = None if time_shift is not None: if selected_gain_channel is not None: time_shift = time_shift[selected_gain_channel, pixel_index] @@ -274,7 +275,7 @@ def _calibrate_dl1(self, event, tel_id): ) # correct non-integer remainder of the shift if given - if self.apply_peak_time_shift.tel[tel_id] and time_shift is not None: + if self.apply_peak_time_shift.tel[tel_id] and remaining_shift is not None: dl1.peak_time -= remaining_shift # Calibrate extracted charge From 4ad3edb4d7b33223e7202b6d2f059e7560a072c7 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 11 Apr 2024 17:29:48 +0200 Subject: [PATCH 37/42] Cache pixel index --- src/ctapipe/calib/camera/calibrator.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/calib/camera/calibrator.py b/src/ctapipe/calib/camera/calibrator.py index 3678e61fb33..baf3d2f1057 100644 --- a/src/ctapipe/calib/camera/calibrator.py +++ b/src/ctapipe/calib/camera/calibrator.py @@ -2,6 +2,7 @@ Definition of the `CameraCalibrator` class, providing all steps needed to apply calibration and image extraction, as well as supporting algorithms. """ +from functools import cache import astropy.units as u import numpy as np @@ -21,10 +22,16 @@ __all__ = ["CameraCalibrator"] +@cache +def _get_pixel_index(n_pixels): + """Cached version of ``np.arange(n_pixels)``""" + return np.arange(n_pixels) + + def _get_invalid_pixels(n_channels, n_pixels, pixel_status, selected_gain_channel): broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool) - index = np.arange(n_pixels) + index = _get_pixel_index(n_pixels) masks = ( pixel_status.hardware_failing_pixels, pixel_status.pedestal_failing_pixels, @@ -224,7 +231,7 @@ def _calibrate_dl1(self, event, tel_id): event.mon.tel[tel_id].pixel_status, selected_gain_channel, ) - pixel_index = np.arange(n_pixels) + pixel_index = _get_pixel_index(n_pixels) dl1_calib = event.calibration.tel[tel_id].dl1 readout = self.subarray.tel[tel_id].camera.readout From 596abb30ed54fde9aa08294f2dfcf907acf1bd57 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 17 Apr 2024 10:34:56 +0200 Subject: [PATCH 38/42] Adressing some minor comments --- src/ctapipe/calib/camera/tests/test_calibrator.py | 2 +- src/ctapipe/calib/camera/tests/test_gainselection.py | 3 +-- src/ctapipe/image/tests/test_extractor.py | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index f9fc4f04f56..55405730f32 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -149,7 +149,7 @@ def test_dl1_charge_calib_1_gain(example_subarray): # Randomize times and create pulses time_offset = random.uniform(-10, +10, (n_channels, n_pixels)) - y = norm.pdf(x, mid + time_offset[:, :, np.newaxis], pulse_sigma).astype("float32") + y = norm.pdf(x, mid + time_offset[..., np.newaxis], pulse_sigma).astype("float32") camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma)[np.newaxis, :] camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate diff --git a/src/ctapipe/calib/camera/tests/test_gainselection.py b/src/ctapipe/calib/camera/tests/test_gainselection.py index 6966791673d..9644887119e 100644 --- a/src/ctapipe/calib/camera/tests/test_gainselection.py +++ b/src/ctapipe/calib/camera/tests/test_gainselection.py @@ -26,8 +26,7 @@ def test_gain_selector(): def test_ndim(): shape = (5, 2, 2048, 128) - waveforms = np.indices(shape)[1] - waveforms[1] *= 2 + waveforms = np.zeros(shape) gain_selector = DummyGainSelector() with pytest.raises(ValueError): diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index a49d2b1324b..72b2e73f3d8 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -678,7 +678,7 @@ def test_waveform_extractor_factory_args(subarray): def test_extractor_tel_param(toymodel): waveforms, subarray, _, _, _, _ = toymodel - n_samples = waveforms.shape[-1] + _, _, n_samples = waveforms.shape config = Config( { From 9697f7ba662bd5791037ea46edcb3984fbb4d066 Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 17 Apr 2024 11:41:39 +0200 Subject: [PATCH 39/42] Merge calibrator test for 1 and 2 gain channels --- .../calib/camera/tests/test_calibrator.py | 283 ++++++------------ src/ctapipe/image/tests/test_extractor.py | 2 +- 2 files changed, 94 insertions(+), 191 deletions(-) diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index 55405730f32..806746728b5 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -130,7 +130,7 @@ def test_check_dl0_empty(example_event, example_subarray): assert (event.dl1.tel[tel_id].image == 2).all() -def test_dl1_charge_calib_1_gain(example_subarray): +def test_dl1_charge_calib(example_subarray): # copy because we mutate the camera, should not affect other tests subarray = deepcopy(example_subarray) camera = subarray.tel[1].camera @@ -139,7 +139,6 @@ def test_dl1_charge_calib_1_gain(example_subarray): sampling_rate = 2 camera.readout.sampling_rate = sampling_rate * u.GHz - n_channels = camera.readout.n_channels n_pixels = camera.geometry.n_pixels n_samples = 96 mid = n_samples // 2 @@ -147,208 +146,112 @@ def test_dl1_charge_calib_1_gain(example_subarray): random = np.random.default_rng(1) x = np.arange(n_samples) - # Randomize times and create pulses - time_offset = random.uniform(-10, +10, (n_channels, n_pixels)) - y = norm.pdf(x, mid + time_offset[..., np.newaxis], pulse_sigma).astype("float32") - - camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma)[np.newaxis, :] - camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate - - # Define absolute calibration coefficients - absolute = random.uniform(100, 1000, (n_channels, n_pixels)).astype("float32") - y *= absolute[..., np.newaxis] - - # Define relative coefficients - relative = random.normal(1, 0.01, (n_channels, n_pixels)) - y /= relative[..., np.newaxis] - - # Define pedestal - pedestal = random.uniform(-4, 4, (n_channels, n_pixels)) - y += pedestal[..., np.newaxis] - - event = ArrayEventContainer() - tel_id = list(subarray.tel.keys())[0] - event.dl0.tel[tel_id].waveform = y - event.dl0.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=int) - event.r1.tel[tel_id].selected_gain_channel = np.zeros(n_pixels, dtype=int) - - # Test default - calibrator = CameraCalibrator( - subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) - ) - calibrator(event) - np.testing.assert_allclose(event.dl1.tel[tel_id].image, y.sum(-1)[0], rtol=1e-4) - - event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal - event.calibration.tel[tel_id].dl1.absolute_factor = absolute - event.calibration.tel[tel_id].dl1.relative_factor = relative - - # Test without timing corrections - calibrator(event) - dl1 = event.dl1.tel[tel_id] - np.testing.assert_allclose(dl1.image, 1, rtol=1e-5) - - # we use only the first gain - expected_peak_time = (mid + time_offset[0]) / sampling_rate - np.testing.assert_allclose(dl1.peak_time, expected_peak_time, rtol=1e-5) - - # test with timing corrections - event.calibration.tel[tel_id].dl1.time_shift = time_offset / sampling_rate - calibrator(event) - - # more rtol since shifting might lead to reduced integral - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-5) - np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 - ) - - # test not applying time shifts - # now we should be back to the result without setting time shift - calibrator.apply_peak_time_shift = False - calibrator.apply_waveform_time_shift = False - calibrator(event) - - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-4) - np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, expected_peak_time, atol=1 - ) - - # We now use GlobalPeakWindowSum to see the effect of missing charge - # due to not correcting time offsets. - calibrator = CameraCalibrator( - subarray=subarray, - image_extractor=GlobalPeakWindowSum(subarray=subarray), - apply_waveform_time_shift=True, - ) - calibrator(event) - # test with timing corrections, should work - # higher rtol because we cannot shift perfectly - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=0.01) - np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 - ) - - # test deactivating timing corrections - calibrator.apply_waveform_time_shift = False - calibrator(event) - - # make sure we chose an example where the time shifts matter - # charges should be quite off due to summing around global shift - assert not np.allclose(event.dl1.tel[tel_id].image, 1, rtol=0.1) - assert not np.allclose(event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1) - - -def test_dl1_charge_calib_2_gain(example_subarray): - # copy because we mutate the camera, should not affect other tests - subarray = deepcopy(example_subarray) - camera = subarray.tel[1].camera - # test with a sampling_rate different than 1 to - # test if we handle time vs. slices correctly - sampling_rate = 2 - camera.readout.sampling_rate = sampling_rate * u.GHz - - n_pixels = camera.geometry.n_pixels - n_samples = 96 - mid = n_samples // 2 - pulse_sigma = 6 - random = np.random.default_rng(1) - x = np.arange(n_samples) - - # Randomize times and create pulses - time_offset_2_gain = random.uniform(-10, +10, (2, n_pixels)) - y_2_gain = np.ones((2, n_pixels, n_samples)) - y_2_gain[:] = norm.pdf( - x, mid + time_offset_2_gain[..., np.newaxis], pulse_sigma - ).astype("float32") - camera.readout.reference_pulse_shape = norm.pdf(x, mid, pulse_sigma) camera.readout.reference_pulse_shape = np.repeat( camera.readout.reference_pulse_shape[np.newaxis, :], 2, axis=0 ) camera.readout.reference_pulse_sample_width = 1 / camera.readout.sampling_rate - # Define absolute calibration coefficients - absolute_2_gain = random.uniform(100, 1000, (2, n_pixels)).astype("float32") - y_2_gain *= absolute_2_gain[..., np.newaxis] - - # Define relative coefficients - relative_2_gain = random.normal(1, 0.01, (2, n_pixels)) - y_2_gain /= relative_2_gain[..., np.newaxis] - - # Define pedestal - pedestal_2_gain = random.uniform(-4, 4, (2, n_pixels)) - y_2_gain += pedestal_2_gain[..., np.newaxis] - - event = ArrayEventContainer() - tel_id = list(subarray.tel.keys())[0] - event.dl0.tel[tel_id].waveform = y_2_gain - event.dl0.tel[tel_id].selected_gain_channel = None - event.r1.tel[tel_id].selected_gain_channel = None - - # Test default - calibrator = CameraCalibrator( - subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) - ) - calibrator(event) - np.testing.assert_allclose( - event.dl1.tel[tel_id].image[0], y_2_gain.sum(2)[0, :], rtol=1e-4 - ) + # test that it works for 1 and 2 gain channel + gain_channel = [1, 2] + for n_channels in gain_channel: + # Randomize times and create pulses + time_offset = random.uniform(-10, +10, (n_channels, n_pixels)) + y = np.zeros((n_channels, n_pixels, n_samples)) + y = norm.pdf(x, mid + time_offset[..., np.newaxis], pulse_sigma).astype( + "float32" + ) + + # Define absolute calibration coefficients + absolute = random.uniform(100, 1000, (n_channels, n_pixels)).astype("float32") + y *= absolute[..., np.newaxis] + + # Define relative coefficients + relative = random.normal(1, 0.01, (n_channels, n_pixels)) + y /= relative[..., np.newaxis] + + # Define pedestal + pedestal = random.uniform(-4, 4, (n_channels, n_pixels)) + y += pedestal[..., np.newaxis] + + event = ArrayEventContainer() + tel_id = list(subarray.tel.keys())[0] + event.dl0.tel[tel_id].waveform = y + selected_gain_channel = None + if n_channels == 1: + selected_gain_channel = np.zeros(n_pixels, dtype=int) + event.dl0.tel[tel_id].selected_gain_channel = selected_gain_channel + event.r1.tel[tel_id].selected_gain_channel = selected_gain_channel + + # Test default + calibrator = CameraCalibrator( + subarray=subarray, image_extractor=FullWaveformSum(subarray=subarray) + ) + calibrator(event) + np.testing.assert_allclose( + event.dl1.tel[tel_id].image, y.sum(-1).squeeze(), rtol=1e-4 + ) - event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal_2_gain - event.calibration.tel[tel_id].dl1.absolute_factor = absolute_2_gain - event.calibration.tel[tel_id].dl1.relative_factor = relative_2_gain + event.calibration.tel[tel_id].dl1.pedestal_offset = pedestal + event.calibration.tel[tel_id].dl1.absolute_factor = absolute + event.calibration.tel[tel_id].dl1.relative_factor = relative - # Test without timing corrections - calibrator(event) - dl1 = event.dl1.tel[tel_id] - np.testing.assert_allclose(dl1.image[0], 1, rtol=1e-5) - expected_peak_time = (mid + time_offset_2_gain) / sampling_rate - np.testing.assert_allclose(dl1.peak_time, expected_peak_time, rtol=1e-5) - - # test with timing corrections - event.calibration.tel[tel_id].dl1.time_shift = time_offset_2_gain / sampling_rate - calibrator(event) + # Test without timing corrections + calibrator(event) + dl1 = event.dl1.tel[tel_id] + np.testing.assert_allclose(dl1.image, 1, rtol=1e-5) - # more rtol since shifting might lead to reduced integral - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-5) - np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 - ) + expected_peak_time = (mid + time_offset) / sampling_rate + np.testing.assert_allclose( + dl1.peak_time, expected_peak_time.squeeze(), rtol=1e-5 + ) - # test not applying time shifts - # now we should be back to the result without setting time shift - calibrator.apply_peak_time_shift = False - calibrator.apply_waveform_time_shift = False - calibrator(event) + # test with timing corrections + event.calibration.tel[tel_id].dl1.time_shift = time_offset / sampling_rate + calibrator(event) - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-4) - np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, expected_peak_time, atol=1 - ) + # more rtol since shifting might lead to reduced integral + np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-5) + np.testing.assert_allclose( + event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + ) - # We now use GlobalPeakWindowSum to see the effect of missing charge - # due to not correcting time offsets. - calibrator = CameraCalibrator( - subarray=subarray, - image_extractor=GlobalPeakWindowSum(subarray=subarray), - apply_waveform_time_shift=True, - ) - calibrator(event) - # test with timing corrections, should work - # higher rtol because we cannot shift perfectly - np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=0.01) - np.testing.assert_allclose( - event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 - ) + # test not applying time shifts + # now we should be back to the result without setting time shift + calibrator.apply_peak_time_shift = False + calibrator.apply_waveform_time_shift = False + calibrator(event) - # test deactivating timing corrections - calibrator.apply_waveform_time_shift = False - calibrator(event) + np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=1e-4) + np.testing.assert_allclose( + event.dl1.tel[tel_id].peak_time, expected_peak_time.squeeze(), atol=1 + ) + + # We now use GlobalPeakWindowSum to see the effect of missing charge + # due to not correcting time offsets. + calibrator = CameraCalibrator( + subarray=subarray, + image_extractor=GlobalPeakWindowSum(subarray=subarray), + apply_waveform_time_shift=True, + ) + calibrator(event) + # test with timing corrections, should work + # higher rtol because we cannot shift perfectly + np.testing.assert_allclose(event.dl1.tel[tel_id].image, 1, rtol=0.01) + np.testing.assert_allclose( + event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + ) + + # test deactivating timing corrections + calibrator.apply_waveform_time_shift = False + calibrator(event) - # make sure we chose an example where the time shifts matter - # charges should be quite off due to summing around global shift - assert not np.allclose(event.dl1.tel[tel_id].image, 1, rtol=0.1) - assert not np.allclose(event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1) + # make sure we chose an example where the time shifts matter + # charges should be quite off due to summing around global shift + assert not np.allclose(event.dl1.tel[tel_id].image, 1, rtol=0.1) + assert not np.allclose( + event.dl1.tel[tel_id].peak_time, mid / sampling_rate, atol=1 + ) def test_shift_waveforms(): diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py index 72b2e73f3d8..a49d2b1324b 100644 --- a/src/ctapipe/image/tests/test_extractor.py +++ b/src/ctapipe/image/tests/test_extractor.py @@ -678,7 +678,7 @@ def test_waveform_extractor_factory_args(subarray): def test_extractor_tel_param(toymodel): waveforms, subarray, _, _, _, _ = toymodel - _, _, n_samples = waveforms.shape + n_samples = waveforms.shape[-1] config = Config( { From f932630619d975de7b7a0b09d0a772afe001407e Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Wed, 17 Apr 2024 13:09:02 +0200 Subject: [PATCH 40/42] Move _apply_correction into the ImageExtractor parent class - _apply_correction is now a static method from ImageExtractors parent class - _calculate_correction is moved to the parent ImageExtractors class to avoid writing the docstrings multiple times. It needs to be overwritten by the child components --- src/ctapipe/image/extractor.py | 144 +++++++++------------------------ 1 file changed, 36 insertions(+), 108 deletions(-) diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py index 8f2c4f07e9a..1d4cfd5881d 100644 --- a/src/ctapipe/image/extractor.py +++ b/src/ctapipe/image/extractor.py @@ -359,15 +359,6 @@ def integration_correction( return correction -def _apply_correction(charge, correction, selected_gain_channel): - """ - Helper function for applying the integration correction for certain `ImageExtractor`s. - """ - if selected_gain_channel is None: - return (charge * correction[:, np.newaxis]).astype(charge.dtype) - return (charge * correction[selected_gain_channel]).astype(charge.dtype) - - class ImageExtractor(TelescopeComponent): def __init__(self, subarray, config=None, parent=None, **kwargs): """ @@ -402,6 +393,37 @@ def __init__(self, subarray, config=None, parent=None, **kwargs): for tel_id, telescope in subarray.tel.items() } + def _calculate_correction(self, tel_id): + """ + Calculate the correction for the extracted charge such that the value + returned would equal 1 for a noise-less unit pulse. `ImageExtractor` types + calculating corrections need to overwrite this method. + + This method should be decorated with @lru_cache to ensure it is only + calculated once per telescope. + + Parameters + ---------- + tel_id : int + + Returns + ------- + correction : ndarray + The correction to apply to an extracted charge using this ImageExtractor + Has size n_channels, as a different correction value might be required + for different gain channels. + """ + pass + + @staticmethod + def _apply_correction(charge, correction, selected_gain_channel): + """ + Helper function for applying the integration correction for certain `ImageExtractor`s. + """ + if selected_gain_channel is None: + return (charge * correction[:, np.newaxis]).astype(charge.dtype) + return (charge * correction[selected_gain_channel]).astype(charge.dtype) + @abstractmethod def __call__( self, waveforms, tel_id, selected_gain_channel, broken_pixels @@ -476,24 +498,6 @@ class FixedWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): - """ - Calculate the correction for the extracted charge such that the value - returned would equal 1 for a noise-less unit pulse. - - This method is decorated with @lru_cache to ensure it is only - calculated once per telescope. - - Parameters - ---------- - tel_id : int - - Returns - ------- - correction : ndarray - The correction to apply to an extracted charge using this ImageExtractor - Has size n_channels, as a different correction value might be required - for different gain channels. - """ readout = self.subarray.tel[tel_id].camera.readout return integration_correction( readout.reference_pulse_shape, @@ -515,7 +519,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - charge = _apply_correction(charge, correction, selected_gain_channel) + charge = self._apply_correction(charge, correction, selected_gain_channel) # reduce dimensions for gain selected data to (n_pixels, ) if selected_gain_channel is not None: @@ -562,24 +566,6 @@ class GlobalPeakWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): - """ - Calculate the correction for the extracted charge such that the value - returned would equal 1 for a noise-less unit pulse. - - This method is decorated with @lru_cache to ensure it is only - calculated once per telescope. - - Parameters - ---------- - tel_id : int - - Returns - ------- - correction : ndarray - The correction to apply to an extracted charge using this ImageExtractor - Has size n_channels, as a different correction value might be required - for different gain channels. - """ readout = self.subarray.tel[tel_id].camera.readout return integration_correction( readout.reference_pulse_shape, @@ -619,7 +605,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - charge = _apply_correction(charge, correction, selected_gain_channel) + charge = self._apply_correction(charge, correction, selected_gain_channel) # reduce dimensions for gain selected data to (n_pixels, ) if selected_gain_channel is not None: @@ -651,24 +637,6 @@ class LocalPeakWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): - """ - Calculate the correction for the extracted charge such that the value - returned would equal 1 for a noise-less unit pulse. - - This method is decorated with @lru_cache to ensure it is only - calculated once per telescope. - - Parameters - ---------- - tel_id : int - - Returns - ------- - correction : ndarray - The correction to apply to an extracted charge using this ImageExtractor - Has size n_channels, as a different correction value might be required - for different gain channels. - """ readout = self.subarray.tel[tel_id].camera.readout return integration_correction( readout.reference_pulse_shape, @@ -691,7 +659,7 @@ def __call__( ) if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - charge = _apply_correction(charge, correction, selected_gain_channel) + charge = self._apply_correction(charge, correction, selected_gain_channel) # reduce dimensions for gain selected data to (n_pixels, ) if selected_gain_channel is not None: @@ -716,28 +684,6 @@ class SlidingWindowMaxSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): - """ - Calculate the correction for the extracted charge such that the value - returned would equal 1 for a noise-less unit pulse. - - This method is decorated with @lru_cache to ensure it is only - calculated once per telescope. - - The same procedure as for the actual SlidingWindowMaxSum extractor is used, but - on the reference pulse_shape (that is also more finely binned) - - Parameters - ---------- - tel_id : int - - Returns - ------- - correction : ndarray - The correction to apply to an extracted charge using this ImageExtractor - Has size n_channels, as a different correction value might be required - for different gain channels. - """ - readout = self.subarray.tel[tel_id].camera.readout # compute the number of slices to integrate in the pulse template @@ -775,7 +721,7 @@ def __call__( if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - charge = _apply_correction(charge, correction, selected_gain_channel) + charge = self._apply_correction(charge, correction, selected_gain_channel) # reduce dimensions for gain selected data to (n_pixels, ) if selected_gain_channel is not None: @@ -813,24 +759,6 @@ class NeighborPeakWindowSum(ImageExtractor): @lru_cache(maxsize=128) def _calculate_correction(self, tel_id): - """ - Calculate the correction for the extracted charge such that the value - returned would equal 1 for a noise-less unit pulse. - - This method is decorated with @lru_cache to ensure it is only - calculated once per telescope. - - Parameters - ---------- - tel_id : int - - Returns - ------- - correction : ndarray - The correction to apply to an extracted charge using this ImageExtractor - Has size n_channels, as a different correction value might be required - for different gain channels. - """ readout = self.subarray.tel[tel_id].camera.readout return integration_correction( readout.reference_pulse_shape, @@ -861,7 +789,7 @@ def __call__( if self.apply_integration_correction.tel[tel_id]: correction = self._calculate_correction(tel_id=tel_id) - charge = _apply_correction(charge, correction, selected_gain_channel) + charge = self._apply_correction(charge, correction, selected_gain_channel) # reduce dimensions for gain selected data to (n_pixels, ) if selected_gain_channel is not None: From 07646b76efb2451158354a491b4ab28e4190362a Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 09:07:19 +0200 Subject: [PATCH 41/42] Remove unnecessary instantiation --- src/ctapipe/calib/camera/tests/test_calibrator.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ctapipe/calib/camera/tests/test_calibrator.py b/src/ctapipe/calib/camera/tests/test_calibrator.py index 806746728b5..b5a3875f62e 100644 --- a/src/ctapipe/calib/camera/tests/test_calibrator.py +++ b/src/ctapipe/calib/camera/tests/test_calibrator.py @@ -157,7 +157,6 @@ def test_dl1_charge_calib(example_subarray): for n_channels in gain_channel: # Randomize times and create pulses time_offset = random.uniform(-10, +10, (n_channels, n_pixels)) - y = np.zeros((n_channels, n_pixels, n_samples)) y = norm.pdf(x, mid + time_offset[..., np.newaxis], pulse_sigma).astype( "float32" ) From e37bb416f72a312ce1f747f96d6d782111ed866f Mon Sep 17 00:00:00 2001 From: Jonas Hackfeld Date: Tue, 23 Apr 2024 09:13:56 +0200 Subject: [PATCH 42/42] Fix wrong merge in table_writer_reader.py - it duplicated the changes made in main, had to remove one of them --- examples/core/table_writer_reader.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/examples/core/table_writer_reader.py b/examples/core/table_writer_reader.py index 7c8e54cbc2f..0bb5197901b 100644 --- a/examples/core/table_writer_reader.py +++ b/examples/core/table_writer_reader.py @@ -236,14 +236,6 @@ def create_stream(n_event): h5.close() -###################################################################### -# close the file afterwards -# - - -h5.close() - - ###################################################################### # Reading one-row-at-a-time: # ~~~~~~~~~~~~~~~~~~~~~~~~~~