switch statistical calculation from RooStats to scipy

kratsg · Jul 18, 2019 · 3acfd8b · 3acfd8b
1 parent ed6410c
commit 3acfd8b
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 2 deletions.
diff --git a/root_optimize/utils.py b/root_optimize/utils.py
@@ -15,6 +15,7 @@
 import itertools
 import numpy as np
 import numexpr as ne
+import scipy
 import os
 import sys
 from time import clock
@@ -196,6 +197,18 @@ def get_scaleFactor(weights, did):
   logger.info( "‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾")
   return scaleFactor
 
+def significance(signalExp, backgroundExp, relativeBkgUncert):
+    """ Numpy/Scipy port of the RooStats function `BinomialExpZ'
+
+    See: https://root.cern.ch/doc/master/NumberCountingUtils_8cxx_source.html
+    """
+    # pylint: disable=invalid-name
+    mainInf = signalExp + backgroundExp
+    tau = 1.0 / backgroundExp / (relativeBkgUncert * relativeBkgUncert)
+    auxiliaryInf = backgroundExp * tau
+    P_Bi = scipy.special.betainc(mainInf, auxiliaryInf + 1, 1.0 / (1.0 + tau))
+    return - scipy.special.ndtri(P_Bi)
+
 #@echo(write=logger.debug)
 def get_significance(signal, bkgd, insignificanceThreshold, bkgdUncertainty, bkgdStatUncertainty, rawBkgd):
   # if not enough events, return string of which one did not have enough
@@ -210,7 +223,7 @@ def get_significance(signal, bkgd, insignificanceThreshold, bkgdUncertainty, bkg
     sig = -3
   else:
     # otherwise, calculate!
-    sig = ROOT.RooStats.NumberCountingUtils.BinomialExpZ(signal, bkgd, bkgdUncertainty)
+    sig = significance(signal, bkgd, bkgdUncertainty)
   return sig
 
 #@echo(write=logger.debug)

diff --git a/root_optimize/version.py b/root_optimize/version.py
@@ -1 +1 @@
-__version__ = '0.6.0'
+__version__ = '0.6.1'