From 8d666a3fec6c4af0b0a846f232c4ce3138f3e0f5 Mon Sep 17 00:00:00 2001
From: Jack Christopher Hutchinson Rolph <jack.rolph@desy.de>
Date: Wed, 11 Jan 2023 20:05:37 +0100
Subject: [PATCH] Update PeakOTron.py
---
PeakOTron.py | 139 +++++++++++++++++++--------------------------------
1 file changed, 52 insertions(+), 87 deletions(-)
diff --git a/PeakOTron.py b/PeakOTron.py
index d6f33b8..53a2f36 100644
--- a/PeakOTron.py
+++ b/PeakOTron.py
@@ -64,6 +64,7 @@ class PeakOTron:
"bin_centres":None,
"counts_bg":None,
"counts_bgsub":None,
+ "truncate_nsigma0_do_min":None,
"bw":None,
"peaks":{
"peak_position":None,
@@ -102,7 +103,7 @@ class PeakOTron:
self._verbose=verbose
self._n_call_minuit = n_call_minuit
self._n_iterations_minuit = n_iterations_minuit
-
+
self._max_n_peaks = None
self._max_n_peaks_dcr = 6
@@ -350,107 +351,68 @@ class PeakOTron:
##
## Calculates the Chi2 of the resultant fit relative to the original density estimate.
########################################################################################################
-
-
-
-# def __CalcChi2(self, x, y, y_hat, y_err):
-# chi2 = np.nansum(((y_hat - y)/y_err)**2)
-# ndof = len(x) - 10
-# return chi2, ndof
-
-
-
-# def GetChi2(self, prefit=False):
-
-# if(self._m_DRM is None):
-# raise Exception ("Fit has not yet been performed. Please first perform a fit.")
-# else:
-
-# x_all = self._hist_data["bin_numbers"]
-# y_all = self._hist_data["counts"]
-# y_err_all = self._hist_data["counts_error"]
-
-# if(self._fit_kwargs["truncate_nsigma0_do"] is not None):
-# lim_nsigma0 = self._prefit_values_bin["Q_0"] - self._fit_kwargs["truncate_nsigma0_do"]*self._prefit_values_bin["sigma_0"]
-# cond_sel_nsigma0 = (self._hist_data["bin_numbers"]>lim_nsigma0)
-
-# else:
-# cond_sel_nsigma0 = np.full(len(self._hist_data["bin_numbers"]), True).astype(bool)
-
-# x = x_all[cond_sel_nsigma0]
-# y = y_all[cond_sel_nsigma0]
-# y_err = y_err_all[cond_sel_nsigma0]
-
-# y_hat = self.GetModel(x, prefit=prefit, bin_units=True)
-# y_hat_err = np.sqrt(y_hat)
-
-
-# plt.figure()
-# plt.plot(x, (y-y_hat)/y_err)
-# plt.show()
-
-# chi2, ndof = self.__CalcChi2(x, y, y_hat, y_err)
-
-# return chi2, ndof
-
+
- def GetChi2(self, pars=None, prefit=False, full=False, n_sigma_lims=[None, None]):
+ def GetChi2(self, pars=None, prefit=False, n_sigma_lims=[None, None]):
- x_all = self._hist_data["bin_numbers"]
- y_all = self._hist_data["counts"]
- y_err_all = self._hist_data["counts_error"]
+ _x_all = self._hist_data["bin_numbers"]
+ _y_all = self._hist_data["counts"]
+ _y_err_all = self._hist_data["counts_error"]
- Q_0_prefit = self._prefit_values_bin["Q_0"]
- sigma_0_prefit = self._prefit_values_bin["sigma_0"]
+ _Q_0_prefit = self._prefit_values_bin["Q_0"]
+ _sigma_0_prefit = self._prefit_values_bin["sigma_0"]
+
+
if(pars is None):
- y_hat_all = self.GetModel(x_all, prefit=prefit, bin_units=True)
+ _y_hat_all = self.GetModel(_x_all, prefit=prefit, bin_units=True)
else:
- y_hat_all = DRM(x_all, **pars)
+ _y_hat_all = DRM(_x_all, **pars)
- conds = []
-
-
- cond_count = (y_all>1) & (y_hat_all>1)
+ cond_count = (_y_all>1)
- conds.append(cond_count)
-# print("Check nsigma = {:3.3f} for chi2".format(self._fit_kwargs["truncate_nsigma0_do"]))
-
+
if(n_sigma_lims[0] is None):
- if(self._fit_kwargs["truncate_nsigma0_do"] is not None and not full):
- n_sigma_lims[0] = self._fit_kwargs["truncate_nsigma0_do"]
-
+ if(self._hist_data["truncate_nsigma0_do_min"] is not None):
+ _lim_nsigma0_low = _Q_0_prefit -self._hist_data["truncate_nsigma0_do_min"]*_sigma_0_prefit
+ else:
+ _lim_nsigma0_low = -np.inf
+ else:
+ _lim_nsigma0_low = _Q_0_prefit - n_sigma_lims[0]*_sigma_0_prefit
+
+ if(n_sigma_lims[1] is None):
+ _lim_nsigma0_hi = np.inf
+ else:
+ _lim_nsigma0_hi = _Q_0_prefit + n_sigma_lims[1]*_sigma_0_prefit
- if(n_sigma_lims[0] is not None and not full):
- lim_nsigma0_low = Q_0_prefit - n_sigma_lims[0]*sigma_0_prefit
- cond_nsigma0_low = (x_all > lim_nsigma0_low)
- conds.append(cond_nsigma0_low)
-
- if(n_sigma_lims[1] is not None and not full):
- lim_nsigma0_hi = Q_0_prefit + n_sigma_lims[1]*sigma_0_prefit
- cond_nsigma0_hi = (x_all < lim_nsigma0_hi)
- conds.append(cond_nsigma0_hi)
+
+
+ cond_nsigma0_low = (_x_all > _lim_nsigma0_low)
+ cond_nsigma0_hi = (_x_all < _lim_nsigma0_hi)
- conds = np.vstack(conds).all(axis=0)
+ conds = cond_count & cond_nsigma0_low & cond_nsigma0_hi
+
- x = x_all[conds]
- y = y_all[conds]
- y_err = y_err_all[conds]
- y_hat = y_hat_all[conds]
- y_hat_err = np.sqrt(y_hat_all[conds])
+ _x = _x_all[conds]
+ _y = _y_all[conds]
+ _y_err = _y_err_all[conds]
+ _y_hat = _y_hat_all[conds]
+ _y_hat_err = np.sqrt(_y_hat_all[conds])
+
- chi2 = np.nansum(((y_hat - y)/y_hat_err)**2)
- ndof = len(x)-10
-
+ chi2 = np.sum(((_y_hat - _y)/_y_hat_err)**2)
+ ndof = len(_x)-10
+
+
return chi2, ndof
@@ -1076,8 +1038,8 @@ class PeakOTron:
sigma_0_prefit = self._prefit_values_bin["sigma_0"]
- if(self._fit_kwargs["truncate_nsigma0_do"] is not None):
- lim_nsigma0 = Q_0_prefit - self._fit_kwargs["truncate_nsigma0_do"]*sigma_0_prefit
+ if(self._hist_data["truncate_nsigma0_do_min"] is not None):
+ lim_nsigma0 = Q_0_prefit - self._hist_data["truncate_nsigma0_do_min"]*sigma_0_prefit
cond_sel_nsigma0 = (self._hist_data["bin_numbers"]>lim_nsigma0)
else:
@@ -1089,10 +1051,11 @@ class PeakOTron:
y_hat = self.GetModel(x, prefit=prefit, bin_units=plot_in_bins)*bw
y_hat_err = np.sqrt(y_hat)
-
-
+
chi2, ndf = self.GetChi2(prefit=prefit)
+
+
if(xlabel is None):
if(plot_in_bins):
xlabel="Bin"
@@ -1680,14 +1643,14 @@ class PeakOTron:
if(self._fit_kwargs["truncate_nsigma0_chi2scanlim"] is not None):
- if(red_chi2_full<self._fit_kwargs["truncate_nsigma0_chi2scanlim"]):
+ if(red_chi2_range<self._fit_kwargs["truncate_nsigma0_chi2scanlim"]):
if(self._verbose):
print("Fit chi2/NDF less than limit. (Fit: {:3.3f} < Limit {:3.3f}. Scan complete".format(red_chi2,
self._fit_kwargs["truncate_nsigma0_chi2scanlim"]))
self._m_DRM = deepcopy(m_DRM)
self._time_fit = np.sum(time_fits)
- self._fit_kwargs["truncate_nsigma0_do"] = _nsigma0low
+ self._hist_data["truncate_nsigma0_do_min"] = _nsigma0low
break
@@ -1715,8 +1678,10 @@ class PeakOTron:
self._m_DRM = deepcopy(m_DRMs[fit_i_best])
self._time_fit = np.sum(time_fits)
- self._fit_kwargs["truncate_nsigma0_do"] = nsigma_0_scan_range[fit_i_best]
+ #print(self._hist_data["truncate_nsigma0_do_min"] )
+ self._hist_data["truncate_nsigma0_do_min"] = nsigma_0_scan_range[fit_i_best]
+ #print(self._hist_data["truncate_nsigma0_do_min"] )
else:
f_DRM = BinnedLH(DRM,
--
GitLab