diff --git a/_PeakOTronAP2.py b/_PeakOTronAP2.py
deleted file mode 100644
index dd2c6601caf8df7acb06ca5f8ae4c121b792d8a7..0000000000000000000000000000000000000000
--- a/_PeakOTronAP2.py
+++ /dev/null
@@ -1,2388 +0,0 @@
-## code starts###
-
-import os
-import numpy as np
-import pandas as pd
-from itertools import chain, zip_longest
-from iminuit import Minuit
-from iminuit.util import describe
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-from matplotlib import cm
-from numba import njit
-from difflib import ndiff
-import re
-from KDEpy import FFTKDE
-import scipy.special as sc
-from scipy.integrate import quad
-from scipy.stats import poisson, skew, norm, binom
-from scipy.interpolate import interp1d
-from scipy.signal import find_peaks, peak_widths
-from scipy.signal import fftconvolve as convolve
-from astropy.stats import knuth_bin_width, freedman_bin_width, scott_bin_width, bootstrap
-from sklearn.preprocessing import RobustScaler
-from AdditionalPDFs import gpoisson, borel
-
-
-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:]))
-
-
-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):
-        try:
-            x_kde, y_kde = FFTKDE(kernel = self.kernel, bw=bw).fit(self.data).evaluate(self.n_kde_samples)
-        
-            loss = -np.sum(np.log(y_kde))
-            
-        except:
-            
-            loss = self.eps_inv
-            
-        return loss
-        
-    def _PPF(self, data):
-        """ Compute ECDF """
-        x = np.sort(data)
-        n = x.size
-        y = np.arange(1, n+1) / n
-        _ppf = interp1d(y, x)
-        return _ppf
-
-        
-
-    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
-        self.PPF = self._PPF(self.data)
-        self.data = self.data[(self.data<self.PPF(alpha))]
-        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
-
-
-               
-                                              
-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)
-        self.log_eps_inv = np.log(self.eps_inv)
-        self.loss = []
-        
-        
-    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)
-        logL = self.y*(self.Logify(y_hat) - self.Logify(self.y)) + (self.y - y_hat)
-        nlogL = -np.sum(logL)
-            
-        self.loss.append(nlogL)
-        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.35):
-        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)
-        
-
-
-@njit
-def _SelectRangeNumba(array, low_lim, hi_lim):
-    index = (array >=  low_lim) & (array <= hi_lim)
-    return array[index]
-        
-
-class PeakOTron:
-    
-    def __init__(self,
-                 verbose=False,
-                 n_call_minuit=2000,
-                 n_iterations_minuit = 20
-                ):
-        
-        
-        self.Init()
-        self.greek_letters = ["Alpha",
-                      "alpha",
-                      "Beta",
-                      "beta",
-                      "Gamma",
-                      "gamma",
-                      "Delta",
-                      "delta",
-                      "Epsilon",
-                      "epsilon",
-                      "Zeta",
-                      "zeta",
-                      "Eta",
-                      "eta",
-                      "Theta",
-                      "theta",
-                      "Iota",
-                      "iota",
-                      "Kappa",
-                      "kappa",
-                      "Lbda",
-                      "lbda",
-                      "Mu",
-                      "mu",
-                      "Nu",
-                      "nu",
-                      "Xi",
-                      "xi",
-                      "Omicron",
-                      "omicron",
-                      "Pi",
-                      "pi",
-                      "Rho",
-                      "rho",
-                      "Sigma",
-                      "sigma",
-                      "Tau",
-                      "tau",
-                      "Upsilon",
-                      "upsilon",
-                      "Phi",
-                      "phi",
-                      "Chi",
-                      "chi",
-                      "Psi",
-                      "psi",
-                      "Omega",
-                      "omega"]
-        
-        self._eps = np.finfo(np.float64).eps * 10
-        self._eps_kde = 1e-5
-        self._FWHM2Sigma = 1/(2*np.sqrt(2*np.log(2)))
-        
-        self._plot_figsize= (10,10)
-        self._plot_fontsize= 20
-        self._cmap = cm.get_cmap('viridis')
-        
-        self._max_n_peaks = 12
-        self._max_n_dcr_peaks = 6
-        self._len_DCR_pad = int(100)
-        
-        self._verbose=verbose
-        self._n_call_minuit = n_call_minuit
-        self._n_iterations_minuit = n_iterations_minuit
-        self._is_histogram = False
-        
-  
-        self._default_fit_kwargs={
-                "tau":None,
-                "t_0":None,
-                "r_fast":None,
-                "tgate":None,
-                "pedestal":None,
-                "truncate":True,
-                "n_bootstrap": 1000,
-                "n_bootstrap_kde":100,
-                "bw":None,
-                "peak_dist_factor":0.8,
-                "peak_width_factor":0.5,
-                "peak_ppf":0.95,
-                "peak_prominence": 0.0,
-                "peak_min_n_evts": 100,
-                "kernel_kde":"gaussian",
-                "bandwidth_kde":"ISJ",
-                "ppf_mainfit":1 - 1e-6,
-                "bin_method":"knuth"
-        }
-        
-
-    
-
-        
-    def LevDist(self, str1, str2):
-        counter = {"+": 0, "-": 0}
-        distance = 0
-        for edit_code, *_ in ndiff(str1, str2):
-            if edit_code == " ":
-                distance += max(counter.values())
-                counter = {"+": 0, "-": 0}
-            else:
-                counter[edit_code] += 1
-        distance += max(counter.values())
-        return distance
-    
-
-    ###HELPER FUNCTIONS
-        
-    def Init(self):
-        
-        self._peak_data = {}
-        self._GD_data = {}
-        self._DCR_data={}
-        
-        
-        self._fit_values= {}
-        self._fit_errors = {}
-        self._fit = None
-        self._failed = False
-        
-        
-        
-    def LatexFormat(self, f, scirange=[0.01,1000]):
-        if(np.abs(f)>scirange[0] and np.abs(f)<scirange[1]):
-            float_str = r"${:3.3f}$".format(f)
-        else:
-            try:
-                float_str = "{:3.3E}".format(f)
-                base, exponent = float_str.split("E")
-                float_str = r"${0} \times 10^{{{1}}}$".format(base, int(exponent))
-            except:
-                float_str=str(f)
-        return float_str
-
-
-        
-
-    def SelectRange(self, array, low_lim, hi_lim):
-        return _SelectRangeNumba(array, low_lim, hi_lim)
-    
-    
-    
-    def Logify(self, 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 SetMaxNPeaks(self, max_n_peaks):
-        max_n_peaks = int(max_n_peaks)
-        if(max_n_peaks<1):
-            raise Exception("Maximum number of peaks must be greater than 2.")
-        
-
-        self._max_n_peaks =  max_n_peaks
-        print("Set maximum number of peaks to {:d}.".format(self._max_n_peaks))
-        
-
-                
-        
-    def SetMaxNDCRPeaks(self, max_n_dcr_peaks):
-        max_n = int( max_n_dcr_peaks)
-        if(max_n_dcr_peaks<2):
-            raise Exception("Maximum number of peaks must be greater than 2.")
-        
-
-        self._max_n_dcr_peaks =  max_n_dcr_peaks
-        print("Set maximum number of dark peaks to {:d}.".format(self._max_n_dcr_peaks))
-
-
-    def DummyLinear(self, x, m, c):
-        return m*x + c
-
-    def DummyLogGPois(self, k, mu, lbda, A):
-        
-        
-        return A + gpoisson.logpmf(k, mu, lbda)
-
-    def SigmaK(self, k, sigma_0, sigma_1):
-            return np.sqrt(sigma_0**2 + k*sigma_1**2)
-    
-
-    def PoisPPF(self, q, mu):
-        vals = np.ceil(sc.pdtrik(q, mu))
-        vals1 = np.maximum(vals - 1, 0)
-        temp = sc.pdtr(vals1, mu)
-        return np.where(temp >= q, vals1, vals)
-    
-    def GPoisPPF(self, q, mu, lbda, k_max):
-        _sum = 0
-        _count = 0
-        while _count<=k_max:
-            _sum +=  gpoisson.pmf(_count, mu, lbda)
-            if(_sum<q):
-                _count+=1
-            else:
-                break
-        
-        return _count
-    
-
-
-    def EstLbda(self, data_sub_x_0):
-        _skew = skew(data_sub_x_0)
-        _mean = np.mean(data_sub_x_0)
-        _std = np.std(data_sub_x_0)
-
-        _lbda = 0.5*(((_mean*_skew)/(_std))- 1)
-        return _lbda
-
-    def EstGain(self, data_sub_x_0):
-        _skew = skew(data_sub_x_0)
-        _mean = np.mean(data_sub_x_0)
-        _std = np.std(data_sub_x_0)
-
-        _lbda = 0.5*(((_mean*_skew)/(_std))- 1)
-        _gain = (_std**2/_mean)*((1 - _lbda)**2)
-        return _gain
-
-    def EstMu(self, data_sub_x_0):
-        _skew = skew(data_sub_x_0)
-        _mean = np.mean(data_sub_x_0)
-        _std = np.std(data_sub_x_0)
-        _lbda = 0.5*(((_mean*_skew)/(_std))- 1)
-        _mu = (1/(1-_lbda))*(_mean**2/_std**2)
-        return _mu
-    
-  
-
-    def Bootstrap(self, data, statistic, alpha=0.95, n_samples=10000):
-        if not (0 < alpha < 1):
-            raise ValueError("confidence level must be in (0, 1)")
-
-        n = n_samples
-        if n <= 0:
-            raise ValueError("data must contain at least one measurement.")
-
-
-
-        boot_stat = bootstrap(data, n_samples, 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(self, data, kernel = "gaussian", bw="optimise",
-                     n_kde_samples=2**13,
-                     alpha=0.95,
-                     n_samples=1000,
-                     limits=None
-                    ):
-
-        if not (0 < alpha < 1):
-            raise ValueError("Bootstrap confidence level, alpha, must be in (0, 1).")
-
-        n = n_samples
-        if n <= 0:
-            raise ValueError("Bootstrap data must contain at least one measurement.")
-
-        kernels = list(FFTKDE._available_kernels.keys())
-        if kernel not in kernels:
-            raise ValueError("Selected KDE kernel not available. Please choose from these options: {:s}.".format(", ".join(kernels)))
-
-        
-        if(not isinstance(bw, float)):
-            bws = list(FFTKDE._bw_methods.keys()) + ["optimise"]
-            if bw not in bws:
-                raise ValueError("Selected KDE bandwidth estimation method not available. Please choose from these options: {:s}.".format(", ".join(bws)))
-
-        
-        
-        if(bw=="optimise"):
-            
-            try:
-                BWO= BandWidthOptimiser(data, kernel=kernel)
-                bw_scott = scott_bin_width(data)
-                minuit_kde = Minuit(BWO, bw=bw_scott/2)
-                minuit_kde.limits["bw"]=(self._eps_kde, 3*bw_scott)
-                minuit_kde.strategy=2
-                minuit_kde.migrad(ncall= self._n_call_minuit,
-                           iterate =self._n_iterations_minuit)
-
-                bw = minuit_kde.values["bw"]
-                if(self._verbose):
-                    print("KDE Bandwidth optimised to: {0}".format(bw))
-            except:
-                if(self._verbose):
-                    print("KDE 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_samples)
-        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 DRM_basic(self, 
-                  x, 
-                  mu, 
-                  x_0, 
-                  G, 
-                  lbda, 
-                  sigma_0, 
-                  sigma_1, 
-                  k_low, 
-                  k_hi):
-        
-        k = np.arange(k_low, k_hi)
-        k = sorted(list(k[k!=0]))
-        k_is = [k[0:int(_k)] for _k in k]
-
-        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,  self.SigmaK(_k, sigma_0, sigma_1))
-                for _k in k
-        ])
-
-
-        f = np.sum(np.vstack([f0, f1s]), axis=0)
-
-        return f
-    
-    
-    def dPApdPH(self, 
-               x, 
-               k, 
-               x_0, 
-               G,
-               tau,
-               r_Ap, 
-               tgate,
-               pAp):
-    
-        PH = (x - x_0)/G - k
-        tauAp = r_Ap*tau
-
-        PH_max = 1 + np.exp(-tgate/tau) - 2*np.exp(-0.5*tgate/tau)
-
-        cond_PH_max = (PH<PH_max)
-
-        PH_temp = np.where(cond_PH_max, PH, PH_max)
-
-        t_b0 = 0.5*tgate - tau*np.arccosh(0.5*((1 - PH_temp)*np.exp(tgate/(2*tau)) + np.exp(-tgate/(2*tau))))
-        t_b1 = tgate-t_b0
-
-        dpApdt_b0 = pAp*(1-np.exp(-t_b0/tau))*np.exp(-t_b0/tauAp)/tauAp
-        dpApdt_b1 = pAp*(1-np.exp(-t_b1/tau))*np.exp(-t_b1/tauAp)/tauAp
-
-        dPHdt_b0 = -2*np.sinh((t_b0-0.5*tgate)/tau)*np.exp(-0.5*tgate/tau)/tau
-        dPHdt_b1 = -2*np.sinh((t_b1-0.5*tgate)/tau)*np.exp(-0.5*tgate/tau)/tau
-
-        dpApdPH_b0 = dpApdt_b0/((np.where(np.abs(dPHdt_b0)>1e-5, np.abs(dPHdt_b0), 1e-5)))
-        dpApdPH_b1 = dpApdt_b1/((np.where(np.abs(dPHdt_b1)>1e-5, np.abs(dPHdt_b1), 1e-5)))
-
-
-
-        dpApdPH_b0 = np.where(cond_PH_max, dpApdPH_b0, 0)
-        dpApdPH_b1 = np.where(cond_PH_max, dpApdPH_b1, 0)
-
-
-        dpApdPH_b0 = np.where(dpApdPH_b0>0, dpApdPH_b0, 0)
-        dpApdPH_b1 = np.where(dpApdPH_b0>0, dpApdPH_b1, 0)
-
-
-        return PH_max*(dpApdPH_b0)/G
-
-
-    
-    def ApplyAp(self, 
-                x, 
-                y_light, 
-                mu, 
-                x_0, 
-                G, 
-                lbda, 
-                DCR, 
-                tgate, 
-                t_0, 
-                tau, 
-                r_Ap, 
-                pAp, 
-                k_low, 
-                k_hi, 
-                k_dcr_low, 
-                k_dcr_hi):
-        
-        k = np.arange(k_low, k_hi)
-        k = sorted(list(k[k!=0]))
-        
-        mu_dcr = DCR*(tgate+abs(t_0))
-        
-        k_dcr = np.arange(k_dcr_low, k_dcr_hi)
-        k_dcr = sorted(list(k_dcr[k_dcr!=0]))
-        
-        
-        ps_light = [gpoisson.pmf(_k, mu, lbda) for _k in k]
-        ps_dark = [gpoisson.pmf(_k, mu_dcr, lbda) for _k in k_dcr]
-        
-        f0s = np.asarray([_p_l + _p_d for _p_l, _p_d in zip_longest(ps_light, ps_dark, fillvalue=0)]).reshape(-1,1)
-        
-            
-        f1s = np.asarray([
-                _k*self.dPApdPH(x, _k, x_0, G, tau, r_Ap, tgate, pAp)
-                for _k in k
-        ])
-        
-        f = np.sum(f0s*f1s, axis=0)
-
-
-        return y_light+f
-
-
-
-  
-
-    def PH_range(self, t, t_0, r_fast, tau, tgate):
-
-        if((t>t_0) and (t<0)):
-            PH_min = (1-r_fast)*np.exp(t_0/tau)*(1 - np.exp(-tgate/tau))
-            PH_max = (1-r_fast)*(1 - np.exp(-tgate/tau))  
-        elif((t>0) and (t<tgate)):
-            PH_min = r_fast
-            PH_max = (1 - (1-r_fast)*np.exp(-tgate/tau))
-        else:
-            PH_min = 0
-            PH_max = 0
-
-        return PH_min, PH_max
-
-
-    def dpDCRdPH(self, x, x_0, G, tau, tgate, t_0, r_fast):
-
-
-        PH_bfg_low, PH_bfg_hi = self.PH_range(-abs(t_0)/2, -abs(t_0), r_fast, tau, tgate)
-        PH_dg_low, PH_dg_hi = self.PH_range(tgate/2, -abs(t_0), r_fast, tau, tgate)
-
-        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(self, 
-                 x, 
-                 y_light, 
-                 x_0, 
-                 G, 
-                 lbda, 
-                 DCR, 
-                 tgate, 
-                 t_0, 
-                 tau, 
-                 r_fast,  
-                 dx, 
-                 k_dcr_low, 
-                 k_dcr_hi):
-    
-        mu_dcr = DCR*(abs(t_0) + tgate)
-        
-        x_lin = np.arange(0, len(x))
-        x_pad_lin = np.arange(-self._len_DCR_pad, len(x)+self._len_DCR_pad)
-        f_lin = interp1d(x, x_lin, fill_value='extrapolate')
-        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))
-
-        
-        n_dcr = np.arange(k_dcr_low, max(2, k_dcr_hi) )
-
-        fs = []
-        pfs = []
-        hs = []
-        phs = []
-
-        for _n_dcr in n_dcr:
-
-            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(len(fs)<=1):
-                    f = self.dpDCRdPH(x_pad_lin, x_0_lin, G_lin, tau, tgate, t_0, r_fast)
-                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)
-
-                h = self.dpDCRdPH(x_pad_lin, x_0_lin, G_lin*(_n_dcr+1), tau, tgate, t_0, r_fast)
-                h = h/ np.trapz(h, dx = 1)
-                ph = poisson.pmf(1, mu_dcr)*borel.pmf(_n_dcr, lbda)/(_n_dcr+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_model = convolve(y_dark,
-                              np.pad(y_light,
-                                     (self._len_DCR_pad, self._len_DCR_pad),
-                                     "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(
-            self,
-            x,
-            mu, 
-            lbda,
-            x_0, 
-            G, 
-            sigma_0, 
-            sigma_1, 
-            pAp, 
-            k_low, 
-            k_hi,
-            DCR,
-            tau,
-            r_Ap,
-            tgate,
-            t_0,
-            r_fast,
-            dx,
-            k_dcr_low,
-            k_dcr_hi,
-           ):
-        
-        
-        
-        y_light = self.DRM_basic(x, 
-                                 mu, 
-                                 x_0, 
-                                 G, 
-                                 lbda, 
-                                 sigma_0, 
-                                 sigma_1, 
-                                 k_low, 
-                                 k_hi)
-        
-        y_lightap = self.ApplyAp(x, 
-                               y_light, 
-                               mu, 
-                               x_0, 
-                               G, 
-                               lbda, 
-                               DCR, 
-                               tgate, 
-                               t_0, 
-                               tau, 
-                               r_Ap, 
-                               pAp, 
-                               k_low, 
-                               k_hi, 
-                               k_dcr_low, 
-                               k_dcr_hi)
-
-
-
-        y_model = self.ApplyDCR(x, 
-                                 y_lightap, 
-                                 x_0, 
-                                 G, 
-                                 lbda, 
-                                 DCR, 
-                                 tgate, 
-                                 t_0, 
-                                 tau, 
-                                 r_fast,  
-                                 dx, 
-                                 k_dcr_low, 
-                                 k_dcr_hi)
-            
-
-
-            
-        return y_model
-    
-    
-
-    ###PLOTTING
-    
-
-    
-    def PlotOriginalHistogram(self, ax):
-        
-        ax.plot(self._GD_data["x_s"],
-                    self._GD_data["y_s"],
-                    lw=5,
-                    label="Data")
-        ax.set_title("Normalised \n Input \n Distribution",fontsize=self._plot_fontsize)
-        ax.set_xlabel("IQR Normalised Charge Units [Q.U.]", fontsize=self._plot_fontsize)
-        ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-        ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-        ax.set_yscale("log")
-        ax.grid(which="both")
-        ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
-        
-    
-    def PlotKDE(self, ax):
-        ax.plot(self._GD_data["x_s"],
-                    self._GD_data["y_s"],
-                    label="Data",
-                    lw=5
-                    )
-        ax.plot(self._GD_data["x_kde_s"],
-                self._GD_data["y_kde_s"],
-                lw=2.5,
-                label="KDE")
-
-        ax.set_title("Kernel Density Estimate", fontsize=self._plot_fontsize)
-        ax.set_xlabel("IQR Normalised Charge Units [Q.U.]", fontsize=self._plot_fontsize)
-        ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-        ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-        ax.set_yscale("log")
-        ax.grid(which="both")
-        ax.set_ylim([1e-5, None])
-        ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
-        
-      
-
-
-        
-        
-    def PlotPeaks(self, ax):
-    
-            ax.plot(self._GD_data["x_s"],
-                     self._GD_data["y_s"],
-                     label="Input Pulse \n Height \n Spectrum",
-                     lw=5
-                    )
-            
-            color_vals = [self._cmap(_v) for _v in np.linspace(0,0.75, len(self._peak_data["x_peak_s"]))]
-            
-            for _i_x, (_x, _x_w, _c_v) in enumerate(zip(
-                                self._peak_data["x_peak_s"],
-                                self._peak_data["x_width_s"],
-                                color_vals
-                            )):
-                _x_peak =_x
-                _x_width_low = _x - _x_w/2
-                _x_width_hi = _x + _x_w/2
-                _y_peak = self._GD_data["y_s"][np.argmin(abs(self._GD_data["x_s"]-_x_peak))]
-                _y_width_low = self._GD_data["y_s"][np.argmin(abs(self._GD_data["x_s"]-_x_width_low))]
-                _y_width_hi = self._GD_data["y_s"][np.argmin(abs(self._GD_data["x_s"]-_x_width_hi))]
-
-
-                if(_i_x == 0):
-                    annotation_text = "Pedestal"
-                    color = "red"
-                else:
-                    annotation_text = "Peak {:d}".format(_i_x)
-                    color = _c_v
-                    
-                ax.vlines(x=[_x],
-                          ymin=0,
-                          ymax=_y_peak,
-                          color=color,
-                          linestyle="-",
-                          lw=2.5,
-                          label = '\n{:s}: \n $x_{{pos}}$ = {:3.3f} Q.U \n $\Delta x_{{pos}}$ {:3.3f} Q.U \n'.format(
-                                                                                     annotation_text,
-                                                                                     _x_peak,
-                                                                                     _x_width_hi - _x_width_low)
-                        
-                         )
-                
-                ax.vlines(x=[_x_width_low],
-                          ymin=0,
-                          ymax=_y_width_low,
-                          color=color,
-                          lw=2.5,
-                          linestyle="--")
-                
-                ax.vlines(x=[_x_width_hi],
-                          ymin=0,
-                          ymax=_y_width_hi,
-                          color=color,
-                          lw=2.5,
-                          linestyle="--")
-
-        
-            if(self._GD_data["truncated"]):
-                ax.axvline(x=[self._GD_data["limit_s"]],
-                           color = "purple",
-                           lw=5,
-                           label="Pedestal Threshold ($PH_{{0}} - \\frac{G_{est}}{2}$)")
-
-            box = ax.get_position()
-            ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
-
-            # Put a legend to the right of the current axis
-            ax.legend(loc='center left',
-                      bbox_to_anchor=(1, 0.5),
-                      fontsize=max(10, self._plot_fontsize//2),
-                      ncol=len(self._peak_data["x_peak_s"])//6 + 1
-                     )
-            ax.set_title("Peak Finding \n # Peaks = {:d}".format(len(self._peak_data["x_peak_s"])) ,
-                         fontsize=self._plot_fontsize)
-
-            ax.set_xlabel("IQR Normalised Charge Units [Q.U.]", fontsize=self._plot_fontsize)
-            ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-            ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-            ax.set_yscale("log")
-            ax.grid(which="both")
-            
-            
-    def PlotGenPois(self, ax):
-    
-            ax.scatter(self._peak_data["n_peak_s"],
-                        self.Logify(self._peak_data["y_peak_s"]),
-                        label="Extracted Peak Heights"
-                       )
-
-
-            y_lgp = self.DummyLogGPois(self._peak_data["n_peak_s"],
-                                        self._GD_data["fit_GP_mu"],
-                                        self._GD_data["fit_GP_lbda"],
-                                        self._GD_data["fit_GP_A"]
-                                      )
-            
-            ax.plot(self._peak_data["n_peak_s"],
-                     y_lgp,
-                     linestyle="--",
-                     lw=2.5,
-                     color="red",
-                     label="Fit")
-            
-            ax.set_title(
-                "Generalised Poisson Fit \n $\\mu$ = {0} $\pm$ {1} \n $\\lambda$ = {2} $\pm$ {3} \n $A$ = {4} $\pm$ {5}".format(
-                self.LatexFormat(self._GD_data["fit_GP_mu"]),
-                self.LatexFormat(self._GD_data["fit_GP_mu_err"]),
-                self.LatexFormat(self._GD_data["fit_GP_lbda"]),
-                self.LatexFormat(self._GD_data["fit_GP_lbda_err"]),
-                self.LatexFormat(self._GD_data["fit_GP_A"]),
-                self.LatexFormat(self._GD_data["fit_GP_A_err"])
-                
-            ), fontsize=self._plot_fontsize)
-
-
-            ax.set_xlabel("Peak", fontsize=self._plot_fontsize)
-            ax.set_ylabel("$\\log{\\frac{dP}{dPH}}$", fontsize=self._plot_fontsize)
-            ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-            ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-            ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
-            ax.grid(which="both")
-            
-            
-            
-    def PlotSigmaEst(self, ax):
-        
-        c = self._peak_data["cond_valid_peak"]
-        
-        ax.set_title("Estimated Peak Variance \n $\\sigma^{{2}}_{{0}}$={0} $\pm$ {1} \n $\\sigma^{{2}}_{{1}}$={2} $\pm$ {3}".format(
-                                                                                        self.LatexFormat(self._GD_data["fit_var0"]),
-                                                                                        self.LatexFormat(self._GD_data["fit_var0_err"]),
-                                                                                        self.LatexFormat(self._GD_data["fit_var1"]),
-                                                                                        self.LatexFormat(self._GD_data["fit_var1_err"]),
-                                                                                       ), fontsize=self._plot_fontsize)
-
-
-        ax.errorbar(self._peak_data["n_peak_s"][c],
-         self._peak_data["peak_var_s"][c],
-         yerr=self._peak_data["peak_var_err_s"][c],
-         label="Extracted Variances",
-         fmt="o")
-
-
-        ax.plot(self._peak_data["n_peak_s"][c],
-                self.DummyLinear(self._peak_data["n_peak_s"][c],
-                                 self._GD_data["fit_var1"],
-                                 self._GD_data["fit_var0"]
-                                 ),
-                 linestyle="--",
-                 lw=2.5,
-                 color="red",
-                 label="Fit")
-
-        ax.set_xlabel("Peak", fontsize=self._plot_fontsize)
-        ax.set_ylabel("$\sigma^{2}_{peak}$", fontsize=self._plot_fontsize)
-        ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-        ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-        ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
-        ax.grid(which="both")
-        
-        
-    
-    def PlotGainEst(self, ax):
-
-        c = self._peak_data["cond_valid_peak"]
-        
-        ax.set_title("Estimated Gain \n $G$={0} $\pm$ {1} \n $PH_{{0}}$={2} $\pm$ {3}".format(
-                                                                                        self.LatexFormat(self._GD_data["fit_G"]),
-                                                                                        self.LatexFormat(self._GD_data["fit_G_err"]),
-                                                                                        self.LatexFormat(self._GD_data["fit_x_0"]),
-                                                                                        self.LatexFormat(self._GD_data["fit_x_0_err"]),
-                                                                                       ), fontsize=self._plot_fontsize)
-
-
-        ax.errorbar(self._peak_data["n_peak_s"][c],
-         self._peak_data["peak_mean_s"][c],
-         yerr=self._peak_data["peak_mean_err_s"][c],
-         label="Extracted Peak Means",
-         fmt="o")
-
-
-        ax.plot(self._peak_data["n_peak_s"][c],
-                self.DummyLinear(self._peak_data["n_peak_s"][c],
-                                 self._GD_data["fit_G"],
-                                 self._GD_data["fit_x_0"]
-                                 ),
-                 linestyle="--",
-                 lw=2.5,
-                 color="red",
-                 label="Fit")
-
-        ax.set_xlabel("Peak", fontsize=self._plot_fontsize)
-        ax.set_ylabel("$<Peak>$", fontsize=self._plot_fontsize)
-        ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-        ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-        ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
-        ax.grid(which="both")
-        
-            
-        
-        
-
-            
-    def PlotDCREst(self, ax):
-        
-        
-        
-        x_dc_kde = self._DCR_data["x_dc_kde"]
-        y_dc_kde = self._DCR_data["y_dc_kde"]
-        y_dc_kde_err = self._DCR_data["y_dc_kde_err"]
-        
-        cond_p_0toLM  = self._DCR_data["cond_p_0toLM"]
-        cond_p_LMto2LM = self._DCR_data["cond_p_LMtoL2M"]
-        cond_interpeak = self._DCR_data["cond_interpeak"]
-        
-        
-        int_p_0toLM = self._DCR_data["int_p_0toLM"]
-        dint_p_0toLM = self._DCR_data["dint_p_0toLM"]
-        int_p_LMto2LM = self._DCR_data["int_p_LMto2LM"]
-        dint_p_LMto2LM = self._DCR_data["dint_p_LMto2LM"]
-        mean_p_interpeak = self._DCR_data["mean_p_interpeak"]
-        dmean_p_interpeak = self._DCR_data["dmean_p_interpeak"]
-        
-  
-        ax.fill_between(
-                     x_dc_kde[cond_p_0toLM],
-                     y_dc_kde[cond_p_0toLM],
-                     color='none',
-                     hatch="///",
-                     edgecolor="b",
-                     label = r"$I_{{0 \rightarrow x_{{min}}}}$ = {:s} $\pm$ {:s}".format(self.LatexFormat(int_p_0toLM),
-                                                                                   self.LatexFormat(dint_p_0toLM)                                                                                  )
-                    )
-        
-        ax.fill_between(
-                     x_dc_kde[cond_p_LMto2LM],
-                     y_dc_kde[cond_p_LMto2LM],
-                     color='none',
-                     hatch="///",
-                     edgecolor="r",
-                     label = r"$I_{{x_{{min}} \rightarrow 2x_{{min}}}}$ = {:s} $\pm$ {:s}".format(self.LatexFormat(int_p_LMto2LM),
-                                                                                   self.LatexFormat(dint_p_LMto2LM)                                                                                  )
-                    )
-
-
-        ax.fill_between(
-                     x_dc_kde[cond_interpeak],
-                     y_dc_kde[cond_interpeak],
-                     color='none',
-                     hatch="xxx",
-                     edgecolor="g",
-                     label = r"$<I_{{x_{{min}} - x_{{lim}} \rightarrow x_{{min}} + x_{{lim}} }}>$ = {:s} $\pm$ {:s}".format(
-                                                                                  self.LatexFormat(mean_p_interpeak),
-                                                                                  self.LatexFormat(dmean_p_interpeak)                                                                                  )
-                    )
-        ax.plot(x_dc_kde,
-                y_dc_kde,
-                linestyle="-",
-                color="purple",
-                label="KDE, 95% Confidence"
-               )
-        
-        ax.fill_between(x_dc_kde,
-                         y1=y_dc_kde - 1.96*y_dc_kde_err,
-                         y2=y_dc_kde + 1.96*y_dc_kde_err,
-                         alpha=0.2,
-                         color="purple",
-                         label="KDE, 95% Confidence"
-                       )
-        
-        
-        ax.set_title("Dark Count Rate \n DCR = {0} $\pm$ {1}".format(
-                      self.LatexFormat(self._DCR_data["fit_DCR"]),
-                      self.LatexFormat(self._DCR_data["fit_DCR_err"])),
-                      fontsize=self._plot_fontsize)
-        
-        ax.set_xlabel("$PH_{norm, est}$", fontsize=self._plot_fontsize)
-        ax.grid(which="both")
-        ax.set_ylabel("$\\frac{dp_{DCR}}{PH_{norm, est}}$", fontsize=self._plot_fontsize)
-        
-        ax.tick_params(axis="x", labelsize=self._plot_fontsize)
-        ax.tick_params(axis="y", labelsize=self._plot_fontsize)
-        ax.legend(fontsize=max(10, self._plot_fontsize//1.5))
-        ax.set_yscale("log")
-        
-        
-
-    
-      
-    def PlotFit(self,
-                figsize=(10.5,10.5),
-                labelsize=None,
-                ticksize=None,
-                titlesize=None,
-                legendsize=None,
-                title=None,
-                xlabel = None,
-                scaled=False,
-                save_directory=None,
-                y_limits = [None, None],
-                residual_limits = [-5, 5],
-                linewidth_main = 5,  
-                x_limits = [None, None],
-                display=True
-               ):
-    
-    
-        if(labelsize is None):
-            labelsize=self._plot_fontsize
-        
-        if(ticksize is None):
-            ticksize=self._plot_fontsize
-            
-        if(titlesize is None):
-            titlesize = self._plot_fontsize
-            
-        if(legendsize is None):
-            legendsize = max(10, self._plot_fontsize//1.5)
-            
-
-        if(scaled):
-            x = self._GD_data["x_s"]
-            y = self._GD_data["y_s"]
-            y_err = self._GD_data["y_err_s"]
-            y_hat = self.GetModelIQR(self._GD_data["x_s"])
-            if(xlabel is None):
-                xlabel = "IQR Normalised Charge Units [Q.U.]"
-            
-        else:
-            x = self._GD_data["x"]
-            y = self._GD_data["y"]
-            y_err = self._GD_data["y_err"]
-            y_hat = self.GetModel(self._GD_data["x"])
-            if(xlabel is None):
-                xlabel = "Input Charge Units [Q.U.]"
-            
-        
-    
-        fig = plt.figure(figsize=figsize)
-        gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
-        ax0 = plt.subplot(gs[0])    
-        ax0.plot(x, y, label="Data", lw=5, color="C0")
-        ax0.plot(x, y_hat, linestyle="--", label="Model", lw=5, color="C1")
-        ax0.set_yscale("log")
-        ax0.set_xlabel("$\\frac{dp}{dPH}$", fontsize=labelsize)
-        ax0.tick_params(axis="x", labelsize=ticksize)
-        ax0.tick_params(axis="y", labelsize=ticksize)
-        ax0.set_ylim(y_limits)
-        ax0.set_xlim(x_limits)
-        if(title is not None):
-            ax0.set_title(title, fontsize=titlesize)
-            
-        ax0.grid(which="both")
-        ax0.legend(fontsize=legendsize)
-
-
-        ax1 = plt.subplot(gs[1], sharex = ax0)
-        ax1.scatter(x, (y - y_hat)/(y_err), color="C1")    
-
-
-        ax1.tick_params(axis="x", labelsize=ticksize)
-        ax1.tick_params(axis="y", labelsize=ticksize)
-        ax1.set_xlabel(xlabel, fontsize=labelsize)
-        ax1.set_ylim(residual_limits)
-        ax1.axhline(y=-1.0, color="purple", lw=2.5, linestyle="--")
-        ax1.axhline(y=1.0, color="purple", lw=2.5, linestyle="--")
-        ax1.tick_params(axis="x", labelsize=self._plot_fontsize)
-        ax1.tick_params(axis="y", labelsize=self._plot_fontsize)
-        ax1.grid(which="both")
-
-        fig.tight_layout()
-
-        fig.subplots_adjust(hspace=.0)
-        if(save_directory is not None):
-            print("Saving figure to {0}...".format(save_directory))
-            fig.savefig(save_directory)
-        if(display):
-            plt.pause(0.01)
-            fig.show()
-        else:
-            plt.close(fig)
-        
-    
-  
-    def PlotSummary(self, display=True, save_directory=None):
-    
-            fig = plt.figure(figsize=(20,40))
-            gs = gridspec.GridSpec(4, 2)
-            gs.update(wspace=0.5, hspace=0.5)
-              
-            ax0 = fig.add_subplot(gs[0,0])
-            self.PlotOriginalHistogram(ax0)
-            
-            ax1 = fig.add_subplot(gs[0,1])
-            self.PlotKDE(ax1)
-            
-            ax2 = fig.add_subplot(gs[1,:])
-            self.PlotPeaks(ax2)
-                        
-            ax3 = fig.add_subplot(gs[2,0])
-            self.PlotGenPois(ax3)
-            
-            ax4 = fig.add_subplot(gs[2,1])
-            self.PlotSigmaEst(ax4)
-            
-            ax5 = fig.add_subplot(gs[3,0])
-            self.PlotGainEst(ax5)
-          
-            ax6 = fig.add_subplot(gs[3,1])
-            self.PlotDCREst(ax6)
-
-            if(save_directory is not None):
-                print("Saving figure to {0}...".format(save_directory))
-                fig.savefig(save_directory)
-            if(display):
-                plt.pause(0.01)
-                fig.show()
-            else:
-                plt.close(fig)
- 
-            
-    def GetModel(self, x):
-        return self.DRM(x, **self._fit_values)
-        
-    def GetModelIQR(self, x):
-        return self.DRM(x, **self._fit_values_s)
-        
-        
-    def GetOriginalDistribution(self):
-        return _self._GD_data["x"], self._GD_data["y"]
-  
-    def GetFitResults(self):
-        return self._fit_values, self._fit_errors
-    
-    def Failed(self):
-        return self._failed
-    
-    
-    ## PREFIT
-    
-            
-    
-    def GetHistAndPeaks(self,
-                        data,
-                        bw,
-                        peak_dist_factor,
-                        peak_width_factor,
-                        peak_ppf,
-                        peak_prominence,
-                        peak_min_n_evts,
-                        truncate,
-                        pedestal,
-                        bin_method,
-                        n_bootstrap,
-                        n_bootstrap_kde,
-                        bandwidth_kde,
-                        kernel_kde):
-        
- 
-        self.scaler = RobustScaler()
-        data_s = np.squeeze(self.scaler.fit_transform(data.reshape(-1,1)))
-        pedestal_s = (pedestal - self.scaler.center_[0])/self.scaler.scale_[0]
-        
-        if(bw is None):
-            if(bin_method == "knuth"):
-                bw_s, bins_s = knuth_bin_width(data_s, return_bins=True)
-            elif(bin_method=="freedman"):
-                bw_s, bins_s = freedman_bin_width(data_s, return_bins=True)
-            elif(bin_method=="scott"):
-                bw_s, bins_s = scott_bin_width(data_s, return_bins=True)
-            else:
-                raise Exception ("Please specify a bin width by setting 'bw' or choose a method from 'knuth', 'freedman' or 'scott'")
-    
-            nbins_s = len(bins_s)
-        
-            bw = bw_s*self.scaler.scale_[0]
-            nbins = np.ceil((data.max() - data.min()) / bw)
-            nbins = max(1, nbins)
-            bins = data.min() + bw * np.arange(nbins + 1)
-        
-
-        else:
-            
-            bw_s = bw/self.scaler.scale_[0]
-            
-            nbins = np.ceil((data.max() - data.min()) / bw)
-            nbins = max(1, nbins)
-            bins = data.min() + bw * np.arange(nbins + 1)
-                        
-            nbins_s = np.ceil((data_s.max() - data_s.min()) / bw_s)
-            nbins_s = max(1, nbins_s)
-            bins_s = data_s.min() + bw_s * np.arange(nbins_s + 1)
-            
-    
-    
-                       
-        if(nbins<50):
-            print("Binning with bw = {:3.3E} produced only {:d} bins, less than limit of 50 bins. Setting fit fail status and continuing...".format(bw, nbins))
-            self._failed = True
-            return
-            
-        #Calculate histogram information
-            
-        y_N, x = np.histogram(data, bins=bins, density=False)
-        
-        N = np.sum(y_N)
-        
-        y = y_N/N/bw
-        
-        y_err = y*np.sqrt(1/y_N + 1/N)
-        
-        x = (x[:-1] + x[1:])/2.
-        
-
-            
-        
-        y_N_s, x_s = np.histogram(data_s, bins=bins_s, density=False)
-        
-        N_s = np.sum(y_N_s)
-        
-        y_s = y_N_s/N_s/bw_s
-        
-        y_err_s = y_s*np.sqrt(1/y_N_s + 1/N_s)
-        
-        x_s = (x_s[:-1] + x_s[1:])/2.
-
-        
-        
-        
-        try:
-            
-            x_kde_s, y_kde_s, y_kde_err_s, _, = self.BootstrapKDE(data_s,
-                                                       kernel=kernel_kde,
-                                                       bw=bandwidth_kde,
-                                                       n_samples=n_bootstrap_kde)
-
-            f_lin_kde_s = interp1d(x_kde_s, np.arange(len(x_kde_s)), fill_value='extrapolate')
-
-
-            
-        except:
-            print("Initial KDE failed. Setting fit fail status and continuing...")
-            self._failed = True
-            return
-          
-          
-
-        
-        est_mu, std_err_mu, _ = self.Bootstrap(data_s-pedestal_s,
-                                       self.EstMu,
-                                       n_samples=n_bootstrap)
-        
-        est_lbda, std_err_lbda, _ = self.Bootstrap(data_s - pedestal_s,
-                                   self.EstLbda,
-                                   n_samples=n_bootstrap)
-        
-        
-        est_gain, std_err_gain, _ = self.Bootstrap(data_s - pedestal_s,
-                                           self.EstGain,
-                                           n_samples=n_bootstrap)
-        
-        
-
-        
-        est_gain_lin = abs(f_lin_kde_s(est_gain) - f_lin_kde_s(0))
-        
-        
-        limit_s = pedestal_s - 0.5*est_gain
-        limit = limit_s*self.scaler.scale_[0] + self.scaler.center_[0]
-        
-        if(truncate):
-            
-            cond_limit = (x>limit)
-            cond_limit_s = (x_s>limit_s)
-            cond_limit_kde_s = (x_kde_s>limit_s)
-        
-            x = x[cond_limit]
-            y = y[cond_limit]
-            y_err = y_err[cond_limit]
-            
-            x_s = x_s[cond_limit_s]
-            y_s = y_s[cond_limit_s]
-            y_err_s = y_err_s[cond_limit_s]
-            
-            x_kde_s = x_kde_s[cond_limit_kde_s]
-            y_kde_s = y_kde_s[cond_limit_kde_s]
-            y_kde_err_s =  y_kde_err_s[cond_limit_kde_s]
-            
-            if(self._verbose):
-                print("Truncated distributions to only fit PH > pedestal - statistical gain/2 ...")
-                
-            
-        
-        
-        peaks, properties = find_peaks(y_kde_s,
-                                       distance=peak_dist_factor*est_gain_lin,
-                                       prominence=peak_prominence
-                                      )
-        
-    
-       
-        
-        if(len(peaks)<3):
-            print("Not enough peaks found with current distance criterion to fit. Setting fit fail status and continuing..")
-            self._failed=True
-            return
-        
-                  
-  
-        
-        x_peaks_s = x_kde_s[peaks]
-        y_peaks_s = y_kde_s[peaks]
-        y_peaks_err_s = y_kde_err_s[peaks]
-        
-            
-        temp_widths = np.asarray(peak_widths(y_kde_s, peaks, rel_height=peak_width_factor))
-        temp_widths[2:] = x_kde_s[temp_widths[2:].astype(int)]
-        x_widths_s = temp_widths[3] - temp_widths[2]
-        
-
-        self._peak_data["x_peak_s"] = []
-        self._peak_data["y_peak_s"] = []
-        self._peak_data["y_peak_err_s"] = []
-        self._peak_data["n_peak_s"] = []
-        self._peak_data["x_width_s"] = []
-
-        
-        if(peak_ppf is not None):
-            peak_threshold = self.GPoisPPF(peak_ppf, est_mu, est_lbda, self._max_n_peaks)
-            if(self._verbose):
-                print("For {:3.3}% of distribution, accepting all peaks below {:d}".format(peak_ppf*100.0, peak_threshold))
-        else:
-            peak_threshold = 0
-            print("Accepting all peaks.")    
-
-        n_valid_peaks = 0
-        for x_i, (y_peak, y_peak_err, x_peak, x_width) in enumerate(zip(y_peaks_s, y_peaks_err_s, x_peaks_s, x_widths_s)):
-            
-            N_evt_peak = len(self.SelectRange(data_s, x_peak-x_width/2, x_peak+x_width/2))
-      
-            if(x_i>peak_threshold and N_evt_peak<peak_min_n_evts):
-                if(self._verbose):
-                    print("Peak {:d} has {:d} events. Less than event threshold {:d}. Ignoring...".format(x_i,
-                                                                                                N_evt_peak,
-                                                                                                peak_min_n_evts))
-            else:
-                if(self._verbose):
-                    print("Peak {:d} has {:d} events. Adding...".format(x_i,
-                                                                              N_evt_peak,
-                                                                              peak_min_n_evts))
-                
-                self._peak_data["x_peak_s"].append(x_peak)
-                self._peak_data["y_peak_s"].append(y_peak)
-                self._peak_data["y_peak_err_s"].append(y_peak_err)
-                self._peak_data["n_peak_s"].append(x_i)
-                self._peak_data["x_width_s"].append(x_width)
-                n_valid_peaks+=1
-          
-        
-                    
-        self._peak_data["x_peak_s"] = np.asarray(self._peak_data["x_peak_s"])
-        self._peak_data["y_peak_s"] = np.asarray(self._peak_data["y_peak_s"])
-        self._peak_data["y_peak_err_s"] = np.asarray(self._peak_data["y_peak_err_s"])
-        self._peak_data["n_peak_s"] = np.asarray(self._peak_data["n_peak_s"])
-        self._peak_data["x_width_s"] = np.asarray(self._peak_data["x_width_s"])
-        
-        if(n_valid_peaks<3):
-            print("Minimum 3 peaks covering a threshold {0} events required. Too few events in peaks to proceed. Setting fit fail status and continuing...".format(peak_min_n_evts))
-            self._failed = True
-            return
-        
-        elif(self._peak_data["n_peak_s"][0]>0):
-            if(self._verbose):
-                print("Estimated Pedestal had fewer than threshold {0} events. Shifting the pedestal forward...")
-            self._peak_data["n_peak_s"] = self._peak_data["n_peak_s"] - self._peak_data["n_peak_s"][0]
-                        
-                
-        self._GD_data["x"] = x
-        self._GD_data["y"] = y
-        self._GD_data["y_err"] = y_err
-        self._GD_data["bw"] = bw
-        self._GD_data["N"] = N
-        self._GD_data["nbins"] = nbins
-
-                
-
-        self._GD_data["x_s"] = x_s
-        self._GD_data["y_s"] = y_s
-        self._GD_data["y_err_s"] = y_err_s
-        self._GD_data["x_kde_s"] = x_kde_s
-        self._GD_data["y_kde_s"] = y_kde_s
-        self._GD_data["y_kde_err_s"] = y_kde_err_s
-        self._GD_data["bw_s"] = bw_s
-        self._GD_data["N_s"] = N_s
-        self._GD_data["nbins_s"] = nbins_s
-        
-        self._GD_data["fit_GP_mu"] = max(est_mu, self._eps)
-        self._GD_data["fit_GP_mu_err"] = std_err_mu
-        self._GD_data["fit_GP_lbda"] = max(est_lbda, self._eps)
-        self._GD_data["fit_GP_lbda_err"] = std_err_lbda
-        self._GD_data["fit_G"] = est_gain
-        self._GD_data["fit_G_err"] = std_err_gain
-        self._GD_data["fit_x_0"] = self._peak_data["x_peak_s"][0]
-        self._GD_data["fit_x_0_err"] = self._peak_data["x_width_s"][0]
-        self._GD_data["scaler"] = self.scaler
-        self._GD_data["limit_s"] = limit_s
-        self._GD_data["truncated"] = truncate
-        
-        
-        
-        
-    
-        
-
-    
-    def PreFitGenPoisToPeaks(self, ppf_mainfit, n_bootstrap):
-        
-
-        
-        ########
-        #STEP 1
-        ########
-
-        
-        fit_lbda = self._GD_data["fit_GP_lbda"]
-        fit_dlbda = self._GD_data["fit_GP_lbda_err"]
-        fit_x_0 = self._GD_data["fit_x_0"]
-        
-        fit_mu = self._GD_data["fit_GP_mu"]
-        fit_dmu = self._GD_data["fit_GP_mu_err"]
-        
-        x_peak = self._peak_data["x_peak_s"]
-        n_peak = self._peak_data["n_peak_s"]
-        y_peak = self._peak_data["y_peak_s"]
-        y_peak_err = self._peak_data["y_peak_err_s"]
-
-        
-        
-        log_y_peak = self.Logify(y_peak)
-        dlog_y_peak = y_peak_err/y_peak
-        
-        fit_A_estspec = gpoisson.logpmf(n_peak, fit_mu, fit_lbda)
-        fit_A = np.nanmean(log_y_peak - fit_A_estspec)
-        fit_dA = np.nanmean(2*dlog_y_peak)
-        
-        
-        
-        
-        fit_dict_pois = {
-                "A": fit_A,
-                "mu":fit_mu,
-                "lbda":fit_lbda,
-        }
-
-        fit_pois = Chi2Regression(
-                            self.DummyLogGPois,
-                            n_peak,
-                            log_y_peak,
-                            dlog_y_peak
-                           )
-
-
-        minuit_pois = Minuit(fit_pois, **fit_dict_pois)
-        minuit_pois.errordef=Minuit.LIKELIHOOD
-        minuit_pois.strategy=2
-
-        
-        
-        minuit_pois.errors["mu"] = fit_dmu
-        minuit_pois.errors["lbda"] = fit_dlbda
-        
-        minuit_pois.limits["lbda"] = (self._eps, 1-self._eps)
-        minuit_pois.limits["mu"] = (self._eps, None)
-        
-        minuit_pois.migrad(ncall= self._n_call_minuit,
-                           iterate =self._n_iterations_minuit)
-        minuit_pois.hesse()
-        
-        
-    
-        
-        
-        
-        if(not minuit_pois.valid):
-            if(self._verbose):
-                print('Minuit returned invalid for Generalised Poisson fit to peaks. Using basic estimated parameters instead...')
-            self._GD_data["fit_GP_mu"] = fit_mu
-            self._GD_data["fit_GP_mu_err"] = fit_dmu
-            self._GD_data["fit_GP_lbda"] = fit_lbda
-            self._GD_data["fit_GP_lbda_err"] = fit_dlbda
-            self._GD_data["fit_GP_A"] = fit_A
-            self._GD_data["fit_GP_A_err"] = fit_dA
-            
-  
-            
-        else:
-            self._GD_data["fit_GP_mu"] = minuit_pois.values["mu"]
-            self._GD_data["fit_GP_mu_err"] = minuit_pois.errors["mu"]
-            self._GD_data["fit_GP_lbda"] = minuit_pois.values["lbda"]
-            self._GD_data["fit_GP_lbda_err"] = minuit_pois.errors["lbda"]
-            self._GD_data["fit_GP_A"] =  minuit_pois.values["A"]
-            self._GD_data["fit_GP_A_err"] =  minuit_pois.errors["A"]
-            
-            
-            
-        k_low = 0
-        k_hi = self.GPoisPPF(ppf_mainfit,
-                             self._GD_data["fit_GP_mu"],
-                             self._GD_data["fit_GP_lbda"],
-                             self._max_n_peaks
-                            )
-
-        self._GD_data["fit_GP_k_low"] = k_low
-        self._GD_data["fit_GP_k_hi"] = k_hi
-
-            
-    
-    
-        if(self._GD_data["fit_GP_k_hi"]<2):
-            if(self._verbose):
-                print("Too few primary Geiger peaks to fit. Artificially increasing max number of peaks from {:d} to 2...".format(self._GD_data["fit_GP_k_hi"]))
-            self._GD_data["fit_GP_k_hi"] = max(self._GD_data["fit_GP_k_hi"], 2)
-            
-        elif(self._GD_data["fit_GP_k_hi"]>self._max_n_peaks):
-            if(self._verbose):
-                print("Too many primary Geiger peaks to fit. Artificially decreasing max number of peaks from {:d} to {:d}...".format(self._GD_data["fit_GP_k_hi"],
-                                                                                                                   self._max_n_peaks))
-            self._GD_data["fit_GP_k_hi"] = min(self._GD_data["fit_GP_k_hi"], self._max_n_peaks)
-
-        
-    
-    
-    
-    def PreFitSigmaAndGain(self, data, n_bootstrap):
-    
-        
-        peak_val = self._peak_data["n_peak_s"]
-        
-        data_s = np.squeeze(self.scaler.transform(data.reshape(-1,1)))
-        
-        self._peak_data["peak_mean_s"]=[]
-        self._peak_data["peak_mean_err_s"]=[]
-        self._peak_data["peak_var_s"]=[]
-        self._peak_data["peak_var_err_s"]=[]
-            
-        cond_peak = np.ones(len(peak_val), bool)
-        for x_i, y_peak, x_peak, x_width in zip(self._peak_data["n_peak_s"],
-                                                           self._peak_data["y_peak_s"],
-                                                           self._peak_data["x_peak_s"],
-                                                           self._peak_data["x_width_s"]):
-        
-            data_range = self.SelectRange(data_s, x_peak - x_width/2, x_peak + x_width/2)
-
-            try:
-                est_mean, std_err_mean, _ = self.Bootstrap(data_range,
-                                                               np.mean,
-                                                               n_samples=n_bootstrap)
-
-                est_var, std_err_var, _ = self.Bootstrap(data_range,
-                                                             np.var,
-                                                             n_samples=n_bootstrap)
-
-                self._peak_data["peak_mean_s"].append(est_mean)
-                self._peak_data["peak_mean_err_s"].append(std_err_mean)
-                self._peak_data["peak_var_s"].append(est_var)
-                self._peak_data["peak_var_err_s"].append(std_err_var)
-                cond_peak[x_i]=True
-            except:
-                cond_peak[x_i]=False
-                
-                
-                self._peak_data["peak_mean_s"].append(None)
-                self._peak_data["peak_mean_err_s"].append(None)
-                self._peak_data["peak_var_s"].append(None)
-                self._peak_data["peak_var_err_s"].append(None)
-                
-        if(sum(cond_peak)<3):
-            print("Fewer than 3 peaks had adequate statistics to determine peak width. Setting fit fail status...")
-            self._failed=True
-            return
-        
-        self._peak_data["cond_valid_peak"] = cond_peak
-        self._peak_data["peak_mean_s"] = np.asarray(self._peak_data["peak_mean_s"])
-        self._peak_data["peak_mean_err_s"] = np.asarray(self._peak_data["peak_mean_err_s"])
-        self._peak_data["peak_var_s"] = np.asarray(self._peak_data["peak_var_s"])
-        self._peak_data["peak_var_err_s"] = np.asarray(self._peak_data["peak_var_err_s"])
-        
-        n_peaks_select = self._peak_data["n_peak_s"][cond_peak]
-        peaks_var_select = self._peak_data["peak_var_s"][cond_peak]
-        peaks_var_err_select = self._peak_data["peak_var_err_s"][cond_peak]
-        peaks_mean_select = self._peak_data["peak_mean_s"][cond_peak]
-        peaks_mean_err_select = self._peak_data["peak_mean_err_s"][cond_peak]
-
-        fit_lin_var = HuberRegression(self.DummyLinear,
-                                 n_peaks_select,
-                                 peaks_var_select,
-                                 peaks_var_err_select
-                       )
-                
-        
-        m_var_temp = np.median(np.gradient(peaks_var_select, n_peaks_select))
-        c_var_temp = peaks_var_select[np.argmin(peaks_var_select)]
-        
-        fit_dict_lin_var = {
-                    "m": m_var_temp,
-                    "c": c_var_temp,
-                   }
-        
-       
-        minuit_lin_var = Minuit(fit_lin_var, **fit_dict_lin_var)
-        minuit_lin_var.limits["m"] = (self._eps, None)
-        
-        
-        
-        minuit_lin_var.strategy=2
-        
-        minuit_lin_var.migrad(ncall= self._n_call_minuit,
-                              iterate =self._n_iterations_minuit)
-        minuit_lin_var.hesse()
-        
-        
-        
-        fit_lin_G = HuberRegression(self.DummyLinear,
-                                n_peaks_select,
-                                peaks_mean_select,
-                                peaks_mean_err_select
-                               )
-
-        
-        m_G_temp = self._GD_data["fit_G"]
-        c_G_temp = self._GD_data["fit_x_0"]
-        
-
-        
-        fit_dict_lin_G = {
-                    "m": m_G_temp,
-                    "c": c_G_temp,
-                   }
-
-        
-        
-        minuit_lin_G = Minuit(fit_lin_G, **fit_dict_lin_G)
-        
-        minuit_lin_G.limits["m"] = (self._eps, None)
-        minuit_lin_G.strategy=2
-        
-        minuit_lin_G.migrad(ncall= self._n_call_minuit,
-                            iterate =self._n_iterations_minuit)
-        minuit_lin_G.hesse()
-        
-        
-        
-
-            
-        self._GD_data["fit_var0"] = minuit_lin_var.values["c"]
-        self._GD_data["fit_var0_err"] = minuit_lin_var.errors["c"]
-        self._GD_data["fit_var1"] =  minuit_lin_var.values["m"]
-        self._GD_data["fit_var1_err"] =  minuit_lin_var.errors["m"]
-
-
-        self._GD_data["fit_sigma0"] = np.sqrt(np.abs(self._GD_data["fit_var0"]))
-        self._GD_data["fit_sigma0_err"] = (0.5/self._GD_data["fit_sigma0"])*self._GD_data["fit_var0_err"]
-        self._GD_data["fit_sigma1"] =  np.sqrt(np.abs(self._GD_data["fit_var1"]))
-        self._GD_data["fit_sigma1_err"] = (0.5/self._GD_data["fit_sigma1"])*self._GD_data["fit_var1_err"]
-
-
-        self._GD_data["fit_G"] = minuit_lin_G.values["m"]
-        self._GD_data["fit_G_err"] = minuit_lin_G.errors["m"]
-        self._GD_data["fit_x_0"] = minuit_lin_G.values["c"]
-        self._GD_data["fit_x_0_err"] = minuit_lin_G.errors["c"]
-       
-    
-            
-    
-    
-    def PreFitDCR(self,
-                  data,
-                  tau,
-                  tgate,
-                  r_fast,
-                  t_0,
-                  ppf_mainfit,
-                  kernel_kde,
-                  bandwidth_kde,
-                  n_bootstrap_kde,
-                  dcr_lim=0.05):
-        
-
-        G = self._GD_data["fit_G"]
-        mu = self._GD_data["fit_GP_mu"]
-        lbda = self._GD_data["fit_GP_lbda"]
-        x_0 = self._GD_data["fit_x_0"]
-
-        
-        data_norm = (np.squeeze(self.scaler.transform(data.reshape(-1,1))) - x_0)/G
-        
-
-        try:
-            x_dc_kde, y_dc_kde, y_dc_kde_err, _ = self.BootstrapKDE(
-                                               data_norm,
-                                               kernel=kernel_kde,
-                                               bw=bandwidth_kde,
-                                               n_samples=n_bootstrap_kde                      
-                                            )
-        except:
-            print("DCR KDE failed. Setting fit fail status and continuing...")
-            self._failed = True
-            return
-    
-        
-        cond_interpeak = (x_dc_kde>0) & (x_dc_kde<1)
-        x_dc_kde = x_dc_kde[cond_interpeak]
-        y_dc_kde = y_dc_kde[cond_interpeak]
-        y_dc_kde_err = y_dc_kde_err[cond_interpeak]    
-        
-            
-
-            
-        f_lin_kde_dc = interp1d(x_dc_kde, np.arange(0, len(x_dc_kde)), fill_value='extrapolate')
-        f_lin_kde_dc_inv = interp1d(np.arange(0, len(x_dc_kde)), x_dc_kde, fill_value='extrapolate')
-        
-        f_log_kde_dc = interp1d(x_dc_kde, self.Logify(y_dc_kde), fill_value='extrapolate')
-        f_log_kde_upper_dc = interp1d(x_dc_kde, self.Logify(y_dc_kde + y_dc_kde_err), fill_value='extrapolate')
-        f_log_kde_lower_dc = interp1d(x_dc_kde, self.Logify(y_dc_kde - y_dc_kde_err), fill_value='extrapolate')
-      
-            
-        cond_inrange = (x_dc_kde>0) & (x_dc_kde<1)
-        x_dc_kde = x_dc_kde[cond_inrange]
-        y_dc_kde = y_dc_kde[cond_inrange]
-        y_dc_kde_err = y_dc_kde_err[cond_inrange]
-            
-            
-        x_dc_min = x_dc_kde[np.argmin(y_dc_kde)]
-        x_dc_lim_min_low = x_dc_min - dcr_lim
-        x_dc_lim_min_hi =  x_dc_min + dcr_lim
-        x_dc_lim_low = 0
-        x_dc_lim_mid = x_dc_min
-        x_dc_lim_hi = min(1-self._eps, f_lin_kde_dc_inv(f_lin_kde_dc(0) + 2*abs(f_lin_kde_dc(x_dc_min) - f_lin_kde_dc(0))))
-        
-        
-        cond_interpeak = (x_dc_kde> x_dc_lim_min_low) & (x_dc_kde<=x_dc_lim_min_hi)
-        cond_p_0toLM = (x_dc_kde> x_dc_lim_low) & (x_dc_kde<=x_dc_lim_mid)
-        cond_p_LMto2LM = (x_dc_kde> x_dc_lim_mid) & (x_dc_kde<=x_dc_lim_hi)
-        
-        
-        
-        
-        mean_p_interpeak = (quad(lambda x: np.exp(f_log_kde_dc(x)),
-                              x_dc_lim_min_low,
-                              x_dc_lim_min_hi)[0])/(2*dcr_lim)
-        
-        
-        mean_p_interpeak_upper = (quad(lambda x: np.exp(f_log_kde_upper_dc(x)),
-                              x_dc_lim_min_low,
-                              x_dc_lim_min_hi)[0])/(2*dcr_lim)
-        
-        mean_p_interpeak_lower = (quad(lambda x: np.exp(f_log_kde_lower_dc(x)),
-                              x_dc_lim_min_low,
-                              x_dc_lim_min_hi)[0])/(2*dcr_lim)
-        
-        
-        dmean_p_interpeak =  abs(mean_p_interpeak_upper -  mean_p_interpeak_lower)
-        
-        
-        
-        int_p_0toLM = quad(lambda x: np.exp(f_log_kde_dc(x)),
-                              x_dc_lim_low,
-                              x_dc_lim_mid)[0]
-        
-        int_p_0toLM_upper = quad(lambda x: np.exp(f_log_kde_upper_dc(x)),
-                              x_dc_lim_low,
-                              x_dc_lim_mid)[0]
-        int_p_0toLM_lower = quad(lambda x: np.exp(f_log_kde_lower_dc(x)),
-                              x_dc_lim_low,
-                              x_dc_lim_mid)[0]
-        
-        dint_p_0toLM = abs(int_p_0toLM_upper -  int_p_0toLM_lower)
-        
-        
-        
-        int_p_LMto2LM = quad(lambda x: np.exp(f_log_kde_dc(x)),
-                              x_dc_lim_mid,
-                              x_dc_lim_hi)[0]
-        
-        int_p_LMto2LM_upper = quad(lambda x: np.exp(f_log_kde_upper_dc(x)),
-                              x_dc_lim_mid,
-                              x_dc_lim_hi)[0]
-        int_p_LMto2LM_lower = quad(lambda x: np.exp(f_log_kde_lower_dc(x)),
-                              x_dc_lim_mid,
-                              x_dc_lim_hi)[0]
-        
-        dint_p_LMto2LM = abs(int_p_LMto2LM_upper -  int_p_LMto2LM_lower)
-        
-            
-      
-        int_p = int_p_0toLM + 0.5*int_p_LMto2LM
-        dint_p = np.sqrt(dint_p_0toLM**2 + 0.5*dint_p_LMto2LM**2)
-
-        
-              
-        DCR = mean_p_interpeak/(max(int_p, self._eps)*4*tau)
-        dDCR = DCR*np.sqrt((dmean_p_interpeak/max(mean_p_interpeak, self._eps))**2  + (dint_p/max(int_p, self._eps))**2)
-        
-        
-        mu_dcr=  DCR*(abs(t_0) + tgate)
-        k_dcr_low = 0
-        k_dcr_hi = self.GPoisPPF(ppf_mainfit,
-                                 mu_dcr,
-                                 self._GD_data["fit_GP_lbda"],
-                                 self._max_n_dcr_peaks
-                                )
-        
-
-        
-        
-        
-        self._DCR_data["x_dc_kde"] = x_dc_kde
-        self._DCR_data["y_dc_kde"] = y_dc_kde
-        self._DCR_data["y_dc_kde_err"] = y_dc_kde_err
-        
-        
-        self._DCR_data["cond_p_0toLM"] = cond_p_0toLM
-        self._DCR_data["cond_p_LMtoL2M"] = cond_p_LMto2LM
-        self._DCR_data["cond_interpeak"] = cond_interpeak
-
-        self._DCR_data["int_p_0toLM"] = int_p_0toLM
-        self._DCR_data["dint_p_0toLM"] = dint_p_0toLM
-        self._DCR_data["int_p_LMto2LM"] = int_p_LMto2LM
-        self._DCR_data["dint_p_LMto2LM"] = dint_p_LMto2LM
-        self._DCR_data["mean_p_interpeak"] = mean_p_interpeak
-        self._DCR_data["dmean_p_interpeak"] = dmean_p_interpeak
-        
-        
-        self._DCR_data["fit_tau"] = tau
-        self._DCR_data["fit_tgate"] = tgate
-        self._DCR_data["fit_t_0"] = t_0
-        self._DCR_data["fit_DCR"] = DCR
-        self._DCR_data["fit_DCR_err"] = dDCR
-        self._DCR_data["fit_dcr_k_low"] = k_dcr_low
-        self._DCR_data["fit_dcr_k_hi"] = k_dcr_hi
-        if(r_fast is None):
-            self._DCR_data["fit_r_fast"] = 0
-        else:
-            self._DCR_data["fit_r_fast"] = r_fast
-
-        if(self._DCR_data["fit_dcr_k_hi"]<2):
-            print("Too few dark peaks to fit. Artificially increasing max number of peaks from {:d} to 2...".format(self._DCR_data["fit_dcr_k_hi"]))
-            self._DCR_data["fit_dcr_k_hi"] = max(self._DCR_data["fit_dcr_k_hi"], 2)
-            
-        elif(self._DCR_data["fit_dcr_k_hi"]>self._max_n_peaks):
-            print("Too many dark peaks to fit. Artificially decreasing max number of peaks from {:d} to {:d}...".format(self._DCR_data["fit_dcr_k_hi"],
-                                                                                                                   self._max_n_peaks))
-            self._DCR_data["fit_dcr_k_hi"] = min(self._DCR_data["fit_dcr_k_hi"], self._max_n_dcr_peaks)
-
-
-        
-            
-    
-    def InitFit(self,
-                data,
-                tau,
-                t_0,
-                r_fast,
-                tgate,
-                pedestal,
-                n_bootstrap,
-                n_bootstrap_kde,
-                bw,
-                peak_dist_factor,
-                peak_width_factor,
-                peak_ppf,
-                peak_prominence,
-                peak_min_n_evts,
-                truncate,
-                kernel_kde,
-                bandwidth_kde,
-                ppf_mainfit,
-                bin_method
-               ):
-        
- 
-            
-
-        
-        if(self._is_histogram):
-            data = np.array(list(chain.from_iterable([[x]*int(y) for x,y in zip(data[:,0], data[:,1])])))
-            
-        
-        if(not self._failed):
-            if(self._verbose):
-                print("Getting histogram and peak values...")
-                
-            self.GetHistAndPeaks(
-                        data,
-                        bw,
-                        peak_dist_factor,
-                        peak_width_factor,
-                        peak_ppf,
-                        peak_prominence,
-                        peak_min_n_evts,
-                        truncate,
-                        pedestal,
-                        bin_method,
-                        n_bootstrap,
-                        n_bootstrap_kde,
-                        bandwidth_kde,
-                        kernel_kde)
-
-        if(not self._failed):
-            if(self._verbose):
-                print("Fitting Generalised Poisson distribution...")  
-            self.PreFitGenPoisToPeaks(ppf_mainfit, n_bootstrap)
-            
-        if(not self._failed):
-            if(self._verbose):
-                print("Fitting peak widths via bootstrap resampling...")  
-            self.PreFitSigmaAndGain(data, n_bootstrap)
-            
-        if(not self._failed):
-            if(self._verbose):
-                print("Estimating DCR from interpeak region...")
-            
-            self.PreFitDCR(data,
-                           tau,
-                           tgate,
-                           r_fast,
-                           t_0,
-                           ppf_mainfit,
-                           kernel_kde,
-                           bandwidth_kde,
-                           n_bootstrap_kde)
-        if(not self._failed):
-            if(self._verbose):
-                print("Prefitting complete.")
-            
-          
-
-
-
-    # FIT
-        
-        
-    def Fit(self,
-            data,
-            **kwargs):
-        
-        
-        
-        kwargs_fit =  self._default_fit_kwargs.copy()
-        kwargs_fit.update(kwargs)
-        
-        if(not isinstance(data,np.ndarray)):    
-            raise Exception("Data is in {0} format. Please provide data in the format of a NumPy array.".format(type(data)))
-    
-        
-        else:
-            
-            if(data.ndim==1):
-                if(self._verbose):
-                    print("Data is 1D. Assuming list of charges as input...")
-                self._is_histogram = False
-                
-            elif(data.ndim==2):
-                if(self._verbose):
-                    print("Data is 2D. Assuming histogram as input...")
-                self._is_histogram = True
-            else:
-                raise Exception("Please provide data in the format of either a list of charges (1D) or a histogram (2D)")
-                
-                
-        if(kwargs_fit["t_0"] is None):
-            raise Exception("t_0 is not set. Please set t_0 as a keyword argument to be the integration period before gate in nanoseconds.")
-            
-        if(kwargs_fit["tgate"] is None):
-            raise Exception("tgate is not set. Please set tgate as a keyword argument to be the integration gate length in nanoseconds.")
-            
-        if(kwargs_fit["tau"] is None):
-            raise Exception("tau is not set. Please set tau as a keyword argument to be the slow component of the SiPM pulse in nanoseconds.")
-            
-        if(kwargs_fit["r_fast"] is None):
-            raise Exception("r_fast is not set. Please set r_fast as a keyword argument to be the fraction of the fast component of the SiPM pulse in nanoseconds.")
-      
-        if(kwargs_fit["pedestal"] is None):
-            raise Exception("pedestal is not set. Please set pedestal to be the pedestal position in units of the input charge spectrum being fitted.")
-      
-        
-        print("Performing fit initialisation...")
-        self.Init()
-        self.InitFit(data, **kwargs_fit)
-
-        if(not self._failed):
-                        
-            print("Fitting...")
-            
-            self.N = len(data)
-
-
-            G = self._GD_data["fit_G"]
-            dG = self._GD_data["fit_G_err"]
-            mu = self._GD_data["fit_GP_mu"]
-            dmu = self._GD_data["fit_GP_mu_err"]
-            lbda = self._GD_data["fit_GP_lbda"]
-            dlbda = self._GD_data["fit_GP_lbda_err"]
-            k_low = self._GD_data["fit_GP_k_low"]
-            k_hi = self._GD_data["fit_GP_k_hi"]
-            x_0 = self._GD_data["fit_x_0"]
-            dx_0 = self._GD_data["fit_x_0_err"]
-            sigma_0 = self._GD_data["fit_sigma0"]
-            sigma_1 = self._GD_data["fit_sigma1"]
-            dsigma_0 = self._GD_data["fit_sigma0_err"]
-            dsigma_1 = self._GD_data["fit_sigma1_err"]
-            dx = self._GD_data["bw_s"]
-
-            tau = self._DCR_data["fit_tau"]
-            tgate = self._DCR_data["fit_tgate"]
-            t_0 = self._DCR_data["fit_t_0"]
-            DCR = self._DCR_data["fit_DCR"]
-            dDCR = self._DCR_data["fit_DCR_err"]
-            r_fast = self._DCR_data["fit_r_fast"]
-            k_dcr_low = self._DCR_data["fit_dcr_k_low"]
-            k_dcr_hi =self._DCR_data["fit_dcr_k_hi"]
-
-
-
-            # beta is set to G/2 good guess
-            self.fit_dict = {
-                "mu":mu,
-                "G":G,
-                "x_0": x_0,
-                "sigma_0": sigma_0,
-                "sigma_1": sigma_1,
-                "r_Ap": 0.25,
-                "pAp": 0.25,
-                "lbda":lbda,
-                "tau": tau,
-                "tgate":tgate,
-                "t_0": t_0,
-                "DCR":DCR,
-                "r_fast":r_fast,
-                "k_low":k_low,
-                "k_hi":k_hi,
-                "k_dcr_low":k_dcr_low,
-                "k_dcr_hi":k_dcr_hi,
-                "dx": dx,
-            }
-            
-
-            bl = BinnedLH(self.DRM, self._GD_data["x_s"], self._GD_data["y_s"])
-            minuit = Minuit(bl, **self.fit_dict)
-
-            minuit.errors["pAp"] = 0.05
-            minuit.errors["r_Ap"] = 0.05
-            if(not np.isnan(dmu) and dmu is not None):
-                minuit.errors["mu"] = abs(dmu)
-            if(not np.isnan(dG) and dG is not None):  
-                minuit.errors["G"] = abs(dG)
-            if(not np.isnan(dx_0) and dx_0 is not None):  
-                minuit.errors["x_0"] = abs(dx_0)                
-            if(not np.isnan(dlbda) and dlbda is not None):
-                minuit.errors["lbda"] = abs(dlbda)
-            if(not np.isnan(dsigma_0) and dsigma_0 is not None):
-                minuit.errors["sigma_0"] = abs(dsigma_0)+0.01
-            if(not np.isnan(dsigma_1) and dsigma_1 is not None):
-                minuit.errors["sigma_1"] = abs(dsigma_0)+0.01
-            if(not np.isnan(dDCR) and dDCR is not None):
-                minuit.errors["DCR"] = abs(dDCR)
-            
-            
-            minuit.limits["mu"] = (self._eps, max(1, k_hi))
-            minuit.limits["G"] = (G-5*dG, G+5*dG)
-            minuit.limits["x_0"] = (x_0-G/2, x_0+G/2)
-            minuit.limits["r_Ap"] = (self._eps, 1-self._eps)
-            minuit.limits["pAp"] = (self._eps, 1-self._eps)
-            minuit.limits["lbda"] = (self._eps, 1-self._eps)    
-            minuit.limits["sigma_0"] = (self._eps, G-self._eps)
-            minuit.limits["sigma_1"] = (self._eps, G-self._eps)
-            minuit.limits["DCR"] = (self._eps, 1)
-            
-  
-
-            minuit.fixed["k_low"] = True
-            minuit.fixed["k_hi"] = True
-            minuit.fixed["tgate"] = True
-            minuit.fixed["tau"] = True
-            minuit.fixed["t_0"] = True
-            minuit.fixed["r_fast"] = True
-            minuit.fixed["dx"] = True
-            minuit.fixed["k_dcr_low"] = True
-            minuit.fixed["k_dcr_hi"] = True
-            minuit.fixed["tgate"] = True
-
-
-            minuit.strategy=2
-            minuit.errordef=minuit.LIKELIHOOD
-            minuit.simplex()
-            minuit.tol=1
-            minuit.migrad(ncall= self._n_call_minuit,
-                          iterate =self._n_iterations_minuit)
-         
-            
-            minuit.hesse()
-            
-            print("tauAp = {:3.3f} ns".format(minuit.values["tau"]*minuit.values["r_Ap"]))
-
-            
-            if(self._verbose):
-                print(minuit)
-            
-            self._fit_minuit = minuit
-            self._fit_values_s = minuit.values.to_dict()
-            self._fit_errors_s = minuit.errors.to_dict()
-            
-            self._fit_values = {}
-            self._fit_errors = {}
-            
-    
-            
-            for key, value in self._fit_values_s.items():
-                if(key=="x_0"):
-                    self._fit_values[key] = value*self.scaler.scale_[0] + self.scaler.center_[0]
-                elif(key=="G"):
-                    self._fit_values[key] = value*self.scaler.scale_[0]
-                elif(key=="sigma_0"):
-                    self._fit_values[key] = value*self.scaler.scale_[0]
-                elif(key=="sigma_1"):
-                    self._fit_values[key] = value*self.scaler.scale_[0]
-                elif(key=="dx"):
-                    self._fit_values[key] = self._GD_data["bw"]
-                else:
-                    self._fit_values[key] = value
-                    
-                    
-        
- 
-            for key, value in self._fit_errors_s.items():
-                if(key=="x_0"):
-                    self._fit_errors[key] = value*self.scaler.scale_[0] + self.scaler.center_[0]
-                elif(key=="G"):
-                    self._fit_errors[key] = value*self.scaler.scale_[0]
-                elif(key=="sigma_0"):
-                    self._fit_errors[key] = value*self.scaler.scale_[0]
-                elif(key=="sigma_1"):
-                    self._fit_errors[key] = value*self.scaler.scale_[0]
-                elif(key=="dx"):
-                    self._fit_errors[key] = 0.1*self._GD_data["bw"]
-                else:
-                    self._fit_errors[key] = value
-            
-            
-            
-            data = []
-            prms = []
-            for par in self._fit_minuit.parameters:
-                if(self._fit_minuit.fixed[par]):
-                    continue
-                
-                prms.append(par)
-                data.append([self._fit_values[par], self._fit_errors[par]])
-            
-            data = np.asarray(data)
-            prms = np.asarray(prms)
-            
-            _print = pd.DataFrame(columns = prms,
-                                    index=["Values", "Errors"], data = data.T).T
-
-            _print.style.format({
-              'Values': lambda val: f'{val:,.2E}',
-              'Errors': lambda val: f'{val:,.2E}'
-            })
-            
-