From aac826f5f9f9a9691bafbd33ee9436abfac2b453 Mon Sep 17 00:00:00 2001
From: Jack Christopher Hutchinson Rolph <jack.rolph@desy.de>
Date: Sun, 5 Dec 2021 17:33:35 +0100
Subject: [PATCH] Update PeakOTron.py
---
PeakOTron.py | 92 ++++++++++++++++++++++++----------------------------
1 file changed, 43 insertions(+), 49 deletions(-)
diff --git a/PeakOTron.py b/PeakOTron.py
index d4d6717..72fa552 100644
--- a/PeakOTron.py
+++ b/PeakOTron.py
@@ -101,10 +101,41 @@ class BandWidthOptimiser:
+
+class UnbinnedLH:
+
+
+ def __init__(self, f, data, x):
+ self.f = f
+ self.x = x
+ 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)
+
+
+ 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)
+ 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)
+
-
+ return nlogL
+
class BinnedLH:
@@ -668,13 +699,15 @@ class PeakOTron:
return PHs
- def ApplyDCR(self, x, y_light, x_0, G, DCR, tgate, t_0, tau, r_fast, lbda, dx, k_dcr_low, k_dcr_hi):
+ def ApplyDCR(self, x, y_light, x_0, G, DCR, tgate, t_0, tau, r_fast, lbda, 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')
+ 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)
@@ -808,8 +841,7 @@ class PeakOTron:
tau,
tgate,
t_0,
- r_fast,
- dx
+ r_fast
):
@@ -838,7 +870,6 @@ class PeakOTron:
tau,
r_fast,
lbda,
- dx,
k_dcr_low,
k_dcr_hi)
@@ -1472,30 +1503,7 @@ class PeakOTron:
x_widths_s = temp_widths[3] - temp_widths[2]
pedestal_s = x_kde_s[peaks][0]
- pedestal = pedestal_s*self.scaler.scale_[0] + self.scaler.center_[0]
- est_gain = est_gain_s*self.scaler.scale_[0]
- # cond_ped_kde_s = (x_kde_s>pedestal_s - est_gain_s/2)
- # cond_ped_s = (x_s>pedestal_s - est_gain_s/2)
- # cond_ped = (x>pedestal - est_gain/2)
- # cond_ped_peaks = (x_peaks_s>pedestal_s - est_gain_s/2)
-
- # x_kde_s = x_kde_s[cond_ped_kde_s]
- # y_kde_s = y_kde_s[cond_ped_kde_s]
- # y_kde_err_s = y_kde_err_s[cond_ped_kde_s]
-
- # x_s = x_s[cond_ped_s]
- # y_s = y_s[cond_ped_s]
- # y_err_s = y_err_s[cond_ped_s]
-
- # x_peaks_s = x_peaks_s[cond_ped_peaks]
- # y_peaks_s = y_peaks_s[cond_ped_peaks]
- # y_peaks_err_s = y_peaks_err_s[cond_ped_peaks]
- # x_widths_s = x_widths_s[cond_ped_peaks]
-
- # x = x[cond_ped]
- # y = y[cond_ped]
- # y_err = y_err[cond_ped]
@@ -1553,34 +1561,28 @@ class PeakOTron:
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"])
-
-
+
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["data"] = data
+
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["data_s"] = data_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_s, self._eps)
self._GD_data["fit_GP_mu_err"] = std_err_mu_s
@@ -2133,7 +2135,6 @@ class PeakOTron:
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"]
@@ -2164,12 +2165,11 @@ class PeakOTron:
"k_low":k_low,
"k_hi":k_hi,
"k_dcr_low":k_dcr_low,
- "k_dcr_hi":k_dcr_hi,
- "dx": dx,
+ "k_dcr_hi":k_dcr_hi
}
- bl = BinnedLH(self.DRM, self._GD_data["x_s"], self._GD_data["y_s"])
+ bl = UnbinnedLH(self.DRM, self._GD_data["data_s"], self._GD_data["x_kde_s"])
minuit = Minuit(bl, **self.fit_dict)
minuit.errors["alpha"] = 0.01
@@ -2201,14 +2201,12 @@ class PeakOTron:
minuit.limits["DCR"] = (self._eps, None)
-
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
@@ -2242,8 +2240,6 @@ class PeakOTron:
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
@@ -2259,8 +2255,6 @@ class PeakOTron:
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
--
GitLab