diff --git a/invisible_cities/calib/calib_sensors_functions.py b/invisible_cities/calib/calib_sensors_functions.py index f8db736c3b..e58d2b7f23 100644 --- a/invisible_cities/calib/calib_sensors_functions.py +++ b/invisible_cities/calib/calib_sensors_functions.py @@ -59,8 +59,9 @@ def modes (wfs): return to_col_vector(mode (wfs, axis=1)) def subtract_baseline(wfs, *, bls_mode=BlsMode.mean): """ - Subtract the baseline to all waveforms in the input - with a specific algorithm. + Subtract the baseline from all waveforms in the input + with a specific algorithm, leaving individual samples + that are already 0 untouched. Parameters ---------- @@ -75,16 +76,19 @@ def subtract_baseline(wfs, *, bls_mode=BlsMode.mean): Returns ------- bls: np.ndarray with shape (n, m) - Baseline-subtracted waveforms. + Baseline-subtracted waveforms, with originally-zero + samples left as 0. """ - if bls_mode is BlsMode.mean : return wfs - means (wfs) - elif bls_mode is BlsMode.median : return wfs - medians (wfs) - elif bls_mode is BlsMode.mode : return wfs - modes (wfs) - elif bls_mode is BlsMode.scipymode: return wfs - scipy_mode(wfs, axis=1) + if bls_mode is BlsMode.mean : baseline = means (wfs) + elif bls_mode is BlsMode.median : baseline = medians (wfs) + elif bls_mode is BlsMode.mode : baseline = modes (wfs) + elif bls_mode is BlsMode.scipymode: baseline = scipy_mode(wfs, axis=1) else: raise TypeError(f"Unrecognized baseline subtraction option: {bls_mode}") + return np.where(wfs == 0, 0, wfs - baseline) # subtracts the baseline to non-zero samples only + def calibrate_wfs(wfs, adc_to_pes): """ @@ -136,15 +140,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 ) @@ -159,10 +161,10 @@ def sipm_subtract_median_and_calibrate(sipm_wfs, adc_to_pes): return calibrate_w # Dict of functions for SiPM processing sipm_processing = { - SiPMCalibMode.subtract_mode : subtract_mode ,# For gain extraction - SiPMCalibMode.subtract_median : subtract_median ,# For gain extraction - SiPMCalibMode.subtract_mode_calibrate : sipm_subtract_mode_and_calibrate ,# For PDF calculation - SiPMCalibMode.subtract_mean_calibrate : sipm_subtract_mean_and_calibrate ,# For PDF calculation - SiPMCalibMode.subtract_median_calibrate: sipm_subtract_median_and_calibrate,# For PDF calculation - SiPMCalibMode.subtract_mode_zs : calibrate_sipms # For data processing + SiPMCalibMode.subtract_mode : subtract_mode ,# For gain extraction + SiPMCalibMode.subtract_median : subtract_median ,# For gain extraction + SiPMCalibMode.subtract_mode_calibrate : sipm_subtract_mode_and_calibrate ,# For PDF calculation + SiPMCalibMode.subtract_mean_calibrate : sipm_subtract_mean_and_calibrate ,# For PDF calculation + SiPMCalibMode.subtract_median_calibrate : sipm_subtract_median_and_calibrate,# For PDF calculation + SiPMCalibMode.subtract_baseline_calibrate: calibrate_sipms # For data processing } diff --git a/invisible_cities/calib/calib_sensors_functions_test.py b/invisible_cities/calib/calib_sensors_functions_test.py index 452d22d90a..1b39485619 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 7445c505ae..751b43b8a7 100644 --- a/invisible_cities/cities/components.py +++ b/invisible_cities/cities/components.py @@ -84,6 +84,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 @@ -762,7 +764,8 @@ 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, + sipm_selection_algo): s1_params = dict(time = minmax(min = s1_tmin, max = s1_tmax), length = minmax(min = s1_lmin, @@ -783,8 +786,8 @@ def build_pmap(detector_db, run_number, pmt_samp_wid, sipm_samp_wid, sipm_ids = np.argwhere(datasipm.Active.values==1).flatten() 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, sipm_ids, - pmt_samp_wid, sipm_samp_wid) + s1_params, s2_params, pmt_ids, sipm_ids, + pmt_samp_wid, sipm_samp_wid, sipm_selection_algo) return build_pmap @@ -802,7 +805,7 @@ 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) DataSiPM = DataSiPM.loc[lambda df: df.Active==1] adc_to_pes = np.abs(DataSiPM.adc_to_pes.values) @@ -810,12 +813,77 @@ def calibrate_sipms(dbfile, run_number, thr_sipm): 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 select_cutting_algorithm(algo, **cutting_params): + if algo is CutAlgo.threshold: + return threshold_sipm_selection(**cutting_params) + elif algo is CutAlgo.pyrrha: + return pyrrha_sipm_selection(**cutting_params) + elif algo is CutAlgo.no_cut: + return no_cut_sipm_selection() + else: + raise ValueError(f"Unsupported cutting algorithm: {algo!r}. Expected one of {list(CutAlgo)}") + + +def threshold_sipm_selection( thr_sipm_type + , thr_sipm + , thr_sipm_s2 + , run_number + , detector_db = None): + ''' + Applies energy thresholds to SiPM S2 waveforms: + - thr_sipm: applied per time bin. Waveform samples below this + threshold are set to zero. + - thr_sipm_s2: applied to the integrated waveform charge. If the total + charge is below this threshold, it is set to zero. + ''' + if detector_db is None: + sipm_thr = thr_sipm + else: + 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): + ''' + 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 no_cut_sipm_selection(): + """" + Function that applies no cuts to the SiPM waveforms. + """ + def no_cut_sipm_selection(wfs): + sipm_ids = np.arange(wfs.shape[0]) + return sipm_ids, wfs + return no_cut_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) @@ -1231,8 +1299,8 @@ def integrate_wfs(wfs): # Compound components 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): + s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, + h5out, sipm_selection_algo, sipm_rwf_to_cal=None): # Filter events without signal over threshold indices_pass = fl.map(check_nonempty_indices, @@ -1243,7 +1311,8 @@ 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, + sipm_selection_algo), args = ("ccwfs", "s1_indices", "s2_indices", "sipm"), out = "pmap") diff --git a/invisible_cities/cities/hypathia.py b/invisible_cities/cities/hypathia.py index 91f2abaeb7..d99458039e 100644 --- a/invisible_cities/cities/hypathia.py +++ b/invisible_cities/cities/hypathia.py @@ -28,6 +28,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 @@ -48,36 +49,38 @@ from . components import calibrate_sipms from . components import get_actual_sipm_thr from . components import sensor_masker +from . components import select_cutting_algorithm + +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 + , 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 + , 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) @@ -115,9 +118,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 + sipm_selection_algo = select_cutting_algorithm(cutting_function, **cutting_params) + + event_count_in = fl.spy_count() event_count_out = fl.spy_count() @@ -136,8 +143,8 @@ def hypathia( files_in : OneOrManyFiles compute_pmaps, empty_indices, empty_pmaps = 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) + s2_lmax, s2_lmin, s2_rebin_stride, s2_stride, s2_tmax, s2_tmin, + h5out, sipm_selection_algo, 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 472c0af1e4..e7776e741f 100644 --- a/invisible_cities/cities/irene.py +++ b/invisible_cities/cities/irene.py @@ -25,6 +25,7 @@ from .. io .trigger_io import trigger_writer from .. io .dst_io import df_writer from .. types.symbols import WfType +from .. types.symbols import CutAlgo from .. types.symbols import SiPMThreshold from .. database.load_db import DataPMT @@ -46,36 +47,50 @@ from . components import wf_from_files from . components import get_number_of_pmts from . components import compute_and_write_pmaps -from . components import get_actual_sipm_thr from . components import sensor_masker +from . components import select_cutting_algorithm +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 - , store_db : bool = True +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 + , 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 + , pmt_samp_wid : float, sipm_samp_wid: float + , cutting_function : CutAlgo + , cutting_params : Dict[str, Any] + , store_db : bool = True ): + ''' + `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 @@ -99,9 +114,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 + sipm_selection_algo = select_cutting_algorithm(cutting_function, **cutting_params) + event_count_in = fl.spy_count() event_count_out = fl.spy_count() @@ -124,8 +143,7 @@ def irene( 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, sipm_selection_algo, 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/cities/irene_test.py b/invisible_cities/cities/irene_test.py index d0fe3f4270..a3f7a41513 100644 --- a/invisible_cities/cities/irene_test.py +++ b/invisible_cities/cities/irene_test.py @@ -82,13 +82,15 @@ def test_irene_electrons_40keV(config_tmpdir, ICDATADIR, s12params, nrequired = 2 conf = configure('dummy invisible_cities/config/irene.conf'.split()) - conf.update(dict(detector_db = DetDB.new, - run_number = 0, + cutting_params = conf['cutting_params'].copy() + cutting_params.update(thr_sipm_type = thr_sipm_type, + thr_sipm = thr_sipm_value) + conf.update(dict(detector_db = DetDB.new, + run_number = 0, files_in = PATH_IN, file_out = PATH_OUT, event_range = (0, nrequired), - thr_sipm_type = thr_sipm_type, - thr_sipm = thr_sipm_value, + cutting_params = cutting_params, **unpack_s12params(s12params))) cnt = irene(**conf) @@ -414,20 +416,23 @@ def test_irene_sequential_times(config_tmpdir, ICDATADIR): PATH_OUT = os.path.join(config_tmpdir, 'test_pmaps.h5') conf = configure('dummy invisible_cities/config/irene.conf'.split()) - conf.update(dict(files_in = PATH_IN , - file_out = PATH_OUT , - run_number = 6351 , - n_baseline = 48000 , - thr_sipm = 1 * units.pes, - s1_tmin = 0 * units.mus, - s1_tmax = 640 * units.mus, - s1_lmin = 5 , - s1_lmax = 30 , - s2_tmin = 645 * units.mus, - s2_tmax = 1300 * units.mus, - s2_lmin = 80 , - s2_lmax = 200000 , - thr_sipm_s2 = 5 * units.pes)) + cutting_params = conf['cutting_params'].copy() + cutting_params.update(run_number = 6351, + thr_sipm_s2 = 5 * units.pes, + thr_sipm = 1 * units.pes) + conf.update(dict(files_in = PATH_IN , + file_out = PATH_OUT , + run_number = 6351 , + n_baseline = 48000 , + s1_tmin = 0 * units.mus, + s1_tmax = 640 * units.mus, + s1_lmin = 5 , + s1_lmax = 30 , + s2_tmin = 645 * units.mus, + s2_tmax = 1300 * units.mus, + s2_lmin = 80 , + s2_lmax = 200000 , + cutting_params = cutting_params)) irene(**conf) diff --git a/invisible_cities/config/hypathia.conf b/invisible_cities/config/hypathia.conf index 3aafcd8e41..4d1838bb01 100644 --- a/invisible_cities/config/hypathia.conf +++ b/invisible_cities/config/hypathia.conf @@ -25,11 +25,6 @@ event_range = 0, 2 thr_csum_s1 = 0.5 * pes thr_csum_s2 = 2.0 * pes -# Set thresholds for SiPM -thr_sipm = 1.0 * pes -thr_sipm_type = common - - # Set parameters to search for S1 # Notice that in MC file S1 is in t=100 mus s1_tmin = 99 * mus # position of S1 in MC files at 100 mus @@ -47,8 +42,12 @@ s2_lmin = 80 # 100 x 25 = 2.5 mus s2_lmax = 100000 # maximum value of S2 width s2_rebin_stride = 40 # Rebin by default, 40 25 ns time bins to make one 1us time bin -# Set S2Si parameters -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 = common + , thr_sipm = 1.0 * pes # Threshold for each SiPM time bin + , thr_sipm_s2 = 5 * pes # Threshold for the full sipm waveform + , detector_db = detector_db + , run_number = run_number) diff --git a/invisible_cities/config/irene.conf b/invisible_cities/config/irene.conf index 9b536c7c1c..bc5bdff334 100644 --- a/invisible_cities/config/irene.conf +++ b/invisible_cities/config/irene.conf @@ -26,10 +26,6 @@ thr_maw = 3 * adc thr_csum_s1 = 0.5 * pes thr_csum_s2 = 1.0 * pes -# Set thresholds for SiPM -thr_sipm = 3.5 * pes -thr_sipm_type = common - # Set parameters to search for S1 # Notice that in MC file S1 is in t=100 mus s1_tmin = 99 * mus # position of S1 in MC files at 100 mus @@ -47,8 +43,12 @@ s2_lmin = 100 # 100 x 25 = 2.5 mus s2_lmax = 100000 # maximum value of S2 width s2_rebin_stride = 40 # Rebin by default, 40 25 ns time bins to make one 1us time bin -# Set S2Si parameters -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 = common + , thr_sipm = 3.5 * pes # Threshold for each SiPM time bin + , thr_sipm_s2 = 10 * pes # Threshold for the full sipm waveform + , detector_db = detector_db + , run_number = run_number) diff --git a/invisible_cities/config/irene_pyrrha.conf b/invisible_cities/config/irene_pyrrha.conf new file mode 100644 index 0000000000..8fe2a9d230 --- /dev/null +++ b/invisible_cities/config/irene_pyrrha.conf @@ -0,0 +1,55 @@ +files_in = '$ICDIR/database/test_data/electrons_40keV_z25_RWF.h5' + +# REPLACE /tmp with your output directory +file_out = '/tmp/electrons_40keV_z25_PMP.h5' + +# compression library +compression = 'ZLIB4' + +# run number 0 is for MC +run_number = 0 +detector_db = 'new' + +# How frequently to print events +print_mod = 1 + +# max number of events to run +event_range = 1 + +n_baseline = 28000 # for a window of 800 mus + +# Set MAW for calibrated sum +n_maw = 100 +thr_maw = 3 * adc + +# Set thresholds for calibrated sum +thr_csum_s1 = 0.5 * pes +thr_csum_s2 = 1.0 * pes + +# Set parameters to search for S1 +# Notice that in MC file S1 is in t=100 mus +s1_tmin = 99 * mus # position of S1 in MC files at 100 mus +s1_tmax = 101 * mus # change tmin and tmax if S1 not at 100 mus +s1_stride = 4 # minimum number of 25 ns bins in S1 searches +s1_lmin = 8 # 8 x 25 = 200 ns +s1_lmax = 20 # 20 x 25 = 500 ns +s1_rebin_stride = 1 # Do not rebin S1 by default + +# Set parameters to search for S2 +s2_tmin = 101 * mus # assumes S1 at 100 mus, change if S1 not at 100 mus +s2_tmax = 1199 * mus # end of the window +s2_stride = 40 # 40 x 25 = 1 mus +s2_lmin = 100 # 100 x 25 = 2.5 mus +s2_lmax = 100000 # maximum value of S2 width +s2_rebin_stride = 40 # Rebin by default, 40 25 ns time bins to make one 1us time bin + +pmt_samp_wid = 25 * ns +sipm_samp_wid = 1 * mus + +cutting_function = pyrrha +cutting_params = dict( selection_method = median_std_method + , selection_kwargs = {'nsigma': 3} # method for selecting energetic SiPMs + , proximity_threshold = 25 # energetic SiPMs with no neighbors within this distance (in mm) are discarded + , padding_radius = 50 # amount of padding around selected SiPMs (in mm) + , run_number = run_number + , detector_db = detector_db) diff --git a/invisible_cities/reco/peak_functions.py b/invisible_cities/reco/peak_functions.py index 30ebaba508..3f796176fc 100644 --- a/invisible_cities/reco/peak_functions.py +++ b/invisible_cities/reco/peak_functions.py @@ -10,13 +10,14 @@ import numpy as np -from .. core import system_of_units as units -from .. evm .ic_containers import ZsWf -from .. evm .pmaps import S1 -from .. evm .pmaps import S2 -from .. evm .pmaps import PMap -from .. evm .pmaps import PMTResponses -from .. evm .pmaps import SiPMResponses +from .. core import system_of_units as units +from .. evm .ic_containers import ZsWf +from .. evm .pmaps import S1 +from .. evm .pmaps import S2 +from .. evm .pmaps import PMap +from .. evm .pmaps import PMTResponses +from .. evm .pmaps import SiPMResponses +from .. types .symbols import CutAlgo def indices_and_wf_above_threshold(wf, thr): @@ -28,6 +29,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 @@ -75,14 +77,14 @@ def build_pmt_responses(indices, times, widths, ccwf, return pk_times, pk_widths, PMTResponses(pmt_ids, pmt_wfs) -def build_sipm_responses(indices, times, widths, - sipm_wfs, sipm_ids, rebin_stride, thr_sipm_s2): - _, _, sipm_wfs_ = pick_slice_and_rebin(indices , times, widths, +def build_sipm_responses(indices, times, widths, sipm_wfs, + sipm_ids, rebin_stride, sipm_selection_algo): + _, _, sipm_wfs = pick_slice_and_rebin(indices , times, widths, sipm_wfs, rebin_stride, pad_zeros = False) (sipm_idx, - sipm_wfs) = select_wfs_above_time_integrated_thr(sipm_wfs_, - thr_sipm_s2) + sipm_wfs) = sipm_selection_algo(sipm_wfs) + return SiPMResponses(sipm_ids[sipm_idx], sipm_wfs) @@ -93,7 +95,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): + sipm_selection_algo = CutAlgo.no_cut): sipm_pmt_bin_ratio = int(sipm_samp_wid/pmt_samp_wid) (pk_times , pk_widths, @@ -107,7 +109,7 @@ def build_peak(indices, times, widths * sipm_pmt_bin_ratio, sipm_wfs, sipm_ids, rebin_stride // sipm_pmt_bin_ratio, - thr_sipm_s2) + sipm_selection_algo) else: sipm_r = SiPMResponses.build_empty_instance() @@ -120,7 +122,9 @@ def find_peaks(ccwfs, index, Pk, pmt_ids, sipm_ids=None, pmt_samp_wid = 25*units.ns, sipm_samp_wid = 1*units.mus, - sipm_wfs=None, thr_sipm_s2=0): + sipm_wfs=None, + sipm_selection_algo = CutAlgo.no_cut): + ccwfs = np.array(ccwfs, ndmin=2) peaks = [] @@ -136,20 +140,20 @@ def find_peaks(ccwfs, index, rebin_stride, with_sipms, Pk, pmt_samp_wid, sipm_samp_wid, - sipm_wfs, thr_sipm_s2) + sipm_wfs, sipm_selection_algo) 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, sipm_ids, - pmt_samp_wid, sipm_samp_wid): + s1_params, s2_params, pmt_ids, sipm_ids, + pmt_samp_wid, sipm_samp_wid, sipm_selection_algo = CutAlgo.no_cut): 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_ids=sipm_ids, sipm_wfs = sipm_zs_wf, - thr_sipm_s2 = thr_sipm_s2, + sipm_selection_algo = sipm_selection_algo, 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 dfd2116b91..4f54652e02 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 select_cutting_algorithm 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 @@ -347,14 +350,22 @@ def test_build_pmt_responses(wf_with_indices): def test_build_sipm_responses(wf_with_indices): times, widths, wfs, indices = wf_with_indices ids = np.arange(wfs.shape[0]) - wfs_slice = wfs[:, indices] - peak_integrals = wfs_slice.sum(axis=1) - below_thr_index = np.argmin (peak_integrals) + 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 - sipm_ids = np.arange(len(wfs)) - sipm_r = pf.build_sipm_responses(indices, times, widths, - wfs, sipm_ids, 1, thr) + thr = peak_integrals[below_thr_index] * 1.000001 + sipm_ids = np.arange(len(wfs)) + + cut_params = dict(detector_db = 'None', + thr_sipm_s2 = thr, + thr_sipm = 0, + thr_sipm_type = SiPMThreshold.common, + run_number = 0) + sipm_selection_algo = select_cutting_algorithm(CutAlgo.threshold, **cut_params) + + sipm_r = pf.build_sipm_responses(indices, times, widths, + wfs, sipm_ids, 1, sipm_selection_algo) expected_ids = np.delete( ids, below_thr_index) expected_wfs = np.delete(wfs_slice, below_thr_index, axis=0) @@ -393,11 +404,12 @@ def test_build_peak_development(pmt_and_sipm_wfs_with_indices, peak = pf.build_peak(pmt_indices, times, widths, pmt_wfs, pmt_ids, sipm_ids, - rebin_stride = rebin, - with_sipms = with_sipms, - Pk = Pk, - sipm_wfs = sipm_wfs, - thr_sipm_s2 = -1) + rebin_stride = rebin, + with_sipms = with_sipms, + Pk = Pk, + sipm_wfs = sipm_wfs, + sipm_selection_algo = select_cutting_algorithm(CutAlgo.no_cut) + ) assert_Peak_equality(peak, expected_peak) @@ -474,8 +486,9 @@ def test_find_peaks_s2_style(pmt_and_sipm_wfs_with_indices): time_range, length_range, stride, rebin_stride, S2, pmt_ids, sipm_ids, - sipm_wfs = sipm_wfs, - thr_sipm_s2 = -1) + sipm_wfs = sipm_wfs, + sipm_selection_algo = select_cutting_algorithm(CutAlgo.no_cut) + ) (rebinned_times, rebinned_widths, @@ -503,10 +516,11 @@ 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, sipm_ids = sipm_ids, pmt_samp_wid = pmt_samp_wid , - sipm_samp_wid = sipm_samp_wid) + sipm_samp_wid = sipm_samp_wid, + sipm_selection_algo = select_cutting_algorithm(CutAlgo.no_cut) + ) (rebinned_times , rebinned_widths, diff --git a/invisible_cities/reco/wfm_functions.py b/invisible_cities/reco/wfm_functions.py index 4ecced9e21..d13c15115c 100644 --- a/invisible_cities/reco/wfm_functions.py +++ b/invisible_cities/reco/wfm_functions.py @@ -3,10 +3,16 @@ authors: J.J. Gomez-Cadenas, G. Martinez """ import numpy as np +from typing import Optional +from typing import Tuple from .. core.core_functions import define_window 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 +135,216 @@ 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 : float) -> 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. + + Returns + ------- + 2D array of shape (n_sipms, n_time_bins) containing the input waveforms with entries below threshold set to zero. + """ + if isinstance(zeroing_thr, (int, float)): + thr = np.full((wfs.shape[0], 1), zeroing_thr) + else: + thr = np.reshape(zeroing_thr, (wfs.shape[0], 1)) + + 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 : float, + integration_thr : float) -> 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. + integration_thr : Charge threshold for total SiPM waveform in PE. + + Returns + ------- + Tuple of np arrays including all passing sipm ids and the corresponding waveforms + """ + zwfs = zero_wfs_below_threshold(wfs, zeroing_thr) + + return select_wfs_above_time_integrated_thr(zwfs, integration_thr) + + +def top_n_method(wfs : np.ndarray, + n : int) -> 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. + + 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: + """ + Removes isolated SiPMs from a boolean selection mask of SiPMs. A selected SiPM is considered + isolated if none of the other selected SiPMs lie within the proximity_threshold of it. Isolated + SiPMs are removed from the selection. The output keeps only SiPMs that belong to a cluster + of two or more mutually nearby selected SiPMs. + + Parameters + ---------- + selected_ids : Boolean array of shape (n_sipms,) corresponding to the initial selection of SiPMs. + 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: + """ + Expands a boolean selection mask of SiPMs to include all SiPMs within a given region determined + by the padding_radius. For each selected SiPM, all SiPM within the padding_radius are added to + the selection; the result is the union of all such neighborhoods together with the original selection. + + Parameters + ---------- + selected_ids_no_isolated : Boolean array of shape (n_sipms,). True marks the SiPMs to pad around. + 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 include neighboring SiPMs around each + selected SiPM. + + Returns + ------- + sipm_ids_with_signal : Boolean array of shape (n_sipms,). True for SiPMs that are either in the + original selection or within the "padding_radius" of a selected SiPM. + """ + 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 60c8a9abf8..50ac35b497 100644 --- a/invisible_cities/reco/wfm_functions_test.py +++ b/invisible_cities/reco/wfm_functions_test.py @@ -4,6 +4,7 @@ import tables as tb from pytest import mark +from pytest import fixture from .. database import load_db @@ -48,3 +49,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 793a2e076d..d31b63454e 100644 --- a/invisible_cities/types/symbols.py +++ b/invisible_cities/types/symbols.py @@ -132,12 +132,12 @@ class SensorType(AutoNameEnumBase): class SiPMCalibMode(AutoNameEnumBase): - subtract_mode = auto() - subtract_median = auto() - subtract_mode_calibrate = auto() - subtract_mean_calibrate = auto() - subtract_median_calibrate = auto() - subtract_mode_zs = auto() + subtract_mode = auto() + subtract_median = auto() + subtract_mode_calibrate = auto() + subtract_mean_calibrate = auto() + subtract_median_calibrate = auto() + subtract_baseline_calibrate = auto() class SiPMCharge(AutoNameEnumBase): @@ -155,6 +155,13 @@ class XYReco(AutoNameEnumBase): barycenter = auto() corona = auto() +class CutAlgo(AutoNameEnumBase): + threshold = auto() + pyrrha = auto() + no_cut = auto() + +class SiPMSelectionMethod(AutoNameEnumBase): + median_std_method = auto() class WfType(AutoNameEnumBase): rwf = auto()