From be86ff8ed4dd33e3da65ec141e54afd70f4573c3 Mon Sep 17 00:00:00 2001
From: Jack Christopher Hutchinson Rolph <jack.rolph@desy.de>
Date: Sat, 13 Nov 2021 17:35:57 +0100
Subject: [PATCH] Update PeakOTron.py

---
 PeakOTron.py | 35 +++++++++++++++--------------------
 1 file changed, 15 insertions(+), 20 deletions(-)

diff --git a/PeakOTron.py b/PeakOTron.py
index 456286d..54bda55 100644
--- a/PeakOTron.py
+++ b/PeakOTron.py
@@ -108,12 +108,20 @@ class BinnedLH:
         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
+        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)
         y_hat = np.nan_to_num(y_hat, nan=self.eps_inv)
-        logL = self.y*np.log(y_hat+self.eps) - y_hat
+        logL = self.y*(self.Logify(y_hat) - self.Logify(self.y)) + (self.y - y_hat)
         nlogL = -np.sum(logL)
 
         return nlogL
@@ -187,8 +195,8 @@ class PeakOTron:
     
     def __init__(self,
                  verbose=False,
-                 n_call_minuit=10000,
-                 n_iterations_minuit = 100
+                 n_call_minuit=2000,
+                 n_iterations_minuit = 50
                 ):
         
         
@@ -278,7 +286,6 @@ class PeakOTron:
                 "bin_method":"knuth"
         }
         
-        
 
     
 
@@ -2157,14 +2164,15 @@ class PeakOTron:
             if(not np.isnan(dDCR) and dDCR is not None):
                 minuit.errors["DCR"] = abs(dDCR)
             
+            
             minuit.limits["mu"] = (self._eps, max(1, k_hi))
             minuit.limits["G"] = (self._eps, None)
-            minuit.limits["x_0"] = (x_0-3*G, x_0+3*G)
+            minuit.limits["x_0"] = (x_0-G/2, x_0+G/2)
             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, 3*G)
-            minuit.limits["sigma_1"] = (self._eps, 3*G)
+            minuit.limits["sigma_0"] = (self._eps, G)
+            minuit.limits["sigma_1"] = (self._eps, G)
             minuit.limits["DCR"] = (self._eps, None)
             
   
@@ -2252,16 +2260,3 @@ class PeakOTron:
               'Errors': lambda val: f'{val:,.2E}'
             })
             
-            
-            
-            
-            
-
-            
-
-                
- 
-            
-
-            
-      
-- 
GitLab