Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add protection for empty waveforms #776

Merged
merged 7 commits into from
Mar 2, 2021
Merged

Add protection for empty waveforms #776

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2923a59
Add a test that demonstrates an unwanted behaviour
gondiaz Feb 26, 2021
4636371
Return empty list if no signal is found
gondiaz Feb 26, 2021
ade421f
Add filter for events with no signal
gondiaz Feb 26, 2021
4b5c93b
Rewrite detsim pipe with signal filter implementation
gondiaz Feb 26, 2021
1af6f61
Add empty waveforms test assertions
gondiaz Feb 26, 2021
5cfe87c
Rewrite buffy pipe with signal filter implementation
gondiaz Feb 26, 2021
7ca39e9
Cosmetics
gondiaz Feb 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
53 changes: 26 additions & 27 deletions invisible_cities/cities/buffy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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") )
Expand All @@ -101,20 +101,19 @@ 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 ))
Expand Down
38 changes: 24 additions & 14 deletions invisible_cities/cities/components.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -878,6 +878,13 @@ 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" )
Expand All @@ -891,22 +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 ,
event_times ,
calculate_buffers,
order_sensors ,
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
64 changes: 32 additions & 32 deletions invisible_cities/cities/detsim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'),
Expand Down Expand Up @@ -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"]))

write_nohits_filter = fl.sink(event_filter_writer(h5out, "active_hits") , args=("event_number", "passed"))
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,
fl.branch("event_number" ,
evtnum_collect.sink),
buffer_calculation),
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))

Expand Down
21 changes: 21 additions & 0 deletions invisible_cities/cities/detsim_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,27 @@ 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 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())
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)

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):

PATH_IN = os.path.join(ICDATADIR , "empty_mcfile.h5")
Expand Down
3 changes: 3 additions & 0 deletions invisible_cities/database/test_data/nexus_new_kr83m_fast.newformat.sim.emptywfs.h5
View file
Open in desktop
Git LFS file not shown
11 changes: 6 additions & 5 deletions invisible_cities/detsim/buffer_functions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]:
Expand All @@ -55,9 +55,10 @@ 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
## Just using this and the stand_off for now
## taking first above sum threshold.
## !! To-do: make more robust with min int? or similar
if len(indices) == 0: return []
# 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]

Expand Down