Add Valdrada, the trigger simulation city

invisible_cities/cities/components.py

@@ -125,7 +125,7 @@ def create_timestamp(rate: float) -> float:
 
     Returns
     -------
-    Function to calculate timestamp for the given rate with 
+    Function to calculate timestamp for the given rate with
     event_number as parameter.
     """
 
@@ -142,12 +142,12 @@ def create_timestamp(rate: float) -> float:
     def create_timestamp_(event_number: Union[int, float]) -> float:
         """
         Calculates timestamp for a given Event Number and Rate.
-    
+
         Parameters
         ----------
         event_number : Union[int, float]
                        ID value of the current event.
-    
+
         Returns
         -------
         Calculated timestamp : float
@@ -326,6 +326,33 @@ def deconv_pmt(RWF):
     return deconv_pmt
 
 
+def deconv_pmt_fpga(dbfile     : str
+                   ,run_number : int
+                   ,selection  : List[int] = None
+                   ) -> Callable:
+    '''
+    Apply deconvolution as done in the FPGA.
+
+    Parameters
+    ----------
+    dbfile     : Database to be used
+    run_number : Run number of the database
+    selection  : List of PMT IDs to apply deconvolution to.
+
+    Returns
+    ----------
+    deconv_pmt : Function that will apply the deconvolution.
+    '''
+    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)[pmt_active]
+    coeff_blr  = DataPMT.coeff_blr.values.astype(np.double)[pmt_active]
+    def deconv_pmt(RWF):
+        return zip(*map(blr.deconvolve_signal_fpga, RWF[pmt_active],
+                        coeff_c                   , coeff_blr      ))
+    return deconv_pmt
+
+
 def get_run_number(h5in):
     if   "runInfo" in h5in.root.Run: return h5in.root.Run.runInfo[0]['run_number']
     elif "RunInfo" in h5in.root.Run: return h5in.root.Run.RunInfo[0]['run_number']

invisible_cities/sierpe/blr.pxd

@@ -32,4 +32,10 @@ cpdef deconvolve_signal(double [:] signal_daq,
                         double     coeff_clean     = *,
                         double     coeff_blr       = *,
                         double     thr_trigger     = *,
-                        int accum_discharge_length = *)
+                        int accum_discharge_length = *)
+
+cpdef deconvolve_signal_fpga(short [:] signal_daq
+                            ,double    coeff_clean = *
+                            ,double    coeff_blr   = *
+                            ,double    thr_trigger = *
+                            ,size_t    base_window = *)

invisible_cities/sierpe/blr.pyx

@@ -1,3 +1,4 @@
+#cython: language_level=3
 import  numpy as np
 cimport numpy as np
 from scipy import signal as SGN
@@ -59,7 +60,6 @@ cpdef deconvolve_signal(double [:] signal_daq,
 
         if (signal_daq[k] < trigger_line) and (acum[k-1] < thr_acum):
             # discharge accumulator
-
             if acum[k-1] > 1:
                 acum[k] = acum[k-1] * (1 - coef)
                 if j < accum_discharge_length - 1:
@@ -70,4 +70,138 @@ cpdef deconvolve_signal(double [:] signal_daq,
                 acum[k] = 0
                 j = 0
     # return recovered signal
-    return np.asarray(signal_r)
+    return np.asarray(signal_r)
+
+
+cpdef deconvolve_signal_fpga( short [:] signal_daq
+                            , double    coeff_clean = 2.905447E-06
+                            , double    coeff_blr   = 1.632411E-03
+                            , double    thr_trigger = 5
+                            , size_t    base_window = 1024        ):
+
+    """
+    Simulate the deconvolution process as in the daq, differences compared to
+    usual offline deconvolution:
+      - Baseline is calculated as a moving average of 1024 counts (FPGA).
+      - Flips result between raw and deconvolved signal when outside of both
+        signal and discharge regions.
+      - Discharge made at fixed value (0.995)
+      - Result is truncated to integers.
+      - ADC threshold acting as absolute threshold.
+
+    Parameters
+    ----------
+    signal_daq  : short array
+         PMT raw waveform
+    coeff_clean : double
+         Characteristic parameter of the high pass filter
+    coeff_blr   : double
+         Characteristic parameter of BLR
+    thr_trigger : double
+         Threshold in ADCs to activate BLR
+    base_window : size
+         Moving average window for baseline calculation
+
+    Returns
+    ----------
+    Tuple of arrays:
+    - Deconvolved waveform (int16  [:])
+    - Baseline             (int16  [:])
+    """
+    cdef double coef = coeff_blr
+    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 [:] signal_f = np.zeros(len_signal_daq, dtype=np.double)
+    cdef double [:] acum     = np.zeros(len_signal_daq, dtype=np.double)
+
+    cdef int j
+
+    # compute noise
+    cdef double noise =  0
+    cdef int nn = 400 # fixed at 10 mus
+
+    for j in range(nn):
+        noise += signal_daq[j] * signal_daq[j]
+    noise /= nn
+    cdef double noise_rms = np.sqrt(noise)
+
+    # trigger line
+    cdef double trigger_line = thr_trigger #* noise_rms
+    # cleaning signal
+    cdef double [:]  b_cf
+    cdef double [:]  a_cf
+
+    ### Baseline related variables
+    cdef double [:]  top   = np.zeros(len_signal_daq, dtype=np.double)
+    cdef short  [:]  aux   = np.zeros(len_signal_daq, dtype=np.int16)
+    cdef long ped          = np.sum  (signal_daq[0:base_window], dtype=np.int32)
+    cdef size_t delay      = 2 # Delay in the FPGA in the bin used for baseline substraction
+    cdef unsigned int iaux = base_window
+
+    top[0:base_window] = ped/base_window
+    aux[0:base_window] = signal_daq[0:base_window]
+
+    b_cf, a_cf = SGN.butter(1, coeff_clean, 'high', analog=False);
+    g, a1 = b_cf[0], -a_cf[-1]
+
+    ### Initiate filt
+    filt_hr   = 0
+    #1st
+    filt_x    = g  * -(signal_daq[0] - top[0])
+    filt_a1hr = a1 * filt_hr
+
+    filt_h    = filt_x + filt_a1hr
+    signal_f[0] = filt_h - filt_hr
+    filt_hr     = filt_h
+
+    #2nd
+    filt_x      = g  * -(signal_daq[1] - top[0])
+    filt_a1hr   = a1 * filt_hr
+
+    filt_h      = filt_x + filt_a1hr
+    signal_f[1] = filt_h - filt_hr
+    filt_hr     = filt_h
+
+
+    cdef int k
+
+    #Initiate BLR
+    for k in range(0, delay):
+        signal_r[k] = signal_daq[k]
+
+    for k in range(2, len_signal_daq):
+        # always update signal and accumulator
+        current     = signal_daq[k]
+        baseline    = top[k-delay]
+
+        ### High-pass filter
+        filt_x      = g  * -(current - baseline)
+        filt_a1hr   = a1 * filt_hr
+
+        filt_h      = filt_x + filt_a1hr
+        signal_f[k] = filt_h - filt_hr
+
+        ### BLR restoration
+        blr_val     = (signal_f[k] + signal_f[k]*(coef / 2) +
+                       coef * acum[k-1]) + baseline
+
+        acum[k] = acum[k-1] + signal_f[k]
+
+        if (signal_f[k] < trigger_line) and (acum[k-1] < thr_acum):
+            # discharge accumulator
+            if acum[k-1] > 1:
+                acum[k] = acum[k-1] * .995 #Fixed discharge in FPGA code
+            else:
+                acum [k] = 0
+                blr_val = current # When outside of BLR/discharge, flip to raw signal
+                if k>=base_window:
+                    ped       = ped + current - aux[iaux-base_window]
+                    aux[iaux] = current
+                    iaux     += 1
+
+        signal_r[k] = blr_val
+        top[k]      = ped/base_window
+
+    return np.asarray(signal_r, dtype='int16'), np.asarray(top, dtype='int16')