Merge pull request #41 from rkansal47/master

`include_threshold` option for AutoMCStats

src/rhalphalib/model.py

@@ -365,13 +365,19 @@ def renderCard(self, outputFilename, workspaceName):
                     if effect != "-":
                         fout.write(effect + "\n")
 
-    def autoMCStats(self, epsilon=0, threshold=0, include_signal=0, channel_name=None):
+    def autoMCStats(self, epsilon=0, threshold=0, include_signal=0, channel_name=None, include_threshold=0):
         """
         Barlow-Beeston-lite method i.e. single stats parameter for all processes per bin.
         Same general algorithm as described in
         https://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part2/bin-wise-stats/
         but *without the analytic minimisation*.
-        `include_signal` only refers to whether signal stats are included in the *decision* to use bb-lite or not.
+
+        epsilon: changes the minimum possible value from 0 -> epsilon.
+        threshold: the threshold used to decide whether or not bblite is used.
+        include_signal: whether signal stats are included in the *decision* only to use bb-lite or not.
+        channel_name: in case a custom name for the channel needs to be specified, e.g. for tying together mcstats parameters in different bins.
+        include_threshold: the *uncertainty* threshold used to decide whether an mcstats parameter is used at all.
+            i.e. if `threshold = 0.01`, then bins with less than a 1% mc stats uncertainty will not be included.
         """
         if not len(self._samples):
             raise RuntimeError("Channel %r has no samples for which to run autoMCStats" % (self))
@@ -402,16 +408,19 @@ def autoMCStats(self, epsilon=0, threshold=0, include_signal=0, channel_name=Non
                 for sample in self._samples.values():
                     if sample._sampletype == Sample.SIGNAL:
                         sample_name = None if channel_name is None else channel_name + "_" + sample._name[sample._name.find("_") + 1 :]
-                        sample.autoMCStats(epsilon=epsilon, sample_name=sample_name, bini=i)
+                        sample.autoMCStats(epsilon=epsilon, threshold=include_threshold, sample_name=sample_name, bini=i)
 
                 continue
 
             neff_bb = ntot_bb**2 / etot2_bb
             if neff_bb <= threshold:
                 for sample in self._samples.values():
                     sample_name = None if channel_name is None else channel_name + "_" + sample._name[sample._name.find("_") + 1 :]
-                    sample.autoMCStats(epsilon=epsilon, sample_name=sample_name, bini=i)
+                    sample.autoMCStats(epsilon=epsilon, threshold=include_threshold, sample_name=sample_name, bini=i)
             else:
+                if (np.sqrt(etot2_bb) / (ntot_bb + 1e-12)) < include_threshold:
+                    continue
+
                 effect_up = np.ones_like(first_sample._nominal)
                 effect_down = np.ones_like(first_sample._nominal)