WIP: First implementation of ATLAS-style global impacts

python/tool_base/Impacts.py

@@ -81,6 +81,12 @@ def attach_args(self, group):
             help="""write output json to a
             file""",
         )
+        group.add_argument(
+            "--globalImpacts",
+            "-g",
+            action="store_true",
+            help="""Run the global impacts calculation""",
+        )
         group.add_argument("--approx", default=None, choices=["hesse", "robust"], help="""Calculate impacts using the covariance matrix instead""")
         group.add_argument("--noInitialFit", action="store_true", default=False, help="""Do not look for results from the initial Fit""")
 
@@ -118,7 +124,7 @@ def run_method(self):
             sys.exit(0)
         elif self.args.approx == "robust" and self.args.doFits:
             self.job_queue.append(
-                "combine -M MultiDimFit -n _approxFit_%(name)s --algo none --redefineSignalPOIs %(poistr)s --floatOtherPOIs 1 --saveInactivePOI 1 --robustHesse 1 %(pass_str)s"
+                "combine -M MultiDimFit -n _approxFit_%(name)s --algo none --redefineSignalPOIs %(poistr)s --floatOtherPOIs 1 --saveInactivePOI 1 --robustHesse 1 --saveFitResult %(pass_str)s"
                 % {"name": name, "poistr": poistr, "pass_str": pass_str}
             )
             self.flush_queue()
@@ -134,17 +140,19 @@ def run_method(self):
             if self.args.splitInitial:
                 for poi in poiList:
                     self.job_queue.append(
-                        "combine -M MultiDimFit -n _initialFit_%(name)s_POI_%(poi)s --algo singles --redefineSignalPOIs %(poistr)s --floatOtherPOIs 1 --saveInactivePOI 1 -P %(poi)s %(pass_str)s"
+                        "combine -M MultiDimFit -n _initialFit_%(name)s_POI_%(poi)s --algo singles --redefineSignalPOIs %(poistr)s --floatOtherPOIs 1 --saveInactivePOI 1 -P %(poi)s %(pass_str)s --saveFitResult"
                         % {"name": name, "poi": poi, "poistr": poistr, "pass_str": pass_str}
                     )
             else:
                 self.job_queue.append(
-                    "combine -M MultiDimFit -n _initialFit_%(name)s --algo singles --redefineSignalPOIs %(poistr)s %(pass_str)s"
+                    "combine -M MultiDimFit -n _initialFit_%(name)s --algo singles --redefineSignalPOIs %(poistr)s %(pass_str)s --saveFitResult"
                     % {"name": name, "poistr": poistr, "pass_str": pass_str}
                 )
             self.flush_queue()
             sys.exit(0)
 
+        constVarValues = {} # The values of constant vars (i.e. global obs) in the fit
+
         # Read the initial fit results
         if not self.args.noInitialFit:
             initialRes = {}
@@ -154,20 +162,27 @@ def run_method(self):
                     rfr = fResult.Get("fit_mdf")
                     fResult.Close()
                     initialRes = utils.get_roofitresult(rfr, poiList, poiList)
+                    constVarValues = utils.get_rfr_constvars(f"multidimfit_approxFit_{name}.root", "fit_mdf")
+
                 elif self.args.approx == "robust":
                     fResult = ROOT.TFile(f"robustHesse_approxFit_{name}.root")
                     floatParams = fResult.Get("floatParsFinal")
                     rfr = fResult.Get("h_correlation")
                     rfr.SetDirectory(0)
                     fResult.Close()
                     initialRes = utils.get_robusthesse(floatParams, rfr, poiList, poiList)
+                    constVarValues = utils.get_rfr_constvars(f"multidimfit_approxFit_{name}.root", "fit_mdf")
+
             elif self.args.splitInitial:
-                for poi in poiList:
+                for idx, poi in enumerate(poiList):
                     initialRes.update(
                         utils.get_singles_results("higgsCombine_initialFit_%(name)s_POI_%(poi)s.MultiDimFit.mH%(mh)s.root" % vars(), [poi], poiList)
                     )
+                    if idx == 0: # We only need to get this once, it will be the same in each file
+                        constVarValues = utils.get_rfr_constvars(f"multidimfit_initialFit_{name}_POI_{poi}.root", "fit_mdf")
             else:
                 initialRes = utils.get_singles_results("higgsCombine_initialFit_%(name)s.MultiDimFit.mH%(mh)s.root" % vars(), poiList, poiList)
+                constVarValues = utils.get_rfr_constvars(f"multidimfit_initialFit_{name}.root", "fit_mdf")
 
         ################################################
         # Build the parameter list
@@ -216,7 +231,10 @@ def run_method(self):
                     set_parameters_str += setParam + ","
             self.args.setParameters = set_parameters_str.rstrip(",")
 
-        prefit = utils.prefit_from_workspace(ws, "w", paramList, self.args.setParameters)
+        # TODO: Now that we extract the const values from the actual fit result, we probably don't need
+        # to parse --setParameters here - the only ones of interest are those which correspond to global
+        # observables ("X_In") and these should always be in constVarValues
+        prefit = utils.prefit_from_workspace(ws, "w", paramList, self.args.setParameters, constVarValues)
         res = {}
         if not self.args.noInitialFit:
             res["POIs"] = []
@@ -235,6 +253,13 @@ def run_method(self):
                     "combine -M MultiDimFit -n _paramFit_%(name)s_%(param)s --algo impact --redefineSignalPOIs %(poistr)s -P %(param)s --floatOtherPOIs 1 --saveInactivePOI 1 %(pass_str)s"
                     % vars()
                 )
+                if self.args.globalImpacts and "prefit" in pres and pres["type"] != "Unconstrained":
+                    gobsHi = pres["prefit"][2]
+                    gobsLo = pres["prefit"][0]
+                    self.job_queue.append(
+                        f"combine -M MultiDimFit -n _globalFit_{name}_{param}_hi --algo fixed --redefineSignalPOIs {poistr} -P {param}_In --floatOtherPOIs 1 --saveInactivePOI 1 {pass_str} --fixedPointPOIs {param}_In={gobsHi}")
+                    self.job_queue.append(
+                        f"combine -M MultiDimFit -n _globalFit_{name}_{param}_lo --algo fixed --redefineSignalPOIs {poistr} -P {param}_In --floatOtherPOIs 1 --saveInactivePOI 1 {pass_str} --fixedPointPOIs {param}_In={gobsLo}")
             else:
                 if self.args.approx == "hesse":
                     paramScanRes = utils.get_roofitresult(rfr, [param], poiList + [param])
@@ -247,6 +272,11 @@ def run_method(self):
                     paramScanRes = utils.get_singles_results(
                         "higgsCombine_paramFit_%(name)s_%(param)s.MultiDimFit.mH%(mh)s.root" % vars(), [param], poiList + [param]
                     )
+                    if self.args.globalImpacts and pres["type"] != "Unconstrained":
+                        globalFitHiRes = utils.get_fixed_results(
+                            f"higgsCombine_globalFit_{name}_{param}_hi.MultiDimFit.mH{mh}.root", poiList)
+                        globalFitLoRes = utils.get_fixed_results(
+                            f"higgsCombine_globalFit_{name}_{param}_lo.MultiDimFit.mH{mh}.root", poiList)
                 if paramScanRes is None:
                     missing.append(param)
                     continue
@@ -258,6 +288,18 @@ def run_method(self):
                             "impact_" + p: max(list(map(abs, (x - paramScanRes[param][p][1] for x in (paramScanRes[param][p][2], paramScanRes[param][p][0]))))),
                         }
                     )
+                    if self.args.globalImpacts and pres["type"] != "Unconstrained":
+                        if self.args.approx is not None:
+                            # print(param)
+                            # print(prefit)
+                            symm_prefit = (prefit[param]["prefit"][2] - prefit[param]["prefit"][0]) / 2.
+                            symm_postfit = (pres["fit"][2] - pres["fit"][0]) / 2.
+                            red_factor = symm_postfit / symm_prefit
+                            imp_hi = ((paramScanRes[param][p][2] - paramScanRes[param][p][1]) * red_factor) + paramScanRes[param][p][1]
+                            imp_lo = ((paramScanRes[param][p][0] - paramScanRes[param][p][1]) * red_factor) + paramScanRes[param][p][1]
+                            pres.update({f"global_{p}": [imp_lo, paramScanRes[param][p][1], imp_hi]})
+                        else:
+                            pres.update({f"global_{p}": [globalFitLoRes["fixedpoint"][p], paramScanRes[param][p][1], globalFitHiRes["fixedpoint"][p]]})
             res["params"].append(pres)
         self.flush_queue()
 

python/tool_base/utils.py

@@ -38,25 +38,33 @@ def list_from_workspace(file, workspace, set):
     return res
 
 
-def prefit_from_workspace(file, workspace, params, setPars=None):
+def prefit_from_workspace(file, workspace, params, setPars=None, constPars=None):
     """Given a list of params, return a dictionary of [-1sig, nominal, +1sig]"""
     res = {}
     wsFile = ROOT.TFile(file)
     ws = wsFile.Get(workspace)
     ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
+    updatePars = dict()
     if setPars is not None:
         parsToSet = [tuple(x.split("=")) for x in setPars.split(",")]
-        allParams = ws.allVars()
-        allParams.add(ws.allCats())
         for par, val in parsToSet:
-            tmp = allParams.find(par)
-            isrvar = tmp.IsA().InheritsFrom(ROOT.RooRealVar.Class())
-            if isrvar:
+            updatePars[par] = (float(val), True)
+    if constPars is not None:
+        for par, val in constPars.items():
+            updatePars[par] = (float(val), False)
+    allParams = ws.allVars()
+    allParams.add(ws.allCats())
+    for par, (val, verbose) in updatePars.items():
+        tmp = allParams.find(par)
+        isrvar = tmp.IsA().InheritsFrom(ROOT.RooRealVar.Class())
+        if isrvar:
+            if verbose:
                 print(f"Setting parameter {par} to {float(val):g}")
-                tmp.setVal(float(val))
-            else:
+            tmp.setVal(float(val))
+        else:
+            if verbose:
                 print(f"Setting index {par} to {float(val):g}")
-                tmp.setIndex(int(val))
+            tmp.setIndex(int(val))
 
     for p in params:
         res[p] = {}
@@ -87,6 +95,29 @@ def prefit_from_workspace(file, workspace, params, setPars=None):
             errlo = -1 * var.getErrorLo()
             errhi = +1 * var.getErrorHi()
             res[p]["prefit"] = [val - errlo, val, val + errhi]
+
+            # For non-Gaussian, the best fit and uncertainties of the gobs (given x),
+            # may not be the same as the best fit and uncertainties on x.
+            # Let's calculate these here in case we want them later
+            var.setConstant(True)
+            gobs.setConstant(False)
+            nll2 = ROOT.RooConstraintSum("NLL", "", ROOT.RooArgSet(pdf), ROOT.RooArgSet(gobs))
+            minim = ROOT.RooMinimizer(nll2)
+            minim.setEps(0.001)  # Might as well get some better precision...
+            minim.setErrorLevel(0.5)  # Unlike for a RooNLLVar we must set this explicitly
+            minim.setPrintLevel(-1)
+            minim.setVerbose(False)
+            # Run the fit then run minos for the error
+            minim.minimize("Minuit2", "migrad")
+            minim.minos(ROOT.RooArgSet(gobs))
+            # Should really have checked that these converged ok...
+            # var.Print()
+            # pdf.Print()
+            val = gobs.getVal()
+            errlo = -1 * gobs.getErrorLo()
+            errhi = +1 * gobs.getErrorHi()
+            res[p]["globalobs"] = [val - errlo, val, val + errhi]
+
             if pdf.IsA().InheritsFrom(ROOT.RooGaussian.Class()):
                 res[p]["type"] = "Gaussian"
             elif pdf.IsA().InheritsFrom(ROOT.RooPoisson.Class()):
@@ -183,3 +214,13 @@ def get_fixed_results(file, params):
     res["deltaNLL"] = getattr(t, "deltaNLL")
     res["pvalue"] = getattr(t, "quantileExpected")
     return res
+
+def get_rfr_constvars(filename, rfr_name):
+    f = ROOT.TFile(filename)
+    rfr = f.Get(rfr_name)
+    res = dict()
+    constpars = rfr.constPars()
+    for v in constpars:
+        res[v.GetName()] = v.getVal()
+    f.Close()
+    return res

scripts/groupGlobalImpacts.py

@@ -0,0 +1,92 @@
+#!/usr/bin/env python3
+import math
+import json
+import argparse
+
+def IsConstrained(param_info):
+    return param_info["type"] != "Unconstrained"
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--input", "-i", help="input json file")
+args = parser.parse_args()
+
+# Load the json output of combineTool.py -M Impacts
+data = {}
+with open(args.input) as jsonfile:
+    data = json.load(jsonfile)
+
+# # Set the global plotting style
+# plot.ModTDRStyle(l=args.left_margin, b=0.10, width=(900 if args.checkboxes else 700), height=args.height)
+
+# # We will assume the first POI is the one to plot
+POIs = [ele["name"] for ele in data["POIs"]]
+POIdata = {ele["name"]: ele for ele in data["POIs"]}
+proto = {POI: 0. for POI in POIs}
+fit_val = {}
+err_tot_hi = {}
+err_tot_lo = {}
+
+# Get the list of groups
+groups = set(["TOTAL"])
+for ele in data["params"]:
+    groups.update(ele["groups"])
+
+err_hi = {grp : dict(proto) for grp in groups}
+err_lo = {grp : dict(proto) for grp in groups}
+
+
+for POI in POIs:
+    fit_val[POI] = POIdata[POI]["fit"][1]
+    err_tot_hi[POI] = POIdata[POI]["fit"][2] - POIdata[POI]["fit"][1]
+    err_tot_lo[POI] = POIdata[POI]["fit"][1] - POIdata[POI]["fit"][0]
+
+
+for ele in data["params"]:
+    if IsConstrained(ele):
+        name = ele["name"]
+        # We expect to find "global_[POI]"
+        for POI in POIs:
+            impact = ele[f"global_{POI}"]
+            impact_hi = impact[2] - impact[1] # impact from shifting global obs up
+            impact_lo = impact[0] - impact[1] # impact from shifting global obs down
+            # Check different situations
+            contrib_hi = 0.
+            contrib_lo = 0.
+            if impact_hi >= 0. and impact_lo <= 0.:
+                contrib_hi = impact_hi
+                contrib_lo = impact_lo
+            elif impact_hi <=0. and impact_lo >= 0.:
+                contrib_lo = impact_hi
+                contrib_hi = impact_lo
+            elif impact_hi >= 0. and impact_lo >= 0.:
+                print(f"Warning, parameter {name} has global impact on {POI} that are both positive: ({impact_hi},{impact_lo}), we will take the max of the two")
+                contrib_hi = max(impact_hi, impact_lo)
+            elif impact_lo <= 0. and impact_lo <= 0.:
+                print(f"Warning, parameter {name} has global impact on {POI} that are both negative: ({impact_hi},{impact_lo}), we will take the min of the two")
+
+                contrib_lo = min(impact_hi, impact_lo)
+            for grp in ["TOTAL"] + ele["groups"]:
+                err_hi[grp][POI] += pow(contrib_hi, 2)
+                err_lo[grp][POI] += pow(contrib_lo, 2)
+            # print(ele["name"], POI, impact_hi, impact_lo, contrib_hi, contrib_lo)
+
+err_hi["STAT"] = dict(proto)
+err_lo["STAT"] = dict(proto)
+for POI in POIs:
+    err_hi["STAT"][POI] = math.sqrt(pow(err_tot_hi[POI], 2) - err_hi["TOTAL"][POI])
+    err_lo["STAT"][POI] = math.sqrt(pow(err_tot_lo[POI], 2) - err_lo["TOTAL"][POI])
+for grp in groups:
+    for POI in POIs:
+        err_hi[grp][POI] = math.sqrt(err_hi[grp][POI])
+        err_lo[grp][POI] = math.sqrt(err_lo[grp][POI])
+
+for POI in POIs:
+    print(f'POI: {POI:<20} {"Best-fit":>20} : {fit_val[POI]:<10.3f}')
+    print(f'POI: {POI:<20} {"Total uncertainty":>20} : +{err_tot_hi[POI]:<10.3f} -{err_tot_lo[POI]:<10.3f}')
+    print(f'POI: {POI:<20} {"Stat":>20} : +{err_hi["STAT"][POI]:<10.3f} -{err_lo["STAT"][POI]:<10.3f}')
+    print(f'POI: {POI:<20} {"Syst":>20} : +{err_hi["TOTAL"][POI]:<10.3f} -{err_lo["TOTAL"][POI]:<10.3f}')
+    for grp in groups:
+        if grp in ["TOTAL", "STAT"]:
+            continue
+        print(f'POI: {POI:<20} {grp:>20} : +{err_hi[grp][POI]:<10.3f} -{err_lo[grp][POI]:<10.3f}')
+