Update PeakOTron.py

PeakOTron.py

@@ -96,7 +96,6 @@ class BandWidthOptimiser:
 
                
                                               
-                                
 class BinnedLH:
     
 
@@ -108,10 +107,12 @@ 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)
 
@@ -120,6 +121,7 @@ class BinnedLH:
     
 
     
+    
 class Chi2Regression:
     
 
@@ -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)