From 41a15ab3bb3cbda97abeb59032845d002efd3d8c Mon Sep 17 00:00:00 2001 From: Ander Date: Fri, 28 May 2021 11:15:36 +0200 Subject: [PATCH 01/14] Add test files --- invisible_cities/database/test_data/blr_fpga_examples.h5 | 3 +++ .../database/test_data/exact_result_multipeak_valdrada.h5 | 3 +++ invisible_cities/database/test_data/exact_result_valdrada.h5 | 3 +++ 3 files changed, 9 insertions(+) create mode 100755 invisible_cities/database/test_data/blr_fpga_examples.h5 create mode 100644 invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5 create mode 100644 invisible_cities/database/test_data/exact_result_valdrada.h5 diff --git a/invisible_cities/database/test_data/blr_fpga_examples.h5 b/invisible_cities/database/test_data/blr_fpga_examples.h5 new file mode 100755 index 000000000..37be52aec --- /dev/null +++ b/invisible_cities/database/test_data/blr_fpga_examples.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5acb7051d8dd3f7f0433baa616937089fd66d54d5ddd8bd470ad711d9be7cc0f +size 15111797 diff --git a/invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5 b/invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5 new file mode 100644 index 000000000..fe38b0535 --- /dev/null +++ b/invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0f11a28f2437e57f243594ee48aec241f821543f8d5fb3dd94545b4d6e2dbb5d +size 21870 diff --git a/invisible_cities/database/test_data/exact_result_valdrada.h5 b/invisible_cities/database/test_data/exact_result_valdrada.h5 new file mode 100644 index 000000000..54ce6b3c6 --- /dev/null +++ b/invisible_cities/database/test_data/exact_result_valdrada.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:599e156701fc2e8b2cf35d3d421990feecb4448523afcae07c261eecadb769c5 +size 21863 From aa6f470743cced512d2996039d4a7e1366174e62 Mon Sep 17 00:00:00 2001 From: Ander Date: Fri, 28 May 2021 11:17:20 +0200 Subject: [PATCH 02/14] Add config file --- invisible_cities/config/valdrada.conf | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 invisible_cities/config/valdrada.conf diff --git a/invisible_cities/config/valdrada.conf b/invisible_cities/config/valdrada.conf new file mode 100644 index 000000000..cfd0424cf --- /dev/null +++ b/invisible_cities/config/valdrada.conf @@ -0,0 +1,19 @@ +files_in = '$ICDIR/database/test_data/blr_fpga_examples.h5' +file_out = '$ICDIR/database/test_data/valdrada.h5' + +compression = 'ZLIB4' +run_number = 8093 +detector_db = 'new' +event_range = all + +trigger_config = {'coincidence_window':64, 'discard_width':40, + 'multipeak':{'q_min':100000, 'time_min':2*mus, 'time_after':800*mus}} + +channel_config = {0: {'q_min' :5000, 'q_max' :50000, + 'time_min' :2*mus, 'time_max': 40*mus, + 'baseline_dev': 10, 'amp_max' : 1000, + 'pulse_valid_ext':50*ns}, + 2: {'q_min' :5000, 'q_max' :50000, + 'time_min' :2*mus, 'time_max': 40*mus, + 'baseline_dev': 10, 'amp_max' : 1000, + 'pulse_valid_ext':50*ns}} From aac25b036dcb8eb81574cfcb43087f7784da3dc4 Mon Sep 17 00:00:00 2001 From: Ander Date: Fri, 28 May 2021 11:27:56 +0200 Subject: [PATCH 03/14] Add related tests --- invisible_cities/cities/valdrada_test.py | 53 +++ invisible_cities/conftest.py | 43 +++ invisible_cities/io/trigger_io_test.py | 56 +++ .../reco/trigger_functions_test.py | 330 ++++++++++++++++++ invisible_cities/sierpe/blr_test.py | 84 ++++- 5 files changed, 563 insertions(+), 3 deletions(-) create mode 100644 invisible_cities/cities/valdrada_test.py create mode 100644 invisible_cities/io/trigger_io_test.py create mode 100644 invisible_cities/reco/trigger_functions_test.py diff --git a/invisible_cities/cities/valdrada_test.py b/invisible_cities/cities/valdrada_test.py new file mode 100644 index 000000000..586643fa9 --- /dev/null +++ b/invisible_cities/cities/valdrada_test.py @@ -0,0 +1,53 @@ +import os +import tables as tb + +from . valdrada import valdrada +from .. core.testing_utils import assert_tables_equality + + +def test_valdrada_contains_all_tables(trigger_config): + conf, PATH_OUT = trigger_config + valdrada(**conf) + with tb.open_file(PATH_OUT) as h5out: + assert "Run" in h5out.root + assert "Run/events" in h5out.root + assert "Run/runInfo" in h5out.root + assert "Trigger" in h5out.root + assert "Trigger/events" in h5out.root + assert "Trigger/trigger" in h5out.root + assert "Trigger/DST" in h5out.root + + +def test_valdrada_exact_result_multipeak(ICDATADIR, trigger_config): + true_out = os.path.join(ICDATADIR, "test_trigger_fpga_multipeak.h5") + conf, PATH_OUT = trigger_config + valdrada(**conf) + + tables = ("Trigger/events", "Trigger/trigger", "Trigger/DST", + "Run/events" , "Run/runInfo") + + with tb.open_file(true_out) as true_output_file: + with tb.open_file(PATH_OUT) as output_file: + for table in tables: + assert hasattr(output_file.root, table) + got = getattr( output_file.root, table) + expected = getattr(true_output_file.root, table) + assert_tables_equality(got, expected) + + +def test_valdrada_exact_result(ICDATADIR, trigger_config): + true_out = os.path.join(ICDATADIR, "test_trigger_fpga.h5") + conf, PATH_OUT = trigger_config + conf['trigger_config']['multipeak'] = None + valdrada(**conf) + + tables = ("Trigger/events", "Trigger/trigger", "Trigger/DST", + "Run/events" , "Run/runInfo") + + with tb.open_file(true_out) as true_output_file: + with tb.open_file(PATH_OUT) as output_file: + for table in tables: + assert hasattr(output_file.root, table) + got = getattr( output_file.root, table) + expected = getattr(true_output_file.root, table) + assert_tables_equality(got, expected) diff --git a/invisible_cities/conftest.py b/invisible_cities/conftest.py index 32bfbd1ae..1f99e9f21 100644 --- a/invisible_cities/conftest.py +++ b/invisible_cities/conftest.py @@ -50,6 +50,11 @@ def example_blr_wfs_filename(ICDATADIR): return os.path.join(ICDATADIR, "blr_examples.h5") +@pytest.fixture(scope='session') +def example_blr_fpga_wfs_filename(ICDATADIR): + return os.path.join(ICDATADIR, "blr_fpga_examples.h5") + + @pytest.fixture(scope = 'session', params = ['electrons_40keV_z250_MCRD.h5']) def electron_MCRD_file(request, ICDATADIR): @@ -699,6 +704,44 @@ def deconvolution_config(ICDIR, ICDATADIR, PSFDIR, config_tmpdir): return conf, PATH_OUT +<<<<<<< HEAD +======= + +@pytest.fixture(scope='function') # Needs to be function as the config dict is modified when running +def trigger_config(ICDIR, ICDATADIR, config_tmpdir): + PATH_IN = os.path.join(ICDATADIR , "blr_fpga_examples.h5") + PATH_OUT = os.path.join(config_tmpdir, "trigger_city_out.h5") + nevt_req = 10 + conf = dict(files_in = PATH_IN , + file_out = PATH_OUT, + event_range = nevt_req, + compression = 'ZLIB4', + print_mod = 1000, + run_number = 8093, + trigger_config = dict(coincidence_window = 64 + ,discard_width = 40 + ,multipeak = dict(q_min = 100000 + ,time_min = 2 *units.mus + ,time_after = 800*units.mus)), + channel_config = {0 : dict(q_min = 5000 + ,q_max = 50000 + ,time_min = 2 *units.mus + ,time_max = 40*units.mus + ,baseline_dev = 10 + ,amp_max = 1000 + ,pulse_valid_ext = 50*units.ns) + ,2 : dict(q_min = 5000 + ,q_max = 50000 + ,time_min = 2 *units.mus + ,time_max = 40*units.mus + ,baseline_dev = 10 + ,amp_max = 1000 + ,pulse_valid_ext = 50*units.ns)}) + + return conf, PATH_OUT + + +>>>>>>> f8fa2343... fixup! Add related tests ## To make very slow tests only run with specific option def pytest_addoption(parser): parser.addoption( diff --git a/invisible_cities/io/trigger_io_test.py b/invisible_cities/io/trigger_io_test.py new file mode 100644 index 000000000..635f1ac6c --- /dev/null +++ b/invisible_cities/io/trigger_io_test.py @@ -0,0 +1,56 @@ +import os + +import numpy as np +import tables as tb + +from numpy.testing import assert_allclose +from hypothesis import given +from hypothesis.strategies import integers +from hypothesis.strategies import lists +from hypothesis.strategies import booleans +from hypothesis.strategies import composite + +from ..core.testing_utils import assert_dataframes_equal +from ..io.dst_io import load_dst +from . trigger_io import trigger_dst_writer + + +@composite +def trigger_values(draw): + size = draw(integers(min_value=1, max_value=10)) + elements = [] + for _ in range(6): + elements.append(draw(lists(integers(min_value=0, max_value=np.iinfo(np.uint32).max), + min_size=size, max_size=size))) + for _ in range(4): + elements.append(draw(lists(booleans(), min_size=size, max_size=size))) + for _ in range(5): + elements.append(draw(lists(integers(min_value=np.iinfo(np.int32).min, max_value=np.iinfo(np.int32).max), + min_size=size, max_size=size))) + return elements + + +@given(triggers=trigger_values()) +def test_trigger_dst_writer(config_tmpdir, triggers): + output_file = os.path.join(config_tmpdir, "test_trigger_dst.h5") + + trigger_save = [list(a) for a in zip(*triggers)] + with tb.open_file(output_file, 'w') as h5out: + write = trigger_dst_writer(h5out) + write(trigger_save) + + trigger_dst = load_dst(output_file, group='Trigger', node='DST') + + col_names = ["event" , "pmt" , "trigger_time" , "q" + ,"width" , "height" , "valid_q" , "valid_w" + ,"valid_h" , "valid_peak", "valid_all" , "baseline" + ,"max_height", "n_coinc" , "closest_ttime", "closest_pmt"] + + assert np.all([tname == cname for tname, cname in zip(trigger_dst.columns.values, col_names)]) + + # Check values stored, the if is needed because the valid_all condition is done in the writer, not on the input trigger + for i, colname in enumerate(col_names): + if colname == 'valid_all': + assert_allclose(np.all(triggers[6:10], axis=0), trigger_dst.loc[:, colname]) + elif i>10: assert_allclose(triggers[i-1] , trigger_dst.loc[:, colname]) + else : assert_allclose(triggers[i ] , trigger_dst.loc[:, colname]) diff --git a/invisible_cities/reco/trigger_functions_test.py b/invisible_cities/reco/trigger_functions_test.py new file mode 100644 index 000000000..74d9893f4 --- /dev/null +++ b/invisible_cities/reco/trigger_functions_test.py @@ -0,0 +1,330 @@ +import numpy as np + +from pytest import fixture +from pytest import mark +from hypothesis import given +from numpy .testing import assert_allclose +from hypothesis.strategies import integers +from hypothesis.strategies import lists +from hypothesis.strategies import booleans +from hypothesis.strategies import composite + +from .. reco.trigger_functions import get_trigger_candidates +from .. reco.trigger_functions import retrieve_trigger_information +from .. reco.trigger_functions import check_trigger_coincidence +from .. core.core_functions import in_range + + +@fixture(scope="session") +def channel_conf(): + channel_dict = {'q_min' : 5000, 'q_max' : 50000, + 'time_min' : 2000, 'time_max': 10000, + 'baseline_dev' : 5, 'amp_max' : 1000, + 'pulse_valid_ext': 15} + return channel_dict + + +@fixture(scope="session") +def trigger_conf(): + trigger_dict = {'coincidence_window':64, 'discard_width':40, 'multipeak':None} + return trigger_dict + + +@fixture(scope="session") +def double_square_wf(): + # Fake signal with two-nearby square pulses and a smaller additional one at the beginning. + n_baseline = 64000 + base_window = 1024 + thr = 5 + wf = np.zeros(n_baseline) + start1 = np.random.randint(2*base_window, n_baseline // 2) + length = np.random.randint(n_baseline // 50, n_baseline // 4) + stop1 = start1 + length + start2 = stop1 + 10 + stop2 = start2 + length + + wf[start1:stop1] = np.full(length, np.random.randint(thr+1, 50)) + wf[start2:stop2] = np.full(length, np.random.randint(thr+1, 50)) + wf[1500 : 1520] = np.full(20 , np.random.randint(thr+1, 50)) # Extra pulse to check discard_width + return wf, thr, length, start1, start2 + + +@composite +def trigger_time_and_validity(draw): + # Set of trigger times an validity flags, to test trigger coincidence + size = draw (integers(min_value=1, max_value=100)) + flag = lists(integers(min_value=0, max_value=1), min_size=4, max_size=4) + + i = draw(lists(integers(min_value=0), min_size=size, max_size=size)) + flags = draw(lists(flag, min_size=size, max_size=size)) + return (i, flags) + + +@given(evt_id=integers(1, 10), pmt_id=integers(1, 10)) +def test_get_trigger_candidates_single(double_square_wf, evt_id, pmt_id, channel_conf, trigger_conf): + # Given a square pulse check that is triggering where the pulse starts and values match the expected ones. + channel_config = channel_conf.copy() + trigger_config = trigger_conf.copy() + + wf = double_square_wf[0] + thr = double_square_wf[1] + length = double_square_wf[2] + start1 = double_square_wf[3] + start2 = double_square_wf[4] + + peak = wf[start1:start2+length+channel_config['pulse_valid_ext']] + t_true = start1 + q_true = peak[:-2].sum() + val_q = in_range(q_true + , channel_config['q_min'], channel_config['q_max'] + , left_closed=False, right_closed=False) + w_true = len(peak) + val_w = in_range(w_true + , channel_config['time_min'], channel_config['time_max'] + , left_closed=True , right_closed=False) + h_true = peak.max() + val_h = h_true < channel_config['amp_max'] + val_p = True + + + triggers = get_trigger_candidates(wf, evt_id, pmt_id + ,channel_config, trigger_config + ,np.zeros(wf.shape), wf.max()) + + elements = [evt_id , pmt_id + ,t_true , q_true , w_true, h_true + ,val_q , val_w , val_h , val_p + ,wf[start1-1], wf.max() + ,-1 , -1 , -1] + + assert len(triggers) == 1 + for i, element in enumerate(elements): + assert element == triggers[0][i] + + +@given(evt_id=integers(1, 10), pmt_id=integers(1, 10), pulse_valid=integers(0, 9)) +def test_get_trigger_candidates_double(double_square_wf, evt_id, pmt_id + ,channel_conf , trigger_conf , pulse_valid): + # Check that the square pulses are identified as separated when pulse valid ext is lower than their separation. + channel_config = channel_conf.copy() + trigger_config = trigger_conf.copy() + + wf = double_square_wf[0] + thr = double_square_wf[1] + length = double_square_wf[2] + start1 = double_square_wf[3] + start2 = double_square_wf[4] + channel_config = channel_conf.copy() + channel_config['pulse_valid_ext'] = pulse_valid + triggers = get_trigger_candidates(wf, evt_id, pmt_id + ,channel_config, trigger_config + ,np.zeros(wf.shape), wf.max()) + assert len(triggers) == 2 + for j, start in enumerate([start1, start2]): + peak = wf[start:start+length+pulse_valid] + t_true = start + q_true = peak[:-2].sum() + val_q = in_range(q_true + , channel_config['q_min'], channel_config['q_max'] + , left_closed=False, right_closed=False) + w_true = len(peak) + val_w = in_range(w_true + , channel_config['time_min'], channel_config['time_max'] + , left_closed=True , right_closed=False) + h_true = peak.max() + val_h = h_true < channel_config['amp_max'] + val_p = True + + elements = [evt_id , pmt_id + ,t_true , q_true , w_true, h_true + ,val_q , val_w , val_h , val_p + ,wf[start1-1], wf.max() + ,-1 , -1 , -1] + + for i, element in enumerate(elements): + assert element == triggers[j][i] + + +@given(evt_id=integers(1, 10), pmt_id=integers(1, 10), discard=integers(0, 30)) +def test_get_trigger_candidates_discard(double_square_wf, evt_id, pmt_id + ,channel_conf , trigger_conf , discard): + # Check that the small pulse is kept when discard_width is low enough + channel_config = channel_conf.copy() + trigger_config = trigger_conf.copy() + + wf = double_square_wf[0] + thr = double_square_wf[1] + length = double_square_wf[2] + start1 = double_square_wf[3] + start2 = double_square_wf[4] + + trigger_config['discard_width'] = discard + triggers = get_trigger_candidates(wf, evt_id, pmt_id + ,channel_config, trigger_config + ,np.zeros(wf.shape), wf.max()) + assert len(triggers) == 2 + for j in range(2): + if j == 0: + peak = wf[1500:1520+channel_config['pulse_valid_ext']] + t_true = 1500 + else: + peak = wf[start1:start2+length+channel_config['pulse_valid_ext']] + t_true = start1 + q_true = peak[:-2].sum() + val_q = in_range(q_true + , channel_config['q_min'], channel_config['q_max'] + , left_closed=False, right_closed=False) + w_true = len(peak) + val_w = in_range(w_true + , channel_config['time_min'], channel_config['time_max'] + , left_closed=True , right_closed=False) + h_true = peak.max() + val_h = h_true < channel_config['amp_max'] + val_p = True + + elements = [evt_id , pmt_id + ,t_true , q_true , w_true, h_true + ,val_q , val_w , val_h , val_p + ,wf[start1-1], wf.max() + ,-1 , -1 , -1] + + for i, element in enumerate(elements): + assert element == triggers[j][i] + + +@given(evt_id=integers(1, 10), pmt_id=integers(1, 10)) +def test_get_trigger_candidates_multipeak(double_square_wf, evt_id, pmt_id + ,channel_conf , trigger_conf): + # Check that the multipeak protection discards the trigger if conditions are met. + channel_config = channel_conf.copy() + trigger_config = trigger_conf.copy() + + wf = double_square_wf[0] + thr = double_square_wf[1] + length = double_square_wf[2] + start1 = double_square_wf[3] + start2 = double_square_wf[4] + channel_config = channel_conf.copy() + channel_config['pulse_valid_ext'] = 0 + triggers = get_trigger_candidates(wf, evt_id, pmt_id + ,channel_config, trigger_config + ,np.zeros(wf.shape), wf.max()) + for t in triggers: + assert t[-6] == True + + ### Second peak characteristics + peak = wf[start2:start2+length+channel_config['pulse_valid_ext']] + t_true = start2 + q_true = peak[:-2].sum() + val_q = in_range(q_true + , channel_config['q_min'], channel_config['q_max'] + , left_closed=False, right_closed=False) + w_true = len(peak) + val_w = in_range(w_true + , channel_config['time_min'], channel_config['time_max'] + , left_closed=True , right_closed=False) + h_true = peak.max() + val_h = h_true < channel_config['amp_max'] + val_p = True + + # Multipeak condition metretrieve_trigger_information + trigger_config['multipeak'] = {'q_min':q_true//2, 'time_min':w_true//2, 'time_after':length+10} + triggers = get_trigger_candidates(wf, evt_id, pmt_id + ,channel_config, trigger_config + ,np.zeros(wf.shape), wf.max()) + + assert triggers[0][-6] == False + assert triggers[1][-6] == True + + # Multipeak ondition no longer met + trigger_config['multipeak'] = {'q_min':q_true*2, 'time_min':w_true*2, 'time_after':length+10} + triggers = get_trigger_candidates(wf, evt_id, pmt_id + ,channel_config, trigger_config + ,np.zeros(wf.shape), wf.max()) + + for t in triggers: + assert t[-6] == True + + +@given(evt_id=integers(0, 3)) +def test_retrieve_trigger_information(double_square_wf, evt_id + ,channel_conf , trigger_conf): + channel_config = channel_conf.copy() + trigger_config = trigger_conf.copy() + + retriever = retrieve_trigger_information({0:channel_config, 1:channel_config}, trigger_config) + + wf = double_square_wf[0] + wfs = np.array([i*wf for i in range(1, 4)]) + thr = double_square_wf[1] + length = double_square_wf[2] + start1 = double_square_wf[3] + start2 = double_square_wf[4] + + triggers = retriever(-wfs, wfs[:2], np.zeros(wfs.shape), evt_id) + + assert len(triggers) == 2 + for j in range(2): + peak = wfs[j][start1:start2+length+channel_config['pulse_valid_ext']] + t_true = start1 + q_true = peak[:-2].sum() + val_q = in_range(q_true + , channel_config['q_min'], channel_config['q_max'] + , left_closed=False, right_closed=False) + w_true = len(peak) + val_w = in_range(w_true + , channel_config['time_min'], channel_config['time_max'] + , left_closed=True , right_closed=False) + h_true = peak.max() + val_h = h_true < channel_config['amp_max'] + val_p = True + + elements = [evt_id , j + ,t_true , q_true , w_true, h_true + ,val_q , val_w , val_h , val_p + ,wf[start1-1], -wfs[j].max() + ,-1 , -1 , -1] + + for i, element in enumerate(elements): + assert element == triggers[j][i] + + + +@given(trigger_time_and_validity(), integers()) +def test_check_trigger_coincidence(trigger_info, coincidence_window): + trigger_time = trigger_info[0] + validity = trigger_info[1] + trigg = np.zeros((len(validity), 15)) + + # Initiate sim. trigger values + trigg[:, 1 ] = list(range(len(validity))) + trigg[:, 2 ] = trigger_time + trigg[:, 6:10] = 1#validity + trigg[:, -3: ] = np.full((len(validity), 3), -1) + + order_idx = np.lexsort((trigg[:,1], trigger_time)) + old_order = np.argsort(order_idx) + + ordered_times = trigg[:, 2 ][order_idx] + ordered_val = trigg[:, 6:10][order_idx] + ordered_ids = trigg[:, 1 ][order_idx] + trigger_back = trigg[order_idx].copy() + + for i, time1 in enumerate(ordered_times): + if not np.all(ordered_val[i]): continue + trigger_back[i, -3] = 1 + nearest = False + for val2, id2, time2 in zip(ordered_val[i+1:], ordered_ids[i+1:], ordered_times[i+1:]): + if not np.all(val2): continue + if not nearest: + trigger_back[i, -2] = time2-time1 + trigger_back[i, -1] = id2 + nearest = True + if (time2-time1) Date: Fri, 28 May 2021 11:37:29 +0200 Subject: [PATCH 04/14] Add blr function imitating the full procedure that happens in FPGA A BLR function is added that computes the baseline as a moving average. It applies the high-pass filter after substracting such baseline (and thus it's done within the loop) and applies deconvolution. As happens in the FPGA, the returned signal is equal to the raw signal whenever the BLR is not active. --- invisible_cities/cities/components.py | 33 +++++- invisible_cities/sierpe/blr.pxd | 8 +- invisible_cities/sierpe/blr.pyx | 138 +++++++++++++++++++++++++- 3 files changed, 173 insertions(+), 6 deletions(-) diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py index 3d18490bc..8a55f4c8f 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -125,7 +125,7 @@ def create_timestamp(rate: float) -> float: Returns ------- - Function to calculate timestamp for the given rate with + Function to calculate timestamp for the given rate with event_number as parameter. """ @@ -142,12 +142,12 @@ def create_timestamp(rate: float) -> float: def create_timestamp_(event_number: Union[int, float]) -> float: """ Calculates timestamp for a given Event Number and Rate. - + Parameters ---------- event_number : Union[int, float] ID value of the current event. - + Returns ------- Calculated timestamp : float @@ -326,6 +326,33 @@ def deconv_pmt(RWF): return deconv_pmt +def deconv_pmt_fpga(dbfile : str + ,run_number : int + ,selection : List[int] = None + ) -> Callable: + ''' + Apply deconvolution as done in the FPGA. + + Parameters + ---------- + dbfile : Database to be used + run_number : Run number of the database + selection : List of PMT IDs to apply deconvolution to. + + Returns + ---------- + deconv_pmt : Function that will apply the deconvolution. + ''' + DataPMT = load_db.DataPMT(dbfile, run_number = run_number) + pmt_active = np.nonzero(DataPMT.Active.values)[0].tolist() if selection is None else selection + coeff_c = DataPMT.coeff_c .values.astype(np.double)[pmt_active] + coeff_blr = DataPMT.coeff_blr.values.astype(np.double)[pmt_active] + def deconv_pmt(RWF): + return zip(*map(blr.deconvolve_signal_fpga, RWF[pmt_active], + coeff_c , coeff_blr )) + return deconv_pmt + + def get_run_number(h5in): if "runInfo" in h5in.root.Run: return h5in.root.Run.runInfo[0]['run_number'] elif "RunInfo" in h5in.root.Run: return h5in.root.Run.RunInfo[0]['run_number'] diff --git a/invisible_cities/sierpe/blr.pxd b/invisible_cities/sierpe/blr.pxd index ad6fbd15c..76710c354 100644 --- a/invisible_cities/sierpe/blr.pxd +++ b/invisible_cities/sierpe/blr.pxd @@ -32,4 +32,10 @@ cpdef deconvolve_signal(double [:] signal_daq, double coeff_clean = *, double coeff_blr = *, double thr_trigger = *, - int accum_discharge_length = *) \ No newline at end of file + int accum_discharge_length = *) + +cpdef deconvolve_signal_fpga(short [:] signal_daq + ,double coeff_clean = * + ,double coeff_blr = * + ,double thr_trigger = * + ,size_t base_window = *) diff --git a/invisible_cities/sierpe/blr.pyx b/invisible_cities/sierpe/blr.pyx index 3db786312..279da44f9 100644 --- a/invisible_cities/sierpe/blr.pyx +++ b/invisible_cities/sierpe/blr.pyx @@ -1,3 +1,4 @@ +#cython: language_level=3 import numpy as np cimport numpy as np from scipy import signal as SGN @@ -59,7 +60,6 @@ cpdef deconvolve_signal(double [:] signal_daq, if (signal_daq[k] < trigger_line) and (acum[k-1] < thr_acum): # discharge accumulator - if acum[k-1] > 1: acum[k] = acum[k-1] * (1 - coef) if j < accum_discharge_length - 1: @@ -70,4 +70,138 @@ cpdef deconvolve_signal(double [:] signal_daq, acum[k] = 0 j = 0 # return recovered signal - return np.asarray(signal_r) \ No newline at end of file + return np.asarray(signal_r) + + +cpdef deconvolve_signal_fpga( short [:] signal_daq + , double coeff_clean = 2.905447E-06 + , double coeff_blr = 1.632411E-03 + , double thr_trigger = 5 + , size_t base_window = 1024 ): + + """ + Simulate the deconvolution process as in the daq, differences compared to + usual offline deconvolution: + - Baseline is calculated as a moving average of 1024 counts (FPGA). + - Flips result between raw and deconvolved signal when outside of both + signal and discharge regions. + - Discharge made at fixed value (0.995) + - Result is truncated to integers. + - ADC threshold acting as absolute threshold. + + Parameters + ---------- + signal_daq : short array + PMT raw waveform + coeff_clean : double + Characteristic parameter of the high pass filter + coeff_blr : double + Characteristic parameter of BLR + thr_trigger : double + Threshold in ADCs to activate BLR + base_window : size + Moving average window for baseline calculation + + Returns + ---------- + Tuple of arrays: + - Deconvolved waveform (int16 [:]) + - Baseline (int16 [:]) + """ + cdef double coef = coeff_blr + cdef double thr_acum = thr_trigger / coef + cdef int len_signal_daq = len(signal_daq) + + cdef double [:] signal_r = np.zeros(len_signal_daq, dtype=np.double) + cdef double [:] signal_f = np.zeros(len_signal_daq, dtype=np.double) + cdef double [:] acum = np.zeros(len_signal_daq, dtype=np.double) + + cdef int j + + # compute noise + cdef double noise = 0 + cdef int nn = 400 # fixed at 10 mus + + for j in range(nn): + noise += signal_daq[j] * signal_daq[j] + noise /= nn + cdef double noise_rms = np.sqrt(noise) + + # trigger line + cdef double trigger_line = thr_trigger #* noise_rms + # cleaning signal + cdef double [:] b_cf + cdef double [:] a_cf + + ### Baseline related variables + cdef double [:] top = np.zeros(len_signal_daq, dtype=np.double) + cdef short [:] aux = np.zeros(len_signal_daq, dtype=np.int16) + cdef long ped = np.sum (signal_daq[0:base_window], dtype=np.int32) + cdef size_t delay = 2 # Delay in the FPGA in the bin used for baseline substraction + cdef unsigned int iaux = base_window + + top[0:base_window] = ped/base_window + aux[0:base_window] = signal_daq[0:base_window] + + b_cf, a_cf = SGN.butter(1, coeff_clean, 'high', analog=False); + g, a1 = b_cf[0], -a_cf[-1] + + ### Initiate filt + filt_hr = 0 + #1st + filt_x = g * -(signal_daq[0] - top[0]) + filt_a1hr = a1 * filt_hr + + filt_h = filt_x + filt_a1hr + signal_f[0] = filt_h - filt_hr + filt_hr = filt_h + + #2nd + filt_x = g * -(signal_daq[1] - top[0]) + filt_a1hr = a1 * filt_hr + + filt_h = filt_x + filt_a1hr + signal_f[1] = filt_h - filt_hr + filt_hr = filt_h + + + cdef int k + + #Initiate BLR + for k in range(0, delay): + signal_r[k] = signal_daq[k] + + for k in range(2, len_signal_daq): + # always update signal and accumulator + current = signal_daq[k] + baseline = top[k-delay] + + ### High-pass filter + filt_x = g * -(current - baseline) + filt_a1hr = a1 * filt_hr + + filt_h = filt_x + filt_a1hr + signal_f[k] = filt_h - filt_hr + + ### BLR restoration + blr_val = (signal_f[k] + signal_f[k]*(coef / 2) + + coef * acum[k-1]) + baseline + + acum[k] = acum[k-1] + signal_f[k] + + if (signal_f[k] < trigger_line) and (acum[k-1] < thr_acum): + # discharge accumulator + if acum[k-1] > 1: + acum[k] = acum[k-1] * .995 #Fixed discharge in FPGA code + else: + acum [k] = 0 + blr_val = current # When outside of BLR/discharge, flip to raw signal + if k>=base_window: + ped = ped + current - aux[iaux-base_window] + aux[iaux] = current + iaux += 1 + + signal_r[k] = blr_val + top[k] = ped/base_window + + return np.asarray(signal_r, dtype='int16'), np.asarray(top, dtype='int16') From 6920f75bb4405b6147fe8cf7d897c0d63af3f42a Mon Sep 17 00:00:00 2001 From: Ander Date: Fri, 28 May 2021 11:40:55 +0200 Subject: [PATCH 05/14] Add trigger functions. These functions extract the trigger parameters and check for coincidences between trigger candidates. --- invisible_cities/reco/trigger_functions.py | 174 +++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 invisible_cities/reco/trigger_functions.py diff --git a/invisible_cities/reco/trigger_functions.py b/invisible_cities/reco/trigger_functions.py new file mode 100644 index 000000000..ad79c471b --- /dev/null +++ b/invisible_cities/reco/trigger_functions.py @@ -0,0 +1,174 @@ +import numpy as np + +from typing import List +from typing import Dict +from typing import Callable + +from .. core.core_functions import in_range + +def get_trigger_candidates(signal : np.ndarray + ,event_id : int , pmt_id : int + ,channel_conf : Dict[int, int], trigger_conf : Dict + ,baseline : np.ndarray , wvf_max_height : int + ) -> List: + """ + Calculates the trigger information of a given signal. + + Parameters + ---------- + signal : Deconvolved waveform. + event_id : Event ID associated to the waveform. + pmt_id : PMT ID of the waveform's PMT. + channel_conf : PMT specific trigger configuration. + trigger_conf : General trigger configuration. + baseline : Array with the baseline value at each time bin. + wvf_max_height : Maximum height of the waveform. + + Returns + ---------- + triggers : List with all the trigger related information. + """ + # Trigger parameters + qrange = [channel_conf['q_min' ], channel_conf['q_max' ]] + wrange = [channel_conf['time_min' ], channel_conf['time_max']] + hrange = [channel_conf['baseline_dev' ], channel_conf['amp_max' ]] + pulse_ext = channel_conf['pulse_valid_ext'] + + fMulti = trigger_conf['multipeak' ] # Protection for later peaks + discard_width = trigger_conf['discard_width'] # Don't store events below certain width + + # Multipeak protection parameters + if fMulti: + multi_q = fMulti['q_min' ] # Minimum charge of later peaks to discard a trigger + multi_w = int(fMulti['time_min' ]) # Minimum time width of later peaks to discard a trigger + multi_t = int(fMulti['time_after']) # Time window for a later peak to be considered + + # Signal above threshold, allowing for pulse_ext counts below threshold + triggers = [] + indexes = np.where(signal > hrange[0])[0] + if len(indexes) == 0: return triggers # Protection in case min. threshold is not surpassed + index_sel = np.diff(indexes)>(pulse_ext+1) + candidates = np.split(indexes, np.where(index_sel)[0] + 1) + + for p in candidates: + peak_slice = slice(p[0]-1, p[-1] + pulse_ext + 1) # DAQ gets pulse_ext counts after peaks; also, for Q it takes a bin before TODO Fix for N100 (remove -1) + peak = signal[peak_slice] + if not (len(peak) > max(discard_width, 1)): continue + width = len(peak) - 1 # TODO Remove -1 in N100 + height = peak[1:].max() # TODO Full peak for N100 + peak = peak[:-2][peak[:-2]>0] # Only counts above 0 are used for Q, last bin not considered TODO ONLY -1 for N100 + charge = peak.sum() + if len(peak)==0: continue + start_time = p[0] + + # Trigger validity + trigger_w = in_range(width , *wrange, left_closed=True , right_closed=False) + trigger_q = in_range(charge, *qrange, left_closed=False, right_closed=False) + trigger_h = height < hrange[1] + + # Baseline at the start of the pulse + baseline_t0 = baseline[start_time] + triggers .append([event_id , pmt_id + ,start_time , charge , width , height + ,trigger_q , trigger_w , trigger_h, True # By default multi flag is set as true + ,baseline_t0, wvf_max_height + ,-1 , -1 , -1]) # Filled later in the coincidence module + + # Extra-peaks protection + if fMulti is not None: + # Evaluate if later peaks (must be over discard_width) are candidates for the multipeak protection + for i, t0 in zip(range(len(triggers[:-1])), triggers[:-1]): + for t1 in triggers[i+1:]: + multi_candidate = ((t1[4] >= multi_w ) & + (t1[3] > multi_q ) & + (t1[2] <= t0[2] + t0[4] + multi_t+2)) # Start of trigger + peak width + time window + FPGA delay (2) + if not multi_candidate: + continue + # In case the later peak ends after the multipeak protection window ends + if (t1[2]+t1[4]) > (t0[2] + t0[4] + multi_t+2): # Event partially contained in the multipeak window + peak_slice = slice(t1[2]-1, t0[2]+t0[4]+multi_t+2+1) + peak = signal[peak_slice] + width = len(peak) - 1 # TODO Remove -1 in N100 + charge = peak[peak>0].sum() + + triggers[i][-6] = (width >= multi_w) & (charge > multi_q) + else: + triggers[i][-6] = False + if not triggers[i][-6]: break + + return triggers + + +def retrieve_trigger_information(channel_config : Dict[int, int] + ,trigger_config : Dict + ) -> Callable: + """ + Calculates the trigger information of deconvolved waveforms. + + Parameters + ---------- + channel_config : Channel-specific trigger configuration. + trigger_config : General trigger configuration. + + Returns + ---------- + trigger_on_channels : Function that returns each channel's trigger candidates. + """ + def trigger_on_channels(rwfs : np.ndarray + ,cwfs : np.ndarray + ,baselines : np.ndarray + ,event_number : int): + triggers = [] + for i, dict_items in enumerate(channel_config.items()): + pmt_id, pmt_conf = dict_items + baseline = baselines[i] + signal = cwfs[i] - baseline + max_height = rwfs[pmt_id].min() + triggers.extend(get_trigger_candidates(signal + ,event_number, pmt_id + ,pmt_conf , trigger_config + ,baseline , max_height )) + return np.array(triggers) + return trigger_on_channels + + +def check_trigger_coincidence(coinc_window : int) -> Callable: + """ + Checks if there is time coincidence between trigger candidates and modifies + the input list of trigger candidates accordingly. + + Parameters + ---------- + coinc_window : Time window to considerate the coincidence condition as valid. + + Returns + ---------- + pmt_coincidence : Function that checks for coincidence and modififes triggers. + """ + + def pmt_coincidence(triggers : np.ndarray) -> np.ndarray: + valid_sel = np.all(triggers[:, 6:10], axis=1).astype('bool') # Use only valid triggers + times = triggers[:, 2][valid_sel] # Times of valid triggers + ids = triggers[:, 1][valid_sel] # Ids of valid triggers + + # Order by time, index so we can use this to order both times and ids + order_idx = np.lexsort((ids, times)) + order_time = times[order_idx] + pmts = np.full(ids.shape, -1) + # Calculate all time differences between each trigger and all later ones + time_diffs = np.array([order_time[i+1:] - t for i, t in enumerate(order_time)]) + # Since pmts are ordered by time, closest pmt will be the next one in the array. + pmts[:-1] = ids[order_idx][1:] + # Shortest time is the first value of time_diff for each trigger + shortest = np.array([t[0] if len(t)>0 else -1 for t in time_diffs]) + # Count how many triggers are within the coincidince window of each trigger + n_coinc = np.array([len(t[t Date: Fri, 28 May 2021 11:43:32 +0200 Subject: [PATCH 06/14] Add TriggerTable class Add class that contains the simualted trigger info table. --- invisible_cities/evm/nh5.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/invisible_cities/evm/nh5.py b/invisible_cities/evm/nh5.py index 8d1be863f..4c502f30d 100644 --- a/invisible_cities/evm/nh5.py +++ b/invisible_cities/evm/nh5.py @@ -186,3 +186,22 @@ class VoxelsTable(tb.IsDescription): class EventPassedFilter(tb.IsDescription): event = tb.Int32Col(pos=0) passed = tb. BoolCol(pos=1) + + +class TriggerTable(tb.IsDescription): + event = tb.UInt32Col(pos= 0) + pmt = tb.UInt32Col(pos= 1) + trigger_time = tb.UInt32Col(pos= 2) + q = tb.UInt32Col(pos= 3) + width = tb.UInt32Col(pos= 4) + height = tb.UInt32Col(pos= 5) + valid_q = tb.BoolCol (pos= 6) + valid_w = tb.BoolCol (pos= 7) + valid_h = tb.BoolCol (pos= 8) + valid_peak = tb.BoolCol (pos= 9) + valid_all = tb.BoolCol (pos=10) + baseline = tb.Int32Col (pos=11) + max_height = tb.Int32Col (pos=12) + n_coinc = tb.Int32Col (pos=13) + closest_ttime = tb.Int32Col (pos=14) + closest_pmt = tb.Int32Col (pos=15) From bf1a78c4d75e7aeea22b4fa4674436fdf18e3403 Mon Sep 17 00:00:00 2001 From: Ander Date: Fri, 28 May 2021 11:45:50 +0200 Subject: [PATCH 07/14] Add trigger info writer Adds the writer for the simulated trigger dataframe. --- invisible_cities/io/trigger_io.py | 39 +++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/invisible_cities/io/trigger_io.py b/invisible_cities/io/trigger_io.py index 8f077e191..752eb0963 100644 --- a/invisible_cities/io/trigger_io.py +++ b/invisible_cities/io/trigger_io.py @@ -4,6 +4,7 @@ from .. evm import nh5 as table_formats from .. reco.tbl_functions import filters as tbl_filters +from .. io .table_io import make_table def store_trigger(tables, trg_type, trg_channels): @@ -41,3 +42,41 @@ def _make_tables(hdf5_file, n_sensors, compression="ZLIB4"): trg_tables = trg_type, trg_channels return trg_tables + + +def trigger_dst_writer(hdf5_file, **kwargs):#{ + trigger_table = make_table(hdf5_file, + group = "Trigger", + name = "DST" , + fformat = table_formats.TriggerTable, + description = "Simulated trigger data", + compression = 'ZLIB4') + def write_trigger(trigger_info):#{ + [event , pmt , + trigger_time , charge , width , height , + valid_q , valid_w , valid_h , valid_peak, + mean_baseline, max_height , + n_coinc , closest_ttime, closest_pmt] = trigger_info + row = trigger_table.row + row["event" ] = event + row["pmt" ] = pmt + row["trigger_time" ] = trigger_time + row["q" ] = charge + row["width" ] = width + row["height" ] = height + row["valid_q" ] = valid_q + row["valid_w" ] = valid_w + row["valid_h" ] = valid_h + row["valid_peak" ] = valid_peak + row["valid_all" ] = (valid_q + valid_w + valid_h + valid_peak) == 4 + row["baseline" ] = mean_baseline + row["max_height" ] = max_height + row["n_coinc" ] = n_coinc + row["closest_ttime"] = closest_ttime + row["closest_pmt" ] = closest_pmt + row.append() + + def write_triggers(triggers): + for t in triggers: write_trigger(t) + + return write_triggers From 0c46d2622758033f8f1e4e3e0662bc6178881957 Mon Sep 17 00:00:00 2001 From: Ander Date: Fri, 28 May 2021 11:47:19 +0200 Subject: [PATCH 08/14] Add valdrada city Adds valdrada, this city takes rwf and produces a dataframe with the simulated trigger information. It takes a configurable set of parameters that imitate the detector acquisition ones and their naming convention. It provides all possible trigger candidates above a given width and flags them as valid or not based on the configuration. --- invisible_cities/cities/valdrada.py | 172 +++++++++++++++++++++++ invisible_cities/cities/valdrada_test.py | 4 +- invisible_cities/conftest.py | 3 - 3 files changed, 174 insertions(+), 5 deletions(-) create mode 100644 invisible_cities/cities/valdrada.py diff --git a/invisible_cities/cities/valdrada.py b/invisible_cities/cities/valdrada.py new file mode 100644 index 000000000..396a35bcc --- /dev/null +++ b/invisible_cities/cities/valdrada.py @@ -0,0 +1,172 @@ +""" +----------------------------------------------------------------------- + valdrada +----------------------------------------------------------------------- + +This city finds simulates the trigger procedure at the FPGA over PMT waveforms. +This includes a number of tasks: + - Use a moving average to obtain the baseline + - Truncate to integers the deconvolution output + - Find trigger candidates and compute their characteristics, taking into + account the unique circumstances derived from working on a continuous scheme + (for example, delay in signals) + - Evaluate coincidences between trigger candidates. + +A number of general configuration parameters are needed (input as a dict, trigger_config): +- coincidence_window : Time window (in time bins) in which valid trigger coincidences are counted +- discard_width : Any trigger with width less than this parameter is discarded. +- multipeak : Dictionary with multipeak protection parameters, otherwise None. + - q_min : Integrated charge threshold of a post-trigger peak to discard + a trigger due to multipeak protection. In ADC counts. + - time_min : Minimum width of a post-trigger peak to discard a trigger due to + multipeak protection. In time bins. + - time_after : For how long is the multipeak protection evaluated after a + valid trigger. In mus. + +A individual trigger configuration can be given per channel, through a dict +with keys equal to PMT IDs, which marks the validity range of the peak +characteristics: +- q_min, q_max : Range for the integrated charge of the peak (q_min < q < q_max). + In ADC counts. +- time_min, time_max : Range for the peak width (time_min <= width < time_max). + In time mus. +- baseline_dev, amp_max : Range for peak height (baseline_dev < height < amp_max). + In ADC counts. +- pulse_valid_ext : Time allowed for a pulse to go below baseline_dev. In ns. + +The result of the city is a dataframe containing the event ID and PMT ID of each +trigger candidate. For each trigger candidate a number of parameters is computed: + - trigger_time : time bin at which the trigger candidate starts. + - q : integrated ADC counts within the peak. + - width : Length of the peak in time bins. + - height : Maximum ADC values in a given time bien within the peak. + - baseline : Value of the baseline at trigger_time. + - max_height : Maximum (minimum due to PMT negative signals) height in the wvf. + +Additionally, a set of flags are assigned depending on wether the parameters are +within range of the trigger configuration: + - valid_q : If peak is within q range. + - valid_w : If peak is within width range. + - valid_h : If peak is within height range. + - valid_peak : Only if 'multipeak' is active, True if there isn't a post-trigger + candidate with the configuration parametersself. + - valid_all : boolean and of previous valid flags. + + Finally, a series of coincidence-related values are also given: + - n_coinc : Number of valid triggers within the coincidence_window, + starting at the trigger trigger_time and including the trigger itself. + -1 indicates no valid trigger (including the trigger itself). + - closest_ttime : Time difference to closest valid trigger. + -1 if there are none aside from the trigger itself. + - closest_pmt : PMT ID of the closest valid trigger. + -1 if there are none aside from the trigger itself. +""" + +import tables as tb +import numpy as np + +from .. reco import tbl_functions as tbl +from .. io .run_and_event_io import run_and_event_writer +from .. io .trigger_io import trigger_writer, trigger_dst_writer + +from .. dataflow import dataflow as fl +from .. dataflow.dataflow import push +from .. dataflow.dataflow import pipe +from .. dataflow.dataflow import sink + +from . components import city +from . components import print_every +from . components import collect +from . components import copy_mc_info +from . components import deconv_pmt_fpga +from . components import WfType +from . components import wf_from_files +from . components import get_number_of_active_pmts + +from .. reco.trigger_functions import retrieve_trigger_information +from .. reco.trigger_functions import get_trigger_candidates +from .. reco.trigger_functions import check_trigger_coincidence + +from .. core import system_of_units as units + + +def check_empty_trigger(triggers) -> bool: + """ + Filter for valdrada flow. + The flow stops if there are no candidate triggers/ + """ + return len(triggers) > 0 + +@city +def valdrada(files_in, file_out, compression, event_range, print_mod, detector_db, run_number, + trigger_config = dict(), channel_config = dict()): + + # Change time units into number of waveform bins. + if trigger_config['multipeak'] is not None: + trigger_config['multipeak']['time_min' ] /= 25*units.ns + trigger_config['multipeak']['time_after'] /= 25*units.ns + for k in channel_config.keys(): + channel_config[k]['time_min'] /= 25*units.ns + channel_config[k]['time_max'] /= 25*units.ns + channel_config[k]['pulse_valid_ext'] = int(channel_config[k]['pulse_valid_ext']/25*units.ns) + + #### Define data transformations + # Raw WaveForm to deconvolved waveforms + rwf_to_cwf = fl.map(deconv_pmt_fpga(detector_db, run_number, list(channel_config.keys())), + args = "pmt", + out = ("cwf", "baseline")) + + # Extract all possible trigger candidates of each PMT + trigger_on_channels = fl.map(retrieve_trigger_information(channel_config, trigger_config), + args = ("pmt", "cwf", "baseline", "event_number"), + out = "triggers") + + # Add coincidence between channels + check_coincidences = fl.map(check_trigger_coincidence(trigger_config['coincidence_window']), + item = "triggers") + + # Filter events with zero triggers + filter_empty_trigger = fl.map(check_empty_trigger, + args = "triggers", + out = "empty_trigger") + + event_count_in = fl.spy_count() + event_count_out = fl.spy_count() + events_no_triggers = fl.count_filter(bool, args = "empty_trigger") + evtnum_collect = collect() + + with tb.open_file(file_out, "w", filters = tbl.filters(compression)) as h5out: + # Define writers... + write_event_info_ = run_and_event_writer(h5out) + write_trigger_info_ = trigger_writer (h5out, get_number_of_active_pmts(detector_db, run_number)) + write_trigger_dst_ = trigger_dst_writer (h5out) + # ... and make them sinks + + write_event_info = sink(write_event_info_ , args=( "run_number" , "event_number" , "timestamp" )) + write_trigger_info = sink(write_trigger_info_, args=( "trigger_type", "trigger_channels" )) + write_trigger_dst = sink(write_trigger_dst_ , args= "triggers" ) + + result = push(source = wf_from_files(files_in, WfType.rwf), + pipe = pipe(fl.slice(*event_range, close_all=True), + print_every(print_mod), + event_count_in.spy, + rwf_to_cwf, + trigger_on_channels, + filter_empty_trigger, + events_no_triggers.filter, + check_coincidences, + event_count_out.spy, + fl.branch("event_number", evtnum_collect.sink), + fl.fork(write_trigger_dst , + write_event_info , + write_trigger_info)), + result = dict(events_in = event_count_in .future, + events_out = event_count_out .future, + evtnum_list = evtnum_collect .future, + events_pass = events_no_triggers.future)) + + if run_number <= 0: + copy_mc_info(files_in, h5out, result.evtnum_list, + detector_db, run_number) + + return result diff --git a/invisible_cities/cities/valdrada_test.py b/invisible_cities/cities/valdrada_test.py index 586643fa9..d6da771cb 100644 --- a/invisible_cities/cities/valdrada_test.py +++ b/invisible_cities/cities/valdrada_test.py @@ -19,7 +19,7 @@ def test_valdrada_contains_all_tables(trigger_config): def test_valdrada_exact_result_multipeak(ICDATADIR, trigger_config): - true_out = os.path.join(ICDATADIR, "test_trigger_fpga_multipeak.h5") + true_out = os.path.join(ICDATADIR, "exact_result_multipeak_valdrada.h5") conf, PATH_OUT = trigger_config valdrada(**conf) @@ -36,7 +36,7 @@ def test_valdrada_exact_result_multipeak(ICDATADIR, trigger_config): def test_valdrada_exact_result(ICDATADIR, trigger_config): - true_out = os.path.join(ICDATADIR, "test_trigger_fpga.h5") + true_out = os.path.join(ICDATADIR, "exact_result_valdrada.h5") conf, PATH_OUT = trigger_config conf['trigger_config']['multipeak'] = None valdrada(**conf) diff --git a/invisible_cities/conftest.py b/invisible_cities/conftest.py index 1f99e9f21..fd81ed404 100644 --- a/invisible_cities/conftest.py +++ b/invisible_cities/conftest.py @@ -704,8 +704,6 @@ def deconvolution_config(ICDIR, ICDATADIR, PSFDIR, config_tmpdir): return conf, PATH_OUT -<<<<<<< HEAD -======= @pytest.fixture(scope='function') # Needs to be function as the config dict is modified when running def trigger_config(ICDIR, ICDATADIR, config_tmpdir): @@ -741,7 +739,6 @@ def trigger_config(ICDIR, ICDATADIR, config_tmpdir): return conf, PATH_OUT ->>>>>>> f8fa2343... fixup! Add related tests ## To make very slow tests only run with specific option def pytest_addoption(parser): parser.addoption( From 1f029d995890e4493ef96dde99ab4f40c84a4bc5 Mon Sep 17 00:00:00 2001 From: Ander Date: Tue, 1 Jun 2021 13:34:38 +0200 Subject: [PATCH 09/14] Description correction --- invisible_cities/cities/valdrada.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/invisible_cities/cities/valdrada.py b/invisible_cities/cities/valdrada.py index 396a35bcc..10e8ff842 100644 --- a/invisible_cities/cities/valdrada.py +++ b/invisible_cities/cities/valdrada.py @@ -49,7 +49,7 @@ - valid_w : If peak is within width range. - valid_h : If peak is within height range. - valid_peak : Only if 'multipeak' is active, True if there isn't a post-trigger - candidate with the configuration parametersself. + candidate with the configuration parameters. - valid_all : boolean and of previous valid flags. Finally, a series of coincidence-related values are also given: From 450d88fa1b0c05d5647c52f4f3a97e24e715cce2 Mon Sep 17 00:00:00 2001 From: Ander Date: Tue, 22 Jun 2021 17:53:51 +0200 Subject: [PATCH 10/14] Modified tests according to PR --- invisible_cities/sierpe/blr_test.py | 35 +++++++++++++++-------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/invisible_cities/sierpe/blr_test.py b/invisible_cities/sierpe/blr_test.py index 641573a6f..5bc6c9351 100644 --- a/invisible_cities/sierpe/blr_test.py +++ b/invisible_cities/sierpe/blr_test.py @@ -4,17 +4,20 @@ import tables as tb import pandas as pd -from pytest import fixture -from pytest import mark -from flaky import flaky +from pytest import fixture +from pytest import mark +from flaky import flaky +from hypothesis import given +from hypothesis import settings +from hypothesis.strategies import integers from .. reco import calib_sensors_functions as csf from . import blr deconv_params = namedtuple("deconv_params", - "coeff_clean coeff_blr " - "thr_trigger accum_discharge_length") + "coeff_clean coeff_blr " + "thr_trigger accum_discharge_length") deconv_fpga_params = namedtuple("deconv_params", "coeff_clean coeff_blr " "thr_trigger base_window") @@ -163,13 +166,12 @@ def test_deconv_pmt_ad_hoc_signals_dead_sensors(ad_hoc_blr_signals): np.allclose(blr_wfs, evt_true_blr_wfs[pmt_active]) +@settings(max_examples=1, deadline=1000) +@given (evt_no=integers(0, 9)) @mark.slow -def test_deconv_pmt_fpga_ad_hoc_signals_all(ad_hoc_blr_fpga_signals): +def test_deconv_pmt_fpga_ad_hoc_signals_all(evt_no, ad_hoc_blr_fpga_signals): all_rwfs, all_true_blr_wfs, all_true_baseline, params = ad_hoc_blr_fpga_signals - # This test takes long, so we pick a random event. - # Its exhaustiveness relies on repeated test runs. - evt_no = np.random.choice(all_rwfs.shape[0]) evt_rwfs = all_rwfs [evt_no] evt_true_blr_wfs = all_true_blr_wfs [evt_no] evt_true_baseline = all_true_baseline[evt_no] @@ -184,8 +186,10 @@ def test_deconv_pmt_fpga_ad_hoc_signals_all(ad_hoc_blr_fpga_signals): assert np.allclose(baseline, evt_true_baseline) +@settings(max_examples=1, deadline=1000) +@given (evt_no=integers(0, 9)) @mark.slow -def test_deconv_pmt_fpga_ad_hoc_signals_dead_sensors(ad_hoc_blr_fpga_signals): +def test_deconv_pmt_fpga_ad_hoc_signals_dead_sensors(evt_no, ad_hoc_blr_fpga_signals): all_rwfs, all_true_blr_wfs, all_true_baseline, params = ad_hoc_blr_fpga_signals n_evts, n_pmts, _ = all_rwfs.shape @@ -193,9 +197,6 @@ def test_deconv_pmt_fpga_ad_hoc_signals_dead_sensors(ad_hoc_blr_fpga_signals): n_alive = np.random.randint(1, n_pmts - 1) pmt_active = np.random.choice(pmt_active, size=n_alive, replace=False) - # This test takes long, so we pick a random event. - # Its exhaustiveness relies on repeated test runs. - evt_no = np.random.choice(all_rwfs.shape[0]) evt_rwfs = all_rwfs [evt_no] evt_true_blr_wfs = all_true_blr_wfs [evt_no] evt_true_baseline = all_true_baseline[evt_no] @@ -210,11 +211,11 @@ def test_deconv_pmt_fpga_ad_hoc_signals_dead_sensors(ad_hoc_blr_fpga_signals): assert np.allclose(baseline, evt_true_baseline[pmt_active]) -def test_moving_average_baseline(): +@given(thr=integers(2, 100), window=integers(10, 5000)) +def test_moving_average_baseline(thr, window): samples = 10000 - thr = np.random.randint(2, 100) - window = np.random.randint(10, samples//2) - wvf = np.random.randint(0, thr-1, samples, dtype=np.int16) + np.random.seed(1234) + wvf = np.random.randint(0 , thr-1, samples, dtype=np.int16) _, baseline = blr.deconvolve_signal_fpga(wvf, 1e-6, 1e-3, thr, window) From 1945ed1d786337d0381f1008bec04f67564e2bd2 Mon Sep 17 00:00:00 2001 From: Ander Date: Tue, 22 Jun 2021 17:54:27 +0200 Subject: [PATCH 11/14] Modified tests to accommodate for PR review --- invisible_cities/io/trigger_io_test.py | 72 ++++++++++++++++++-------- 1 file changed, 49 insertions(+), 23 deletions(-) diff --git a/invisible_cities/io/trigger_io_test.py b/invisible_cities/io/trigger_io_test.py index 635f1ac6c..d30efdf11 100644 --- a/invisible_cities/io/trigger_io_test.py +++ b/invisible_cities/io/trigger_io_test.py @@ -3,54 +3,80 @@ import numpy as np import tables as tb -from numpy.testing import assert_allclose +from numpy .testing import assert_allclose from hypothesis import given from hypothesis.strategies import integers from hypothesis.strategies import lists from hypothesis.strategies import booleans from hypothesis.strategies import composite -from ..core.testing_utils import assert_dataframes_equal -from ..io.dst_io import load_dst -from . trigger_io import trigger_dst_writer +from ..core.testing_utils import assert_dataframes_equal +from ..io.dst_io import load_dst +from . trigger_io import trigger_dst_writer +from .. reco.trigger_functions import TriggerInfo + + +@composite +def trigger_value(draw): + uints = integers(min_value=0 , max_value=np.iinfo(np.uint32).max) + ints = integers(min_value=np.iinfo(np.int32).min, max_value=np.iinfo(np.int32).max) + + # Event and PMT ID + evt_id = draw(uints) + pmt_id = draw(uints) + + # Peak time, q, width and height + time = draw(uints) + q = draw(uints) + width = draw(uints) + height = draw(uints) + + # Trigger validity flags + valid_q = draw(booleans()) + valid_w = draw(booleans()) + valid_h = draw(booleans()) + valid_peak = draw(booleans()) + valid_all = all([valid_q, valid_w, valid_h, valid_peak]) + + # Baseline and wvf max. height + baseline = draw(ints) + max_height = draw(ints) + + # Number of coincidences, closest pmt and tbin + n_coinc = draw(ints) + closest_time = draw(ints) + closest_pmt = draw(ints) + + trigger = [evt_id , pmt_id , time , q , width , height + ,valid_q , valid_w , valid_h, valid_peak , valid_all + ,baseline, max_height, n_coinc, closest_time, closest_pmt ] + + return TriggerInfo(*trigger) @composite def trigger_values(draw): size = draw(integers(min_value=1, max_value=10)) - elements = [] - for _ in range(6): - elements.append(draw(lists(integers(min_value=0, max_value=np.iinfo(np.uint32).max), - min_size=size, max_size=size))) - for _ in range(4): - elements.append(draw(lists(booleans(), min_size=size, max_size=size))) - for _ in range(5): - elements.append(draw(lists(integers(min_value=np.iinfo(np.int32).min, max_value=np.iinfo(np.int32).max), - min_size=size, max_size=size))) - return elements + return draw(lists(trigger_value(), min_size=size, max_size=size)) @given(triggers=trigger_values()) def test_trigger_dst_writer(config_tmpdir, triggers): output_file = os.path.join(config_tmpdir, "test_trigger_dst.h5") - trigger_save = [list(a) for a in zip(*triggers)] with tb.open_file(output_file, 'w') as h5out: write = trigger_dst_writer(h5out) - write(trigger_save) + write(triggers) trigger_dst = load_dst(output_file, group='Trigger', node='DST') col_names = ["event" , "pmt" , "trigger_time" , "q" - ,"width" , "height" , "valid_q" , "valid_w" + ,"width" , "height" , "valid_q" , "valid_w" ,"valid_h" , "valid_peak", "valid_all" , "baseline" ,"max_height", "n_coinc" , "closest_ttime", "closest_pmt"] assert np.all([tname == cname for tname, cname in zip(trigger_dst.columns.values, col_names)]) - # Check values stored, the if is needed because the valid_all condition is done in the writer, not on the input trigger - for i, colname in enumerate(col_names): - if colname == 'valid_all': - assert_allclose(np.all(triggers[6:10], axis=0), trigger_dst.loc[:, colname]) - elif i>10: assert_allclose(triggers[i-1] , trigger_dst.loc[:, colname]) - else : assert_allclose(triggers[i ] , trigger_dst.loc[:, colname]) + # Check values stored + for i, row in trigger_dst.iterrows(): + assert_allclose(triggers[i], list(row.values)) From 48f212c3f60bcbd89dd65a3b8cb6602117b09469 Mon Sep 17 00:00:00 2001 From: Ander Date: Tue, 22 Jun 2021 17:54:58 +0200 Subject: [PATCH 12/14] Modified tests to accommodate for PR review --- .../reco/trigger_functions_test.py | 133 +++++++++--------- 1 file changed, 67 insertions(+), 66 deletions(-) diff --git a/invisible_cities/reco/trigger_functions_test.py b/invisible_cities/reco/trigger_functions_test.py index 74d9893f4..35e68bd7f 100644 --- a/invisible_cities/reco/trigger_functions_test.py +++ b/invisible_cities/reco/trigger_functions_test.py @@ -12,40 +12,43 @@ from .. reco.trigger_functions import get_trigger_candidates from .. reco.trigger_functions import retrieve_trigger_information from .. reco.trigger_functions import check_trigger_coincidence +from .. reco.trigger_functions import TriggerInfo from .. core.core_functions import in_range -@fixture(scope="session") -def channel_conf(): +def channel_conf(**kwargs): channel_dict = {'q_min' : 5000, 'q_max' : 50000, 'time_min' : 2000, 'time_max': 10000, 'baseline_dev' : 5, 'amp_max' : 1000, 'pulse_valid_ext': 15} + channel_dict.update(kwargs) return channel_dict -@fixture(scope="session") -def trigger_conf(): +def trigger_conf(**kwargs): trigger_dict = {'coincidence_window':64, 'discard_width':40, 'multipeak':None} + trigger_dict.update(kwargs) return trigger_dict -@fixture(scope="session") -def double_square_wf(): +@composite +def double_square_wf(draw): # Fake signal with two-nearby square pulses and a smaller additional one at the beginning. n_baseline = 64000 base_window = 1024 thr = 5 wf = np.zeros(n_baseline) - start1 = np.random.randint(2*base_window, n_baseline // 2) - length = np.random.randint(n_baseline // 50, n_baseline // 4) + start1 = draw(integers(2*base_window , n_baseline // 2)) + length = draw(integers(n_baseline // 50, n_baseline // 4)) stop1 = start1 + length start2 = stop1 + 10 stop2 = start2 + length - wf[start1:stop1] = np.full(length, np.random.randint(thr+1, 50)) - wf[start2:stop2] = np.full(length, np.random.randint(thr+1, 50)) - wf[1500 : 1520] = np.full(20 , np.random.randint(thr+1, 50)) # Extra pulse to check discard_width + heights = draw(lists(integers(thr+1, 50), min_size=3, max_size=3)) + wf[start1:stop1] = np.full(length, heights[0]) + wf[start2:stop2] = np.full(length, heights[1]) + wf[1500 : 1520] = np.full(20 , heights[2]) # Extra pulse to check discard_width + return wf, thr, length, start1, start2 @@ -53,18 +56,18 @@ def double_square_wf(): def trigger_time_and_validity(draw): # Set of trigger times an validity flags, to test trigger coincidence size = draw (integers(min_value=1, max_value=100)) - flag = lists(integers(min_value=0, max_value=1), min_size=4, max_size=4) + flag = lists(integers(min_value=0, max_value=1 ), min_size=4, max_size=4) i = draw(lists(integers(min_value=0), min_size=size, max_size=size)) - flags = draw(lists(flag, min_size=size, max_size=size)) + flags = draw(lists(flag , min_size=size, max_size=size)) return (i, flags) -@given(evt_id=integers(1, 10), pmt_id=integers(1, 10)) -def test_get_trigger_candidates_single(double_square_wf, evt_id, pmt_id, channel_conf, trigger_conf): +@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10)) +def test_get_trigger_candidates_single(double_square_wf, evt_id, pmt_id): # Given a square pulse check that is triggering where the pulse starts and values match the expected ones. - channel_config = channel_conf.copy() - trigger_config = trigger_conf.copy() + channel_config = channel_conf() + trigger_config = trigger_conf() wf = double_square_wf[0] thr = double_square_wf[1] @@ -85,7 +88,7 @@ def test_get_trigger_candidates_single(double_square_wf, evt_id, pmt_id, channel h_true = peak.max() val_h = h_true < channel_config['amp_max'] val_p = True - + val_t = all([val_q, val_w, val_h, val_p]) triggers = get_trigger_candidates(wf, evt_id, pmt_id ,channel_config, trigger_config @@ -93,29 +96,26 @@ def test_get_trigger_candidates_single(double_square_wf, evt_id, pmt_id, channel elements = [evt_id , pmt_id ,t_true , q_true , w_true, h_true - ,val_q , val_w , val_h , val_p + ,val_q , val_w , val_h , val_p , val_t ,wf[start1-1], wf.max() ,-1 , -1 , -1] assert len(triggers) == 1 - for i, element in enumerate(elements): - assert element == triggers[0][i] + assert elements == list(triggers[0]) -@given(evt_id=integers(1, 10), pmt_id=integers(1, 10), pulse_valid=integers(0, 9)) -def test_get_trigger_candidates_double(double_square_wf, evt_id, pmt_id - ,channel_conf , trigger_conf , pulse_valid): +@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10), pulse_valid=integers(0, 9)) +def test_get_trigger_candidates_double(double_square_wf, evt_id, pmt_id, pulse_valid): # Check that the square pulses are identified as separated when pulse valid ext is lower than their separation. - channel_config = channel_conf.copy() - trigger_config = trigger_conf.copy() + channel_config = channel_conf(pulse_valid_ext=pulse_valid) + trigger_config = trigger_conf() wf = double_square_wf[0] thr = double_square_wf[1] length = double_square_wf[2] start1 = double_square_wf[3] start2 = double_square_wf[4] - channel_config = channel_conf.copy() - channel_config['pulse_valid_ext'] = pulse_valid + triggers = get_trigger_candidates(wf, evt_id, pmt_id ,channel_config, trigger_config ,np.zeros(wf.shape), wf.max()) @@ -134,23 +134,22 @@ def test_get_trigger_candidates_double(double_square_wf, evt_id, pmt_id h_true = peak.max() val_h = h_true < channel_config['amp_max'] val_p = True + val_t = all([val_q, val_w, val_h, val_p]) elements = [evt_id , pmt_id ,t_true , q_true , w_true, h_true - ,val_q , val_w , val_h , val_p + ,val_q , val_w , val_h , val_p , val_t ,wf[start1-1], wf.max() ,-1 , -1 , -1] - for i, element in enumerate(elements): - assert element == triggers[j][i] + assert elements == list(triggers[j]) -@given(evt_id=integers(1, 10), pmt_id=integers(1, 10), discard=integers(0, 30)) -def test_get_trigger_candidates_discard(double_square_wf, evt_id, pmt_id - ,channel_conf , trigger_conf , discard): +@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10), discard=integers(0, 30)) +def test_get_trigger_candidates_discard(double_square_wf, evt_id, pmt_id, discard): # Check that the small pulse is kept when discard_width is low enough - channel_config = channel_conf.copy() - trigger_config = trigger_conf.copy() + channel_config = channel_conf() + trigger_config = trigger_conf(discard_width=discard) wf = double_square_wf[0] thr = double_square_wf[1] @@ -158,7 +157,6 @@ def test_get_trigger_candidates_discard(double_square_wf, evt_id, pmt_id start1 = double_square_wf[3] start2 = double_square_wf[4] - trigger_config['discard_width'] = discard triggers = get_trigger_candidates(wf, evt_id, pmt_id ,channel_config, trigger_config ,np.zeros(wf.shape), wf.max()) @@ -181,36 +179,36 @@ def test_get_trigger_candidates_discard(double_square_wf, evt_id, pmt_id h_true = peak.max() val_h = h_true < channel_config['amp_max'] val_p = True + val_t = all([val_q, val_w, val_h, val_p]) elements = [evt_id , pmt_id ,t_true , q_true , w_true, h_true - ,val_q , val_w , val_h , val_p + ,val_q , val_w , val_h , val_p , val_t ,wf[start1-1], wf.max() ,-1 , -1 , -1] - for i, element in enumerate(elements): - assert element == triggers[j][i] + assert elements == list(triggers[j]) -@given(evt_id=integers(1, 10), pmt_id=integers(1, 10)) -def test_get_trigger_candidates_multipeak(double_square_wf, evt_id, pmt_id - ,channel_conf , trigger_conf): +@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10)) +def test_get_trigger_candidates_multipeak(double_square_wf, evt_id, pmt_id): # Check that the multipeak protection discards the trigger if conditions are met. - channel_config = channel_conf.copy() - trigger_config = trigger_conf.copy() + channel_config = channel_conf(pulse_valid_ext=0) + trigger_config = trigger_conf() wf = double_square_wf[0] thr = double_square_wf[1] length = double_square_wf[2] start1 = double_square_wf[3] start2 = double_square_wf[4] - channel_config = channel_conf.copy() - channel_config['pulse_valid_ext'] = 0 + triggers = get_trigger_candidates(wf, evt_id, pmt_id ,channel_config, trigger_config ,np.zeros(wf.shape), wf.max()) for t in triggers: - assert t[-6] == True + assert t.valid_all == all(t[6:10]) + assert t.valid_peak == True + ### Second peak characteristics peak = wf[start2:start2+length+channel_config['pulse_valid_ext']] @@ -226,6 +224,7 @@ def test_get_trigger_candidates_multipeak(double_square_wf, evt_id, pmt_id h_true = peak.max() val_h = h_true < channel_config['amp_max'] val_p = True + val_t = all([val_q, val_w, val_h, val_p]) # Multipeak condition metretrieve_trigger_information trigger_config['multipeak'] = {'q_min':q_true//2, 'time_min':w_true//2, 'time_after':length+10} @@ -233,24 +232,26 @@ def test_get_trigger_candidates_multipeak(double_square_wf, evt_id, pmt_id ,channel_config, trigger_config ,np.zeros(wf.shape), wf.max()) - assert triggers[0][-6] == False - assert triggers[1][-6] == True + assert triggers[0].valid_peak == False + assert triggers[1].valid_peak == True + for t in triggers: + assert t.valid_all == all(t[6:10]) - # Multipeak ondition no longer met + # Multipeak condition no longer met trigger_config['multipeak'] = {'q_min':q_true*2, 'time_min':w_true*2, 'time_after':length+10} triggers = get_trigger_candidates(wf, evt_id, pmt_id ,channel_config, trigger_config ,np.zeros(wf.shape), wf.max()) for t in triggers: - assert t[-6] == True + assert t.valid_all == all(t[6:10]) + assert t.valid_peak == True -@given(evt_id=integers(0, 3)) -def test_retrieve_trigger_information(double_square_wf, evt_id - ,channel_conf , trigger_conf): - channel_config = channel_conf.copy() - trigger_config = trigger_conf.copy() +@given(double_square_wf=double_square_wf(), evt_id=integers(0, 3)) +def test_retrieve_trigger_information(double_square_wf, evt_id): + channel_config = channel_conf() + trigger_config = trigger_conf() retriever = retrieve_trigger_information({0:channel_config, 1:channel_config}, trigger_config) @@ -278,35 +279,34 @@ def test_retrieve_trigger_information(double_square_wf, evt_id h_true = peak.max() val_h = h_true < channel_config['amp_max'] val_p = True + val_t = all([val_q, val_w, val_h, val_p]) elements = [evt_id , j ,t_true , q_true , w_true, h_true - ,val_q , val_w , val_h , val_p + ,val_q , val_w , val_h , val_p , val_t ,wf[start1-1], -wfs[j].max() ,-1 , -1 , -1] - for i, element in enumerate(elements): - assert element == triggers[j][i] - + assert np.allclose(elements, list(triggers[j])) @given(trigger_time_and_validity(), integers()) def test_check_trigger_coincidence(trigger_info, coincidence_window): trigger_time = trigger_info[0] validity = trigger_info[1] - trigg = np.zeros((len(validity), 15)) + trigg = np.zeros((len(validity), 16), dtype='object') # Initiate sim. trigger values trigg[:, 1 ] = list(range(len(validity))) trigg[:, 2 ] = trigger_time - trigg[:, 6:10] = 1#validity + trigg[:, 6:11] = True#validity trigg[:, -3: ] = np.full((len(validity), 3), -1) order_idx = np.lexsort((trigg[:,1], trigger_time)) old_order = np.argsort(order_idx) ordered_times = trigg[:, 2 ][order_idx] - ordered_val = trigg[:, 6:10][order_idx] + ordered_val = trigg[:, 6:11][order_idx] ordered_ids = trigg[:, 1 ][order_idx] trigger_back = trigg[order_idx].copy() @@ -325,6 +325,7 @@ def test_check_trigger_coincidence(trigger_info, coincidence_window): else: break checks_coinc = check_trigger_coincidence(coincidence_window) - new_t = checks_coinc(trigg.copy()) + new_t = checks_coinc([TriggerInfo(*t) for t in trigg]) - assert np.allclose(trigger_back[old_order], new_t) + for t_old, t_new in zip(trigger_back[old_order], new_t): + assert list(t_old) == list(t_new) From 03fdce8228f80f80988334b36e3e8029897a98fd Mon Sep 17 00:00:00 2001 From: Ander Date: Tue, 22 Jun 2021 18:02:20 +0200 Subject: [PATCH 13/14] Changes according to PR review --- invisible_cities/io/trigger_io.py | 47 ++++++++++++++----------------- 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/invisible_cities/io/trigger_io.py b/invisible_cities/io/trigger_io.py index 752eb0963..236a2910b 100644 --- a/invisible_cities/io/trigger_io.py +++ b/invisible_cities/io/trigger_io.py @@ -46,34 +46,29 @@ def _make_tables(hdf5_file, n_sensors, compression="ZLIB4"): def trigger_dst_writer(hdf5_file, **kwargs):#{ trigger_table = make_table(hdf5_file, - group = "Trigger", - name = "DST" , - fformat = table_formats.TriggerTable, - description = "Simulated trigger data", - compression = 'ZLIB4') + group = "Trigger", + name = "DST" , + fformat = table_formats.TriggerTable, + description = "Simulated trigger data", + compression = 'ZLIB4') def write_trigger(trigger_info):#{ - [event , pmt , - trigger_time , charge , width , height , - valid_q , valid_w , valid_h , valid_peak, - mean_baseline, max_height , - n_coinc , closest_ttime, closest_pmt] = trigger_info row = trigger_table.row - row["event" ] = event - row["pmt" ] = pmt - row["trigger_time" ] = trigger_time - row["q" ] = charge - row["width" ] = width - row["height" ] = height - row["valid_q" ] = valid_q - row["valid_w" ] = valid_w - row["valid_h" ] = valid_h - row["valid_peak" ] = valid_peak - row["valid_all" ] = (valid_q + valid_w + valid_h + valid_peak) == 4 - row["baseline" ] = mean_baseline - row["max_height" ] = max_height - row["n_coinc" ] = n_coinc - row["closest_ttime"] = closest_ttime - row["closest_pmt" ] = closest_pmt + row["event" ] = trigger_info.event + row["pmt" ] = trigger_info.pmt + row["trigger_time" ] = trigger_info.trigger_time + row["q" ] = trigger_info.q + row["width" ] = trigger_info.width + row["height" ] = trigger_info.height + row["valid_q" ] = trigger_info.valid_q + row["valid_w" ] = trigger_info.valid_w + row["valid_h" ] = trigger_info.valid_h + row["valid_peak" ] = trigger_info.valid_peak + row["valid_all" ] = trigger_info.valid_all + row["baseline" ] = trigger_info.baseline + row["max_height" ] = trigger_info.max_height + row["n_coinc" ] = trigger_info.n_coinc + row["closest_ttime"] = trigger_info.closest_ttime + row["closest_pmt" ] = trigger_info.closest_pmt row.append() def write_triggers(triggers): From 7f0ce4cdac8956f2229a9022959f6f50719a3d41 Mon Sep 17 00:00:00 2001 From: Ander Date: Tue, 22 Jun 2021 18:03:34 +0200 Subject: [PATCH 14/14] Modifications following PR comments --- invisible_cities/reco/trigger_functions.py | 98 ++++++++++++++-------- 1 file changed, 61 insertions(+), 37 deletions(-) diff --git a/invisible_cities/reco/trigger_functions.py b/invisible_cities/reco/trigger_functions.py index ad79c471b..2e8c68201 100644 --- a/invisible_cities/reco/trigger_functions.py +++ b/invisible_cities/reco/trigger_functions.py @@ -1,11 +1,19 @@ import numpy as np -from typing import List -from typing import Dict -from typing import Callable +from typing import List +from typing import Dict +from typing import Callable + +from collections import namedtuple from .. core.core_functions import in_range +TriggerInfo = namedtuple('TriggerInfo', + 'event pmt ' + 'trigger_time q width height ' + 'valid_q valid_w valid_h valid_peak valid_all ' + 'baseline max_height n_coinc closest_ttime closest_pmt') + def get_trigger_candidates(signal : np.ndarray ,event_id : int , pmt_id : int ,channel_conf : Dict[int, int], trigger_conf : Dict @@ -34,14 +42,14 @@ def get_trigger_candidates(signal : np.ndarray hrange = [channel_conf['baseline_dev' ], channel_conf['amp_max' ]] pulse_ext = channel_conf['pulse_valid_ext'] - fMulti = trigger_conf['multipeak' ] # Protection for later peaks + multipk = trigger_conf['multipeak' ] # Protection for later peaks discard_width = trigger_conf['discard_width'] # Don't store events below certain width # Multipeak protection parameters - if fMulti: - multi_q = fMulti['q_min' ] # Minimum charge of later peaks to discard a trigger - multi_w = int(fMulti['time_min' ]) # Minimum time width of later peaks to discard a trigger - multi_t = int(fMulti['time_after']) # Time window for a later peak to be considered + if multipk: + multi_q = multipk['q_min' ] # Minimum charge of later peaks to discard a trigger + multi_w = int(multipk['time_min' ]) # Minimum time width of later peaks to discard a trigger + multi_t = int(multipk['time_after']) # Time window for a later peak to be considered # Signal above threshold, allowing for pulse_ext counts below threshold triggers = [] @@ -54,47 +62,52 @@ def get_trigger_candidates(signal : np.ndarray peak_slice = slice(p[0]-1, p[-1] + pulse_ext + 1) # DAQ gets pulse_ext counts after peaks; also, for Q it takes a bin before TODO Fix for N100 (remove -1) peak = signal[peak_slice] if not (len(peak) > max(discard_width, 1)): continue - width = len(peak) - 1 # TODO Remove -1 in N100 - height = peak[1:].max() # TODO Full peak for N100 - peak = peak[:-2][peak[:-2]>0] # Only counts above 0 are used for Q, last bin not considered TODO ONLY -1 for N100 + width = len(peak) - 1 # TODO Remove -1 in N100 + height = peak[1: ].max() # TODO Full peak for N100 + peak = peak[ :-2][peak[:-2]>0] # Only counts above 0 are used for Q, last bin not considered TODO ONLY -1 for N100 charge = peak.sum() if len(peak)==0: continue start_time = p[0] # Trigger validity - trigger_w = in_range(width , *wrange, left_closed=True , right_closed=False) - trigger_q = in_range(charge, *qrange, left_closed=False, right_closed=False) - trigger_h = height < hrange[1] + trigger_w = in_range(width , *wrange, left_closed=True , right_closed=False) + trigger_q = in_range(charge, *qrange, left_closed=False, right_closed=False) + trigger_h = height < hrange[1] + trigger_all = all([trigger_w, trigger_q, trigger_h]) # Baseline at the start of the pulse baseline_t0 = baseline[start_time] - triggers .append([event_id , pmt_id - ,start_time , charge , width , height - ,trigger_q , trigger_w , trigger_h, True # By default multi flag is set as true - ,baseline_t0, wvf_max_height - ,-1 , -1 , -1]) # Filled later in the coincidence module + + trigger = TriggerInfo(event_id , pmt_id + ,start_time , charge , width , height + ,trigger_q , trigger_w , trigger_h, True, trigger_all # By default multiflag is set as true + ,baseline_t0, wvf_max_height + ,-1 , -1 , -1) # Filled later in the coincidence module + triggers .append(trigger) # Extra-peaks protection - if fMulti is not None: + if multipk is not None: # Evaluate if later peaks (must be over discard_width) are candidates for the multipeak protection for i, t0 in zip(range(len(triggers[:-1])), triggers[:-1]): + new_val = t0.valid_peak for t1 in triggers[i+1:]: - multi_candidate = ((t1[4] >= multi_w ) & - (t1[3] > multi_q ) & - (t1[2] <= t0[2] + t0[4] + multi_t+2)) # Start of trigger + peak width + time window + FPGA delay (2) + multi_candidate = ((t1.width >= multi_w ) & + (t1.q > multi_q ) & + (t1.trigger_time <= t0.trigger_time + t0.width + multi_t+2)) # Start of trigger + peak width + time window + FPGA delay (2) if not multi_candidate: continue # In case the later peak ends after the multipeak protection window ends - if (t1[2]+t1[4]) > (t0[2] + t0[4] + multi_t+2): # Event partially contained in the multipeak window - peak_slice = slice(t1[2]-1, t0[2]+t0[4]+multi_t+2+1) + if (t1.trigger_time + t1.width) > (t0.trigger_time + t0.width + multi_t + 2): # Event partially contained in the multipeak window + peak_slice = slice(t1.trigger_time - 1, t0.trigger_time + t0.width + multi_t + 2 + 1) peak = signal[peak_slice] width = len(peak) - 1 # TODO Remove -1 in N100 charge = peak[peak>0].sum() - triggers[i][-6] = (width >= multi_w) & (charge > multi_q) + new_val = (width >= multi_w) & (charge > multi_q) else: - triggers[i][-6] = False - if not triggers[i][-6]: break + new_val = False + if not new_val: break + triggers[i] = triggers[i]._replace(valid_peak=new_val, valid_all=all([triggers[i].valid_all, new_val])) return triggers @@ -128,7 +141,7 @@ def trigger_on_channels(rwfs : np.ndarray ,event_number, pmt_id ,pmt_conf , trigger_config ,baseline , max_height )) - return np.array(triggers) + return triggers return trigger_on_channels @@ -147,28 +160,39 @@ def check_trigger_coincidence(coinc_window : int) -> Callable: """ def pmt_coincidence(triggers : np.ndarray) -> np.ndarray: - valid_sel = np.all(triggers[:, 6:10], axis=1).astype('bool') # Use only valid triggers - times = triggers[:, 2][valid_sel] # Times of valid triggers - ids = triggers[:, 1][valid_sel] # Ids of valid triggers + valid_sel = np.array([t.valid_all for t in triggers]) # Use only valid triggers + times = np.array([t.trigger_time for t in triggers])[valid_sel] # Times of valid triggers + ids = np.array([t.pmt for t in triggers])[valid_sel] # Ids of valid triggers # Order by time, index so we can use this to order both times and ids order_idx = np.lexsort((ids, times)) order_time = times[order_idx] pmts = np.full(ids.shape, -1) + # Calculate all time differences between each trigger and all later ones - time_diffs = np.array([order_time[i+1:] - t for i, t in enumerate(order_time)]) + time_diffs = [order_time[i+1:] - t for i, t in enumerate(order_time)] + # Since pmts are ordered by time, closest pmt will be the next one in the array. pmts[:-1] = ids[order_idx][1:] + # Shortest time is the first value of time_diff for each trigger shortest = np.array([t[0] if len(t)>0 else -1 for t in time_diffs]) + # Count how many triggers are within the coincidince window of each trigger n_coinc = np.array([len(t[t