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,31 @@ 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:
-                print(f"Setting parameter {par} to {float(val):g}")
-                tmp.setVal(float(val))
-            else:
-                print(f"Setting index {par} to {float(val):g}")
-                tmp.setIndex(int(val))
+            updatePars[par] = float(val)
+    if constPars is not None:
+        updatePars.update(constPars)
+    allParams = ws.allVars()
+    allParams.add(ws.allCats())
+    print(updatePars)
+    for par, val in updatePars.items():
+        tmp = allParams.find(par)
+        isrvar = tmp.IsA().InheritsFrom(ROOT.RooRealVar.Class())
+        if isrvar:
+            print(f"Setting parameter {par} to {float(val):g}")
+            tmp.setVal(float(val))
+        else:
+            print(f"Setting index {par} to {float(val):g}")
+            tmp.setIndex(int(val))
 
     for p in params:
         res[p] = {}
@@ -87,6 +93,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 +212,12 @@ 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()
+    return res

scripts/plotImpacts.py

@@ -50,6 +50,7 @@ def IsConstrained(param_info):
 parser.add_argument("--POI", default=None, help="Specify a POI to draw")
 parser.add_argument("--sort", "-s", choices=["impact", "constraint", "pull"], default="impact", help="The metric to sort the list of parameters")
 parser.add_argument("--relative", "-r", action="store_true", help="Show impacts relative to the uncertainty on the POI")
+parser.add_argument("--show-global", "-g", action="store_true", help="Additionally show the global impacts (if ran in the Impacts tool)")
 parser.add_argument("--summary", action="store_true", help="Produce additional summary page, named [output]_summary.pdf")
 args = parser.parse_args()
 
@@ -82,25 +83,50 @@ def IsConstrained(param_info):
 
 POI_fit = POI_info["fit"]
 
-# Pre-compute info for each parameter
-params = data["params"]
-for ele in params:
-    # Calculate impacts and relative impacts. Note that here the impacts are signed.
-    ele["impact_hi"] = ele[POI][2] - ele[POI][1]
-    ele["impact_lo"] = ele[POI][0] - ele[POI][1]
-    # Some care needed with the relative ones, since we don't know the signs of hi and lo.
-    # We want to divide any positive impact by the positive uncert. on the POI, and similar for negative.
-    # We also need to be careful in case the uncertainties on the POI came out as zero (shouldn't happen...)
+def ComputeImpact(ele, POI, POI_fit, prefix="", checkRelative=False):
+    entry = f"{prefix}{POI}"
+    if entry not in ele:
+        return
+    ele[f"{prefix}impact_hi"] = ele[entry][2] - ele[entry][1]
+    ele[f"{prefix}impact_lo"] = ele[entry][0] - ele[entry][1]
     if (POI_fit[2] - POI_fit[1]) > 0.0 and (POI_fit[1] - POI_fit[0]) > 0.0:
-        ele["impact_rel_hi"] = ele["impact_hi"] / ((POI_fit[2] - POI_fit[1]) if ele["impact_hi"] >= 0 else (POI_fit[1] - POI_fit[0]))
-        ele["impact_rel_lo"] = ele["impact_lo"] / ((POI_fit[2] - POI_fit[1]) if ele["impact_lo"] >= 0 else (POI_fit[1] - POI_fit[0]))
+        ele[f"{prefix}impact_rel_hi"] = ele[f"{prefix}impact_hi"] / ((POI_fit[2] - POI_fit[1]) if ele[f"{prefix}impact_hi"] >= 0 else (POI_fit[1] - POI_fit[0]))
+        ele[f"{prefix}impact_rel_lo"] = ele[f"{prefix}impact_lo"] / ((POI_fit[2] - POI_fit[1]) if ele[f"{prefix}impact_lo"] >= 0 else (POI_fit[1] - POI_fit[0]))
     else:
-        ele["impact_rel_hi"] = 0.0
-        ele["impact_rel_lo"] = 0.0
-        if args.relative:
+        ele[f"{prefix}impact_rel_hi"] = 0.0
+        ele[f"{prefix}impact_rel_lo"] = 0.0
+        if checkRelative:
             # Now we have a real problem, best throw:
             raise RuntimeError("Relative impacts requested (--relative), but uncertainty on the POI is zero")
 
+
+# Pre-compute info for each parameter
+params = data["params"]
+for ele in params:
+    ComputeImpact(ele, POI, POI_fit, prefix="", checkRelative=args.relative)
+    if args.show_global:
+        ComputeImpact(ele, POI, POI_fit, prefix="global_", checkRelative=args.relative)
+    # # Calculate impacts and relative impacts. Note that here the impacts are signed.
+    # ele["impact_hi"] = ele[POI][2] - ele[POI][1]
+    # ele["impact_lo"] = ele[POI][0] - ele[POI][1]
+    # g_entry = f"global_{POI}"
+    # if g_entry in ele:
+    #     # The global impact may or may not be available
+    #     ele["global_impact_hi"] = ele[g_entry][2] - ele[g_entry][1]
+    #     ele["global_impact_lo"] = ele[g_entry][0] - ele[g_entry][1]
+    # # Some care needed with the relative ones, since we don't know the signs of hi and lo.
+    # # We want to divide any positive impact by the positive uncert. on the POI, and similar for negative.
+    # # We also need to be careful in case the uncertainties on the POI came out as zero (shouldn't happen...)
+    # if (POI_fit[2] - POI_fit[1]) > 0.0 and (POI_fit[1] - POI_fit[0]) > 0.0:
+    #     ele["impact_rel_hi"] = ele["impact_hi"] / ((POI_fit[2] - POI_fit[1]) if ele["impact_hi"] >= 0 else (POI_fit[1] - POI_fit[0]))
+    #     ele["impact_rel_lo"] = ele["impact_lo"] / ((POI_fit[2] - POI_fit[1]) if ele["impact_lo"] >= 0 else (POI_fit[1] - POI_fit[0]))
+    # else:
+    #     ele["impact_rel_hi"] = 0.0
+    #     ele["impact_rel_lo"] = 0.0
+    #     if args.relative:
+    #         # Now we have a real problem, best throw:
+    #         raise RuntimeError("Relative impacts requested (--relative), but uncertainty on the POI is zero")
+
     if IsConstrained(ele):
         pre = ele["prefit"]
         fit = ele["fit"]
@@ -362,7 +388,6 @@ def DrawBoxes(color):
 
     canv.Print(f"{args.output}_summary.pdf")
 
-
 if args.summary:
     MakeSummaryPage()
 
@@ -406,6 +431,8 @@ def DrawBoxes(color):
     g_pull = ROOT.TGraph(n_params)
     g_impacts_hi = ROOT.TGraphAsymmErrors(n_params)
     g_impacts_lo = ROOT.TGraphAsymmErrors(n_params)
+    g_glob_impacts_hi = ROOT.TGraphAsymmErrors(n_params)
+    g_glob_impacts_lo = ROOT.TGraphAsymmErrors(n_params)
     g_check = ROOT.TGraphAsymmErrors()
     g_check_i = 0
 
@@ -457,6 +484,8 @@ def DrawBoxes(color):
             redo_boxes.append(i)
         g_impacts_hi.SetPoint(i, 0, float(i) + 0.5)
         g_impacts_lo.SetPoint(i, 0, float(i) + 0.5)
+        g_glob_impacts_hi.SetPoint(i, 0, float(i) + 0.5)
+        g_glob_impacts_lo.SetPoint(i, 0, float(i) + 0.5)
         if args.checkboxes:
             pboxes = pdata[p]["checkboxes"]
             for pbox in pboxes:
@@ -466,6 +495,9 @@ def DrawBoxes(color):
         imp = pdata[p][POI]
         g_impacts_hi.SetPointError(i, 0, par[impt_prefix + "_hi"], 0.5, 0.5)
         g_impacts_lo.SetPointError(i, -1.0 * par[impt_prefix + "_lo"], 0, 0.5, 0.5)
+        if "global_impact_hi" in par:
+            g_glob_impacts_hi.SetPointError(i, 0, par[f"global_{impt_prefix}_hi"], 0.30, 0.30)
+            g_glob_impacts_lo.SetPointError(i, -1.0 * par[f"global_{impt_prefix}_lo"], 0, 0.30, 0.30)
         max_impact = max(max_impact, abs(par[impt_prefix + "_hi"]), abs(par[impt_prefix + "_lo"]))
         col = colors.get(tp, 2)
         if args.color_groups is not None and len(pdata[p]["groups"]) >= 1:
@@ -544,7 +576,7 @@ def DrawBoxes(color):
     g_pull.Draw("PSAME")
     # And back to the second pad to draw the impacts graphs
     pads[1].cd()
-    alpha = 0.7
+    alpha = 0.7 if not args.show_global else 0.3
 
     lo_color = {"default": 38, "hesse": ROOT.kOrange - 3, "robust": ROOT.kGreen + 1}
     hi_color = {"default": 46, "hesse": ROOT.kBlue, "robust": ROOT.kAzure - 5}
@@ -555,14 +587,27 @@ def DrawBoxes(color):
     g_impacts_hi.Draw("2SAME")
     g_impacts_lo.SetFillColor(plot.CreateTransparentColor(lo_color[method], alpha))
     g_impacts_lo.Draw("2SAME")
+    g_glob_impacts_hi.SetLineColor(hi_color[method])
+    g_glob_impacts_hi.SetFillStyle(0)
+    g_glob_impacts_hi.SetLineWidth(2)
+    g_glob_impacts_lo.SetLineColor(lo_color[method])
+    g_glob_impacts_lo.SetLineWidth(2)
+    g_glob_impacts_lo.SetFillStyle(0)
+    g_glob_impacts_hi.Draw("5SAME")
+    g_glob_impacts_lo.Draw("5SAME")
     pads[1].RedrawAxis()
 
+    showGlobalImpacts = True
     legend = ROOT.TLegend(0.02, 0.02, 0.35, 0.09, "", "NBNDC")
-    legend.SetNColumns(2)
+    legend.SetNColumns(3 if showGlobalImpacts else 2)
     legend.AddEntry(g_fit, "Fit", "LP")
     legend.AddEntry(g_impacts_hi, "+1#sigma Impact", "F")
+    if showGlobalImpacts:
+        legend.AddEntry(g_glob_impacts_hi, "+1#sigma Global", "F")
     legend.AddEntry(g_pull, "Pull", "P")
     legend.AddEntry(g_impacts_lo, "-1#sigma Impact", "F")
+    if showGlobalImpacts:
+        legend.AddEntry(g_glob_impacts_lo, "-1#sigma Global", "F")
     legend.Draw()
 
     leg_width = pads[0].GetLeftMargin() - 0.01