Adding CYTHON functions and LSC notebooks

CYTHON:
1) deleted cBLR.pyx from top level ICYTHON directory (it lives in
SIERPE)
2) Added BLR.pyx (pxd) with the latest version of the BLR algorithm
3) Added peakFuncitions to CORE, with pmt calibrated sum

Notebooks
1) Added draft notebooks in LSC directory

ICython/Core/peakFunctions.pxd

@@ -0,0 +1,21 @@
+"""
+Cython version of some peak functions
+JJGC December, 2016
+
+calibrated_pmt_sum computes the ZS calibrated sum of the PMTs
+after correcting the baseline with a MAU to suppress low frequency noise.
+input:
+CWF:    Corrected waveform (passed by BLR)
+adc_to_pes: a vector with calibration constants
+n_MAU:  length of the MAU window
+thr_MAU: treshold above MAU to select sample
+
+"""
+cimport numpy as np
+import numpy as np
+from scipy import signal
+
+cpdef calibrated_pmt_sum(double [:, :] CWF,
+                         double [:] adc_to_pes,
+                         int n_MAU=*,
+                         double thr_MAU=*)

ICython/Core/peakFunctions.pyx

@@ -0,0 +1,47 @@
+"""
+Cython version of some peak functions
+JJGC December, 2016
+"""
+cimport numpy as np
+import numpy as np
+from scipy import signal
+
+cpdef calibrated_pmt_sum(double [:, :] CWF,
+                         double [:] adc_to_pes,
+                         int n_MAU=200,
+                         double thr_MAU=5):
+    """
+    Computes the ZS calibrated sum of the PMTs
+    after correcting the baseline with a MAU to suppress low frequency noise.
+    input:
+    CWF:    Corrected waveform (passed by BLR)
+    adc_to_pes: a vector with calibration constants
+    n_MAU:  length of the MAU window
+    thr_MAU: treshold above MAU to select sample
+
+    """
+
+    cdef int j, k
+    cdef int NPMT = CWF.shape[0]
+    cdef int NWF = CWF.shape[1]
+    cdef double [:] MAU = np.array(np.ones(n_MAU),
+                                   dtype=np.double)*(1./float(n_MAU))
+
+
+    # CWF if above MAU threshold
+    cdef double [:, :] pmt_thr = np.zeros((NPMT,NWF), dtype=np.double)
+    cdef double [:] csum = np.zeros(NWF, dtype=np.double)
+    cdef double [:] MAU_pmt = np.zeros(NWF, dtype=np.double)
+
+    for j in range(NPMT):
+        # MAU for each of the PMTs, following the waveform
+        MAU_pmt = signal.lfilter(MAU,1,CWF[j,:])
+
+        for k in range(NWF):
+            if CWF[j,k] > MAU_pmt[k] + thr_MAU:
+                pmt_thr[j,k] = CWF[j,k]
+
+    for j in range(NPMT):
+        for k in range(NWF):
+            csum[k] += pmt_thr[j, k]*1./adc_to_pes[j]
+    return csum

ICython/Sierpe/BLR.pxd

@@ -0,0 +1,40 @@
+"""
+Definition file for BLR algorithm
+JJGC December 2016
+
+deconvolve_signal takes the raw signal output from a PMT FEE
+and returns the deconvolved (BLR) waveform.
+input:
+signal_daq: the raw signal (in doubles)
+n_baseline: length of the baseline to compute the mean. It tipically takes
+            most of the baseline, since the mean over the whole raw signal is
+            zero (thus the calculation including the full window gives
+                  the pedestal accurately)
+coeff_clean: the coefficient for the cleaning part of the BLR algorithm
+coef_blr: the coefficient for the accumulator part of the BLR algorithm
+thr_trigger: threshold to decide whether signal is on
+acum_discharge_length: length to dicharge the accumulator in the absence
+                       of signal.
+
+
+deconv_pmt performs the deconvolution for the PMTs of the energy plane
+input:
+pmtrwf:    the raw waveform for all the PMTs (shorts)
+coeff_c:   a vector of deconvolution coefficients (cleaning)
+coeff_blr: a vector of deconvolution coefficients (blr)
+n_baseline, thr_trigger as described above
+
+"""
+import numpy as np
+cimport numpy as np
+
+cpdef deconvolve_signal(double [:] signal_daq,
+                        int n_baseline=*,
+                        double coef_clean=*,
+                        double coef_blr=*,
+                        double thr_trigger=*,
+                        int acum_discharge_length=*)
+
+cpdef deconv_pmt(np.ndarray[np.int16_t, ndim=2] pmtrwf,
+                 double [:] coeff_c, double [:] coeff_blr,
+                 int n_baseline=*, double thr_trigger=*)

ICython/Sierpe/BLR.pyx

@@ -0,0 +1,113 @@
+import numpy as np
+cimport numpy as np
+from scipy import signal as SGN
+
+cpdef deconvolve_signal(double [:] signal_daq,
+                        int n_baseline=28000,
+                        double coef_clean=2.905447E-06,
+                        double coef_blr=1.632411E-03,
+                        double thr_trigger=5,
+                        int acum_discharge_length = 5000):
+
+    """
+    The accumulator approach by Master VHB
+    decorated and cythonized  by JJGC
+    Current version using memory views
+
+    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 = coef_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.
+    cdef int nn = 400 # fixed at 10 mus
+
+    for j in range(0,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
+
+    b_cf, a_cf = SGN.butter(1, coef_clean, 'high', analog=False);
+    signal_daq = SGN.lfilter(b_cf,a_cf,signal_daq)
+
+    cdef int k
+    j=0
+    signal_r[0] = signal_daq[0]
+    for k in range(1,len_signal_daq):
+
+        # always update signal and accumulator
+        signal_r[k] = signal_daq[k] + signal_daq[k]*(coef/2.0) +\
+                      coef * acum[k-1]
+
+        acum[k] = acum[k-1] + signal_daq[k]
+
+        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 < acum_discharge_length - 1:
+                    j = j + 1
+                else:
+                    j = acum_discharge_length - 1
+            else:
+                acum[k]=0
+                j=0
+    # return recovered signal
+    return signal_r
+
+
+cpdef deconv_pmt(np.ndarray[np.int16_t, ndim=2] pmtrwf,
+                 double [:] coeff_c, double [:] coeff_blr,
+                 int n_baseline=28000, double thr_trigger=5):
+    """
+    Deconvolution of all the PMTs in the event cython function
+    """
+
+    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 int pmt
+    for pmt in range(NPMT):
+        signal_r = deconvolve_signal(signal_i[pmt],
+                                     n_baseline=n_baseline,
+                                     coef_clean=coeff_c[pmt],
+                                     coef_blr=coeff_blr[pmt],
+                                     thr_trigger=thr_trigger)
+
+        CWF.append(signal_r)
+
+
+    return np.array(CWF)