diff --git a/PeakOTron.py b/PeakOTron.py index fcbd0e846c01f45accd46a25c5ac09fb8759e141..456286d7ad733b0f2b806cf185d0fbf088744913 100644 --- a/PeakOTron.py +++ b/PeakOTron.py @@ -95,8 +95,7 @@ class BandWidthOptimiser: - - + class BinnedLH: @@ -108,15 +107,18 @@ class BinnedLH: 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 def __call__(self, *arg): self.last_arg = arg y_hat = self.f(self.x, *arg) + y_hat = np.nan_to_num(y_hat, nan=self.eps_inv) logL = self.y*np.log(y_hat+self.eps) - y_hat nlogL = -np.sum(logL) return nlogL + @@ -186,7 +188,7 @@ class PeakOTron: def __init__(self, verbose=False, n_call_minuit=10000, - n_iterations_minuit = 50 + n_iterations_minuit = 100 ): @@ -267,7 +269,7 @@ class PeakOTron: "n_bootstrap_kde":100, "bw":None, "peak_dist_factor":0.8, - "peak_width_factor":0.5, + "peak_width_factor":0.65, "prominence": 0.0, "min_n_peak_evts": 100, "kernel_kde":"gaussian", @@ -1681,7 +1683,6 @@ class PeakOTron: minuit_lin_var = Minuit(fit_lin_var, **fit_dict_lin_var) - minuit_lin_var.limits["m"] = (self._eps, None) minuit_lin_var.strategy=2 @@ -1713,7 +1714,6 @@ class PeakOTron: minuit_lin_G = Minuit(fit_lin_G, **fit_dict_lin_G) - minuit_lin_G.limits["m"] = (self._eps, None) minuit_lin_G.strategy=2 @@ -1725,7 +1725,7 @@ class PeakOTron: - if(minuit_lin_var.valid): + if((minuit_lin_var.values["m"]>self._eps) and minuit_lin_var.valid): self._GD_data["fit_var0"] = minuit_lin_var.values["c"] self._GD_data["fit_var0_err"] = minuit_lin_var.errors["c"] @@ -1753,7 +1753,7 @@ class PeakOTron: - if(minuit_lin_G.valid): + if((minuit_lin_G.values["m"]>self._eps) and minuit_lin_G.valid): 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"] @@ -2140,8 +2140,8 @@ class PeakOTron: bl = BinnedLH(self.DRM, self._GD_data["x_s"], self._GD_data["y_s"]) minuit = Minuit(bl, **self.fit_dict) - minuit.errors["alpha"] = 0.025 - minuit.errors["r_beta"] = 0.025 + minuit.errors["alpha"] = 0.01 + minuit.errors["r_beta"] = 0.01 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): @@ -2163,8 +2163,8 @@ class PeakOTron: 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) - minuit.limits["sigma_0"] = (self._eps, None) - minuit.limits["sigma_1"] = (self._eps, None) + minuit.limits["sigma_0"] = (self._eps, 3*G) + minuit.limits["sigma_1"] = (self._eps, 3*G) minuit.limits["DCR"] = (self._eps, None)