From 88c398c9fca86b2f3ea669e4f066d2965d28ef04 Mon Sep 17 00:00:00 2001 From: Luigi Pertoldi Date: Wed, 8 May 2024 16:31:09 +0200 Subject: [PATCH 1/4] [evt.modules.spms] remove outdated hacks for old data prod --- src/pygama/evt/modules/larveto.py | 3 +- src/pygama/evt/modules/spms.py | 37 +++++------------ tests/evt/test_build_evt.py | 67 ++++++++++++++++--------------- 3 files changed, 44 insertions(+), 63 deletions(-) diff --git a/src/pygama/evt/modules/larveto.py b/src/pygama/evt/modules/larveto.py index f26ae3b8e..1e6196f2d 100644 --- a/src/pygama/evt/modules/larveto.py +++ b/src/pygama/evt/modules/larveto.py @@ -45,8 +45,7 @@ def l200_combined_test_stat( amp = ak.flatten(amp, axis=-1) # subtract the HPGe t0 from the SiPM pulse t0s - # HACK: remove 16 when units will be fixed - rel_t0 = 16 * t0 - geds_t0 + rel_t0 = t0 - geds_t0 return l200_test_stat(rel_t0, amp, ts_bkg_prob, rc_density) diff --git a/src/pygama/evt/modules/spms.py b/src/pygama/evt/modules/spms.py index d29981a7e..9a43c1100 100644 --- a/src/pygama/evt/modules/spms.py +++ b/src/pygama/evt/modules/spms.py @@ -254,13 +254,6 @@ def make_pulse_data_mask( drop_empty=False, ) - # HACK: handle units - # HACK: remove me once units are fixed in the dsp tier - if "units" in pulse_t0.attrs and pulse_t0.attrs["units"] == "ns": - pulse_t0_ns = pulse_t0.view_as("ak") - else: - pulse_t0_ns = pulse_t0.view_as("ak") * 16 - pulse_amp = gather_pulse_data( datainfo, tcm, @@ -283,6 +276,7 @@ def make_pulse_data_mask( t_loc_ns = ak.fill_none(ak.nan_to_none(t_loc_ns), t_loc_default_ns) # start with all-true mask + pulse_t0_ns = pulse_t0.view_as("ak") mask = pulse_t0_ns == pulse_t0_ns # apply p.e. threshold @@ -334,7 +328,14 @@ def geds_coincidence_classifier( ) # we'll need to remove pulses below noise threshold later - is_good_pulse = gather_is_valid_hit(datainfo, tcm, table_names).view_as("ak") + is_good_pulse = gather_pulse_data( + datainfo, + tcm, + table_names, + observable="hit.is_valid_hit", + pulse_mask=None, + drop_empty=False, + ).view_as("ak") # load the data data = {} @@ -369,23 +370,3 @@ def geds_coincidence_classifier( ) return types.Array(ts_data) - - -# REMOVE ME: not needed anymore with VectorOfVectors DSP outputs -def gather_is_valid_hit(datainfo, tcm, table_names): - data = {} - for field in ("is_valid_hit", "trigger_pos"): - data[field] = gather_pulse_data( - datainfo, - tcm, - table_names, - observable=f"hit.{field}", - pulse_mask=None, - drop_empty=False, - ).view_as("ak") - - return types.VectorOfVectors( - data["is_valid_hit"][ - ak.local_index(data["is_valid_hit"]) < ak.num(data["trigger_pos"], axis=-1) - ] - ) diff --git a/tests/evt/test_build_evt.py b/tests/evt/test_build_evt.py index 99bf66d6f..add7c7fdd 100644 --- a/tests/evt/test_build_evt.py +++ b/tests/evt/test_build_evt.py @@ -123,50 +123,51 @@ def test_field_nesting(lgnd_test_data, files_config): assert sorted(evt.sub2.keys()) == ["dummy", "multiplicity"] -def test_spms_module(lgnd_test_data, files_config): - build_evt( - files_config, - config=f"{config_dir}/spms-module-config.yaml", - wo_mode="of", - ) +# FIXME: this can't be properly tested until proper testdata is available +# def test_spms_module(lgnd_test_data, files_config): +# build_evt( +# files_config, +# config=f"{config_dir}/spms-module-config.yaml", +# wo_mode="of", +# ) - outfile = files_config["evt"][0] +# outfile = files_config["evt"][0] - evt = lh5.read("/evt", outfile) +# evt = lh5.read("/evt", outfile) - t0 = ak.fill_none(ak.nan_to_none(evt.t0.view_as("ak")), 48_000) - tr_pos = evt.trigger_pos.view_as("ak") * 16 - assert ak.all(tr_pos > t0 - 30_000) - assert ak.all(tr_pos < t0 + 30_000) +# t0 = ak.fill_none(ak.nan_to_none(evt.t0.view_as("ak")), 48_000) +# tr_pos = evt.trigger_pos.view_as("ak") * 16 +# assert ak.all(tr_pos > t0 - 30_000) +# assert ak.all(tr_pos < t0 + 30_000) - mask = evt._pulse_mask - assert isinstance(mask, VectorOfVectors) - assert len(mask) == 10 - assert mask.ndim == 3 +# mask = evt._pulse_mask +# assert isinstance(mask, VectorOfVectors) +# assert len(mask) == 10 +# assert mask.ndim == 3 - full = evt.spms_amp_full.view_as("ak") - amp = evt.spms_amp.view_as("ak") - assert ak.all(amp > 0.1) +# full = evt.spms_amp_full.view_as("ak") +# amp = evt.spms_amp.view_as("ak") +# assert ak.all(amp > 0.1) - assert ak.all(full[mask.view_as("ak")] == amp) +# assert ak.all(full[mask.view_as("ak")] == amp) - wo_empty = evt.spms_amp_wo_empty.view_as("ak") - assert ak.all(wo_empty == amp[ak.count(amp, axis=-1) > 0]) +# wo_empty = evt.spms_amp_wo_empty.view_as("ak") +# assert ak.all(wo_empty == amp[ak.count(amp, axis=-1) > 0]) - rawids = evt.rawid.view_as("ak") - assert rawids.ndim == 2 - assert ak.count(rawids) == 30 +# rawids = evt.rawid.view_as("ak") +# assert rawids.ndim == 2 +# assert ak.count(rawids) == 30 - idx = evt.hit_idx.view_as("ak") - assert idx.ndim == 2 - assert ak.count(idx) == 30 +# idx = evt.hit_idx.view_as("ak") +# assert idx.ndim == 2 +# assert ak.count(idx) == 30 - rawids_wo_empty = evt.rawid_wo_empty.view_as("ak") - assert ak.count(rawids_wo_empty) == 7 +# rawids_wo_empty = evt.rawid_wo_empty.view_as("ak") +# assert ak.count(rawids_wo_empty) == 7 - vhit = evt.is_valid_hit.view_as("ak") - vhit.show() - assert ak.all(ak.num(vhit, axis=-1) == ak.num(full, axis=-1)) +# vhit = evt.is_valid_hit.view_as("ak") +# vhit.show() +# assert ak.all(ak.num(vhit, axis=-1) == ak.num(full, axis=-1)) def test_vov(lgnd_test_data, files_config): From a3a83dbbd8e63412fe8a75317aaa5b843aa43c23 Mon Sep 17 00:00:00 2001 From: Luigi Pertoldi Date: Fri, 10 May 2024 12:58:22 +0200 Subject: [PATCH 2/4] [evt.modules.spms] allow user pulse mask before computing classifier --- src/pygama/evt/modules/larveto.py | 4 +++- src/pygama/evt/modules/spms.py | 32 ++++++++++++++++++------------- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/pygama/evt/modules/larveto.py b/src/pygama/evt/modules/larveto.py index 1e6196f2d..fcbbc8a3a 100644 --- a/src/pygama/evt/modules/larveto.py +++ b/src/pygama/evt/modules/larveto.py @@ -37,7 +37,9 @@ def l200_combined_test_stat( probability for a pulse coming from some uncorrelated physics (uniform distribution). needed for the LAr scintillation time pdf. rc_density - array of densities (probabilities) of uncorrelated number of photoelectrons in a 6µs window. + density array of the random coincidence LAr energy distribution (total + energy summed over all channels, in p.e.). Derived from forced trigger + data. """ # flatten the data in the last axis (i.e. merge all channels together) # TODO: implement channel distinction diff --git a/src/pygama/evt/modules/spms.py b/src/pygama/evt/modules/spms.py index 9a43c1100..a036da417 100644 --- a/src/pygama/evt/modules/spms.py +++ b/src/pygama/evt/modules/spms.py @@ -305,6 +305,7 @@ def geds_coincidence_classifier( geds_t0_ns: types.Array, ts_bkg_prob: float, rc_density: Sequence[float] | None = None, + pulse_mask: types.VectorOfVectors | None = None, ) -> types.Array: """Calculate the HPGe / SiPMs coincidence classifier. @@ -315,9 +316,20 @@ def geds_coincidence_classifier( ---------- datainfo, tcm, table_names positional arguments automatically supplied by :func:`.build_evt`. + geds_t0_ns + t0 (ns) of the HPGe signal. + ts_bkg_prob + probability for a pulse coming from some uncorrelated physics (uniform + distribution). needed for the LAr scintillation time pdf. + rc_density + density array of the random coincidence LAr energy distribution (total + energy summed over all channels, in p.e.). Derived from forced trigger + data. + pulse_mask + optional user mask to filter pulses before classifier computation. """ # mask for windowing data around the HPGe t0 - pulse_mask = make_pulse_data_mask( + geds_t0_mask = make_pulse_data_mask( datainfo, tcm, table_names, @@ -327,16 +339,6 @@ def geds_coincidence_classifier( t_loc_default_ns=48_000, ) - # we'll need to remove pulses below noise threshold later - is_good_pulse = gather_pulse_data( - datainfo, - tcm, - table_names, - observable="hit.is_valid_hit", - pulse_mask=None, - drop_empty=False, - ).view_as("ak") - # load the data data = {} for k, obs in {"amp": "hit.energy_in_pe", "t0": "hit.trigger_pos"}.items(): @@ -349,8 +351,12 @@ def geds_coincidence_classifier( drop_empty=False, ).view_as("ak") - # remove pulses below noise threshold and outside the HPGe trigger window - data[k] = all_data[is_good_pulse & pulse_mask.view_as("ak")] + # remove pulses outside the HPGe trigger window + data[k] = all_data[geds_t0_mask.view_as("ak")] + + # remove pulses as required by user + if pulse_mask is not None: + data[k] = data[pulse_mask.view_as("ak")] # load the channel info # rawids = spms.gather_tcm_id_data( From 4355d57f8ece413a289eba01b65e9f982142a474 Mon Sep 17 00:00:00 2001 From: Luigi Pertoldi Date: Fri, 10 May 2024 15:38:15 +0200 Subject: [PATCH 3/4] [evt.modules.spms] change classifier function signature --- src/pygama/evt/modules/larveto.py | 2 +- src/pygama/evt/modules/spms.py | 48 +++++++++---------------------- 2 files changed, 14 insertions(+), 36 deletions(-) diff --git a/src/pygama/evt/modules/larveto.py b/src/pygama/evt/modules/larveto.py index fcbbc8a3a..8ec0c9625 100644 --- a/src/pygama/evt/modules/larveto.py +++ b/src/pygama/evt/modules/larveto.py @@ -6,7 +6,7 @@ import awkward as ak import numpy as np -import scipy +import scipy.stats from numpy.typing import ArrayLike diff --git a/src/pygama/evt/modules/spms.py b/src/pygama/evt/modules/spms.py index a036da417..f8f4e6e7d 100644 --- a/src/pygama/evt/modules/spms.py +++ b/src/pygama/evt/modules/spms.py @@ -302,10 +302,11 @@ def geds_coincidence_classifier( tcm: utils.TCMData, table_names: Sequence[str], *, + spms_t0: types.VectorOfVectors, + spms_amp: types.VectorOfVectors, geds_t0_ns: types.Array, ts_bkg_prob: float, rc_density: Sequence[float] | None = None, - pulse_mask: types.VectorOfVectors | None = None, ) -> types.Array: """Calculate the HPGe / SiPMs coincidence classifier. @@ -316,6 +317,10 @@ def geds_coincidence_classifier( ---------- datainfo, tcm, table_names positional arguments automatically supplied by :func:`.build_evt`. + t0 + arrival times of pulses in ns, split by channel. + amp + amplitude of pulses in p.e., split by channel. geds_t0_ns t0 (ns) of the HPGe signal. ts_bkg_prob @@ -325,8 +330,6 @@ def geds_coincidence_classifier( density array of the random coincidence LAr energy distribution (total energy summed over all channels, in p.e.). Derived from forced trigger data. - pulse_mask - optional user mask to filter pulses before classifier computation. """ # mask for windowing data around the HPGe t0 geds_t0_mask = make_pulse_data_mask( @@ -337,42 +340,17 @@ def geds_coincidence_classifier( t_loc_ns=geds_t0_ns, dt_range_ns=(-1_000, 5_000), t_loc_default_ns=48_000, - ) - - # load the data - data = {} - for k, obs in {"amp": "hit.energy_in_pe", "t0": "hit.trigger_pos"}.items(): - all_data = gather_pulse_data( - datainfo, - tcm, - table_names, - observable=obs, - pulse_mask=None, - drop_empty=False, - ).view_as("ak") - - # remove pulses outside the HPGe trigger window - data[k] = all_data[geds_t0_mask.view_as("ak")] - - # remove pulses as required by user - if pulse_mask is not None: - data[k] = data[pulse_mask.view_as("ak")] + ).view_as("ak") - # load the channel info - # rawids = spms.gather_tcm_id_data( - # datainfo, - # tcm, - # table_names, - # pulse_mask=pulse_mask, - # drop_empty=True, - # ) + # remove pulses outside the HPGe trigger window + spms_t0 = spms_t0.view_as("ak")[geds_t0_mask] + spms_amp = spms_amp.view_as("ak")[geds_t0_mask] - # (HPGe) trigger position can vary among events! - if isinstance(geds_t0_ns, types.Array): - geds_t0_ns = geds_t0_ns.view_as("ak") + # get the (HPGe) trigger position + geds_t0_ns = geds_t0_ns.view_as("ak") ts_data = larveto.l200_combined_test_stat( - data["t0"], data["amp"], geds_t0_ns, ts_bkg_prob, rc_density + spms_amp, spms_amp, geds_t0_ns, ts_bkg_prob, rc_density ) return types.Array(ts_data) From 9185cfe5eacbe49dfeecfb554318942486cdfd2e Mon Sep 17 00:00:00 2001 From: Luigi Pertoldi Date: Fri, 10 May 2024 16:00:43 +0200 Subject: [PATCH 4/4] [evt.modules.spms] also remove HPGe window masking --- src/pygama/evt/modules/spms.py | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/src/pygama/evt/modules/spms.py b/src/pygama/evt/modules/spms.py index f8f4e6e7d..5954471ec 100644 --- a/src/pygama/evt/modules/spms.py +++ b/src/pygama/evt/modules/spms.py @@ -331,26 +331,12 @@ def geds_coincidence_classifier( energy summed over all channels, in p.e.). Derived from forced trigger data. """ - # mask for windowing data around the HPGe t0 - geds_t0_mask = make_pulse_data_mask( - datainfo, - tcm, - table_names, - a_thr_pe=None, - t_loc_ns=geds_t0_ns, - dt_range_ns=(-1_000, 5_000), - t_loc_default_ns=48_000, - ).view_as("ak") - - # remove pulses outside the HPGe trigger window - spms_t0 = spms_t0.view_as("ak")[geds_t0_mask] - spms_amp = spms_amp.view_as("ak")[geds_t0_mask] - - # get the (HPGe) trigger position - geds_t0_ns = geds_t0_ns.view_as("ak") - - ts_data = larveto.l200_combined_test_stat( - spms_amp, spms_amp, geds_t0_ns, ts_bkg_prob, rc_density + return types.Array( + larveto.l200_combined_test_stat( + spms_t0.view_as("ak"), + spms_amp.view_as("ak"), + geds_t0_ns.view_as("ak"), + ts_bkg_prob, + rc_density, + ) ) - - return types.Array(ts_data)