diff --git a/invisible_cities/calib/calib_sensors_functions.py b/invisible_cities/calib/calib_sensors_functions.py index f8db736c3..38d522f20 100644 --- a/invisible_cities/calib/calib_sensors_functions.py +++ b/invisible_cities/calib/calib_sensors_functions.py @@ -136,15 +136,13 @@ def pmt_subtract_maw(cwfs, n_maw=100): return cwfs - maw -def calibrate_sipms(sipm_wfs, adc_to_pes, thr, *, bls_mode=BlsMode.mode): +def calibrate_sipms(sipm_wfs, adc_to_pes, *, bls_mode=BlsMode.mode): """ - Subtracts the baseline, calibrates waveforms to pes - and suppresses values below `thr` (in pes). + Subtract baseline, and calibrates waveforms to pes. """ - thr = to_col_vector(np.full(sipm_wfs.shape[0], thr)) bls = subtract_baseline(sipm_wfs, bls_mode=bls_mode) cwfs = calibrate_wfs(bls, adc_to_pes) - return np.where(cwfs > thr, cwfs, 0) + return cwfs def subtract_mean (wfs): return subtract_baseline(wfs, bls_mode=BlsMode.mean ) diff --git a/invisible_cities/calib/calib_sensors_functions_test.py b/invisible_cities/calib/calib_sensors_functions_test.py index 452d22d90..1b3948561 100644 --- a/invisible_cities/calib/calib_sensors_functions_test.py +++ b/invisible_cities/calib/calib_sensors_functions_test.py @@ -13,6 +13,7 @@ from . import calib_sensors_functions as csf from .. core import core_functions as cf +from .. reco import wfm_functions as wfm from .. sierpe import fee as FE from .. sierpe import waveform_generator as wfg @@ -207,10 +208,11 @@ def test_calibrate_sipms_stat(oscillating_waveform_with_baseline, baseline ) = oscillating_waveform_with_baseline #n_maw = n_samples // 500 - ccwfs = csf.calibrate_sipms(wfs, adc_to_pes, nsigma * noise_sigma, bls_mode=BlsMode.mode) + ccwfs = csf.calibrate_sipms(wfs, adc_to_pes, bls_mode=BlsMode.mode) + zeroed_ccwfs = wfm.zero_wfs_below_threshold(ccwfs, zeroing_thr=nsigma * noise_sigma) - number_of_zeros = np.count_nonzero(ccwfs == 0) - assert number_of_zeros > fraction * ccwfs.size + number_of_zeros = np.count_nonzero(zeroed_ccwfs == 0) + assert number_of_zeros > fraction * zeroed_ccwfs.size def test_calibrate_sipms_common_threshold(toy_sipm_signal): @@ -218,8 +220,8 @@ def test_calibrate_sipms_common_threshold(toy_sipm_signal): signal_zs_common_threshold, _, common_threshold, _) = toy_sipm_signal - zs_wf = csf.calibrate_sipms(signal_adc, adc_to_pes, - common_threshold, bls_mode=BlsMode.mode) + ccwf = csf.calibrate_sipms(signal_adc, adc_to_pes, bls_mode=BlsMode.mode) + zs_wf = wfm.zero_wfs_below_threshold(ccwf, zeroing_thr=common_threshold) for actual, expected in zip(zs_wf, signal_zs_common_threshold): assert actual == approx(expected) @@ -230,10 +232,9 @@ def test_calibrate_sipms_individual_thresholds(toy_sipm_signal): _, signal_zs_individual_thresholds, _, individual_thresholds) = toy_sipm_signal + ccwf = csf.calibrate_sipms(signal_adc, adc_to_pes, bls_mode=BlsMode.mode) + zs_wf = wfm.zero_wfs_below_threshold(ccwf, zeroing_thr=individual_thresholds) - zs_wf = csf.calibrate_sipms(signal_adc, adc_to_pes, - individual_thresholds, - bls_mode=BlsMode.mode) for actual, expected in zip(zs_wf, signal_zs_individual_thresholds): assert actual == approx(expected) @@ -309,11 +310,12 @@ def test_area_of_sum_equals_sum_of_areas_pmts(square_pmt_and_sipm_waveforms): def test_area_of_sum_equals_sum_of_areas_sipms(square_pmt_and_sipm_waveforms): _, nsensors, _, _, _, sipms_wfm, _ = square_pmt_and_sipm_waveforms - adc_to_pes = np.full(nsensors, 100, dtype=float) - cwfs = csf.calibrate_sipms(sipms_wfm, adc_to_pes, thr=10, bls_mode=BlsMode.mode) - stot = np.sum(cwfs[0]) * nsensors - sums = np.sum(cwfs, axis=1) - stot2 = reduce(add, sums) + adc_to_pes = np.full(nsensors, 100, dtype=float) + cwfs = csf.calibrate_sipms(sipms_wfm, adc_to_pes, bls_mode=BlsMode.mode) + zeroed_cwfs = wfm.zero_wfs_below_threshold(cwfs, zeroing_thr=10) + stot = np.sum(zeroed_cwfs[0]) * nsensors + sums = np.sum(zeroed_cwfs, axis=1) + stot2 = reduce(add, sums) assert stot == approx(stot2, rel=1e-3) diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py index e5945bd7d..6a69741e1 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -83,6 +83,8 @@ from .. types .ic_types import types_dict_summary from .. types .ic_types import types_dict_tracks from .. types .symbols import WfType +from .. types .symbols import CutAlgo +from .. types .symbols import SiPMSelectionMethod from .. types .symbols import RebinMethod from .. types .symbols import SiPMCharge from .. types .symbols import BlsMode @@ -746,7 +748,7 @@ def sensor_data(path, wf_type): def build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, - s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2): + s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, apply_cut): s1_params = dict(time = minmax(min = s1_tmin, max = s1_tmax), length = minmax(min = s1_lmin, @@ -766,8 +768,8 @@ def build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, def build_pmap(ccwf, s1_indx, s2_indx, sipmzs): # -> PMap return pkf.get_pmap(ccwf, s1_indx, s2_indx, sipmzs, - s1_params, s2_params, thr_sipm_s2, pmt_ids, - pmt_samp_wid, sipm_samp_wid) + s1_params, s2_params, pmt_ids, + pmt_samp_wid, sipm_samp_wid, apply_cut) return build_pmap @@ -785,19 +787,77 @@ def calibrate_pmts(cwf):# -> CCwfs: return calibrate_pmts -def calibrate_sipms(dbfile, run_number, thr_sipm): +def calibrate_sipms(dbfile, run_number): DataSiPM = load_db.DataSiPM(dbfile, run_number) adc_to_pes = np.abs(DataSiPM.adc_to_pes.values) def calibrate_sipms(rwf): return csf.calibrate_sipms(rwf, adc_to_pes = adc_to_pes, - thr = thr_sipm, bls_mode = BlsMode.mode) return calibrate_sipms +def apply_cutting_function(algo, **cutting_params): + + if algo is CutAlgo.threshold: + func = threshold_sipm_selection(**cutting_params) + elif algo is CutAlgo.pyrrha: + func = pyrrha_sipm_selection(**cutting_params) + else: + raise ValueError(f"Unsupported cutting algorithm: {algo!r}. Expected one of {list(CutAlgo)}") + + def apply_cutting_function(wfs): + return func(wfs) + + return apply_cutting_function + + +def threshold_sipm_selection(thr_sipm_type + , thr_sipm + , thr_sipm_s2 + , run_number + , detector_db = None): + ''' + Function that applies thresholding to the sipms in standard irene manner, + by zeroing all waveform values below a threshold. + ''' + # assume that if the detector_db is None, you return sipm threshold as the number provided + if detector_db is None: + sipm_thr = thr_sipm + else: + # extract sipm threshold + sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number) + + def threshold_sipm_selection(wfs): + return wfm.charge_threshold_method(wfs, zeroing_thr = sipm_thr, integration_thr=thr_sipm_s2) + + return threshold_sipm_selection + + +def pyrrha_sipm_selection(selection_method : SiPMSelectionMethod + , selection_kwargs : dict + , proximity_threshold : float + , padding_radius : float + , run_number : int + , detector_db : str): + ''' + Function that applies a generic selection function to the sipms, which can be used to + implement a spatial SiPM selection method (called Pyrrha). + ''' + def pyrrha_sipm_selection(wfs): + return wfm.spatial_selection_method(wfs, + selection_method, + selection_kwargs, + proximity_threshold, + padding_radius, + run_number, + detector_db) + + return pyrrha_sipm_selection + + def calibrate_with_mean(dbfile, run_number): DataSiPM = load_db.DataSiPM(dbfile, run_number) adc_to_pes = np.abs(DataSiPM.adc_to_pes.values) @@ -1214,7 +1274,7 @@ def integrate_wfs(wfs): def compute_and_write_pmaps(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2, - h5out, sipm_rwf_to_cal=None): + h5out, apply_cut, sipm_rwf_to_cal=None): # Filter events without signal over threshold indices_pass = fl.map(check_nonempty_indices, @@ -1225,7 +1285,7 @@ def compute_and_write_pmaps(detector_db, run_number, pmt_samp_wid, sipm_samp_wid # Build the PMap compute_pmap = fl.map(build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, - s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2), + s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, apply_cut), args = ("ccwfs", "s1_indices", "s2_indices", "sipm"), out = "pmap") diff --git a/invisible_cities/cities/hypathia.py b/invisible_cities/cities/hypathia.py index 6a125b542..68f038c6c 100644 --- a/invisible_cities/cities/hypathia.py +++ b/invisible_cities/cities/hypathia.py @@ -30,6 +30,7 @@ from .. io .run_and_event_io import run_and_event_writer from .. io . trigger_io import trigger_writer from .. types.symbols import WfType +from .. types.symbols import CutAlgo from .. types.symbols import SiPMThreshold from .. dataflow import dataflow as fl @@ -49,36 +50,40 @@ from . components import simulate_sipm_response from . components import calibrate_sipms from . components import get_actual_sipm_thr +from . components import apply_cutting_function + +from typing import Dict +from typing import Any @city -def hypathia( files_in : OneOrManyFiles - , file_out : str - , compression : str - , event_range : EventRangeType - , print_mod : int - , detector_db : str - , run_number : int - , sipm_noise_cut : float - , filter_padding : int - , thr_sipm : float - , thr_sipm_type : SiPMThreshold - , pmt_wfs_rebin : int - , pmt_pe_rms : float - , s1_lmin : int , s1_lmax : int - , s1_tmin : float, s1_tmax : float - , s1_rebin_stride : int , s1_stride : int - , thr_csum_s1 : float - , s2_lmin : int , s2_lmax : int - , s2_tmin : float, s2_tmax : float - , s2_rebin_stride : int , s2_stride : int - , thr_csum_s2 : float, thr_sipm_s2 : float - , pmt_samp_wid : float - , sipm_samp_wid : float +def hypathia( files_in : OneOrManyFiles + , file_out : str + , compression : str + , event_range : EventRangeType + , print_mod : int + , detector_db : str + , run_number : int + , sipm_noise_cut : float + , filter_padding : int + , thr_sipm : float + , thr_sipm_type : SiPMThreshold + , pmt_wfs_rebin : int + , pmt_pe_rms : float + , s1_lmin : int , s1_lmax : int + , s1_tmin : float, s1_tmax : float + , s1_rebin_stride : int , s1_stride : int + , thr_csum_s1 : float + , s2_lmin : int , s2_lmax : int + , s2_tmin : float, s2_tmax : float + , s2_rebin_stride : int , s2_stride : int + , thr_csum_s2 : float, thr_sipm_s2 : float + , pmt_samp_wid : float + , sipm_samp_wid : float + , cutting_function : CutAlgo + , cutting_params : Dict[str, Any] ): - sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number) - #### Define data transformations sd = sensor_data(files_in[0], WfType.mcrd) @@ -112,9 +117,13 @@ def hypathia( files_in : OneOrManyFiles item="sipm") # SiPMs calibration - sipm_rwf_to_cal = fl.map(calibrate_sipms(detector_db, run_number, sipm_thr), + sipm_rwf_to_cal = fl.map(calibrate_sipms(detector_db, run_number), item = "sipm") + # apply function depending on user input, from provided list of functions + apply_cut = apply_cutting_function(cutting_function, **cutting_params) + + event_count_in = fl.spy_count() event_count_out = fl.spy_count() @@ -134,7 +143,7 @@ def hypathia( files_in : OneOrManyFiles detector_db, run_number, pmt_samp_wid, sipm_samp_wid, s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2, - h5out, sipm_rwf_to_cal) + h5out, apply_cut, sipm_rwf_to_cal) result = push(source = wf_from_files(files_in, WfType.mcrd), pipe = pipe(fl.slice(*event_range, close_all=True), diff --git a/invisible_cities/cities/irene.py b/invisible_cities/cities/irene.py index 1c3e91ec8..31a6b8041 100644 --- a/invisible_cities/cities/irene.py +++ b/invisible_cities/cities/irene.py @@ -25,6 +25,7 @@ from .. io .run_and_event_io import run_and_event_writer from .. io .trigger_io import trigger_writer from .. types.symbols import WfType +from .. types.symbols import CutAlgo from .. types.symbols import SiPMThreshold from .. dataflow import dataflow as fl @@ -43,34 +44,49 @@ from . components import wf_from_files from . components import get_number_of_active_pmts from . components import compute_and_write_pmaps -from . components import get_actual_sipm_thr +from . components import apply_cutting_function +from typing import Dict +from typing import Any @city -def irene( files_in : OneOrManyFiles - , file_out : str - , compression : str - , event_range : EventRangeType - , print_mod : int - , detector_db : str - , run_number : int - , n_baseline : int - , n_maw : int - , thr_maw : float - , thr_sipm : float - , thr_sipm_type : SiPMThreshold - , s1_lmin : int , s1_lmax : int - , s1_tmin : float, s1_tmax : float - , s1_rebin_stride : int , s1_stride : int - , thr_csum_s1 : float - , s2_lmin : int , s2_lmax : int - , s2_tmin : float, s2_tmax : float - , s2_rebin_stride : int , s2_stride : int - , thr_csum_s2 : float, thr_sipm_s2 : float - , pmt_samp_wid : float, sipm_samp_wid: float +def irene( files_in : OneOrManyFiles + , file_out : str + , compression : str + , event_range : EventRangeType + , print_mod : int + , detector_db : str + , run_number : int + , n_baseline : int + , n_maw : int + , thr_maw : float + , thr_sipm : float + , thr_sipm_type : SiPMThreshold + , s1_lmin : int , s1_lmax : int + , s1_tmin : float, s1_tmax : float + , s1_rebin_stride : int , s1_stride : int + , thr_csum_s1 : float + , s2_lmin : int , s2_lmax : int + , s2_tmin : float, s2_tmax : float + , s2_rebin_stride : int , s2_stride : int + , thr_csum_s2 : float, thr_sipm_s2 : float + , pmt_samp_wid : float, sipm_samp_wid: float + , cutting_function : CutAlgo + , cutting_params : Dict[str, Any] ): + ''' + `cutting_function` is defined within components.py, and can vary, resulting + in the need for `cutting_params`, which are defined as such to allow for + the prior function to run. Currently implemented are `threshold_sipm_selection`, + with differing selections methods added soon. - sipm_thr = get_actual_sipm_thr(thr_sipm_type, thr_sipm, detector_db, run_number) + + params for `threshold_sipm_selection`: + thr_sipm_type : SiPMThreshold + thr_sipm : float + detector_db : str + run_number : int + ''' #### Define data transformations @@ -89,9 +105,13 @@ def irene( files_in : OneOrManyFiles args = ("cwf_sum", "cwf_sum_maw"), out = ("s1_indices", "s2_indices", "s2_energies")) + # Remove baseline and calibrate SiPMs - sipm_rwf_to_cal = fl.map(calibrate_sipms(detector_db, run_number, sipm_thr), + sipm_rwf_to_cal = fl.map(calibrate_sipms(detector_db, run_number), item = "sipm") + # apply function depending on user input, from provided list of functions + apply_cut = apply_cutting_function(cutting_function, **cutting_params) + event_count_in = fl.spy_count() event_count_out = fl.spy_count() @@ -115,7 +135,7 @@ def irene( files_in : OneOrManyFiles s1_lmax, s1_lmin, s1_rebin_stride, s1_stride, s1_tmax, s1_tmin, s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, thr_sipm_s2, - h5out, sipm_rwf_to_cal) + h5out, apply_cut, sipm_rwf_to_cal) result = push(source = wf_from_files(files_in, WfType.rwf), pipe = pipe(fl.slice(*event_range, close_all=True), diff --git a/invisible_cities/config/hypathia.conf b/invisible_cities/config/hypathia.conf index 3aafcd8e4..d3029f25f 100644 --- a/invisible_cities/config/hypathia.conf +++ b/invisible_cities/config/hypathia.conf @@ -52,3 +52,10 @@ thr_sipm_s2 = 5 * pes # Threshold for the full sipm waveform pmt_samp_wid = 25 * ns sipm_samp_wid = 1 * mus + +cutting_function = threshold +cutting_params = dict( thr_sipm_type = thr_sipm_type + , thr_sipm = thr_sipm + , thr_sipm_s2 = thr_sipm_s2 + , detector_db = detector_db + , run_number = run_number) diff --git a/invisible_cities/config/irene.conf b/invisible_cities/config/irene.conf index 9b536c7c1..1f2312a78 100644 --- a/invisible_cities/config/irene.conf +++ b/invisible_cities/config/irene.conf @@ -52,3 +52,10 @@ thr_sipm_s2 = 10 * pes # Threshold for the full sipm waveform pmt_samp_wid = 25 * ns sipm_samp_wid = 1 * mus + +cutting_function = threshold +cutting_params = dict( thr_sipm_type = thr_sipm_type + , thr_sipm = thr_sipm + , thr_sipm_s2 = thr_sipm_s2 + , detector_db = detector_db + , run_number = run_number) diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py index 686b16c87..c440cda96 100644 --- a/invisible_cities/reco/peak_functions.py +++ b/invisible_cities/reco/peak_functions.py @@ -28,6 +28,7 @@ def indices_and_wf_above_threshold(wf, thr): def select_wfs_above_time_integrated_thr(wfs, thr): selected_ids = np.where(np.sum(wfs, axis=1) >= thr)[0] selected_wfs = wfs[selected_ids] + return selected_ids, selected_wfs @@ -76,13 +77,18 @@ def build_pmt_responses(indices, times, widths, ccwf, def build_sipm_responses(indices, times, widths, - sipm_wfs, rebin_stride, thr_sipm_s2): - _, _, sipm_wfs_ = pick_slice_and_rebin(indices , times, widths, + sipm_wfs, rebin_stride, apply_cut): + _, _, sipm_wfs = pick_slice_and_rebin(indices , times, widths, sipm_wfs, rebin_stride, pad_zeros = False) - (sipm_ids, - sipm_wfs) = select_wfs_above_time_integrated_thr(sipm_wfs_, - thr_sipm_s2) + if apply_cut is not None: + # apply cut before slicing and rebinning + (sipm_ids, + sipm_wfs) = apply_cut(sipm_wfs) + else: + # give all sipm ids as index if no cut is applied + sipm_ids = np.arange(sipm_wfs.shape[0]) + return SiPMResponses(sipm_ids, sipm_wfs) @@ -93,7 +99,7 @@ def build_peak(indices, times, pmt_samp_wid = 25 * units.ns, sipm_samp_wid = 1 * units.mus, sipm_wfs = None, - thr_sipm_s2 = 0): + apply_cut = None): sipm_pmt_bin_ratio = int(sipm_samp_wid/pmt_samp_wid) (pk_times , pk_widths, @@ -107,7 +113,7 @@ def build_peak(indices, times, widths * sipm_pmt_bin_ratio, sipm_wfs, rebin_stride // sipm_pmt_bin_ratio, - thr_sipm_s2) + apply_cut) else: sipm_r = SiPMResponses.build_empty_instance() @@ -120,7 +126,8 @@ def find_peaks(ccwfs, index, Pk, pmt_ids, pmt_samp_wid = 25*units.ns, sipm_samp_wid = 1*units.mus, - sipm_wfs=None, thr_sipm_s2=0): + sipm_wfs=None, apply_cut = None): + ccwfs = np.array(ccwfs, ndmin=2) peaks = [] @@ -136,20 +143,20 @@ def find_peaks(ccwfs, index, rebin_stride, with_sipms, Pk, pmt_samp_wid, sipm_samp_wid, - sipm_wfs, thr_sipm_s2) + sipm_wfs, apply_cut) peaks.append(pk) return peaks def get_pmap(ccwf, s1_indx, s2_indx, sipm_zs_wf, - s1_params, s2_params, thr_sipm_s2, pmt_ids, - pmt_samp_wid, sipm_samp_wid): + s1_params, s2_params, pmt_ids, + pmt_samp_wid, sipm_samp_wid, apply_cut = None): return PMap(find_peaks(ccwf, s1_indx, Pk=S1, pmt_ids=pmt_ids, pmt_samp_wid=pmt_samp_wid, **s1_params), find_peaks(ccwf, s2_indx, Pk=S2, pmt_ids=pmt_ids, sipm_wfs = sipm_zs_wf, - thr_sipm_s2 = thr_sipm_s2, + apply_cut = apply_cut, pmt_samp_wid = pmt_samp_wid, sipm_samp_wid = sipm_samp_wid, **s2_params)) diff --git a/invisible_cities/reco/peak_functions_test.py b/invisible_cities/reco/peak_functions_test.py index 16efa9bcf..bdde0a2b3 100644 --- a/invisible_cities/reco/peak_functions_test.py +++ b/invisible_cities/reco/peak_functions_test.py @@ -13,6 +13,7 @@ from hypothesis.strategies import integers from hypothesis.extra.numpy import arrays +from ..cities.components import apply_cutting_function from ..core.testing_utils import exactly from ..core.testing_utils import previous_float from ..core.testing_utils import assert_Peak_equality @@ -26,6 +27,8 @@ from ..evm .pmaps import PMap from ..io .pmaps_io import load_pmaps from ..types.ic_types import minmax +from ..types.symbols import SiPMThreshold +from ..types.symbols import CutAlgo from . import peak_functions as pf @@ -350,10 +353,18 @@ def test_build_sipm_responses(wf_with_indices): wfs_slice = wfs[:, indices] peak_integrals = wfs_slice.sum(axis=1) below_thr_index = np.argmin (peak_integrals) - # next_float doesn't work here thr = peak_integrals[below_thr_index] * 1.000001 + + # next_float doesn't work here + cut_params = dict(detector_db = 'None', + thr_sipm_s2 = thr, + thr_sipm = 0, + thr_sipm_type = SiPMThreshold.common, + run_number = 0) + apply_cut = apply_cutting_function(CutAlgo.threshold, **cut_params) + sipm_r = pf.build_sipm_responses(indices, times, widths, - wfs, 1, thr) + wfs, 1, apply_cut) expected_ids = np.delete( ids, below_thr_index) expected_wfs = np.delete(wfs_slice, below_thr_index, axis=0) @@ -396,7 +407,7 @@ def test_build_peak_development(pmt_and_sipm_wfs_with_indices, with_sipms = with_sipms, Pk = Pk, sipm_wfs = sipm_wfs, - thr_sipm_s2 = -1) + ) assert_Peak_equality(peak, expected_peak) @@ -473,8 +484,7 @@ def test_find_peaks_s2_style(pmt_and_sipm_wfs_with_indices): time_range, length_range, stride, rebin_stride, S2, pmt_ids, - sipm_wfs = sipm_wfs, - thr_sipm_s2 = -1) + sipm_wfs = sipm_wfs) (rebinned_times, rebinned_widths, @@ -502,7 +512,6 @@ def test_get_pmap(s1_and_s2_with_indices): sipm_samp_wid = 1 * units.mus pmap = pf.get_pmap(pmt_wfs, s1_indx, s2_indx, sipm_wfs, s1_params, s2_params, - thr_sipm_s2 = -1, pmt_ids = pmt_ids, pmt_samp_wid = pmt_samp_wid , sipm_samp_wid = sipm_samp_wid) diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py index 4ecced9e2..0f2feba26 100644 --- a/invisible_cities/reco/wfm_functions.py +++ b/invisible_cities/reco/wfm_functions.py @@ -3,10 +3,17 @@ authors: J.J. Gomez-Cadenas, G. Martinez """ import numpy as np +from typing import Optional +from typing import Callable +from typing import Tuple -from .. core.core_functions import define_window +from .. core.core_functions import define_window, to_col_vector from .. calib import calib_sensors_functions as csf from .. sierpe import blr +from .. database import load_db +from .. reco.peak_functions import select_wfs_above_time_integrated_thr + +from .. types .symbols import SiPMSelectionMethod def to_adc(wfs, adc_to_pes): """ @@ -129,3 +136,211 @@ def compare_cwf_blr(cwf, pmtblr, event_list, window_size=500): DIFF.append(diff) return np.array(DIFF) + + +def zero_wfs_below_threshold(wfs : np.ndarray, + zeroing_thr : Optional[float] = 2.) -> np.ndarray: + """ + Zeroes the entries of the input waveforms that are below a given threshold. + + Parameters + ---------- + wfs : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM. + zeroing_thr : Charge threshold for zero suppression in PE, default 2. + + Returns + ------- + 2D array of shape (n_sipms, n_time_bins) containing the input waveforms with entries below threshold set to zero. + """ + thr = to_col_vector(np.full(wfs.shape[0], zeroing_thr)) + return np.where(wfs > thr, wfs, 0) + + +def median_std_method(wfs : np.ndarray, + nsigma : Optional[float] = 3.) -> np.ndarray: + """ + Computes the median and standard deviation of the time summed SiPM + waveforms and selects the SiPMs that are nsigma over the median. + + Parameters + ---------- + wfs : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM. + nsigma : Number of standard deviations above the median, default 3. + + Returns + ------- + Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected. + """ + charges = np.sum(wfs, axis=1) + threshold = np.median(charges) + nsigma * np.std(charges) + return charges >= threshold + + +def charge_threshold_method(wfs : np.ndarray, + zeroing_thr : Optional[float] = 2., + integration_thr : Optional[float] = 5.) -> Tuple[np.ndarray, np.ndarray]: + """ + Selects the SiPMs whose time summed waveforms within the s2 windows are above two thresholds: + - initial zero suprresion threshold (setting values in each waveform below a value to 0) + - threshold over each selected integrated slice of the waveforms + + Parameters + ---------- + wfs : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM. + zeroing_thr : Charge threshold for zero suppression in PE, default 2. + integration_thr : Charge threshold for total SiPM waveform in PE, default 5. + + Returns + ------- + Tuple of np arrays including all passing sipm ids and the corresponding waveforms + """ + # zero entries below threshold + zwfs = zero_wfs_below_threshold(wfs, zeroing_thr) + + # returns selected ids and waveforms above integral + return select_wfs_above_time_integrated_thr(zwfs, integration_thr) + + +def top_n_method(wfs : np.ndarray, + n : Optional[int] = 10) -> np.ndarray: + """ + Selects the SiPMs with the top n highest time summed waveforms. + + Parameters + ---------- + wfs : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM. + n : Number of most energeticSiPMs to select, default 10. + + Returns + ------- + Boolean numpy array of shape (n_sipms,) where True indicates that the SiPM is selected. + """ + charges = np.sum(wfs, axis=1) + idx = np.argsort(charges)[-n:] + + selected_ids = np.zeros_like(charges, dtype=bool) + selected_ids[idx] = True + return selected_ids + + +def kill_isolated_sipms(selected_ids : np.ndarray, + sipm_x : np.ndarray, + sipm_y : np.ndarray, + proximity_threshold : float) -> np.ndarray: + """ + For the SiPMs that have passed the previous selection, scans through the SiPMs to check if they + have neighbouring SiPMs - i.e., within the proximity_threshold - that have also passed the selection. + If no neighbours are found, the SiPMs are classed as isolated, and are removed. + + Parameters + ---------- + selected_ids : Boolean array of shape (n_sipms,) indicating which SiPMs passed the previous selection. + sipm_x : 1D array of shape (n_sipms,) containing the x positions of the SiPMs. + sipm_y : 1D array of shape (n_sipms,) containing the y positions of the SiPMs. + proximity_threshold : Distance threshold in mm used to identify isolated SiPMs. + + Returns + ------- + selected_ids_no_isolated : Boolean array of shape (n_sipms,) where True indicates that the SiPM is selected. + """ + selected_ids_no_isolated = selected_ids.copy() + + for i in np.where(selected_ids)[0]: + x, y = sipm_x[i], sipm_y[i] + + distances = np.sqrt((sipm_x - x)**2 + (sipm_y - y)**2) + + n_neighbors = np.sum((distances < proximity_threshold) & selected_ids) + + if n_neighbors <= 1: + selected_ids_no_isolated[i] = False + + return selected_ids_no_isolated + + +def apply_circular_padding(selected_ids_no_isolated : np.ndarray, + sipm_x : np.ndarray, + sipm_y : np.ndarray, + padding_radius : float) -> np.ndarray: + """ + For the SiPMs that pass the previous selection, creates circular padding of radius padding_radius, + selecting all SiPMs within that radius. Stores the union of all selected SiPMs. + + Parameters + ---------- + selected_ids_no_isolated : Boolean array of shape (n_sipms,) indicating which SiPMs passed the previous selection. + sipm_x : 1D array of shape (n_sipms,) containing the x positions of the SiPMs. + sipm_y : 1D array of shape (n_sipms,) containing the y positions of the SiPMs. + padding_radius : Distance threshold in mm used to create circular padding around selected SiPMs. + + Returns + ------- + sipm_ids_with_signal : Boolean array of shape (n_sipms,) where True indicates that the SiPM is selected. + """ + sipm_ids_with_signal = np.zeros_like(selected_ids_no_isolated, dtype=bool) + + for i in np.where(selected_ids_no_isolated)[0]: + x, y = sipm_x[i], sipm_y[i] + distances = np.sqrt((sipm_x - x)**2 + (sipm_y - y)**2) + sipm_ids_with_signal |= distances <= padding_radius + + return sipm_ids_with_signal + + +def spatial_selection_method(wfs : np.ndarray, + selection_method : SiPMSelectionMethod, + selection_kwargs : dict, + proximity_threshold : float, + padding_radius : float, + run_number : int, + detector_db : str) -> np.ndarray: + """ + SiPM selection function, applies SiPM cuts based on user input. + A first selection of SiPMs is made, isolated SiPMs are removed + and padding is added around the SiPMs that are left. + + Parameters + ---------- + wfs : 2D array of shape (n_sipms, n_time_bins) containing the waveforms of each SiPM. + selection_method : Method used to select SiPMs. + selection_kwargs : Dictionary of arguments passed to the selection function. + proximity_threshold : Threshold used to identify isolated SiPMs. + padding_radius : Radial padding added to each SiPM that passes the selections. + run_number : Run number used to load the detector database. + detector_db : Database used to load the detector geometry. + + Returns + ------- + selected_ids : Array of shape (n_sipms,) containing the indices of the selected SiPMs. + selected_wfs : 2D array of shape (n_selected_sipms, n_time_bins) with the waveforms of the selected SiPMs. + """ + detector_info = load_db.DataSiPM(detector_db, run_number) + sipm_x = np.array(detector_info.X) + sipm_y = np.array(detector_info.Y) + + if selection_method is SiPMSelectionMethod.median_std_method: + starting_ids_ = median_std_method(wfs, **selection_kwargs) + elif selection_method is SiPMSelectionMethod.top_n_method: + starting_ids_ = top_n_method(wfs, **selection_kwargs) + else: + raise ValueError(f"Selection method {selection_method} not recognized.") + + selected_ids_no_isolated_ = kill_isolated_sipms( + starting_ids_, + sipm_x, + sipm_y, + proximity_threshold + ) + + sipm_ids_with_signal_ = apply_circular_padding( + selected_ids_no_isolated_, + sipm_x, + sipm_y, + padding_radius + ) + + selected_ids = np.where(sipm_ids_with_signal_)[0] + selected_wfs = wfs[selected_ids] + + return selected_ids, selected_wfs + diff --git a/invisible_cities/reco/wfm_functions_test.py b/invisible_cities/reco/wfm_functions_test.py index 35ab09a88..4ef8209b3 100644 --- a/invisible_cities/reco/wfm_functions_test.py +++ b/invisible_cities/reco/wfm_functions_test.py @@ -5,6 +5,7 @@ import tables as tb from pytest import mark +from pytest import fixture from .. database import load_db @@ -49,3 +50,172 @@ def test_compare_cwf_blr(dbnew, ICDATADIR): event_list=range(NEVT), window_size=300) assert max(diff) < 0.15 + +@fixture +def sipm_wfs_for_sipm_selection_testing(): + """ + 2D array of SiPM waveforms. All but two SiPMs have 100 PEs integrated charge. + Two outliers: + - SiPM 2: 150 PEs, more than 3 sigma above the median + - SiPM 7: 120 PEs, between 1-2 sigma above the median + """ + n_sipms = 10 + n_time_bins = 100 + + wfs = np.ones((n_sipms, n_time_bins), dtype=np.float32) + + wfs[2, :] = 1.5 + wfs[7, :] = 1.2 + # median ~ 100, std ~ 15.5 + return wfs, [2, 7], [2], [0, 1, 3, 4, 5, 6, 8, 9] + + +def test_zero_wfs_below_threshold(sipm_wfs_for_sipm_selection_testing): + """ + Test function zero_wfs_below_threshold(). The test asserts that the function correctly + sets to zero the entries of the waveforms that are below a specified threshold. + """ + wfs, passing_wfs_11_ids, _, zeroed_wfs_11_ids = sipm_wfs_for_sipm_selection_testing + + zeroed_wfs_11 = wfm.zero_wfs_below_threshold(wfs, zeroing_thr=1.1) + + assert np.all(zeroed_wfs_11[zeroed_wfs_11_ids] == 0) + assert np.all(zeroed_wfs_11[passing_wfs_11_ids] == wfs[passing_wfs_11_ids]) + + +def test_median_std_method(sipm_wfs_for_sipm_selection_testing): + """ + Test function median_std_method(). The test asserts that the function correctly + identifies the outliers based on different standard deviation thresholds. + """ + wfs, expected_outliers_1sigma, expected_outliers_3sigma, _ = sipm_wfs_for_sipm_selection_testing + + passing_sipms_1sigma = np.where(wfm.median_std_method(wfs, nsigma=1))[0].tolist() + passing_sipms_3sigma = np.where(wfm.median_std_method(wfs, nsigma=3))[0].tolist() + + assert passing_sipms_1sigma == expected_outliers_1sigma + assert passing_sipms_3sigma == expected_outliers_3sigma + + +def test_threshold_method(sipm_wfs_for_sipm_selection_testing): + """ + Test function threshold_method(). The test asserts that the function correctly + kills the SiPMs below a certain charge threshold. + """ + wfs, expected_outliers_110pes, expected_outliers_130pes, _ = sipm_wfs_for_sipm_selection_testing + + passing_sipms_110pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=110) + passing_sipms_130pes, _ = wfm.charge_threshold_method(wfs, zeroing_thr=0, integration_thr=130) + + assert passing_sipms_110pes.tolist() == expected_outliers_110pes + assert passing_sipms_130pes.tolist() == expected_outliers_130pes + + +def test_top_n_method(sipm_wfs_for_sipm_selection_testing): + """ + Test function top_n_method(). The test asserts that the function selects the correct + number of SiPMs andcorrectly identifies the top N SiPMs based on their integrated charge. + """ + wfs, expected_outliers_top2, expected_outliers_top1, _ = sipm_wfs_for_sipm_selection_testing + + passing_sipms_top2 = np.where(wfm.top_n_method(wfs, n=2))[0].tolist() + passing_sipms_top1 = np.where(wfm.top_n_method(wfs, n=1))[0].tolist() + + assert len(passing_sipms_top2) == 2 + assert len(passing_sipms_top1) == 1 + + assert passing_sipms_top2 == expected_outliers_top2 + assert passing_sipms_top1 == expected_outliers_top1 + + +@fixture +def sipm_grid_for_isolation_and_padding_testing(): + """ + 5x5 grid of with 10mm spacing containing 4 SiPMs: + - SiPM 7 (x=20, y=10) has nearest neighbours 12 and 13 + - SiPM 12 (x=20, y=20) has nearest neighbours 7 and 13 + - SiPM 13 (x=30, y=20) has nearestneighbours 12 and 7 + - SiPM 24 (x=40, y=40) has no nearest neighbours + With proximity_threshold=15mm SiPMs 7, 12, 13 survive and SiPM 24 is killed. + With proximity_threshold=5mm, all SiPMs are killed. + + With a padding of 12mm around the cluster of SiPMs 7, 12, 13 you also + include SiPMs 2, 6, 8, 11, 14, 17, 18. + """ + spacing = 10 + xs = np.arange(5) * spacing # [0, 10, 20, 30, 40] + ys = np.arange(5) * spacing + + grid_x, grid_y = np.meshgrid(xs, ys) + sipm_x = grid_x.flatten().astype(np.float32) # shape (25,) + sipm_y = grid_y.flatten().astype(np.float32) + + # visual representation of the grid with SiPMs plotted as their IDs: + # X X X X 24 + # X X X X X + # X X 12 13 X + # X X 7 X X + # X X X X X + selected_ids = np.zeros(25, dtype=bool) + cluster_ids = [7, 12, 13] + isolated_id = [24] + for idx in cluster_ids + isolated_id: + selected_ids[idx] = True + + # diagonal SiPM distance ~ 14mm + # only SiPMs with a nearest neighbour pass assuming proximity_threshold=15mm + # SiPM 24 should be killed + expected_survivors_15mm = np.zeros(25, dtype=bool) + for idx in cluster_ids: + expected_survivors_15mm[idx] = True + + # given a SiPM distance of 10mm, there should be no survivors with proximity_threshold=5mm + expected_survivors_5mm = np.zeros(25, dtype=bool) + + # adding a padding of 12mm around the cluster of SiPMs 7, 12, 13 would include + # SiPMs 2, 6, 8, 11, 14, 17, 18 + # X X X X X + # X X 17 18 X + # X 11 12 13 14 + # X 6 7 8 X + # X X 2 X X + padded_cluster_ids = cluster_ids + [2, 6, 8, 11, 14, 17, 18] + expected_survivors_padding12 = np.zeros(25, dtype=bool) + for idx in padded_cluster_ids: + expected_survivors_padding12[idx] = True + + # adding a padding of 0mm around the cluster of SiPMs 7, 12, 13 would include only the cluster itself + expected_survivors_padding0 = np.zeros(25, dtype=bool) + for idx in cluster_ids: + expected_survivors_padding0[idx] = True + + return (sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm, + expected_survivors_padding12, expected_survivors_padding0) + + +def test_kill_isolated_sipms(sipm_grid_for_isolation_and_padding_testing): + """" + Test function kill_isolated_sipms(). The test asserts that the function correctly + identifies and removes isolated SiPMs based on their proximity to other SiPMs. + """ + sipm_x, sipm_y, selected_ids, expected_survivors_15mm, expected_survivors_5mm, _, _ = sipm_grid_for_isolation_and_padding_testing + + surviving_sipms_15mm = wfm.kill_isolated_sipms(selected_ids, sipm_x, sipm_y, proximity_threshold=15.0) + surviving_sipms_5mm = wfm.kill_isolated_sipms(selected_ids, sipm_x, sipm_y, proximity_threshold=5.0) + + assert surviving_sipms_15mm.tolist() == expected_survivors_15mm.tolist() + assert surviving_sipms_5mm.tolist() == expected_survivors_5mm.tolist() + + +def test_apply_circular_padding(sipm_grid_for_isolation_and_padding_testing): + """ + Test function apply_circular_padding(). The test asserts that the function correctly applies + a circular padding around selected SiPMs to include neighboring SiPMs within the specified radius. + """ + sipm_x, sipm_y, _, selected_ids, _, expected_survivors_padding12, expected_survivors_padding0 = sipm_grid_for_isolation_and_padding_testing + + padded_sipms_12mm = wfm.apply_circular_padding(selected_ids, sipm_x, sipm_y, padding_radius=12.0) + padded_sipms_0mm = wfm.apply_circular_padding(selected_ids, sipm_x, sipm_y, padding_radius=0.0) + + assert padded_sipms_12mm.tolist() == expected_survivors_padding12.tolist() + assert padded_sipms_0mm.tolist() == expected_survivors_padding0.tolist() \ No newline at end of file diff --git a/invisible_cities/types/symbols.py b/invisible_cities/types/symbols.py index e5916f326..43a650493 100644 --- a/invisible_cities/types/symbols.py +++ b/invisible_cities/types/symbols.py @@ -130,6 +130,12 @@ class XYReco(AutoNameEnumBase): barycenter = auto() corona = auto() +class CutAlgo(AutoNameEnumBase): + threshold = auto() + pyrrha = auto() + +class SiPMSelectionMethod(AutoNameEnumBase): + median_std_method = auto() class WfType(AutoNameEnumBase): rwf = auto()