From 2923a59d87f9f2c2a01dbe6e83bcd1289c9e074f Mon Sep 17 00:00:00 2001 From: gondiaz Date: Fri, 26 Feb 2021 17:46:40 +0100 Subject: [PATCH 1/7] Add a test that demonstrates an unwanted behaviour When an empty waveform is created in detsim the signal finder does not find any pulse, breaking the city flow. --- invisible_cities/cities/detsim_test.py | 14 ++++++++++++++ .../nexus_new_kr83m_fast.newformat.sim.emptywfs.h5 | 3 +++ 2 files changed, 17 insertions(+) create mode 100644 invisible_cities/database/test_data/nexus_new_kr83m_fast.newformat.sim.emptywfs.h5 diff --git a/invisible_cities/cities/detsim_test.py b/invisible_cities/cities/detsim_test.py index fdc0bdf20..6844bbc02 100644 --- a/invisible_cities/cities/detsim_test.py +++ b/invisible_cities/cities/detsim_test.py @@ -59,6 +59,20 @@ def test_detsim_filter_active_hits(ICDATADIR, output_tmpdir): np.testing.assert_array_equal(filters["passed"], [False, True]) +def test_detsim_filter_empty_waveforms(ICDATADIR, output_tmpdir): + # the first event radious is slighty above NEW active radious of 208.0 mm + PATH_IN = os.path.join(ICDATADIR, "nexus_new_kr83m_fast.newformat.sim.emptywfs.h5") + PATH_OUT = os.path.join(output_tmpdir, "detsim_test.h5") + conf = configure('detsim $ICTDIR/invisible_cities/config/detsim.conf'.split()) + physics_params = conf["physics_params"] + physics_params["transverse_diffusion"] = 0 * units.mm / units.cm**0.5 + conf.update(dict(files_in = PATH_IN, + file_out = PATH_OUT, + run_number = 0, + physics_params = physics_params)) + result = detsim(**conf) + + def test_detsim_empty_input_file(ICDATADIR, output_tmpdir): PATH_IN = os.path.join(ICDATADIR , "empty_mcfile.h5") diff --git a/invisible_cities/database/test_data/nexus_new_kr83m_fast.newformat.sim.emptywfs.h5 b/invisible_cities/database/test_data/nexus_new_kr83m_fast.newformat.sim.emptywfs.h5 new file mode 100644 index 000000000..275a35a5a --- /dev/null +++ b/invisible_cities/database/test_data/nexus_new_kr83m_fast.newformat.sim.emptywfs.h5 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:857a8528be210a21b38e21da4047fd092424fbf97baa129060de7b9a7f423bb5 +size 63458 From 463637137694d8e0181428d0d44dbcd1fece610a Mon Sep 17 00:00:00 2001 From: gondiaz Date: Fri, 26 Feb 2021 19:24:48 +0100 Subject: [PATCH 2/7] Return empty list if no signal is found --- invisible_cities/detsim/buffer_functions.py | 1 + 1 file changed, 1 insertion(+) diff --git a/invisible_cities/detsim/buffer_functions.py b/invisible_cities/detsim/buffer_functions.py index 0653856a4..636835149 100644 --- a/invisible_cities/detsim/buffer_functions.py +++ b/invisible_cities/detsim/buffer_functions.py @@ -55,6 +55,7 @@ def find_signal_start(wfs : Union[pd.Series, np.ndarray], eng_sum = wfs.sum() indices = indices_and_wf_above_threshold(eng_sum, bin_threshold).indices + if len(indices) == 0: return [] ## Just using this and the stand_off for now ## taking first above sum threshold. ## !! To-do: make more robust with min int? or similar From ade421f7256fe8106bb4348493f412a726c8a90e Mon Sep 17 00:00:00 2001 From: gondiaz Date: Fri, 26 Feb 2021 19:25:54 +0100 Subject: [PATCH 3/7] Add filter for events with no signal --- invisible_cities/cities/components.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py index df3063cac..06f502c85 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -369,7 +369,7 @@ def mcsensors_from_file(paths : List[str], """ if rate == 0: - warnings.warn("Zero rate is unphysical, using set period of 1 ms instead", + warnings.warn("Zero rate is unphysical, using set period of 1 ms instead", category=UserWarning) pmt_ids = load_db.DataPMT(db_file, run_number).SensorID @@ -878,6 +878,12 @@ def calculate_and_save_buffers(buffer_length : float , args = "pmt_bin_wfs" , out = "pulses" ) + filter_events_signal = fl.map(lambda x: len(x) > 0, + args= 'pulses', + out = 'passed_signal') + events_passed_signal = fl.count_filter(bool, args='passed_signal') + write_signal_filter = fl.sink(event_filter_writer(h5out, "signal"), args=('event_number', 'passed_signal')) + event_times = fl.map(trigger_times , args = ("pulses", "timestamp", "pmt_bins"), out = "evt_times" ) @@ -902,10 +908,13 @@ def calculate_and_save_buffers(buffer_length : float , saved_buffers )) fn_list = (find_signal , + filter_events_signal, + fl.branch(write_signal_filter), + events_passed_signal.filter, event_times , calculate_buffers, order_sensors , - buffer_writer_ ) + fl.branch(buffer_writer_)) # Filter out order_sensors if it is not set buffer_definition = pipe(*filter(None, fn_list)) From 4b5c93bb336d1231314e1e7310819fa60a765fa2 Mon Sep 17 00:00:00 2001 From: gondiaz Date: Fri, 26 Feb 2021 19:26:47 +0100 Subject: [PATCH 4/7] Rewrite detsim pipe with signal filter implementation Avoid unnecessary branching in detsim --- invisible_cities/cities/detsim.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/invisible_cities/cities/detsim.py b/invisible_cities/cities/detsim.py index 6257af491..306b4fe2c 100644 --- a/invisible_cities/cities/detsim.py +++ b/invisible_cities/cities/detsim.py @@ -144,8 +144,8 @@ def detsim(*, files_in, file_out, event_range, print_mod, compression, filter_events_no_active_hits = fl.map(lambda x:np.any(x), args= 'energy_a', - out = 'passed') - events_passed_active_hits = fl.count_filter(bool, args='passed') + out = 'passed_active') + events_passed_active_hits = fl.count_filter(bool, args='passed_active') simulate_electrons = fl.map(ielectron_simulator(**physics_params_), args = ('x_a', 'y_a', 'z_a', 'time_a', 'energy_a'), @@ -198,7 +198,7 @@ def detsim(*, files_in, file_out, event_range, print_mod, compression, int(buffer_params["length"] / buffer_params["pmt_width"]), int(buffer_params["length"] / buffer_params["sipm_width"])) - write_nohits_filter = fl.sink(event_filter_writer(h5out, "active_hits") , args=("event_number", "passed")) + write_nohits_filter = fl.sink(event_filter_writer(h5out, "active_hits"), args=("event_number", "passed_active")) result = fl.push(source= MC_hits_from_files(files_in), pipe = fl.pipe(fl.slice(*event_range, close_all=True), event_count_in.spy , @@ -213,9 +213,9 @@ def detsim(*, files_in, file_out, event_range, print_mod, compression, create_pmt_waveforms, create_sipm_waveforms, get_bin_edges, - fl.branch("event_number" , - evtnum_collect.sink), - buffer_calculation), + buffer_calculation, + "event_number" , + evtnum_collect.sink), result = dict(events_in = event_count_in.future, evtnum_list = evtnum_collect.future)) From 1af6f614a136da5a9a3fde3461d100b7e5dda4af Mon Sep 17 00:00:00 2001 From: gondiaz Date: Fri, 26 Feb 2021 19:27:17 +0100 Subject: [PATCH 5/7] Add empty waveforms test assertions --- invisible_cities/cities/detsim_test.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/invisible_cities/cities/detsim_test.py b/invisible_cities/cities/detsim_test.py index 6844bbc02..3126bd18d 100644 --- a/invisible_cities/cities/detsim_test.py +++ b/invisible_cities/cities/detsim_test.py @@ -72,6 +72,13 @@ def test_detsim_filter_empty_waveforms(ICDATADIR, output_tmpdir): physics_params = physics_params)) result = detsim(**conf) + assert result.events_in == 2 + assert result.evtnum_list == [1] + + with tb.open_file(PATH_OUT, mode="r") as h5out: + filters = h5out.root.Filters.signal.read() + np.testing.assert_array_equal(filters["passed"], [False, True]) + def test_detsim_empty_input_file(ICDATADIR, output_tmpdir): From 5cfe87c5bf76dd6975498669d94ddbc1406b384d Mon Sep 17 00:00:00 2001 From: gondiaz Date: Fri, 26 Feb 2021 19:28:29 +0100 Subject: [PATCH 6/7] Rewrite buffy pipe with signal filter implementation --- invisible_cities/cities/buffy.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/invisible_cities/cities/buffy.py b/invisible_cities/cities/buffy.py index 7c524e0cf..0ae714c82 100644 --- a/invisible_cities/cities/buffy.py +++ b/invisible_cities/cities/buffy.py @@ -101,20 +101,20 @@ def buffy(files_in , file_out , compression , event_range, result = fl.push(source = mcsensors_from_file(files_in , detector_db, run_number , - rate ) , - pipe = fl.pipe(fl.slice(*event_range , - close_all=True) , - event_count_in.spy , - print_every(print_mod) , - filter_events , - fl.branch(write_filter) , - events_with_resp.filter , - extract_tminmax , - bin_pmt_wf , - bin_sipm_wf , - fl.branch("event_number" , - evtnum_collect.sink), - buffer_calculation ) , + rate ), + pipe = fl.pipe(fl.slice(*event_range , + close_all=True), + event_count_in.spy , + print_every(print_mod) , + filter_events , + fl.branch(write_filter) , + events_with_resp.filter , + extract_tminmax , + bin_pmt_wf , + bin_sipm_wf , + buffer_calculation , + "event_number" , + evtnum_collect.sink ), result = dict(events_in = event_count_in.future , events_resp = events_with_resp.future, evtnum_list = evtnum_collect.future )) From 7ca39e9aa490eed3cf2350f7167be0066898f2cb Mon Sep 17 00:00:00 2001 From: gondiaz Date: Sat, 27 Feb 2021 13:18:11 +0100 Subject: [PATCH 7/7] Cosmetics Fix typo in test comment Rename subpipe Re-align commas in calculate_and_save_buffers Re-align commas in detsim Re-align commas in buffy Rename TODOs in file --- invisible_cities/cities/buffy.py | 51 +++++++++---------- invisible_cities/cities/components.py | 39 +++++++------- invisible_cities/cities/detsim.py | 56 ++++++++++----------- invisible_cities/cities/detsim_test.py | 2 +- invisible_cities/detsim/buffer_functions.py | 10 ++-- 5 files changed, 79 insertions(+), 79 deletions(-) diff --git a/invisible_cities/cities/buffy.py b/invisible_cities/cities/buffy.py index 0ae714c82..52736b201 100644 --- a/invisible_cities/cities/buffy.py +++ b/invisible_cities/cities/buffy.py @@ -77,19 +77,19 @@ def buffy(files_in , file_out , compression , event_range, events_with_resp = fl.count_filter(bool, args="event_passed") with tb.open_file(file_out, "w", filters=tbl.filters(compression)) as h5out: - buffer_calculation = calculate_and_save_buffers(buffer_length , - max_time , - pre_trigger , - pmt_wid , - sipm_wid , - trigger_threshold, - h5out , - run_number , - npmt , - nsipm , - nsamp_pmt , - nsamp_sipm , - order_sensors ) + buffer_calculation = calculate_and_save_buffers( buffer_length + , max_time + , pre_trigger + , pmt_wid + , sipm_wid + , trigger_threshold + , h5out + , run_number + , npmt + , nsipm + , nsamp_pmt + , nsamp_sipm + , order_sensors) write_filter = fl.sink(event_filter_writer(h5out, "detected_events"), args=("event_number", "event_passed") ) @@ -102,19 +102,18 @@ def buffy(files_in , file_out , compression , event_range, detector_db, run_number , rate ), - pipe = fl.pipe(fl.slice(*event_range , - close_all=True), - event_count_in.spy , - print_every(print_mod) , - filter_events , - fl.branch(write_filter) , - events_with_resp.filter , - extract_tminmax , - bin_pmt_wf , - bin_sipm_wf , - buffer_calculation , - "event_number" , - evtnum_collect.sink ), + pipe = fl.pipe( fl.slice(*event_range, close_all=True) + , event_count_in.spy + , print_every(print_mod) + , filter_events + , fl.branch(write_filter) + , events_with_resp.filter + , extract_tminmax + , bin_pmt_wf + , bin_sipm_wf + , buffer_calculation + , "event_number" + , evtnum_collect.sink), result = dict(events_in = event_count_in.future , events_resp = events_with_resp.future, evtnum_list = evtnum_collect.future )) diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py index 06f502c85..16dcef017 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -879,10 +879,11 @@ def calculate_and_save_buffers(buffer_length : float , out = "pulses" ) filter_events_signal = fl.map(lambda x: len(x) > 0, - args= 'pulses', - out = 'passed_signal') + args= 'pulses', + out = 'passed_signal') events_passed_signal = fl.count_filter(bool, args='passed_signal') - write_signal_filter = fl.sink(event_filter_writer(h5out, "signal"), args=('event_number', 'passed_signal')) + write_signal_filter = fl.sink(event_filter_writer(h5out, "signal"), + args=('event_number', 'passed_signal')) event_times = fl.map(trigger_times , args = ("pulses", "timestamp", "pmt_bins"), @@ -897,25 +898,25 @@ def calculate_and_save_buffers(buffer_length : float , saved_buffers = "buffers" if order_sensors is None else "ordered_buffers" max_subevt = max_time // buffer_length + 1 - buffer_writer_ = sink(buffer_writer(h5out , - run_number = run_number, - n_sens_eng = npmt , - n_sens_trk = nsipm , - length_eng = nsamp_pmt , - length_trk = nsamp_sipm, - max_subevt = max_subevt), + buffer_writer_ = sink(buffer_writer( h5out + , run_number = run_number + , n_sens_eng = npmt + , n_sens_trk = nsipm + , length_eng = nsamp_pmt + , length_trk = nsamp_sipm + , max_subevt = max_subevt), args = ("event_number", "evt_times" , saved_buffers )) - fn_list = (find_signal , - filter_events_signal, - fl.branch(write_signal_filter), - events_passed_signal.filter, - event_times , - calculate_buffers, - order_sensors , - fl.branch(buffer_writer_)) + find_signal_and_write_buffers = ( find_signal + , filter_events_signal + , fl.branch(write_signal_filter) + , events_passed_signal.filter + , event_times + , calculate_buffers + , order_sensors + , fl.branch(buffer_writer_)) # Filter out order_sensors if it is not set - buffer_definition = pipe(*filter(None, fn_list)) + buffer_definition = pipe(*filter(None, find_signal_and_write_buffers)) return buffer_definition diff --git a/invisible_cities/cities/detsim.py b/invisible_cities/cities/detsim.py index 306b4fe2c..0bf39dfd0 100644 --- a/invisible_cities/cities/detsim.py +++ b/invisible_cities/cities/detsim.py @@ -185,37 +185,37 @@ def detsim(*, files_in, file_out, event_range, print_mod, compression, evtnum_collect = collect() with tb.open_file(file_out, "w", filters = tbl.filters(compression)) as h5out: - buffer_calculation = calculate_and_save_buffers(buffer_params["length"] , - buffer_params["max_time"] , - buffer_params["pre_trigger"], - buffer_params["pmt_width"] , - buffer_params["sipm_width"] , - buffer_params["trigger_thr"], - h5out , - run_number , - len(datapmt) , - len(datasipm), - int(buffer_params["length"] / buffer_params["pmt_width"]), - int(buffer_params["length"] / buffer_params["sipm_width"])) + buffer_calculation = calculate_and_save_buffers( buffer_params["length"] + , buffer_params["max_time"] + , buffer_params["pre_trigger"] + , buffer_params["pmt_width"] + , buffer_params["sipm_width"] + , buffer_params["trigger_thr"] + , h5out + , run_number + , len(datapmt) + , len(datasipm) + , int(buffer_params["length"] / buffer_params["pmt_width"]) + , int(buffer_params["length"] / buffer_params["sipm_width"])) write_nohits_filter = fl.sink(event_filter_writer(h5out, "active_hits"), args=("event_number", "passed_active")) result = fl.push(source= MC_hits_from_files(files_in), - pipe = fl.pipe(fl.slice(*event_range, close_all=True), - event_count_in.spy , - print_every(print_mod), - select_s1_candidate_hits, - select_active_hits, - filter_events_no_active_hits, - fl.branch(write_nohits_filter), - events_passed_active_hits.filter, - simulate_electrons, - get_buffer_times_and_length, - create_pmt_waveforms, - create_sipm_waveforms, - get_bin_edges, - buffer_calculation, - "event_number" , - evtnum_collect.sink), + pipe = fl.pipe( fl.slice(*event_range, close_all=True) + , event_count_in.spy + , print_every(print_mod) + , select_s1_candidate_hits + , select_active_hits + , filter_events_no_active_hits + , fl.branch(write_nohits_filter) + , events_passed_active_hits.filter + , simulate_electrons + , get_buffer_times_and_length + , create_pmt_waveforms + , create_sipm_waveforms + , get_bin_edges + , buffer_calculation + , "event_number" + , evtnum_collect.sink), result = dict(events_in = event_count_in.future, evtnum_list = evtnum_collect.future)) diff --git a/invisible_cities/cities/detsim_test.py b/invisible_cities/cities/detsim_test.py index 3126bd18d..1153a9110 100644 --- a/invisible_cities/cities/detsim_test.py +++ b/invisible_cities/cities/detsim_test.py @@ -60,7 +60,7 @@ def test_detsim_filter_active_hits(ICDATADIR, output_tmpdir): def test_detsim_filter_empty_waveforms(ICDATADIR, output_tmpdir): - # the first event radious is slighty above NEW active radious of 208.0 mm + # the first event radius is slighty above NEW active radius of 208.0 mm PATH_IN = os.path.join(ICDATADIR, "nexus_new_kr83m_fast.newformat.sim.emptywfs.h5") PATH_OUT = os.path.join(output_tmpdir, "detsim_test.h5") conf = configure('detsim $ICTDIR/invisible_cities/config/detsim.conf'.split()) diff --git a/invisible_cities/detsim/buffer_functions.py b/invisible_cities/detsim/buffer_functions.py index 636835149..6c5443f87 100644 --- a/invisible_cities/detsim/buffer_functions.py +++ b/invisible_cities/detsim/buffer_functions.py @@ -40,8 +40,8 @@ def bin_sensors(sensors : pd.DataFrame, return bins[:-1], bin_sensors -## !! to-do: clarify for non-pmt versions of next -## !! to-do: Check on integral instead of only threshold? +# TODO: clarify for non-pmt versions of next +# TODO: Check on integral instead of only threshold? def find_signal_start(wfs : Union[pd.Series, np.ndarray], bin_threshold: float , stand_off : int ) -> List[int]: @@ -56,9 +56,9 @@ def find_signal_start(wfs : Union[pd.Series, np.ndarray], indices = indices_and_wf_above_threshold(eng_sum, bin_threshold).indices if len(indices) == 0: return [] - ## Just using this and the stand_off for now - ## taking first above sum threshold. - ## !! To-do: make more robust with min int? or similar + # Just using this and the stand_off for now + # taking first above sum threshold. + # TODO: make more robust with min int? or similar all_indx = split_in_peaks(indices, stand_off) return [pulse[0] for pulse in all_indx]