new functionalities added

dpscience · dpscience · commit 57087b89b071 · 2021-03-23T13:40:35.000+01:00
- Gaussian fit of the retrieved distiributions
- results of the Gaussian fits are shown
- IRF consisting of a superposition of Gaussian functions can be used
- mono-decay spectrum is not anymore required
diff --git a/pyDCONTINPALS/pyDCONTINPALS.py b/pyDCONTINPALS/pyDCONTINPALS.py
@@ -28,12 +28,14 @@
 from pathlib import Path
 import matplotlib.pyplot as plt
 import numpy as np
+from scipy.optimize import curve_fit
+from scipy import interpolate
 
 import pyDCONTINPALSSpecSimulator as specSimulator
 import pyDCONTINPALSInput as userInput
 
 def __information__():
-    print("#********************* pyDCONTINPALS 1.01 (20.01.2021) *********************")
+    print("#********************* pyDCONTINPALS 1.02 (16.03.2021) *********************")
     print("#**")
     print("#** Copyright (C) 2020-2021 Danny Petschke")
     print("#**")
@@ -63,6 +65,53 @@ def __licence__():
     print("#**")
     print("#*************************************************************************************************")
 
+def gaussian(x=[],amplitude=1.,loc=0.,scale=1.):
+    return amplitude*np.exp(-0.5*((x-loc)/scale)**2)
+    
+def detectPeaks(x=[],y=[]):
+    peak_list = []
+    peak_list_y = []
+    
+    der_1st = []
+    for i in range(0,len(y)-1):
+        der_1st.append((y[i+1]-y[i])/(x[i+1]-x[i]))
+        
+    #plt.plot(der_1st,'ro')
+    #plt.show()
+    
+    for i in range(0,len(der_1st)-1):
+        p1=der_1st[i]   
+        p2=der_1st[i+1]
+        
+        if np.abs(p1-p2)<=1e-5:
+            continue
+            
+        if (p1 > 0 and p2 < 0) or (p2 > 0 and p1 < 0):
+            peak_list.append(x[i+1])
+            peak_list_y.append(y[i+1])
+                
+    return peak_list,peak_list_y
+    
+def multiPeakFit(x=[],y=[]):
+    peak_list_x,peak_list_y = detectPeaks(x,y)
+    
+    results = []
+    results_uncertainties = []
+    results_curve = []
+    results_area = []
+    
+    for i in range(0,len(peak_list_x)):
+        init_vals = [peak_list_y[i],peak_list_x[i],50.]
+        
+        best_vals,covar = curve_fit(gaussian,x,y,p0=init_vals)
+        
+        results.append(best_vals)
+        results_uncertainties.append(np.sqrt(np.diag(covar)))
+        results_curve.append(gaussian(x,amplitude=best_vals[0],loc=best_vals[1],scale=best_vals[2]))
+        results_area.append(sum(gaussian(x,amplitude=best_vals[0],loc=best_vals[1],scale=best_vals[2])))
+            
+    return results,results_uncertainties,results_curve,results_area
+
 if __name__ == '__main__':
     __information__()
     
@@ -75,15 +124,16 @@ def __licence__():
         
     monoDecayTau   = userInput.__tau_monoDecaySpec_in_ps
     binWidth_in_ps = userInput.__channelResolutionInPs
-    tauGrid_in_ps  = [userInput.__gridTau_start, userInput.__gridTau_stop]
+    tauGrid_in_ps  = [userInput.__gridTau_start,userInput.__gridTau_stop]
     gridPoints     = userInput.__gridPoints
+    binFac         = userInput.__binFactor
         
     # running demo mode ?
     if userInput.__demoMode:
         numberOfBins = userInput.__bkgrd_startIndex + userInput.__bkgrd_count + 10 # set to the maximum for demonstration purposes
         
-        charactLifetimes_in_ps  = [160.0, 400.0, 2000.0]
-        contributionOfLifetimes = [0.90,  0.095, 0.005]
+        charactLifetimes_in_ps  = [180.0, 400.0, 1600.0]
+        contributionOfLifetimes = [0.30,  0.50, 0.20]
         
         specdata_sample = specSimulator.generateCompleteLTSpectrum(numberOfComponents=len(charactLifetimes_in_ps),
                                                                     binWidth_in_ps=binWidth_in_ps, 
@@ -92,49 +142,110 @@ def __licence__():
                                                                     numberOfBins=numberOfBins, 
                                                                     charactLifetimes_in_ps=charactLifetimes_in_ps, 
                                                                     contributionOfLifetimes=contributionOfLifetimes,
-                                                                    irf_tZero_in_ps=3500.0,
+                                                                    irf_tZero_in_ps=userInput.__t_zero*userInput.__channelResolutionInPs,
                                                                     irf_fwhm_in_ps=230.0,
                                                                     noise=True,
                                                                     noiseLevel=1.0)
         
-        specdata_ref = specSimulator.generateCompleteLTSpectrum(numberOfComponents=1,
+        specdata_ref = specSimulator.generateCompleteLTSpectrum(numberOfComponents=3,
                                                                     binWidth_in_ps=binWidth_in_ps, 
                                                                     integralCounts=5000000, 
                                                                     constBkgrdCounts=5, 
                                                                     numberOfBins=numberOfBins, 
-                                                                    charactLifetimes_in_ps=[monoDecayTau], 
-                                                                    contributionOfLifetimes=[1.0],
-                                                                    irf_tZero_in_ps=3500.0,
+                                                                    charactLifetimes_in_ps=[userInput.__tau_monoDecaySpec_in_ps], 
+                                                                    contributionOfLifetimes=[1.],
+                                                                    irf_tZero_in_ps=userInput.__t_zero*userInput.__channelResolutionInPs,
                                                                     irf_fwhm_in_ps=230.0,
                                                                     noise=True,
                                                                     noiseLevel=1.0)
     else:
-        __ ,specdata_sample = np.loadtxt(userInput.__filePathSpec, delimiter=userInput.__specDataDelimiter, skiprows=userInput.__skipRows, unpack=True, dtype='float');
-        __ ,specdata_ref    = np.loadtxt(userInput.__filePathRefOrIRFSpec, delimiter=userInput.__refDataDelimiter, skiprows=userInput.__skipRows, unpack=True, dtype='float');
-
-        numberOfBins = len(specdata_sample)
+        specdata_sample = np.loadtxt(userInput.__filePathSpec, delimiter=userInput.__specDataDelimiter, skiprows=userInput.__skipRows, unpack=True, dtype='float')
+        
+        if userInput.__usingRefSpectrum:
+            specdata_ref = np.loadtxt(userInput.__filePathRefOrIRFSpec, delimiter=userInput.__refDataDelimiter, skiprows=userInput.__skipRows, unpack=True, dtype='float')
+        else:
+            specdata_ref = np.zeros(len(specdata_sample))
+        
+    spec_rebinned = np.zeros(int(np.ceil(len(specdata_sample)/binFac)))
+    
+    i_bin = 0
+    sum_bins = 0.
+    for i in range(0,len(specdata_sample)):
+        sum_bins += specdata_sample[i]
+        if not i%binFac:
+            spec_rebinned[i_bin] = sum_bins
+            i_bin += 1
+            sum_bins = 0.
+            
+    specdata_sample = spec_rebinned
+    
+    if not userInput.__usingRefSpectrum:
+        for i in range(len(specdata_ref)):
+            t = (i-userInput.__t_zero)*binWidth_in_ps
+            
+            a = 0
+            for ii in range(len(userInput.__irf_fwhm)):
+                a += userInput.__irf_intensity[ii]*np.exp(-0.5*((t-userInput.__irf_t0[ii])/(userInput.__irf_fwhm[ii]/2.3548))**2)
+                
+            specdata_ref[i] = a
+        
+    irf_rebinned = np.zeros(int(np.ceil(len(specdata_ref)/binFac)))
+    
+    i_bin = 0
+    sum_bins = 0.
+    for i in range(0,len(specdata_ref)):
+        sum_bins += specdata_ref[i]
+        if not i%binFac:
+            irf_rebinned[i_bin] = sum_bins
+            i_bin += 1
+            sum_bins = 0.
+            
+    specdata_ref = irf_rebinned
+    
+    # adjust values to rebinned data ...
+    binWidth_in_ps *= binFac
+    
+    roi_start = int(np.ceil(userInput.__roi_start/binFac))
+    roi_end   = int(np.ceil(userInput.__roi_end/binFac))
+    
+    bkgrd_startIndex = int(np.ceil(userInput.__bkgrd_startIndex/binFac))
+    bkgrd_count      = int(np.ceil(userInput.__bkgrd_count/binFac))
+    
+    bkgrd_startIndex -= roi_start+1
+    
+    t_zero_chn = int(np.ceil(userInput.__t_zero/binFac))
+    
+    spec_data_roi = specdata_sample[roi_start:roi_end]
+    irf_data_roi  = specdata_ref[roi_start:roi_end]
+    
+    numberOfBins = len(spec_data_roi)
         
     # catch general limitations given by CONTIN-PALS
     assert numberOfBins <= 4000 and numberOfBins >= 10 
     assert gridPoints >= 10 and gridPoints <= 100
     assert binWidth_in_ps >= 10.0
     
-    if monoDecayTau == 0.0:
+    if not userInput.__usingRefSpectrum:
         monoDecayTau = 1E-6
 
     # show the data
     fig, ax = plt.subplots()
-    plt.semilogy(specdata_sample, 'bo', specdata_ref, 'ro')
-    ax.set_ylabel('counts [a.u.]')
-    ax.set_xlabel('channel [a.u.]')
+    
+    if userInput.__usingRefSpectrum:
+        plt.semilogy(spec_data_roi/sum(spec_data_roi), 'bo', irf_data_roi/sum(irf_data_roi), 'r-')
+    else:
+        plt.semilogy(spec_data_roi/sum(spec_data_roi), 'bo')
+        
+    ax.set_ylabel('area normalized counts [a.u.]')
+    ax.set_xlabel('channels [{} ps]'.format(binWidth_in_ps))
     plt.show()
     
     specSamp = (ctypes.c_int*numberOfBins)()
     specRef  = (ctypes.c_int*numberOfBins)()
     
     for i in range(numberOfBins):
-        specSamp[i] = int(specdata_sample[i])
-        specRef[i]  = int(specdata_ref[i])
+        specSamp[i] = int(spec_data_roi[i])
+        specRef[i]  = int(irf_data_roi[i])
     
     program          = __dllPtr.analyseData
     program.restype = ctypes.c_int
@@ -148,8 +259,8 @@ def __licence__():
                      ctypes.c_double(tauGrid_in_ps[0]),
                      ctypes.c_double(tauGrid_in_ps[1]),
                      ctypes.c_int(gridPoints),
-                     ctypes.c_int(userInput.__bkgrd_startIndex),
-                     ctypes.c_int(userInput.__bkgrd_count))
+                     ctypes.c_int(bkgrd_startIndex),
+                     ctypes.c_int(bkgrd_count))
     
     if not result == 1:
         if result == 0:
@@ -187,40 +298,89 @@ def __licence__():
     m_yerr.restype = ctypes.c_double 
     m_res.restype  = ctypes.c_double 
         
-    x     = np.zeros(gridPoints)
-    y     = np.zeros(gridPoints)
-    yerr  = np.zeros(gridPoints)
+    x    = np.zeros(gridPoints)
+    y    = np.zeros(gridPoints)
+    yerr = np.zeros(gridPoints)
     
-    res   = np.zeros(numberOfBins)
+    res  = np.zeros(numberOfBins)
     
     # visualize results
     for i in range(0, gridPoints):
         x[i]    = m_x(ctypes.c_int(i))
         y[i]    = m_y(ctypes.c_int(i))
         yerr[i] = m_yerr(ctypes.c_int(i))
         
+    # fit results to retrieve information ...
+    results,uncertainties,fitData,intensities = multiPeakFit(x,y)    
+    
     for i in range(0, numberOfBins):
         res[i]  = m_res(ctypes.c_int(i))
         
     fig, ax = plt.subplots()
-    plt.errorbar(x, y, yerr=yerr)
+    plt.errorbar(x,y,yerr=yerr,marker='s',mfc='red',mec='blue', ms=2, mew=4,label="CONTIN-PALS results")
+    
+    for i in range(0,len(fitData)):
+        plt.plot(x,fitData[i],'r--',lw=1,label='Gaussian fit at {} ps'.format(results[i][1]))
+        
+        x_v = [results[i][1],results[i][1]]
+        y_v = [0,results[i][0]]
+        
+        plt.plot(x_v,y_v,'r-',lw=2)
     
     # running demo mode ?
     if userInput.__demoMode:
         for xc in charactLifetimes_in_ps:
             plt.axvline(x=xc, color='k', linestyle='--')
         
-    ax.set_ylabel('intensity [a.u.]')
+    ax.set_ylabel('intensity pdf [a.u.]')
     ax.set_xlabel('characteristic lifetimes [ps]')
+    plt.legend(loc='best')
     plt.show()
     
     fig, ax = plt.subplots()
-    plt.plot(res, 'ro')
+    plt.plot(res,'ro',label="error normalized residuals")
+    plt.legend(loc='best')
     
-    hlines = [-2, 0 , 2]    
+    hlines = [-4,-2,0,2,4]    
     for yc in hlines:
         plt.axhline(y=yc, color='k', linestyle='--')
         
     ax.set_ylabel('confidence level [sigma]')
-    ax.set_xlabel('channel [a.u.]')
-    plt.show()
+    ax.set_xlabel('channels [{} ps]'.format(binWidth_in_ps))
+    ax.set_ylim([-6,6])
+    plt.show()
+    
+    print('')
+    print('channel-width:  {} ps (= {} x {} ps)'.format(binWidth_in_ps,binFac,userInput.__channelResolutionInPs))
+    print('')
+    print('ROI:            [{} : {}]'.format(roi_start,roi_end))
+    
+    if not userInput.__usingRefSpectrum:
+        print('t-zero channel: {}'.format(t_zero_chn))
+        
+    print('background:     [{} : {}]'.format(bkgrd_startIndex,bkgrd_startIndex+bkgrd_count))
+    print('')
+    
+    if not userInput.__usingRefSpectrum:
+        for i in range(len(userInput.__irf_fwhm)):
+            print('----- fixed Gaussian IRF components ({}/{}) -----'.format(i+1,len(userInput.__irf_fwhm)))
+            print('')
+            print('intensity:  {} %'.format(100.*userInput.__irf_intensity[i]))
+            print('t:          {} ps'.format(userInput.__irf_t0[i]))
+            print('FWHM:       {} ps'.format(userInput.__irf_fwhm[i]))
+        
+        print('')
+        
+    for i in range(0,len(results)):
+        print('------------ found component ({}/{}) ------------'.format(i+1,len(results)))
+        print('')
+        print('tau-mean:   ({} +/- {}) ps'.format(results[i][1],uncertainties[i][1]))
+        print('tau-sigma:  ({} +/- {}) ps'.format(results[i][2],uncertainties[i][2]))
+        print('')
+        print('intensity:  {} %'.format(100.*intensities[i]/sum(intensities)))
+        print('')
+    
+    
+    
+    
+    
diff --git a/pyDCONTINPALS/pyDCONTINPALSInput.py b/pyDCONTINPALS/pyDCONTINPALSInput.py
@@ -25,7 +25,12 @@
 
 # NOTE: spectrum and IRF (or mono-exponential decay spectrum) data vectors require equal length!
 
+__roi_start                 = 0
+__roi_end                   = 7400
+
 # file path (and name) to the SPECTRUM data:
+__usingRefSpectrum          = True # if FALSE the '__irfXXX' related parameters are considered
+
 __filePathSpec              = 'testData/spectrum_10ps.dat'
 __specDataDelimiter         = '\t'
 
@@ -38,16 +43,22 @@
 
 # fixed mono-decay component in units of picoseconds [ps] (1/lambda = tau):
 # Note: set to something like 1E-6 if you provide numerical IRF data as input
-__tau_monoDecaySpec_in_ps   = 182.0  #[ps]
+__tau_monoDecaySpec_in_ps   = 182.  #[ps]
+
+__t_zero                    = 2000
+__irf_fwhm                  = [270.04,498.63]
+__irf_intensity             = [0.9382,0.0618]
+__irf_t0                    = [0.,6.6]
 
 # grid of characteristic lifetimes with equally distributed grid points defining the resulting intensity spectrum
-__gridTau_start             = 50.0   # [ps]
+__gridTau_start             = 10.0   # [ps]
 __gridTau_stop              = 3000.0 # [ps]
 __gridPoints                = 100    # 10 ... 100 Note: this value is internally limited to 100 by CONTIN
 
 # channel/bin resolution [ps]
-__channelResolutionInPs     = 50.0  # >= 10 ... Note: this value is internally limited by CONTIN
+__channelResolutionInPs     = 5.  # >= 10 ... Note: this value is internally limited by CONTIN
+__binFactor                 = 1
 
 # background estimation:
-__bkgrd_startIndex          = 800;
-__bkgrd_count               = 190; # number of channels with respect to the 'startIndex'
+__bkgrd_startIndex          = 6500;
+__bkgrd_count               = 900; # number of channels with respect to the 'startIndex'
