From bfc5b437bc26d51862e9530e9c2de7cb61444316 Mon Sep 17 00:00:00 2001 From: Andrew Gilbert Date: Mon, 12 Jan 2026 14:52:01 +0100 Subject: [PATCH 1/3] First implementation of global impacts - Added option to combineTool.py -M Impacts to calculate the global impacts both for the standard (asymmetric) impacts (what ATLAS call the shifted global-observable method), and the approximate (symmetric) impacts (what ATLAS call the post-fit nuisance-parameter covariance method) - Enabled by adding `--globalImpacts` or `-g` to all three of the combineTool.py steps. At the moment, runs in addition to (not instead of) the standard impacts calculation. Anticipate support run only global impacts eventually. - plotImpacts.py also has an option `--show-global` or `-g` to draw the global impacts too. --- python/tool_base/Impacts.py | 52 +++++++++++++++++++++--- python/tool_base/utils.py | 60 ++++++++++++++++++++++------ scripts/plotImpacts.py | 79 +++++++++++++++++++++++++++++-------- src/MultiDimFit.cc | 2 +- 4 files changed, 159 insertions(+), 34 deletions(-) diff --git a/python/tool_base/Impacts.py b/python/tool_base/Impacts.py index 51cbf26aefc..798ab72c9b8 100755 --- a/python/tool_base/Impacts.py +++ b/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,6 +162,8 @@ 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") @@ -161,13 +171,18 @@ def run_method(self): 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() diff --git a/python/tool_base/utils.py b/python/tool_base/utils.py index b5885c76db6..d0aa7924bde 100644 --- a/python/tool_base/utils.py +++ b/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 \ No newline at end of file diff --git a/scripts/plotImpacts.py b/scripts/plotImpacts.py index b2cb280f2b6..514b4baa17e 100755 --- a/scripts/plotImpacts.py +++ b/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 diff --git a/src/MultiDimFit.cc b/src/MultiDimFit.cc index 395bdf14100..3b4fdf8cf6b 100644 --- a/src/MultiDimFit.cc +++ b/src/MultiDimFit.cc @@ -197,7 +197,7 @@ bool MultiDimFit::runSpecific(RooWorkspace *w, RooStats::ModelConfig *mc_s, RooS const RooCmdArg &constrainCmdArg = withSystematics ? RooFit::Constrain(*mc_s->GetNuisanceParameters()) : RooCmdArg(); std::unique_ptr res; if (verbose <= 3) RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); - bool doHesse = (algo_ == Singles || algo_ == Impact) || (saveFitResult_) ; + bool doHesse = (algo_ == Singles || algo_ == Impact) || (saveFitResult_); if ( !skipInitialFit_){ std::cout << "Doing initial fit: " << std::endl; res.reset(doFit(pdf, data, (doHesse ? poiList_ : RooArgList()), constrainCmdArg, (saveFitResult_ && !robustHesse_), 1, true, false)); From 623fed41cc05ee58cb96b57c076852fba43c077b Mon Sep 17 00:00:00 2001 From: Andrew Gilbert Date: Mon, 12 Jan 2026 15:49:54 +0100 Subject: [PATCH 2/3] Fix a couple of things found by the CI --- python/tool_base/utils.py | 1 + scripts/plotImpacts.py | 7 +++---- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/tool_base/utils.py b/python/tool_base/utils.py index d0aa7924bde..8624931da84 100644 --- a/python/tool_base/utils.py +++ b/python/tool_base/utils.py @@ -220,4 +220,5 @@ def get_rfr_constvars(filename, rfr_name): constpars = rfr.constPars() for v in constpars: res[v.GetName()] = v.getVal() + f.Close() return res \ No newline at end of file diff --git a/scripts/plotImpacts.py b/scripts/plotImpacts.py index 514b4baa17e..ef16c36d226 100755 --- a/scripts/plotImpacts.py +++ b/scripts/plotImpacts.py @@ -597,16 +597,15 @@ def DrawBoxes(color): 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(3 if showGlobalImpacts else 2) + legend.SetNColumns(3 if args.show_global else 2) legend.AddEntry(g_fit, "Fit", "LP") legend.AddEntry(g_impacts_hi, "+1#sigma Impact", "F") - if showGlobalImpacts: + if args.show_global: 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: + if args.show_global: legend.AddEntry(g_glob_impacts_lo, "-1#sigma Global", "F") legend.Draw() From 4e8a57431963b6ab9d038105f5b88450382a96c6 Mon Sep 17 00:00:00 2001 From: Andrew Gilbert Date: Thu, 22 Jan 2026 14:19:41 +0100 Subject: [PATCH 3/3] Adding small script for combining global impacts into groups --- python/tool_base/utils.py | 14 +++--- scripts/groupGlobalImpacts.py | 92 +++++++++++++++++++++++++++++++++++ scripts/plotImpacts.py | 19 ++++---- 3 files changed, 110 insertions(+), 15 deletions(-) create mode 100755 scripts/groupGlobalImpacts.py diff --git a/python/tool_base/utils.py b/python/tool_base/utils.py index 8624931da84..cce9008dc83 100644 --- a/python/tool_base/utils.py +++ b/python/tool_base/utils.py @@ -48,20 +48,22 @@ def prefit_from_workspace(file, workspace, params, setPars=None, constPars=None) if setPars is not None: parsToSet = [tuple(x.split("=")) for x in setPars.split(",")] for par, val in parsToSet: - updatePars[par] = float(val) + updatePars[par] = (float(val), True) if constPars is not None: - updatePars.update(constPars) + for par, val in constPars.items(): + updatePars[par] = (float(val), False) allParams = ws.allVars() allParams.add(ws.allCats()) - print(updatePars) - for par, val in updatePars.items(): + for par, (val, verbose) 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}") + if verbose: + print(f"Setting parameter {par} to {float(val):g}") tmp.setVal(float(val)) else: - print(f"Setting index {par} to {float(val):g}") + if verbose: + print(f"Setting index {par} to {float(val):g}") tmp.setIndex(int(val)) for p in params: diff --git a/scripts/groupGlobalImpacts.py b/scripts/groupGlobalImpacts.py new file mode 100755 index 00000000000..fdc9f71dc8b --- /dev/null +++ b/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}') + \ No newline at end of file diff --git a/scripts/plotImpacts.py b/scripts/plotImpacts.py index ef16c36d226..9f0a6b92b2d 100755 --- a/scripts/plotImpacts.py +++ b/scripts/plotImpacts.py @@ -495,7 +495,7 @@ 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: + if f"global_{impt_prefix}_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"])) @@ -587,14 +587,15 @@ 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") + if args.show_global: + 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() legend = ROOT.TLegend(0.02, 0.02, 0.35, 0.09, "", "NBNDC")