From 51c703740a40aa8b1833764d62876808d5315cca Mon Sep 17 00:00:00 2001 From: nieves Date: Mon, 15 May 2023 22:47:19 +0200 Subject: [PATCH 1/8] added extra test for timing and solved bug --- ctapipe/image/extractor.py | 4 +- ctapipe/image/tests/test_extractor.py | 64 +++++++++++++++++++++++++-- 2 files changed, 63 insertions(+), 5 deletions(-) diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py index 5893ebf1dd0..87f75bedb18 100644 --- a/ctapipe/image/extractor.py +++ b/ctapipe/image/extractor.py @@ -1681,7 +1681,9 @@ def __call__( d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1) peak_time = adaptive_centroid( d_waveforms, - max(0, min(np.round(peak_index - pz2d).astype(int))), + np.clip( + np.round(peak_index - pz2d).astype(int), 0, np.shape(d_waveforms)[1] + ), leading_edge_rel_descend_limit, ) diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py index 768329bd07c..20764499bda 100644 --- a/ctapipe/image/tests/test_extractor.py +++ b/ctapipe/image/tests/test_extractor.py @@ -10,6 +10,7 @@ from traitlets.traitlets import TraitError from ctapipe.core import non_abstract_children +from ctapipe.image.cleaning import dilate from ctapipe.image.extractor import ( FixedWindowSum, FlashCamExtractor, @@ -109,6 +110,44 @@ def toymodel(subarray): return get_test_toymodel(subarray) +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 + geometry = subarray.tel[tel_id].camera.geometry + n_samples = 96 + readout = subarray.tel[tel_id].camera.readout + + random = np.random.RandomState(1) + charge = random.uniform(minCharge, maxCharge, n_pixels) + tmid = (n_samples // 2) / readout.sampling_rate.to_value(u.GHz) + tmax = (n_samples - 1) / readout.sampling_rate.to_value(u.GHz) + time = random.uniform(0, tmax, n_pixels) + pix_id = 1 + mask = geometry.pix_id == pix_id + + dilated_mask = mask.copy() + for row in range(4): + dilated_mask = dilate(geometry, dilated_mask) + + x_d = subarray.tel[tel_id].camera.geometry.pix_x.value + min_xd = np.min(x_d[dilated_mask]) + diff_xd = x_d[dilated_mask] - min_xd + slope = 15 # ns + intercept = 0 # ns + time[dilated_mask] = tmid + (slope * diff_xd + intercept) + waveform_model = WaveformModel.from_camera_readout(readout) + waveform = waveform_model.get_waveform(charge, time, n_samples) + + selected_gain_channel = np.zeros(charge.size, dtype=np.int64) + + return waveform, subarray, tel_id, selected_gain_channel, charge, time, dilated_mask + + +@pytest.fixture(scope="module") +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) @@ -741,8 +780,26 @@ def test_upsampling(toymodel_mst_fc): ) +def test_FC_timing(toymodel_mst_fc_time): + # Test time on toy model with time gradient (other toy not sensitive to timing bugs!!) + ( + waveforms, + subarray, + tel_id, + selected_gain_channel, + true_charge, + true_time, + mask, + ) = toymodel_mst_fc_time + extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) + broken_pixels = np.zeros(waveforms.shape[0], 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 == True + + def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): - # Test on toy model + # Test charge on standard toy model ( waveforms, subarray, @@ -751,11 +808,10 @@ def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): true_charge, true_time, ) = toymodel_mst_fc - extractor = FlashCamExtractor(subarray=subarray) + extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) broken_pixels = np.zeros(waveforms.shape[0], dtype=bool) dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels) - assert_allclose(dl1.image, true_charge, rtol=0.15) - assert_allclose(dl1.peak_time, true_time, atol=1.0) + assert_allclose(dl1.image, true_charge, rtol=0.1) assert dl1.is_valid == True # Test on prod5 simulations From 68a2495beaecc5f3833f34020e1745d54cdf126a Mon Sep 17 00:00:00 2001 From: nieves Date: Mon, 15 May 2023 22:59:01 +0200 Subject: [PATCH 2/8] modified names --- ctapipe/image/tests/test_extractor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py index 20764499bda..3d7aff91e92 100644 --- a/ctapipe/image/tests/test_extractor.py +++ b/ctapipe/image/tests/test_extractor.py @@ -780,7 +780,7 @@ def test_upsampling(toymodel_mst_fc): ) -def test_FC_timing(toymodel_mst_fc_time): +def test_FC_time(toymodel_mst_fc_time): # Test time on toy model with time gradient (other toy not sensitive to timing bugs!!) ( waveforms, From be38c997a85251761d1f853778c200dc33e9e92f Mon Sep 17 00:00:00 2001 From: nieves Date: Mon, 15 May 2023 23:34:25 +0200 Subject: [PATCH 3/8] added doc --- docs/changes/2333.bug.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/2333.bug.rst diff --git a/docs/changes/2333.bug.rst b/docs/changes/2333.bug.rst new file mode 100644 index 00000000000..021e48bafc3 --- /dev/null +++ b/docs/changes/2333.bug.rst @@ -0,0 +1 @@ +Solve bug in time reconstruction with the FlashCamExtractor \ No newline at end of file From 540fb1ee1dd13c75b003031f2adb5d0470f472a9 Mon Sep 17 00:00:00 2001 From: nieves Date: Mon, 15 May 2023 23:44:33 +0200 Subject: [PATCH 4/8] added test when leading_edge_timing=False --- ctapipe/image/tests/test_extractor.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/ctapipe/image/tests/test_extractor.py b/ctapipe/image/tests/test_extractor.py index 3d7aff91e92..41adff412c3 100644 --- a/ctapipe/image/tests/test_extractor.py +++ b/ctapipe/image/tests/test_extractor.py @@ -791,12 +791,19 @@ def test_FC_time(toymodel_mst_fc_time): true_time, mask, ) = toymodel_mst_fc_time + extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=True) broken_pixels = np.zeros(waveforms.shape[0], 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 == True + extractor = FlashCamExtractor(subarray=subarray, leading_edge_timing=False) + broken_pixels = np.zeros(waveforms.shape[0], 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 == True + def test_flashcam_extractor(toymodel_mst_fc, prod5_gamma_simtel_path): # Test charge on standard toy model From 1d49021afc25990725ecf4e9d7333de969c089dc Mon Sep 17 00:00:00 2001 From: nieves Date: Tue, 16 May 2023 11:42:32 +0200 Subject: [PATCH 5/8] added comment and added out arg in clip --- ctapipe/image/extractor.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py index 87f75bedb18..5787d3a2d41 100644 --- a/ctapipe/image/extractor.py +++ b/ctapipe/image/extractor.py @@ -1679,12 +1679,13 @@ def __call__( if leading_edge_timing: d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1) + peak_index = np.round(peak_index - pz2d).astype( + int + ) # correct the offset between leading edge peak and deconvolved peak + n_samples = d_waveforms.shape[-1] + np.clip(peak_index, 0, n_samples - 1, out=peak_index) peak_time = adaptive_centroid( - d_waveforms, - np.clip( - np.round(peak_index - pz2d).astype(int), 0, np.shape(d_waveforms)[1] - ), - leading_edge_rel_descend_limit, + d_waveforms, peak_index, leading_edge_rel_descend_limit ) if gain != 0: From 149d544c844b1ef9b327ca4382e3feece8e4d21b Mon Sep 17 00:00:00 2001 From: Clara Escanuela Nieves <53437218+clara-escanuela@users.noreply.github.com> Date: Tue, 16 May 2023 11:49:07 +0200 Subject: [PATCH 6/8] code style Co-authored-by: Maximilian Linhoff --- ctapipe/image/extractor.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py index 5787d3a2d41..4a35f14a5ce 100644 --- a/ctapipe/image/extractor.py +++ b/ctapipe/image/extractor.py @@ -1679,9 +1679,8 @@ def __call__( if leading_edge_timing: d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1) - peak_index = np.round(peak_index - pz2d).astype( - int - ) # correct the offset between leading edge peak and deconvolved peak + # correct the offset between leading edge peak and deconvolved peak + peak_index = np.round(peak_index - pz2d).astype(int) n_samples = d_waveforms.shape[-1] np.clip(peak_index, 0, n_samples - 1, out=peak_index) peak_time = adaptive_centroid( From a93f2758d9fc5cdda4a1d723d10027e7293bee6e Mon Sep 17 00:00:00 2001 From: Clara Escanuela Nieves <53437218+clara-escanuela@users.noreply.github.com> Date: Tue, 16 May 2023 11:56:08 +0200 Subject: [PATCH 7/8] Update docs/changes/2333.bug.rst Co-authored-by: Maximilian Linhoff --- docs/changes/2333.bug.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/changes/2333.bug.rst b/docs/changes/2333.bug.rst index 021e48bafc3..0a5d33c3fc4 100644 --- a/docs/changes/2333.bug.rst +++ b/docs/changes/2333.bug.rst @@ -1 +1 @@ -Solve bug in time reconstruction with the FlashCamExtractor \ No newline at end of file +Fix a bug in the peak_time estimation of ``FlashCamExtractor`` (See issue `#2332 `_) \ No newline at end of file From 9ecd144815316b419d09abdc1ba0da3fa507e729 Mon Sep 17 00:00:00 2001 From: nieves Date: Tue, 16 May 2023 11:51:07 +0200 Subject: [PATCH 8/8] code style --- ctapipe/image/extractor.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py index 4a35f14a5ce..8b791a169db 100644 --- a/ctapipe/image/extractor.py +++ b/ctapipe/image/extractor.py @@ -1679,8 +1679,9 @@ def __call__( if leading_edge_timing: d_waveforms = deconvolve(waveforms, 0.0, upsampling, 1) + # correct the offset between leading edge peak and deconvolved peak - peak_index = np.round(peak_index - pz2d).astype(int) + peak_index = np.round(peak_index - pz2d).astype(int) n_samples = d_waveforms.shape[-1] np.clip(peak_index, 0, n_samples - 1, out=peak_index) peak_time = adaptive_centroid(