Blr pedestal and performance

invisible_cities/cities/components.py

@@ -257,18 +257,17 @@ def find_signal(wfs: pd.Series) -> List[int]:
 
 
 # TODO: consider caching database
-def deconv_pmt(dbfile, run_number, n_baseline, selection=None):
+def deconv_pmt(dbfile, run_number, n_baseline,
+               selection=None, pedestal_function=csf.means):
     DataPMT    = load_db.DataPMT(dbfile, run_number = run_number)
     pmt_active = np.nonzero(DataPMT.Active.values)[0].tolist() if selection is None else selection
     coeff_c    = DataPMT.coeff_c  .values.astype(np.double)
     coeff_blr  = DataPMT.coeff_blr.values.astype(np.double)
 
     def deconv_pmt(RWF):
-        return blr.deconv_pmt(RWF,
-                              coeff_c,
-                              coeff_blr,
-                              pmt_active = pmt_active,
-                              n_baseline = n_baseline)
+        CWF = pedestal_function(RWF[:, :n_baseline]) - RWF
+        return np.array(tuple(map(blr.deconvolve_signal, CWF[pmt_active],
+                                  coeff_c              , coeff_blr      )))
     return deconv_pmt
 
 

invisible_cities/reco/sensor_functions_test.py

@@ -11,6 +11,7 @@
 from .                   import wfm_functions as wfm
 from .  sensor_functions import convert_channel_id_to_IC_id
 from .  sensor_functions import simulate_pmt_response
+from .. reco             import calib_sensors_functions as csf
 
 
 def test_cwf_blr(dbnew, electron_MCRD_file):
@@ -23,19 +24,16 @@ def test_cwf_blr(dbnew, electron_MCRD_file):
     coeff_c       = abs(DataPMT.coeff_c .values)
     adc_to_pes    = abs(DataPMT.adc_to_pes.values)
     single_pe_rms = abs(DataPMT.Sigma.values)
+    thr_trigger   = 5
 
     with tb.open_file(electron_MCRD_file, 'r') as h5in:
         event = 0
         pmtrd = h5in.root.pmtrd
         dataPMT, BLR = simulate_pmt_response(event, pmtrd, adc_to_pes, single_pe_rms)
-        RWF = dataPMT.astype(np.int16)
-
-        CWF = blr.deconv_pmt(RWF,
-                             coeff_c,
-                             coeff_blr,
-                             pmt_active,
-                             n_baseline=28000,
-                             thr_trigger=5)
+        ZWF = csf.means(dataPMT[:, :28000]) - dataPMT
+        CWF = np.array(tuple(map(blr.deconvolve_signal, ZWF[pmt_active] ,
+                                 coeff_c              , coeff_blr       ,
+                                 np.repeat(thr_trigger, len(pmt_active)))))
 
         diff = wfm.compare_cwf_blr(cwf         = [CWF],
                                    pmtblr      = [BLR],

invisible_cities/reco/wfm_functions.py

@@ -5,6 +5,7 @@
 import numpy as np
 
 from .. core.core_functions import define_window
+from .. reco                import calib_sensors_functions as csf
 from .. sierpe              import blr
 
 def to_adc(wfs, adc_to_pes):
@@ -102,10 +103,13 @@ def cwf_from_rwf(pmtrwf, event_list, calib_vectors, deconv_params):
 
     CWF=[]
     for event in event_list:
-        CWF.append(blr.deconv_pmt(pmtrwf[event], calib_vectors.coeff_c,
-                             calib_vectors.coeff_blr,
-                             n_baseline=deconv_params.n_baseline,
-                             thr_trigger=deconv_params.thr_trigger))
+        pmt_evt = pmtrwf[event]
+        ZWF     = csf.means(pmt_evt[:, :deconv_params.n_baseline]) - pmt_evt
+        rep_thr = np.repeat(deconv_params.thr_trigger, ZWF.shape[0])
+        CWF.append(np.array(tuple(map(blr.deconvolve_signal, ZWF,
+                                      calib_vectors.coeff_c     ,
+                                      calib_vectors.coeff_blr   ,
+                                      rep_thr                   ))))
     return CWF
 
 

invisible_cities/sierpe/blr.pxd

@@ -29,16 +29,7 @@ import numpy as np
 cimport numpy as np
 
 cpdef deconvolve_signal(double [:] signal_daq,
-                        int        n_baseline      = *,
                         double     coeff_clean     = *,
                         double     coeff_blr       = *,
                         double     thr_trigger     = *,
-                        int accum_discharge_length = *)
-
-cpdef deconv_pmt(np.ndarray[np.int16_t, ndim=2] pmtrwf,
-                 double [:]                     coeff_c,
-                 double [:]                     coeff_blr,
-                 list                           pmt_active             = *,
-                 int                            n_baseline             = *,
-                 double                         thr_trigger            = *,
-                 int                            accum_discharge_length = *)
+                        int accum_discharge_length = *)

invisible_cities/sierpe/blr.pyx

@@ -3,7 +3,6 @@ cimport numpy as np
 from scipy import signal as SGN
 
 cpdef deconvolve_signal(double [:] signal_daq,
-                        int    n_baseline             = 28000,
                         double coeff_clean            = 2.905447E-06,
                         double coeff_blr              = 1.632411E-03,
                         double thr_trigger            =     5,
@@ -17,28 +16,16 @@ cpdef deconvolve_signal(double [:] signal_daq,
     In this version the recovered signal and the accumulator are
     always being charged. At the same time, the accumulator is being
     discharged when there is no signal. This avoids runoffs
-    The baseline is computed using a window of 700 mus (by default)
-    which should be good for Na and Kr
     """
 
     cdef double coef = coeff_blr
-    cdef int nm = n_baseline
     cdef double thr_acum = thr_trigger / coef
     cdef int len_signal_daq = len(signal_daq)
 
     cdef double [:] signal_r = np.zeros(len_signal_daq, dtype=np.double)
     cdef double [:] acum     = np.zeros(len_signal_daq, dtype=np.double)
 
     cdef int j
-    cdef double baseline = 0
-
-    for j in range(0,nm):
-        baseline += signal_daq[j]
-    baseline /= nm
-
-    # reverse sign of signal and subtract baseline
-    for j in range(0,len_signal_daq):
-        signal_daq[j] = baseline - signal_daq[j]
 
     # compute noise
     cdef double noise =  0
@@ -83,49 +70,4 @@ cpdef deconvolve_signal(double [:] signal_daq,
                 acum[k] = 0
                 j = 0
     # return recovered signal
-    return np.asarray(signal_r)
-
-
-cpdef deconv_pmt(np.ndarray[np.int16_t, ndim=2] pmtrwf,
-                 double [:]                     coeff_c,
-                 double [:]                     coeff_blr,
-                 list                           pmt_active             =    [],
-                 int                            n_baseline             = 28000,
-                 double                         thr_trigger            =     5,
-                 int                            accum_discharge_length =  5000):
-    """
-    Deconvolve all the PMTs in the event.
-    :param pmtrwf: array of PMTs holding the raw waveform
-    :param coeff_c:     cleaning coefficient
-    :param coeff_blr:   deconvolution coefficient
-    :param pmt_active:  list of active PMTs (by id number). An empt list
-                        implies that all PMTs are active
-    :param n_baseline:  number of samples taken to compute baseline
-    :param thr_trigger: threshold to start the BLR process
-
-    :returns: an array with deconvoluted PMTs. If PMT is not active
-              wvfs are removed.
-    """
-
-    cdef int NPMT = pmtrwf.shape[0]
-    cdef int NWF  = pmtrwf.shape[1]
-    cdef double [:, :] signal_i = pmtrwf.astype(np.double)
-    cdef double [:]    signal_r = np.zeros(NWF, dtype=np.double)
-    CWF = []
-
-    cdef list PMT = list(range(NPMT))
-    if len(pmt_active) > 0:
-        PMT = pmt_active
-
-    cdef int pmt
-    for pmt in PMT:
-        signal_r = deconvolve_signal(signal_i[pmt],
-                                     n_baseline             = n_baseline,
-                                     coeff_clean            = coeff_c[pmt],
-                                     coeff_blr              = coeff_blr[pmt],
-                                     thr_trigger            = thr_trigger,
-                                     accum_discharge_length = accum_discharge_length)
-
-        CWF.append(signal_r)
-
-    return np.array(CWF)
+    return np.asarray(signal_r)

invisible_cities/sierpe/blr_test.py

@@ -7,11 +7,12 @@
 from pytest import mark
 from flaky  import flaky
 
-from . import blr
+from .. reco import calib_sensors_functions as csf
+from .       import blr
 
 
 deconv_params = namedtuple("deconv_params",
-                           "n_baseline coeff_clean coeff_blr "
+                           "coeff_clean coeff_blr "
                            "thr_trigger accum_discharge_length")
 
 
@@ -22,17 +23,17 @@ def sin_wf_params():
     coeff_blr              = 1e-3
     thr_trigger            = 1e-3
     accum_discharge_length = 10
-    return deconv_params(n_baseline, coeff_clean, coeff_blr,
-                         thr_trigger, accum_discharge_length)
+    return n_baseline, deconv_params(coeff_clean, coeff_blr,
+                                     thr_trigger, accum_discharge_length)
 
 
 @fixture
 def sin_wf(sin_wf_params):
+    n_baseline, _   = sin_wf_params
     baseline        = np.random.uniform(0, 100)
-    wf              = np.full(sin_wf_params.n_baseline, baseline)
-    start           = np.random.choice (sin_wf_params.n_baseline // 2)
-    length          = np.random.randint(sin_wf_params.n_baseline // 50,
-                                        sin_wf_params.n_baseline // 2)
+    wf              = np.full(n_baseline, baseline)
+    start           = np.random.choice (n_baseline // 2)
+    length          = np.random.randint(n_baseline // 50, n_baseline // 2)
     stop            = start + length
 
     # minus sign to simulate inverted pmt response
@@ -46,27 +47,30 @@ def ad_hoc_blr_signals(example_blr_wfs_filename):
         rwf    = file.root.RWF[:]
         blrwf  = file.root.BLR[:]
         attrs  = file.root.BLR.attrs
-        params = deconv_params(attrs.n_baseline            ,
-                               attrs.coeff_c               ,
+        params = deconv_params(attrs.coeff_c               ,
                                attrs.coeff_blr             ,
                                attrs.thr_trigger           ,
                                attrs.accum_discharge_length)
-        return rwf, blrwf, params
+        return rwf, blrwf, attrs.n_baseline, params
 
 
 def test_deconvolve_signal_positive_integral(sin_wf, sin_wf_params):
     # The RWF should have null integral because contains roughly
     # the same number of positive and negative samples.
     # The CWF on the other hand should contain mostly positive
     # samples, therefore the integral should be positive.
-    blr_wf = blr.deconvolve_signal(sin_wf, **sin_wf_params._asdict())
+    n_baseline, params = sin_wf_params
+    sin_wf_noped = np.mean(sin_wf[:n_baseline]) - sin_wf
+    blr_wf = blr.deconvolve_signal(sin_wf_noped, **params._asdict())
     assert np.sum(blr_wf > 0)
 
 
 def test_deconvolve_signal_baseline_is_recovered(sin_wf, sin_wf_params):
     # The RWF contains a baseline. The deconvolution process should
     # remove it. We take the baseline as the most repeated value.
-    blr_wf      = blr.deconvolve_signal(sin_wf, **sin_wf_params._asdict())
+    n_baseline, params = sin_wf_params
+    sin_wf_noped = np.mean(sin_wf[:n_baseline]) - sin_wf
+    blr_wf      = blr.deconvolve_signal(sin_wf_noped, **params._asdict())
     (entries,
      amplitude) = np.histogram(blr_wf, 200)
     baseline    = amplitude[np.argmax(entries)]
@@ -76,7 +80,7 @@ def test_deconvolve_signal_baseline_is_recovered(sin_wf, sin_wf_params):
 @mark.slow
 @flaky(max_runs=3, min_passes=3)
 def test_deconvolve_signal_ad_hoc_signals(ad_hoc_blr_signals):
-    all_rwfs, all_true_blr_wfs, params = ad_hoc_blr_signals
+    all_rwfs, all_true_blr_wfs, n_baseline, params = ad_hoc_blr_signals
 
     # This test takes long, so we pick just one waveform.
     # Its exhaustiveness relies on repeated test runs.
@@ -86,8 +90,9 @@ def test_deconvolve_signal_ad_hoc_signals(ad_hoc_blr_signals):
 
     rwf         = all_rwfs        [evt_no, pmt_no]
     true_blr_wf = all_true_blr_wfs[evt_no, pmt_no]
-    blr_wf = blr.deconvolve_signal(rwf.astype(np.double),
-                                   n_baseline             = params.n_baseline,
+
+    cwf         = np.mean(rwf[:n_baseline]) - rwf
+    blr_wf = blr.deconvolve_signal(cwf,
                                    coeff_clean            = params.coeff_clean[pmt_no],
                                    coeff_blr              = params.coeff_blr  [pmt_no],
                                    thr_trigger            = params.thr_trigger,
@@ -97,29 +102,28 @@ def test_deconvolve_signal_ad_hoc_signals(ad_hoc_blr_signals):
 
 @mark.slow
 def test_deconv_pmt_ad_hoc_signals_all(ad_hoc_blr_signals):
-    all_rwfs, all_true_blr_wfs, params = ad_hoc_blr_signals
-    pmt_active                         = [] # Means all
+    all_rwfs, all_true_blr_wfs, n_baseline, params = ad_hoc_blr_signals
 
     # This test takes long, so we pick a random event.
     # Its exhaustiveness relies on repeated test runs.
     evt_no           = np.random.choice(all_rwfs.shape[0])
     evt_rwfs         = all_rwfs        [evt_no]
     evt_true_blr_wfs = all_true_blr_wfs[evt_no]
 
-    blr_wfs = blr.deconv_pmt(evt_rwfs,
-                             params.coeff_clean,
-                             params.coeff_blr  ,
-                             pmt_active,
-                             params.n_baseline ,
-                             params.thr_trigger,
-                             params.accum_discharge_length)
+    evt_cwfs = csf.means(evt_rwfs[:, :n_baseline]) - evt_rwfs
+
+    rep_thr  = np.repeat(params.thr_trigger           , evt_cwfs.shape[0])
+    rep_acc  = np.repeat(params.accum_discharge_length, evt_cwfs.shape[0])
+    blr_wfs  = np.array(tuple(map(blr.deconvolve_signal, evt_cwfs        ,
+                                  params.coeff_clean   , params.coeff_blr,
+                                  rep_thr              , rep_acc         )))
 
     np.allclose(blr_wfs, evt_true_blr_wfs)
 
 
 @mark.slow
 def test_deconv_pmt_ad_hoc_signals_dead_sensors(ad_hoc_blr_signals):
-    all_rwfs, all_true_blr_wfs, params = ad_hoc_blr_signals
+    all_rwfs, all_true_blr_wfs, n_baseline, params = ad_hoc_blr_signals
 
     n_evts, n_pmts, _ = all_rwfs.shape
     pmt_active        = np.arange(n_pmts)
@@ -132,12 +136,11 @@ def test_deconv_pmt_ad_hoc_signals_dead_sensors(ad_hoc_blr_signals):
     evt_rwfs         = all_rwfs        [evt_no]
     evt_true_blr_wfs = all_true_blr_wfs[evt_no]
 
-    blr_wfs = blr.deconv_pmt(evt_rwfs,
-                             params.coeff_clean,
-                             params.coeff_blr  ,
-                             pmt_active.tolist(),
-                             params.n_baseline ,
-                             params.thr_trigger,
-                             params.accum_discharge_length)
+    evt_cwfs = csf.means(evt_rwfs[:, :n_baseline]) - evt_rwfs
+    rep_thr  = np.repeat(params.thr_trigger           , len(pmt_active))
+    rep_acc  = np.repeat(params.accum_discharge_length, len(pmt_active))
+    blr_wfs  = np.array(tuple(map(blr.deconvolve_signal, evt_cwfs[pmt_active],
+                                  params.coeff_clean   , params.coeff_blr    ,
+                                  rep_thr              , rep_acc             )))
 
     np.allclose(blr_wfs, evt_true_blr_wfs[pmt_active])