Update PeakOTron.py, Bootstrapping.py, Example.py, HelperFunctions.py,...

Update PeakOTron.py, Bootstrapping.py, Example.py, HelperFunctions.py, LossFunctions.py, Model_AP1.py, data/hamamatsu_pcb6/Light/PCB6_MPPC_53.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_53V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_54.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_54V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_55.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_55V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_56.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_56V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_57.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_57V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_58.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_58V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_59.5V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_59V_histo.txt, data/hamamatsu_pcb6/Light/PCB6_MPPC_60V_histo.txt
Deleted SiPMSpectra/03/F1_U28p5_dark_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U28p5_light_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U29p5_dark_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U29p5_light_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U30p5_dark_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U30p5_light_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U31p5_dark_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U31p5_light_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U32p5_dark_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/03/F1_U32p5_light_PM3325WBA0-03_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U28p1_dark_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U28p1_light_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U29p1_dark_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U29p1_light_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U30p6_dark_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U30p6_light_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U31p6_dark_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/04_ZnS/F1_U31p6_light_PM3325WBA0-04_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U28p4_dark_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U28p4_light_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U29p4_dark_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U29p4_light_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U30p4_dark_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U30p4_light_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U31p4_dark_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U31p4_light_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U32p4_dark_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/06_ZnS/F1_U32p4_light_PM3315WBA0-06_ZnS_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U28p1_dark_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U28p1_light_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U30p1_dark_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U30p1_light_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U31p6_dark_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U31p6_light_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U33p1_dark_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U33p1_light_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U35p1_dark_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/07/F1_U35p1_light_PM3315WBA0-07_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U28p1_dark_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U28p1_light_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U29p6_dark_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U29p6_light_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U31p1_dark_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U31p1_light_PM33125WBA0-08_AmpG101_660nm_00001.dat, SiPMSpectra/08/F1_U32p6_dark_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U32p6_light_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U34p1_dark_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/08/F1_U34p1_light_PM33125WBA0-08_AmpG101_660nm_00000.dat, SiPMSpectra/KETEK_spectra.opju, README.txt, RunFit.py, requirements.txt

Bootstrapping.py

+from HelperFunctions import FakeFuncCode, SelectRangeNumba, EmpiricalPPF
+from KDEpy import FFTKDE
+import numpy as np
+from scipy.interpolate import interp1d
+from iminuit import Minuit
+from astropy.stats import bootstrap, scott_bin_width
+import scipy.special as sc
+import logging
+
+            
+class BandWidthOptimiser:
+    
+    #Hobbema, Hermans, and Van den Broeck (1971) and by Duin (1976)
+    #leave-one-out cross-validation approach
+    #https://cran.r-project.org/web/packages/kedd/vignettes/kedd.pdf
+    def _Dummy(self, _, bw):
+        return bw
+    
+    
+    def _KDE(self, bw):
+
+        x_kde, y_kde = FFTKDE(kernel = self.kernel, bw=bw).fit(self.data).evaluate(self.n_kde_samples)        
+        loss = -np.mean(np.log(y_kde))
+        return loss
+
+
+        
+    def __init__(self, data, kernel="gaussian", alpha = 0.95, n_kde_samples=2**13):
+        if(alpha<=0.5):
+            raise ValueError("alpha must be above 0.5")
+        
+        self.data = data
+        PPF = EmpiricalPPF(self.data)
+        x_low, x_hi = PPF(0.5-alpha/2), PPF(0.5+alpha/2)
+        self.data = SelectRangeNumba(self.data, x_low, x_hi)
+        self.N = len(data)
+        self.kernel=kernel
+        self.n_kde_samples=n_kde_samples
+        self.last_arg = None
+        self.func_code = FakeFuncCode(self._Dummy, dock=True)
+        self.eps = np.finfo(np.float64).eps * 10
+        self.eps_inv = 1/self.eps
+        
+
+    def __call__(self, *arg):
+        self.last_arg = arg
+        loss = self._KDE(*arg)
+        return loss
+    
+    
+def Bootstrap(data, statistic, n_bootstrap, alpha=0.95, ):
+    if not (0 < alpha < 1):
+        raise ValueError("confidence level must be in (0, 1)")
+
+
+    if len(data) < 1:
+        raise ValueError("data must contain at least one measurement.")
+
+
+    boot_stat = bootstrap(data, n_bootstrap, bootfunc=statistic)
+
+    stat_data = statistic(data)
+    mean_stat = np.mean(boot_stat)
+    est_stat = 2*stat_data - mean_stat
+    std_err = np.std(boot_stat)
+    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
+    conf_interval = est_stat + z_score*np.array((-std_err, std_err))
+
+
+    return est_stat, std_err, conf_interval
+    
+    
+    
+    
+def BootstrapKDE(data,
+                 n_bootstrap,
+                 n_call=1000,
+                 n_iterations=10,
+                 kernel = "gaussian",
+                 n_kde_samples=2**14,
+                 alpha=0.95,
+                 bw_limits=(1e-6, None),
+                 limits=None,
+
+                ):
+
+    if not (0 < alpha < 1):
+        raise ValueError("Bootstrap confidence level, alpha, must be in (0, 1).")
+
+    if len(data) <= 3:
+        raise ValueError("Bootstrap data must contain at least three measurements.")
+
+    try:
+        logging.debug("Optimising bandwidth...")
+        BWO= BandWidthOptimiser(data, kernel=kernel)
+        bw_scott = scott_bin_width(data)
+        minuit_kde = Minuit(BWO, bw=bw_scott/2)
+        minuit_kde.limits["bw"]=bw_limits
+        minuit_kde.strategy=2
+        minuit_kde.migrad(ncall=n_call,
+                   iterate=n_iterations)
+
+        bw = minuit_kde.values["bw"]
+        logging.debug("Bandwidth optimised to {:3.3E} units.".format(bw))
+    except:
+        logging.warning("Bandwidth optimisation failed. Defaulting to Improved Sheather Jones...")
+        bw = "ISJ"
+
+
+    kde = FFTKDE(kernel = kernel, bw=bw).fit(data)
+
+    x_kde, y_kde_orig = kde.evaluate(n_kde_samples)
+
+
+    boot_data = bootstrap(data, n_bootstrap)
+    y_kde_bs = np.vstack([FFTKDE(kernel = kernel, bw=bw).fit(_data).evaluate(x_kde)
+                               for _data in boot_data])
+
+
+    y_kde_mean = np.mean(y_kde_bs, axis=0)
+    y_kde = 2*y_kde_orig - y_kde_mean
+    y_kde_err = np.std(y_kde_bs, axis=0)
+    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
+    y_kde_conf_interval = y_kde + z_score*np.array((-y_kde_err, y_kde_err))
+
+    if(limits is not None):
+        cond_inrange =(x_kde>np.min(limits)) & (x_kde<np.max(limits))
+        x_kde = x_kde[cond_inrange]
+        y_kde = y_kde[cond_inrange]
+        y_kde_err = y_kde_err[cond_inrange]
+
+    return x_kde, y_kde, y_kde_err, y_kde_conf_interval
+
+
+def Bootstrap(data, statistic, n_bootstrap, alpha=0.95):
+    if not (0 < alpha < 1):
+        raise ValueError("confidence level must be in (0, 1)")
+
+
+    if len(data) < 1:
+        raise ValueError("data must contain at least one measurement.")
+
+
+    boot_stat = bootstrap(data, n_bootstrap, bootfunc=statistic)
+
+    stat_data = statistic(data)
+    mean_stat = np.mean(boot_stat)
+    est_stat = 2*stat_data - mean_stat
+    std_err = np.std(boot_stat)
+    z_score = np.sqrt(2.0)*sc.erfinv(alpha)
+    conf_interval = est_stat + z_score*np.array((-std_err, std_err))
+
+
+    return est_stat, std_err, conf_interval
\ No newline at end of file

RunFit.py → Example.py

@@ -3,6 +3,7 @@ import sys
 import numpy as np
 import pandas as pd
 from PeakOTron import PeakOTron
+from joblib import dump
 
 
 print("--------------------")
@@ -15,11 +16,11 @@ f_data.SetMaxNDCRPeaks(6)
 
 out_dict = {}
 files_to_fit = []
-for root, dirs, files in os.walk("./SiPMSpectra/"):
+for root, dirs, files in os.walk("./data/hamamatsu_pcb6/Light"):
      
         for file in files:
             
-            if file.endswith(".dat") and "light" in file:
+            if file.endswith(".txt"):
                 files_to_fit.append([file, os.path.join(root, file)])
                 
             
@@ -29,37 +30,51 @@ for i, (file, _) in enumerate(files_to_fit):
     
     
     
-for i, (file, path) in enumerate(files_to_fit):
+SiPM = "PM1125NS_SBO"
     
+for i, (file, path) in enumerate(files_to_fit):
     items = file.split('_')
 
-    V = float(items[1].replace('U', '').replace('p', '.'))
-    SiPM = items[3]
+    V = float(items[2].replace('V', '').replace('p', '.'))
+        
+    
+    while (True):
+
+
+       
 
         print("\n Fit {0}: {1} \n".format(i, file))
-    data = np.loadtxt(path)
-    bw = 1.1*(data[:,0][1] -  data[:,0][0])
 
+        data = np.loadtxt(path, skiprows=8)
 
         f_data.Fit(data, 
-          tau=20,  #SLOW PULSE COMPONENT TIME CONSTANT (ns)
-          tgate=100, #GATE LENGTH (ns)
-          r_fast=0.1, #FRACTION OF FAST PULSE COMPONENT 
-          t_0 = 100, #INTEGRATION TIME BEFORE GATE (ns)
-          pedestal = 0,
-          prominence=1e-3,
-          min_n_peak_evts=200,
-          bandwidth_kde="optimise",
-          bw=bw
+              tau=21.953,  #SLOW PULSE COMPONENT TIME CONSTANT (ns)
+              tau_err=5.648,
+              t_gate=104, #GATE LENGTH (ns)
+              t_0 = 64, #INTEGRATION TIME BEFORE GATE (ns)               
 
             ) #BINNING RULE "knuth", "freedman", "scott" - use bw= #### to override. it is still required for DCR calculation. 
 
-    if(not f_data.Failed()):
-        f_data.PlotSummary(save_directory="./Results/{0}_init.png".format(os.path.splitext(file)[0]))
-        f_data.PlotFit(scaled=False, title=file, xlabel="Charge [Vs]", y_limits=[1e7, 1e11], save_directory="./Results/{0}_fit.png".format(os.path.splitext(file)[0]))
+
+        f_data.PlotFit(xlabel="Charge [Vs]", save_directory="./Results/{0}_fit.png".format(os.path.splitext(file)[0]))
+    
+        chi2, ndf = f_data.GetChi2()
+        if(V-51.1>3):
+            if(chi2/ndf<4):
+                break
+        else:
+            if(chi2/ndf<10):
+                break
+    
+    dump(f_data, "./Results/{0}".format(os.path.splitext(file)[0]))
+    
+    
     
     fit_out = {}
     fit_val, fit_err = f_data.GetFitResults()
+    for key, val in fit_val.items():
+        print("{:s} : {:3.3E}".format(key, val))
+    
     fit_out["SiPM"] = SiPM
     fit_out["V"] = V
 
@@ -77,9 +92,5 @@ for i, (file, path) in enumerate(files_to_fit):
 
 
 
-    else:
-        print("Fit failed. \n")
-
-
 df = pd.DataFrame.from_dict(out_dict)
-df.to_csv("./Results/fit_results.csv")
+df.to_csv("./fit_results_{:s}.csv".format(SiPM))

HelperFunctions.py

+from itertools import chain
+from iminuit import Minuit
+from numba import njit
+from iminuit.util import describe
+from scipy.interpolate import interp1d
+import numpy as np
+from scipy.stats import skew as sc_skew
+import matplotlib.pyplot as plt
+
+
+#HELPER FUNCTIONS
+
+
+class FakeFuncCode:
+    def __init__(self, f, prmt=None, dock=0, append=None):
+        #f can either be tuple or function object
+        self.co_varnames = describe(f)
+        self.co_argcount = len(self.co_varnames)
+        self.co_argcount -= dock
+        self.co_varnames = self.co_varnames[dock:]
+
+        if prmt is not None:  #rename parameters from the front
+            for i, p in enumerate(prmt):
+                self.co_varnames[i] = p
+
+        if isinstance(append, str): append = [append]
+
+        if append is not None:
+            old_count = self.co_argcount
+            self.co_argcount += len(append)
+            self.co_varnames = tuple(
+                list(self.co_varnames[:old_count]) +
+                append +
+                list(self.co_varnames[old_count:]))
+
+
+
+            
+def LatexFormat(value, scirange=[0.01,1000]):
+    if(np.abs(value)>scirange[0] and np.abs(value)<scirange[1]):
+        float_str = r"${:3.3f}$".format(value)
+    else:
+        try:
+            float_str = "{:3.3E}".format(value)
+            base, exponent = float_str.split("E")
+            float_str = r"${0} \times 10^{{{1}}}$".format(base, int(exponent))
+        except:
+            float_str=str(value)
+    return float_str
+
+            
+@njit
+def SelectRangeNumba(array, low_lim, hi_lim):
+    index = (array >=  low_lim) & (array <= hi_lim)
+    return array[index]
+
+ 
+def EmpiricalPPF(data):
+    x = np.sort(data)
+    n = x.size
+    y = np.arange(1, n+1) / n
+    ppf = interp1d(y, x, fill_value="extrapolate")
+    return ppf
+
+
+def EmpiricalCDF(data):
+    x = np.sort(data)
+    n = x.size
+    y = np.arange(1, n+1) / n
+    ppf = interp1d(x, y, fill_value="extrapolate")
+    return ppf
+
+
+def Logify(y):
+    log_y = np.log(y)
+    min_log_y = np.min(log_y[y>0])
+
+    return np.where(y>0, log_y, min_log_y)
+
+
+def Linear(x, m, c):
+    return m*x + c
+
+def GetStats(data):
+    return np.mean(data), np.std(data), sc_skew(data)
+
+
+
+def GP_lbda(mu, sigma, gamma):
+    lbda = 0.5*(((mu*gamma)/(sigma))- 1)
+    return lbda
+
+def GP_gain(mu, sigma, gamma):
+    lbda = GP_lbda(mu, sigma, gamma)
+    gain = (sigma**2/mu)*((1 - lbda)**2)
+    return gain
+
+def GP_muGP(mu, sigma, gamma):
+    lbda = GP_lbda(mu, sigma, gamma)
+    mu_gp = (1/(1-lbda))*(mu**2/sigma**2)
+    return mu_gp
+
+
+
+def FormatExponent(ax, axis='y'):
+
+    # Change the ticklabel format to scientific format
+    ax.ticklabel_format(axis=axis, style='sci', scilimits=(-2, 2))
+
+    # Get the appropriate axis
+    if axis == 'y':
+        ax_axis = ax.yaxis
+        x_pos = 0.0
+        y_pos = 1.0
+        horizontalalignment='left'
+        verticalalignment='bottom'
+    else:
+        ax_axis = ax.xaxis
+        x_pos = 1.0
+        y_pos = -0.05
+        horizontalalignment='right'
+        verticalalignment='top'
+
+    # Run plt.tight_layout() because otherwise the offset text doesn't update
+    plt.tight_layout()
+    ##### THIS IS A BUG 
+    ##### Well, at least it's sub-optimal because you might not
+    ##### want to use tight_layout(). If anyone has a better way of 
+    ##### ensuring the offset text is updated appropriately
+    ##### please comment!
+
+    # Get the offset value
+    offset = ax_axis.get_offset_text().get_text()
+
+    if len(offset) > 0:
+        # Get that exponent value and change it into latex format
+        minus_sign = u'\u2212'
+        expo = np.float(offset.replace(minus_sign, '-').split('e')[-1])
+        offset_text = r'x$\mathregular{10^{%d}}$' %expo
+
+        # Turn off the offset text that's calculated automatically
+        ax_axis.offsetText.set_visible(False)
+
+        # Add in a text box at the top of the y axis
+        ax.text(x_pos, y_pos, offset_text, transform=ax.transAxes,
+               horizontalalignment=horizontalalignment,
+               verticalalignment=verticalalignment)
+    return ax
+
+
+
+    
+
+

LossFunctions.py

+import numpy as np
+from scipy.interpolate import interp1d
+from iminuit.util import describe
+from HelperFunctions import FakeFuncCode
+import matplotlib.pyplot as plt
+    
+
+class UnbinnedLH:
+    
+
+    def __init__(self, f, data, bw):
+        self.f = f
+        self.x_min = np.min(data)
+        self.x_max = np.max(data)
+        self.n_bins = np.ceil((np.max(data) - np.min(data))/bw).astype(int)
+        self.x = np.linspace(self.x_min, self.x_max, 10*self.n_bins) 
+        self.data = data
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.ndof = len(self.data) - (self.func_code.co_argcount - 1)
+        self.eps = np.finfo(np.float64).eps * 10
+        self.eps_inv = 1/self.eps
+        self.log_eps = np.log(self.eps)
+        self.n_calls = 0
+        
+        
+    def Logify(self, y):
+        
+        return np.where(y>self.eps, np.log(y), self.log_eps)
+
+          
+
+    def __call__(self, *arg):
+        self.last_arg = arg
+        y_hat = self.f(self.x, *arg)
+        y_hat = np.nan_to_num(y_hat, nan=self.eps)
+        
+        log_y_hat = interp1d(self.x, 
+                             self.Logify(y_hat), 
+                             fill_value=(self.log_eps,self.log_eps), 
+                             bounds_error=False)(self.data)
+
+        
+        
+        nlogL = -np.sum(log_y_hat)
+        
+        self.n_calls+=1
+        
+#         if(self.n_calls%200==0):
+#             print(arg)
+#             plt.figure()
+#             plt.hist(self.data, bins=200, density=True)
+#             plt.plot(self.x, y_hat)
+#             plt.yscale("log")
+#             plt.show()
+       
+        return nlogL
+        
+           
+                                              
+class BinnedLH:
+    
+
+    def __init__(self, f, x, y):
+        self.f = f
+        self.x = x
+        self.y = y
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
+        self.eps = np.finfo(np.float64).eps * 10
+        self.eps_inv = 1/self.eps
+        self.log_eps = np.log(self.eps)
+        
+        
+    def Logify(self, y):
+     
+        return np.where(y>self.eps, np.log(y), self.log_eps)
+
+        
+
+    def __call__(self, *arg):
+        self.last_arg = arg
+        y_hat = self.f(self.x, *arg)
+        # y_hat = np.nan_to_num(y_hat, nan=self.eps_inv)
+        logL = self.y*(self.Logify(y_hat) - self.Logify(self.y)) + (self.y - y_hat)
+        nlogL = -np.sum(logL)
+
+        return nlogL
+    
+    
+    
+    
+class Chi2Regression:
+    
+
+    def __init__(self, f, x, y, y_err, epsilon=1.35):
+        self.f = f
+        self.x = x
+        self.y = y
+        self.y_err = y_err
+        self.eps = np.finfo(np.float64).eps * 10
+        self.y_err[self.y_err<self.eps] = self.eps
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
+      
+
+    def __call__(self, *arg):
+        
+        self.last_arg = arg
+        
+        loss = ((self.f(self.x, *arg) - self.y)/(self.y_err))**2
+        
+        return np.sum(loss)
+    
+    
+    
+    
+class HuberRegression:
+    
+
+    def __init__(self, f, x, y, y_err, delta=1.345):
+        self.f = f
+        self.x = x
+        self.y = y
+        self.y_err = y_err
+        self.delta = delta
+        self.eps = np.finfo(np.float64).eps * 10
+        self.y_err[self.y_err<self.eps] = self.eps
+        self.last_arg = None
+        self.func_code = FakeFuncCode(f, dock=True)
+        self.ndof = len(self.y) - (self.func_code.co_argcount - 1)
+        
+
+    def __call__(self, *arg):
+        
+        self.last_arg = arg
+        
+        a = abs((self.y - self.f(self.x, *arg))/self.y_err)
+        cond_flag = (a>self.delta)
+        
+        loss = np.sum((~cond_flag) * (0.5 * a ** 2) - (cond_flag) * self.delta * (0.5 * self.delta - a), -1)
+        
+        return np.sum(loss)
\ No newline at end of file

Model_AP1.py

+import numpy as np
+from scipy.signal import fftconvolve as convolve
+from scipy.stats import norm, binom, poisson,expon
+from scipy.interpolate import interp1d
+from AdditionalPDFs import borel, gpoisson
+import matplotlib.pyplot as plt
+    
+    
+def SigmaK(k, sigma_0, sigma_1):
+    return np.sqrt(sigma_0**2 + k*sigma_1**2)
+
+
+
+def dPApdPH(x, x_0, G, mu, lbda, tauAp,  pAp, tau, t_gate, k_low, k_hi):
+     
+    x_lin = np.arange(0, len(x))
+    x_pad_lin = np.arange(-100, len(x)+100)
+    f_lin = interp1d(x, x_lin, fill_value='extrapolate')
+    f_lin_inv = interp1d(x_lin, x, fill_value='extrapolate')
+    dx = abs(f_lin_inv(1) - f_lin_inv(0))
+    G_lin = f_lin(x_0+G) - f_lin(x_0)
+    x_0_lin = f_lin(x_0)
+    y_ap = np.zeros_like(x_pad_lin).astype(float)
+
+    idx_orig = np.where(np.in1d(x_pad_lin,
+                                    np.intersect1d(
+                                    x_pad_lin,
+                                    x_lin
+                                )
+                        ), True, False)
+
+    idx_x_0_lin =  np.argmin(abs(x_pad_lin - x_0_lin))
+    
+
+    exp_tgate_2tau = np.exp(t_gate/(2*tau))
+    exp_mtgate_2tau = np.exp(-t_gate/(2*tau))
+    PH_norm = (x_pad_lin - x_0_lin)/G_lin
+    PH_max = (1 - exp_mtgate_2tau)**2
+    norm_factor = PH_max/G_lin 
+                       
+                             
+    for _n_pg in range(1, int(k_hi)+1):
+        _p_pg = gpoisson.pmf(_n_pg, mu, lbda)
+        for _n_ap in range(1, min(5,_n_pg+1)):
+            
+            PH = (PH_norm - _n_pg)/_n_ap
+            cond_PH_max = (PH>0.95*PH_max)
+           
+            
+            PH[cond_PH_max] = 0.95*PH_max
+            
+
+            t_b0 = 0.5*t_gate - tau*np.arccosh(0.5*((1 - PH)*exp_tgate_2tau + exp_mtgate_2tau))
+            
+            cond_gate = (t_b0<0) | (t_b0>t_gate/2)
+            
+            dpApdt_b0 = _p_pg*binom.pmf(_n_ap, _n_pg, pAp*(1-np.exp(-t_b0/tau)))*expon.pdf(t_b0, loc=0, scale=tauAp)
+            dpApdt_b0[cond_gate] = 0
+
+            dPHdt_b0 = abs(-2*np.sinh((t_b0-0.5*t_gate)/tau)*exp_mtgate_2tau/tau)
+            
+            cond_dPHdt = (dPHdt_b0<1e-3)
+            dPHdt_b0[cond_dPHdt] = 1e-3
+            
+            dpApdPH_b0 = dpApdt_b0/dPHdt_b0
+                            
+            dpApdPH_b0[cond_PH_max] = 0
+            dpApdPH_b0[cond_gate] = 0
+            
+
+
+            y_ap+=dpApdPH_b0*(norm_factor/_n_ap)
+    
+    y_ap = y_ap[idx_orig]/dx
+    
+    return y_ap
+
+
+    
+
+
+
+def DRM_basic(x, x_0, G, mu, lbda, sigma_0, sigma_1, DCR, tau, t_gate, t_0, tauAp, pAp,  k_low, k_hi):
+
+        
+        n_pg = np.arange(1, k_hi)
+
+        f0 = gpoisson.pmf(0, mu, lbda)*norm.pdf(x, x_0, sigma_0)
+        f1s = np.asarray([
+                gpoisson.pmf(_k, mu, lbda)*norm.pdf(x, x_0 + _k*G,  SigmaK(_k, sigma_0, sigma_1))
+                for _k in n_pg
+        ])
+
+      
+        f_light = np.sum(np.vstack([f0, f1s]), axis=0)
+        f_ap = dPApdPH(x, x_0, G, mu, lbda, tauAp,  pAp, tau, t_gate, k_low, k_hi)
+        
+
+        return f_light + f_ap
+
+
+
+
+def DC_PH_range(t, t_0, r_fast, tau, t_gate):
+    
+    if((t>t_0) and (t<0)):
+        PH_min = (1-r_fast)*np.exp(t_0/tau)*(1 - np.exp(-t_gate/tau))
+        PH_max = (1-r_fast)*(1 - np.exp(-t_gate/tau))  
+    elif((t>0) and (t<t_gate)):
+        PH_min = r_fast
+        PH_max = (1 - (1-r_fast)*np.exp(-t_gate/tau))
+    else:
+        PH_min = 0
+        PH_max = 0
+
+    return PH_min, PH_max
+
+
+def dpDCRdPH(x, x_0, G, tau, t_gate, t_0):
+
+    PH_bfg_low, PH_bfg_hi = DC_PH_range(-abs(t_0)/2, -abs(t_0), 0, tau, t_gate)
+    PH_dg_low, PH_dg_hi = DC_PH_range(t_gate/2, -abs(t_0), 0, tau, t_gate)
+
+    x_norm = (x-x_0)/(G)
+
+    PHs = np.zeros_like(x_norm)
+
+    cond_bfg = (x_norm>PH_bfg_low) & (x_norm<=PH_bfg_hi)
+    cond_dg = (x_norm>PH_dg_low) & (x_norm<=PH_dg_hi)
+
+    PHs[cond_bfg] += 1/x_norm[cond_bfg]
+    PHs[cond_dg] += 1/(1 - x_norm[cond_dg])  
+
+
+    return PHs
+
+
+
+
+def ApplyDCR(x, y_light, x_0, G, DCR, t_gate, t_0, tau, lbda, k_dcr_low, k_dcr_hi):
+
+    mu_dcr = DCR*(abs(t_0) + t_gate)
+
+    x_lin = np.arange(0, len(x))
+    x_pad_lin = np.arange(-100, len(x)+100)
+    f_lin = interp1d(x, x_lin, fill_value='extrapolate')
+    f_lin_inv = interp1d(x_lin, x, fill_value='extrapolate')
+    dx = abs(f_lin_inv(1) - f_lin_inv(0))
+    G_lin = f_lin(x_0+G) - f_lin(x_0)
+    x_0_lin = f_lin(x_0)
+
+    idx_orig = np.where(np.in1d(x_pad_lin,
+                                    np.intersect1d(
+                                    x_pad_lin,
+                                    x_lin
+                                )
+                        ), True, False)
+
+    idx_x_0_lin =  np.argmin(abs(x_pad_lin - x_0_lin))
+
+
+    fs = []
+    pfs = []
+    hs = []
+    phs = []
+
+
+
+    for _n_dcr in range(0, max(2, int(k_dcr_hi))):
+
+        if(_n_dcr==0):
+
+            f = np.zeros(len(x_pad_lin))
+
+            f[idx_x_0_lin] = 1
+
+            pf = poisson.pmf(0, mu_dcr)
+
+            h = np.zeros(len(x_pad_lin))
+            
+            ph = 0
+            
+            
+        else:
+            if(_n_dcr==1):
+                f = dpDCRdPH(x_pad_lin, x_0_lin, G_lin, tau, t_gate, t_0)
+            else:
+                f = convolve(fs[1], fs[-1])[idx_x_0_lin:len(x_pad_lin)+idx_x_0_lin]
+
+            f = f/ np.trapz(f, dx = 1)
+
+            pf = poisson.pmf(_n_dcr, mu_dcr)*(borel.pmf(0, lbda)**_n_dcr)
+
+
+            
+            if(_n_dcr==1):
+                h = np.zeros(len(x_pad_lin))
+                ph=0
+            else:
+                h = dpDCRdPH(x_pad_lin, x_0_lin, G_lin*((_n_dcr-1)+1), tau, t_gate, t_0)
+                h = h/ np.trapz(h, dx = 1)
+                ph = poisson.pmf(1, mu_dcr)*borel.pmf((_n_dcr-1), lbda)/((_n_dcr-1)+1)
+
+        fs.append(f)
+        pfs.append(pf)
+        hs.append(h)
+        phs.append(ph)
+
+
+
+
+
+    fs = np.asarray(fs)
+    hs = np.asarray(hs)
+    pfs = np.expand_dims(np.asarray(pfs), -1)
+    phs = np.expand_dims(np.asarray(phs), -1)
+    y_dark = np.sum(fs*pfs, axis=0) + np.sum(hs*phs, axis=0)
+    #y_dark = y_dark/np.trapz(h, dx = 1)
+    y_dark = y_dark/np.trapz(y_dark, dx = 1)
+
+
+    y_model = convolve(y_dark,
+                          np.pad(y_light,
+                                 (100, 100),
+                                 "constant",
+                                 constant_values=(0.0,0.0)),
+                        )[idx_x_0_lin:len(x_pad_lin)+idx_x_0_lin]
+
+    y_model = y_model/np.trapz(y_model, dx = 1)/dx
+
+    y_model  = y_model[idx_orig]
+
+    return y_model
+
+
+def DRM(
+        x,
+        mu,
+        x_0,
+        G,
+        sigma_0,
+        sigma_1,
+        tauAp,
+        pAp,
+        lbda,
+        k_low,
+        k_hi,
+        k_dcr_low,
+        k_dcr_hi,
+        DCR,
+        tau,
+        t_gate,
+        t_0
+       ):
+
+
+    y_light = DRM_basic(x, x_0, G, mu, lbda, sigma_0, sigma_1, DCR, tau, t_gate, t_0, tauAp, pAp,  k_low, k_hi)
+
+    y_model = ApplyDCR(x, y_light, x_0, G, DCR, t_gate, t_0, tau, lbda, k_dcr_low, k_dcr_hi)
+
+
+
+    return y_model

README.txt

- - Run "pip install -r requirements.txt" in the terminal to install the necessary packages. 
- - Run python "RunFit.py" to test the fitting class on experimental data. 
- - The output fits will be stored in the Results folder.
- - The output fitted values will be stored in csv format in the main directory, as "fit_results.csv"
- - Please look at the "RunFit.py" example program for ideas on how to input and output from your data.

SiPMSpectra/03/F1_U32p5_dark_PM3325WBA0-03_AmpG101_660nm_00000.dat

--3.6e-010	0
--3.4e-010	0
--3.2e-010	0
--3e-010	0
--2.8e-010	0
--2.6e-010	1
--2.4e-010	7
--2.2e-010	26
--2e-010	55
--1.8e-010	69
--1.6e-010	46
--1.4e-010	57
--1.2e-010	81
--1e-010	136
--8e-011	193
--6e-011	324
--4e-011	1256
--2e-011	4653
-1.221e-024	2860
-2e-011	450
-4e-011	269
-6e-011	319
-8e-011	289
-1e-010	353
-1.2e-010	402
-1.4e-010	629
-1.6e-010	1117
-1.8e-010	1094
-2e-010	379
-2.2e-010	193
-2.4e-010	190
-2.6e-010	194
-2.8e-010	184
-3e-010	227
-3.2e-010	218
-3.4e-010	309
-3.6e-010	395
-3.8e-010	238
-4e-010	110
-4.2e-010	93
-4.4e-010	86
-4.6e-010	93
-4.8e-010	99
-5e-010	88
-5.2e-010	114
-5.4e-010	131
-5.6e-010	103
-5.8e-010	69
-6e-010	50
-6.2e-010	43
-6.4e-010	46
-6.6e-010	37
-6.8e-010	42
-7e-010	43
-7.2e-010	41
-7.4e-010	52
-7.6e-010	46
-7.8e-010	41
-8e-010	30
-8.2e-010	25
-8.4e-010	27
-8.6e-010	20
-8.8e-010	30
-9e-010	25
-9.2e-010	25
-9.4e-010	32
-9.6e-010	27
-9.8e-010	19
-1e-009	26
-1.02e-009	21
-1.04e-009	14
-1.06e-009	17
-1.08e-009	20
-1.1e-009	21
-1.12e-009	27
-1.14e-009	19
-1.16e-009	14
-1.18e-009	15
-1.2e-009	15
-1.22e-009	12
-1.24e-009	8
-1.26e-009	17
-1.28e-009	12
-1.3e-009	13
-1.32e-009	19
-1.34e-009	18
-1.36e-009	11
-1.38e-009	16
-1.4e-009	7
-1.42e-009	13
-1.44e-009	9
-1.46e-009	10
-1.48e-009	11
-1.5e-009	10
-1.52e-009	10
-1.54e-009	13
-1.56e-009	8
-1.58e-009	11
-1.6e-009	8
-1.62e-009	11
-1.64e-009	9
-1.66e-009	8
-1.68e-009	6
-1.7e-009	10
-1.72e-009	9
-1.74e-009	5
-1.76e-009	6
-1.78e-009	10
-1.8e-009	9
-1.82e-009	5
-1.84e-009	6
-1.86e-009	10
-1.88e-009	9
-1.9e-009	6
-1.92e-009	4
-1.94e-009	6
-1.96e-009	4
-1.98e-009	7
-2e-009	9
-2.02e-009	8
-2.04e-009	6
-2.06e-009	3
-2.08e-009	11
-2.1e-009	6
-2.12e-009	6
-2.14e-009	6
-2.16e-009	14
-2.18e-009	8
-2.2e-009	2
-2.22e-009	8
-2.24e-009	6
-2.26e-009	2
-2.28e-009	8
-2.3e-009	6
-2.32e-009	3
-2.34e-009	10
-2.36e-009	3
-2.38e-009	5
-2.4e-009	3
-2.42e-009	6
-2.44e-009	3
-2.46e-009	3
-2.48e-009	2
-2.5e-009	5
-2.52e-009	6
-2.54e-009	4
-2.56e-009	1
-2.58e-009	4
-2.6e-009	2
-2.62e-009	4
-2.64e-009	7
-2.66e-009	5
-2.68e-009	6
-2.7e-009	2
-2.72e-009	4
-2.74e-009	2
-2.76e-009	1
-2.78e-009	2
-2.8e-009	4
-2.82e-009	2
-2.84e-009	3
-2.86e-009	3
-2.88e-009	5
-2.9e-009	5
-2.92e-009	3
-2.94e-009	6
-2.96e-009	1
-2.98e-009	4
-3e-009	2
-3.02e-009	3
-3.04e-009	4
-3.06e-009	1
-3.08e-009	6
-3.1e-009	6
-3.12e-009	1
-3.14e-009	5
-3.16e-009	3
-3.18e-009	3
-3.2e-009	5
-3.22e-009	3
-3.24e-009	2
-3.26e-009	6
-3.28e-009	3
-3.3e-009	3
-3.32e-009	2
-3.34e-009	2
-3.36e-009	3
-3.38e-009	1
-3.4e-009	3
-3.42e-009	4
-3.44e-009	1
-3.46e-009	3
-3.48e-009	2
-3.5e-009	2
-3.52e-009	1
-3.54e-009	5
-3.56e-009	3
-3.58e-009	1
-3.6e-009	0
-3.62e-009	4
-3.64e-009	2
-3.66e-009	0
-3.68e-009	3
-3.7e-009	2
-3.72e-009	1
-3.74e-009	2
-3.76e-009	0
-3.78e-009	2
-3.8e-009	2
-3.82e-009	0
-3.84e-009	3
-3.86e-009	4
-3.88e-009	0
-3.9e-009	2
-3.92e-009	3
-3.94e-009	0
-3.96e-009	4
-3.98e-009	0
-4e-009	2
-4.02e-009	1
-4.04e-009	1
-4.06e-009	1
-4.08e-009	1
-4.1e-009	1
-4.12e-009	1
-4.14e-009	1
-4.16e-009	1
-4.18e-009	3
-4.2e-009	1
-4.22e-009	1
-4.24e-009	4
-4.26e-009	2
-4.28e-009	1
-4.3e-009	1
-4.32e-009	5
-4.34e-009	3
-4.36e-009	3
-4.38e-009	1
-4.4e-009	1
-4.42e-009	2
-4.44e-009	3
-4.46e-009	1
-4.48e-009	3
-4.5e-009	0
-4.52e-009	0
-4.54e-009	2
-4.56e-009	0
-4.58e-009	1
-4.6e-009	0
-4.62e-009	3
-4.64e-009	1
-4.66e-009	1
-4.68e-009	2
-4.7e-009	3
-4.72e-009	1
-4.74e-009	1
-4.76e-009	2
-4.78e-009	0
-4.8e-009	1
-4.82e-009	3
-4.84e-009	2
-4.86e-009	2
-4.88e-009	1
-4.9e-009	0
-4.92e-009	0
-4.94e-009	0
-4.96e-009	2
-4.98e-009	0
-5e-009	1
-5.02e-009	1
-5.04e-009	1
-5.06e-009	0
-5.08e-009	2
-5.1e-009	0
-5.12e-009	0
-5.14e-009	0
-5.16e-009	2
-5.18e-009	0
-5.2e-009	1
-5.22e-009	2
-5.24e-009	0
-5.26e-009	0
-5.28e-009	1
-5.3e-009	2
-5.32e-009	0
-5.34e-009	0
-5.36e-009	0
-5.38e-009	0
-5.4e-009	2
-5.42e-009	0
-5.44e-009	1
-5.46e-009	0
-5.48e-009	0
-5.5e-009	0
-5.52e-009	1
-5.54e-009	0
-5.56e-009	1
-5.58e-009	3
-5.6e-009	1
-5.62e-009	1
-5.64e-009	2
-5.66e-009	0
-5.68e-009	0
-5.7e-009	0
-5.72e-009	0
-5.74e-009	0
-5.76e-009	0
-5.78e-009	0
-5.8e-009	4
-5.82e-009	0
-5.84e-009	0
-5.86e-009	0
-5.88e-009	1
-5.9e-009	2
-5.92e-009	3
-5.94e-009	0
-5.96e-009	0
-5.98e-009	0
-6e-009	0
-6.02e-009	1
-6.04e-009	0
-6.06e-009	0
-6.08e-009	1
-6.1e-009	0
-6.12e-009	0
-6.14e-009	0
-6.16e-009	0
-6.18e-009	0
-6.2e-009	0
-6.22e-009	0
-6.24e-009	1
-6.26e-009	1
-6.28e-009	0
-6.3e-009	1
-6.32e-009	0
-6.34e-009	0
-6.36e-009	1
-6.38e-009	0
-6.4e-009	0
-6.42e-009	1
-6.44e-009	1
-6.46e-009	0
-6.48e-009	0
-6.5e-009	0
-6.52e-009	0
-6.54e-009	1
-6.56e-009	1
-6.58e-009	0
-6.6e-009	2
-6.62e-009	0
-6.64e-009	0
-6.66e-009	0
-6.68e-009	0
-6.7e-009	0
-6.72e-009	0
-6.74e-009	0
-6.76e-009	0
-6.78e-009	0
-6.8e-009	1
-6.82e-009	0
-6.84e-009	0
-6.86e-009	0
-6.88e-009	0
-6.9e-009	0
-6.92e-009	1
-6.94e-009	1
-6.96e-009	0
-6.98e-009	0
-7e-009	1
-7.02e-009	1
-7.04e-009	1
-7.06e-009	0
-7.08e-009	1
-7.1e-009	0
-7.12e-009	0
-7.14e-009	1
-7.16e-009	0
-7.18e-009	1
-7.2e-009	1
-7.22e-009	0
-7.24e-009	0
-7.26e-009	0
-7.28e-009	0
-7.3e-009	1
-7.32e-009	0
-7.34e-009	1
-7.36e-009	0
-7.38e-009	0
-7.4e-009	0
-7.42e-009	0
-7.44e-009	1
-7.46e-009	0
-7.48e-009	0
-7.5e-009	0
-7.52e-009	0
-7.54e-009	0
-7.56e-009	0
-7.58e-009	1
-7.6e-009	0
-7.62e-009	0
-7.64e-009	0
-7.66e-009	1
-7.68e-009	0
-7.7e-009	0
-7.72e-009	0
-7.74e-009	0
-7.76e-009	0
-7.78e-009	0
-7.8e-009	2
-7.82e-009	0
-7.84e-009	0
-7.86e-009	1
-7.88e-009	0
-7.9e-009	1
-7.92e-009	1
-7.94e-009	1
-7.96e-009	0
-7.98e-009	0
-8e-009	1
-8.02e-009	0
-8.04e-009	1
-8.06e-009	0
-8.08e-009	0
-8.1e-009	0
-8.12e-009	0
-8.14e-009	0
-8.16e-009	0
-8.18e-009	0
-8.2e-009	0
-8.22e-009	0
-8.24e-009	1
-8.26e-009	0
-8.28e-009	0
-8.3e-009	0
-8.32e-009	0
-8.34e-009	0
-8.36e-009	0
-8.38e-009	0
-8.4e-009	0
-8.42e-009	2
-8.44e-009	0
-8.46e-009	0
-8.48e-009	0
-8.5e-009	0
-8.52e-009	0
-8.54e-009	0
-8.56e-009	0
-8.58e-009	0
-8.6e-009	0
-8.62e-009	0
-8.64e-009	0
-8.66e-009	0
-8.68e-009	0
-8.7e-009	0
-8.72e-009	0
-8.74e-009	0
-8.76e-009	0
-8.78e-009	0
-8.8e-009	0
-8.82e-009	0
-8.84e-009	0
-8.86e-009	0
-8.88e-009	0
-8.9e-009	0
-8.92e-009	1
-8.94e-009	0
-8.96e-009	0
-8.98e-009	0
-9e-009	0
-9.02e-009	0
-9.04e-009	0
-9.06e-009	0
-9.08e-009	0
-9.1e-009	0
-9.12e-009	0
-9.14e-009	0
-9.16e-009	0
-9.18e-009	0
-9.2e-009	0
-9.22e-009	0
-9.24e-009	1
-9.26e-009	0
-9.28e-009	0
-9.3e-009	0
-9.32e-009	0
-9.34e-009	0
-9.36e-009	0
-9.38e-009	0
-9.4e-009	0
-9.42e-009	1
-9.44e-009	0
-9.46e-009	0
-9.48e-009	0
-9.5e-009	0
-9.52e-009	1
-9.54e-009	0
-9.56e-009	0
-9.58e-009	0
-9.6e-009	0
-9.62e-009	0
-9.64e-009	0