diff --git a/PeakOTron.py b/PeakOTron.py
index 3f1cbaa908292d81a2c183173fbd6987698f149f..a3d724384b6231e5d5fdfa5d292700b5b1a52ff5 100644
--- a/PeakOTron.py
+++ b/PeakOTron.py
@@ -98,6 +98,10 @@ class BandWidthOptimiser:
self.last_arg = arg
loss = self._KDE(*arg)
return loss
+
+
+
+
@@ -278,20 +282,18 @@ class PeakOTron:
"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,
+ "peak_ppf":0.99,
+ "peak_prominence":1e-3,
"kernel_kde":"gaussian",
- "bandwidth_kde":"ISJ",
+ "bandwidth_kde":"optimise",
"ppf_mainfit":1 - 1e-6,
- "bin_method":"knuth"
+ "bin_method":"knuth",
+ "alpha":0.95
}
@@ -388,11 +390,6 @@ class PeakOTron:
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
@@ -406,6 +403,53 @@ class PeakOTron:
return _count
+
+ def EmpiricalPPF(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 GainFFT(self, x_kde, y_kde):
+ n = len(x_kde)
+ dt = abs(x_kde[1]-x_kde[0])
+ fhat = np.fft.fft(y_kde, n) #computes the fft
+ psd = fhat * np.conj(fhat)/n
+ freq = (1/(dt*n)) * np.arange(n) #frequency array
+ idxs_half = np.arange(1, np.floor(n/2), dtype=np.int32)
+
+ y = np.abs(psd[idxs_half])
+ x = freq[idxs_half]
+
+ peaks, _= find_peaks(y)
+ idx_sort = np.argsort(y[peaks])
+ y_peaks = y[peaks][idx_sort]
+ x_peaks = x[peaks][idx_sort]
+
+ x_peak = x_peaks[-1]
+
+ plt.figure(figsize=(10.0, 10.0))
+ plt.title("$G$={:s} Q.U.".format(self.LatexFormat((1/x_peak)*self.scaler.scale_[0])), fontsize=20)
+ plt.plot(x, y, lw=3)
+ plt.ylabel("Energy Spectral Density $|S(f_{PH})|^{2}$", fontsize=25)
+ plt.xlabel("$f_{PH}{$ [Inv. Q.U.]", fontsize=20)
+ plt.axvline(x=x_peak, color="purple", lw=3, linestyle="--", label="Maximum Peak")
+ plt.yscale("log")
+ plt.xscale("log")
+ plt.xticks(fontsize=20)
+ plt.yticks(fontsize=20)
+ plt.grid(which="both")
+ plt.legend(fontsize=15)
+ plt.show()
+
+ gain = 1/x_peak
+ return gain
+
def EstLbda(self, data_sub_x_0):
@@ -905,12 +949,7 @@ class PeakOTron:
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])
@@ -1011,41 +1050,41 @@ class PeakOTron:
- def PlotGainEst(self, ax):
+# def PlotGainEst(self, ax):
- c = self._peak_data["cond_valid_peak"]
+# 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.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.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.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")
+# 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")
@@ -1253,9 +1292,6 @@ class PeakOTron:
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)
@@ -1298,27 +1334,29 @@ class PeakOTron:
peak_width_factor,
peak_ppf,
peak_prominence,
- peak_min_n_evts,
- truncate,
- pedestal,
bin_method,
n_bootstrap,
n_bootstrap_kde,
bandwidth_kde,
- kernel_kde):
+ kernel_kde,
+ alpha
+
+ ):
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):
+
+ ppf = self.EmpiricalPPF(data_s)
+
if(bin_method == "knuth"):
- bw_s, bins_s = knuth_bin_width(data_s, return_bins=True)
+ bw_s, bins_s = knuth_bin_width(data_s[data_s<ppf(alpha)], return_bins=True)
elif(bin_method=="freedman"):
- bw_s, bins_s = freedman_bin_width(data_s, return_bins=True)
+ bw_s, bins_s = freedman_bin_width(data_s[data_s<ppf(alpha)], return_bins=True)
elif(bin_method=="scott"):
- bw_s, bins_s = scott_bin_width(data_s, return_bins=True)
+ bw_s, bins_s = scott_bin_width(data_s[data_s<ppf(alpha)], return_bins=True)
else:
raise Exception ("Please specify a bin width by setting 'bw' or choose a method from 'knuth', 'freedman' or 'scott'")
@@ -1383,7 +1421,9 @@ class PeakOTron:
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)
+ n_samples=n_bootstrap_kde,
+ alpha=alpha
+ )
f_lin_kde_s = interp1d(x_kde_s, np.arange(len(x_kde_s)), fill_value='extrapolate')
@@ -1396,51 +1436,12 @@ class PeakOTron:
-
- 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)
+ est_gain_s = self.GainFFT(x_kde_s, y_kde_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 ...")
-
-
+ est_gain_lin = abs(f_lin_kde_s(est_gain_s) - f_lin_kde_s(0))
peaks, properties = find_peaks(y_kde_s,
@@ -1448,26 +1449,74 @@ class PeakOTron:
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
-
-
+ 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_peaks_s = x_kde_s[peaks]
y_peaks_s = y_kde_s[peaks]
y_peaks_err_s = y_kde_err_s[peaks]
+ 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]
+
+
+ plt.figure()
+ plt.plot(x_kde_s, y_kde_s)
+ for x_peak in x_peaks_s:
+ plt.axvline(x=x_peak)
+ plt.yscale("log")
+ plt.show()
+
+
- 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]
+ est_mu_s, std_err_mu_s, _ = self.Bootstrap(data_s-pedestal_s,
+ self.EstMu,
+ n_samples=n_bootstrap)
+ est_lbda_s, std_err_lbda_s, _ = self.Bootstrap(data_s - pedestal_s,
+ self.EstLbda,
+ n_samples=n_bootstrap)
+
+
+
+ _, std_err_gain_s, _ = self.Bootstrap(data_s - pedestal_s,
+ self.EstGain,
+ n_samples=n_bootstrap)
+
self._peak_data["x_peak_s"] = []
self._peak_data["y_peak_s"] = []
@@ -1475,9 +1524,8 @@ class PeakOTron:
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)
+ peak_threshold = max(3, self.GPoisPPF(peak_ppf, est_mu_s, est_lbda_s, 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:
@@ -1487,18 +1535,19 @@ class PeakOTron:
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))
+ data_peak_s = self.SelectRange(data_s, x_peak-x_width/2, x_peak+x_width/2)
+
+ N_evt_peak = len(data_peak_s)
+
+ if(x_i==0):
+ est_x_0_s, std_err_x_0_s, _ = self.Bootstrap(data_peak_s,
+ np.mean,
+ n_samples=n_bootstrap)
+
- 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(x_i<peak_threshold):
if(self._verbose):
- print("Peak {:d} has {:d} events. Adding...".format(x_i,
- N_evt_peak,
- peak_min_n_evts))
+ print("Peak {:d} has {:d} events. Adding...".format(x_i, N_evt_peak))
self._peak_data["x_peak_s"].append(x_peak)
self._peak_data["y_peak_s"].append(y_peak)
@@ -1515,15 +1564,6 @@ class PeakOTron:
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
@@ -1545,24 +1585,18 @@ class PeakOTron:
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["fit_GP_mu"] = max(est_mu_s, self._eps)
+ self._GD_data["fit_GP_mu_err"] = std_err_mu_s
+ self._GD_data["fit_GP_lbda"] = min(1, max(est_lbda_s, self._eps))
+ self._GD_data["fit_GP_lbda_err"] = std_err_lbda_s
+ self._GD_data["fit_G"] = est_gain_s
+ self._GD_data["fit_G_err"] = std_err_gain_s
+ self._GD_data["fit_x_0"] = est_x_0_s
+ self._GD_data["fit_x_0_err"] = std_err_x_0_s
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):
@@ -1682,15 +1716,13 @@ class PeakOTron:
- def PreFitSigmaAndGain(self, data, n_bootstrap):
+ def PreFitSigma(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"]=[]
@@ -1703,16 +1735,13 @@ class PeakOTron:
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
@@ -1720,8 +1749,6 @@ class PeakOTron:
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)
@@ -1731,16 +1758,12 @@ class PeakOTron:
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,
@@ -1769,38 +1792,7 @@ class PeakOTron:
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"]
@@ -1813,13 +1805,7 @@ class PeakOTron:
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(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"]
-
+
@@ -2009,28 +1995,24 @@ class PeakOTron:
def InitFit(self,
data,
+ bw,
+ bin_method,
tau,
t_0,
r_fast,
tgate,
- pedestal,
n_bootstrap,
n_bootstrap_kde,
- bw,
+ kernel_kde,
+ bandwidth_kde,
peak_dist_factor,
peak_width_factor,
peak_ppf,
peak_prominence,
- peak_min_n_evts,
- truncate,
- kernel_kde,
- bandwidth_kde,
ppf_mainfit,
- bin_method
+ alpha
):
-
-
if(self._is_histogram):
@@ -2041,21 +2023,18 @@ class PeakOTron:
if(self._verbose):
print("Getting histogram and peak values...")
- self.GetHistAndPeaks(
- data,
+ 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)
+ kernel_kde,
+ alpha)
if(not self._failed):
if(self._verbose):
@@ -2065,7 +2044,7 @@ class PeakOTron:
if(not self._failed):
if(self._verbose):
print("Fitting peak widths via bootstrap resampling...")
- self.PreFitSigmaAndGain(data, n_bootstrap)
+ self.PreFitSigma(data, n_bootstrap)
if(not self._failed):
if(self._verbose):
@@ -2131,9 +2110,6 @@ class PeakOTron:
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()
@@ -2219,7 +2195,7 @@ class PeakOTron:
minuit.limits["mu"] = (self._eps, max(1, k_hi))
minuit.limits["G"] = (self._eps, None)
- minuit.limits["x_0"] = (x_0-G/2, x_0+G/2)
+ minuit.limits["x_0"] = (None, None)
minuit.limits["alpha"] = (self._eps, 1-self._eps)
minuit.limits["r_beta"] = (self._eps, 1-self._eps)
minuit.limits["lbda"] = (self._eps, 1-self._eps)
@@ -2243,6 +2219,7 @@ class PeakOTron:
minuit.strategy=2
minuit.errordef=minuit.LIKELIHOOD
+ minuit.simplex()
minuit.migrad(ncall= self._n_call_minuit,
iterate =self._n_iterations_minuit)
minuit.hesse()