From 210fa7bb4fd0ffcda79545bb6774350e866bcd12 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Tue, 21 Mar 2023 12:31:05 +0100 Subject: [PATCH 01/11] First pass at adding nuisance delta nll plot --- data/tutorials/longexercise/diffNuisances.py | 318 ++++++++++++------- 1 file changed, 208 insertions(+), 110 deletions(-) diff --git a/data/tutorials/longexercise/diffNuisances.py b/data/tutorials/longexercise/diffNuisances.py index 8161579db4e..86323da1793 100644 --- a/data/tutorials/longexercise/diffNuisances.py +++ b/data/tutorials/longexercise/diffNuisances.py @@ -111,6 +111,22 @@ type="string", help="Choose the definition of the pull, see python/calculate_pulls.py for options", ) +parser.add_option( + "-w", + "--workspace", + dest="workspace", + default="", + type="string", + help="Workspace to use for evaluating NLL differences." +) +parser.add_option( + "", + "--max-nuis", + dest="max_nuis", + default=65, + type="int", + help="Maximum nuisances for a single plot" +) (options, args) = parser.parse_args() if len(args) == 0: @@ -133,6 +149,11 @@ ) + str(options) file = ROOT.TFile(args[0]) +workspace = None +if options.workspace: + workspace_file = ROOT.TFile( options.workspace, "READ" ) + workspace = workspace_file.Get("w") + if file == None: raise RuntimeError("Cannot open file %s" % args[0]) fit_s = file.Get("fit_s") @@ -155,6 +176,7 @@ fpf_s = fit_s.floatParsFinal() pulls = [] +dnlls = [] nuis_p_i = 0 title = "pull" if options.pullDef else "#theta" @@ -172,25 +194,37 @@ def getGraph(hist,shift): gr.SetPointError(i,float(abs(shift))*0.8,e) return gr """ - -# Also make histograms for pull distributions: -hist_fit_b = ROOT.TH1F("fit_b", "B-only fit Nuisances;;%s " % title, prefit.getSize(), 0, prefit.getSize()) -hist_fit_s = ROOT.TH1F("fit_s", "S+B fit Nuisances ;;%s " % title, prefit.getSize(), 0, prefit.getSize()) -hist_prefit = ROOT.TH1F( - "prefit_nuisancs", - "Prefit Nuisances ;;%s " % title, - prefit.getSize(), - 0, - prefit.getSize(), -) -# Store also the *asymmetric* uncertainties -gr_fit_b = ROOT.TGraphAsymmErrors() -gr_fit_b.SetTitle("fit_b_g") -gr_fit_s = ROOT.TGraphAsymmErrors() -gr_fit_s.SetTitle("fit_b_s") + +n_hists = prefit.getSize() // options.max_nuis +if prefit.getSize() % options.max_nuis == 0: + n_hists -= 1 + +hist_fit_b = [] +hist_fit_s = [] +hist_prefit = [] +gr_fit_b = [] +gr_fit_s = [] +for idx in range( ((prefit.getSize() - 1) // options.max_nuis) + 1 ): + # Also make histograms for pull distributions: + nbins = min(prefit.getSize(), options.max_nuis, prefit.getSize() - options.max_nuis * idx ) + hist_fit_b += [ ROOT.TH1F("fit_b %s" % idx, "B-only fit Nuisances;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_fit_s += [ ROOT.TH1F("fit_s %s" % idx, "S+B fit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_prefit +=[ ROOT.TH1F( + "prefit_nuisancs %s" % idx, + "Prefit Nuisances ;;%s %s" % (title, idx), + nbins, + 0, + nbins, + ) ] + # Store also the *asymmetric* uncertainties + gr_fit_b += [ ROOT.TGraphAsymmErrors() ] + gr_fit_b[idx].SetTitle("fit_b_g % s" % idx) + gr_fit_s += [ ROOT.TGraphAsymmErrors() ] + gr_fit_s[idx].SetTitle("fit_b_s % s" % idx) # loop over all fitted parameters +dnll = [] for i in range(fpf_s.getSize()): nuis_s = fpf_s.at(i) name = nuis_s.GetName() @@ -259,46 +293,58 @@ def getGraph(hist,shift): if nuis_p != None: if options.plotfile: + pdf = workspace.pdf(name + "_Pdf") + var = workspace.var(name) + var_val = var.getVal() if fit_name == "b": nuis_p_i += 1 + hist_idx = (nuis_p_i - 1) // options.max_nuis + bin_idx = (nuis_p_i - 1) % options.max_nuis + 1 + print(hist_idx, bin_idx) if options.pullDef and nuis_p != None: # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) - gr_fit_b.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 + 0.1, nx) - gr_fit_b.SetPointError(nuis_p_i - 1, 0, 0, ned, neu) + gr_fit_b[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 + 0.1, nx) + gr_fit_b[hist_idx].SetPointError(bin_idx - 1, 0, 0, ned, neu) else: - gr_fit_b.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 + 0.1, nuis_x.getVal()) - gr_fit_b.SetPointError( - nuis_p_i - 1, + gr_fit_b[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 + 0.1, nuis_x.getVal()) + gr_fit_b[hist_idx].SetPointError( + bin_idx - 1, 0, 0, abs(nuis_x.getErrorLo()), nuis_x.getErrorHi(), ) - hist_fit_b.SetBinContent(nuis_p_i, nuis_x.getVal()) - hist_fit_b.SetBinError(nuis_p_i, nuis_x.getError()) - hist_fit_b.GetXaxis().SetBinLabel(nuis_p_i, name) - gr_fit_b.GetXaxis().SetBinLabel(nuis_p_i, name) + hist_fit_b[hist_idx].SetBinContent(bin_idx, nuis_x.getVal()) + hist_fit_b[hist_idx].SetBinError(bin_idx, nuis_x.getError()) + hist_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + gr_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + var.setVal(nuis_x.getVal()) + bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) if fit_name == "s": if options.pullDef and nuis_p != None: # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) - gr_fit_s.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 - 0.1, nx) - gr_fit_s.SetPointError(nuis_p_i - 1, 0, 0, ned, neu) + gr_fit_s[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 - 0.1, nx) + gr_fit_s[hist_idx].SetPointError(bin_idx - 1, 0, 0, ned, neu) else: - gr_fit_s.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 - 0.1, nuis_x.getVal()) - gr_fit_s.SetPointError( - nuis_p_i - 1, + gr_fit_s[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 - 0.1, nuis_x.getVal()) + gr_fit_s[hist_idx].SetPointError( + bin_idx - 1, 0, 0, abs(nuis_x.getErrorLo()), nuis_x.getErrorHi(), ) - hist_fit_s.SetBinContent(nuis_p_i, nuis_x.getVal()) - hist_fit_s.SetBinError(nuis_p_i, nuis_x.getError()) - hist_fit_s.GetXaxis().SetBinLabel(nuis_p_i, name) - gr_fit_s.GetXaxis().SetBinLabel(nuis_p_i, name) - hist_prefit.SetBinContent(nuis_p_i, mean_p) - hist_prefit.SetBinError(nuis_p_i, sigma_p) - hist_prefit.GetXaxis().SetBinLabel(nuis_p_i, name) + hist_fit_s[hist_idx].SetBinContent(bin_idx, nuis_x.getVal()) + hist_fit_s[hist_idx].SetBinError(bin_idx, nuis_x.getError()) + hist_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + gr_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + var.setVal(nuis_x.getVal()) + sfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) + hist_prefit[hist_idx].SetBinContent(bin_idx, mean_p) + hist_prefit[hist_idx].SetBinError(bin_idx, sigma_p) + hist_prefit[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) if sigma_p > 0: if options.pullDef: @@ -353,12 +399,32 @@ def getGraph(hist,shift): elif options.show_all_parameters: flag = True + dnll.append( (name, bfit_nll - sfit_nll) ) # end of loop over s and b row += ["%+4.2f" % fit_s.correlation(name, options.poi)] + row += ["%.4f" % (bfit_nll - sfit_nll) ] if flag or options.show_all_parameters: table[name] = row +len_dnll = len(dnll) +hist_dnll = [] +hist_cdnll = [] +for hist_idx in range( ((len_dnll - 1) // options.max_nuis) + 1 ): + nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx ) + hist_dnll += [ ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] + hist_cdnll += [ ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] +dnll.sort( key= lambda x: -x[1] ) +cdnll = 0 +for idx, (nm, val) in enumerate(dnll): + hist_idx = idx // options.max_nuis + bin_idx = idx % options.max_nuis + hist_dnll[hist_idx].SetBinContent(bin_idx+1, val) + hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx+1, nm) + cdnll += val + hist_cdnll[hist_idx].SetBinContent( bin_idx+1, cdnll ) +# fmt_string = "%-40s %.3f" +# print(fmt_string % (nm, val ) ) # end of loop over all fitted parameters # ---------- @@ -369,7 +435,7 @@ def getGraph(hist,shift): print(setUpString) print() -fmtstring = "%-40s %15s %15s %10s" +fmtstring = "%-40s %15s %15s %10s %10s" highlight = "*%s*" morelight = "!%s!" pmsub, sigsub = None, None @@ -381,7 +447,7 @@ def getGraph(hist,shift): fmtstring = "%-40s %15s %30s %30s %10s" print(fmtstring % ("name", "pre fit", "b-only fit", "s+b fit", "rho")) else: - print(fmtstring % ("name", "b-only fit", "s+b fit", "rho")) + print(fmtstring % ("name", "b-only fit", "s+b fit", "rho","dnll")) elif options.format == "latex": pmsub = (r"(\S+) \+/- (\S+)", r"$\1 \\pm \2$") sigsub = ("sig", r"$\\sigma$") @@ -479,13 +545,13 @@ def getGraph(hist,shift): if sigsub != None: v = [re.sub(sigsub[0], sigsub[1], i) for i in v] if (n, "b") in isFlagged: - v[-3] = highlighters[isFlagged[(n, "b")]] % v[-3] + v[-4] = highlighters[isFlagged[(n, "b")]] % v[-4] if (n, "s") in isFlagged: - v[-2] = highlighters[isFlagged[(n, "s")]] % v[-2] + v[-3] = highlighters[isFlagged[(n, "s")]] % v[-3] if options.absolute_values: print(fmtstring % (n, v[0], v[1], v[2], v[3])) else: - print(fmtstring % (n, v[0], v[1], v[2])) + print(fmtstring % (n, v[0], v[1], v[2], v[3])) if options.format == "latex": print(" \\hline\n\\end{tabular}") @@ -499,10 +565,11 @@ def getGraph(hist,shift): fout = ROOT.TFile(options.plotfile, "RECREATE") ROOT.gROOT.SetStyle("Plain") ROOT.gStyle.SetOptFit(1) + ROOT.gStyle.SetPadBottomMargin(0.3) histogram = ROOT.TH1F("pulls", "Pulls", 60, -3, 3) for pull in pulls: histogram.Fill(pull) - canvas = ROOT.TCanvas("asdf", "asdf", 800, 800) + canvas = ROOT.TCanvas("asdf","asdf", 800, 800) if options.pullDef: histogram.GetXaxis().SetTitle("pull") else: @@ -513,70 +580,101 @@ def getGraph(hist,shift): histogram.Draw("pe") fout.WriteTObject(canvas) - canvas_nuis = ROOT.TCanvas("nuisances", "nuisances", 900, 600) - hist_fit_e_s = hist_fit_s.Clone("errors_s") - hist_fit_e_b = hist_fit_b.Clone("errors_b") - # gr_fit_s = getGraph(hist_fit_s,-0.1) - # gr_fit_b = getGraph(hist_fit_b, 0.1) - gr_fit_s.SetLineColor(ROOT.kRed) - gr_fit_s.SetMarkerColor(ROOT.kRed) - gr_fit_b.SetLineColor(ROOT.kBlue) - gr_fit_b.SetMarkerColor(ROOT.kBlue) - gr_fit_b.SetMarkerStyle(20) - gr_fit_s.SetMarkerStyle(20) - gr_fit_b.SetMarkerSize(1.0) - gr_fit_s.SetMarkerSize(1.0) - gr_fit_b.SetLineWidth(2) - gr_fit_s.SetLineWidth(2) - hist_prefit.SetLineWidth(2) - hist_prefit.SetTitle("Nuisance Paramaeters") - hist_prefit.SetLineColor(ROOT.kBlack) - hist_prefit.SetFillColor(ROOT.kGray) - hist_prefit.SetMaximum(3) - hist_prefit.SetMinimum(-3) - hist_prefit.Draw("E2") - hist_prefit.Draw("histsame") - gr_fit_b.Draw("EPsame") - gr_fit_s.Draw("EPsame") - canvas_nuis.SetGridx() - canvas_nuis.RedrawAxis() - canvas_nuis.RedrawAxis("g") - leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg.SetFillColor(0) - leg.SetTextFont(42) - leg.AddEntry(hist_prefit, "Prefit", "FL") - leg.AddEntry(gr_fit_b, "B-only fit", "EPL") - leg.AddEntry(gr_fit_s, "S+B fit", "EPL") - leg.Draw() - fout.WriteTObject(canvas_nuis) - canvas_pferrs = ROOT.TCanvas("post_fit_errs", "post_fit_errs", 900, 600) - for b in range(1, hist_fit_e_s.GetNbinsX() + 1): - hist_fit_e_s.SetBinContent(b, hist_fit_s.GetBinError(b) / hist_prefit.GetBinError(b)) - hist_fit_e_b.SetBinContent(b, hist_fit_b.GetBinError(b) / hist_prefit.GetBinError(b)) - hist_fit_e_s.SetBinError(b, 0) - hist_fit_e_b.SetBinError(b, 0) - hist_fit_e_s.SetFillColor(ROOT.kRed) - hist_fit_e_b.SetFillColor(ROOT.kBlue) - hist_fit_e_s.SetBarWidth(0.4) - hist_fit_e_b.SetBarWidth(0.4) - hist_fit_e_b.SetBarOffset(0.45) - hist_fit_e_b.GetYaxis().SetTitle("#sigma_{#theta}/(#sigma_{#theta} prefit)") - hist_fit_e_b.SetTitle("Nuisance Parameter Uncertainty Reduction") - hist_fit_e_b.SetMaximum(1.5) - hist_fit_e_b.SetMinimum(0) - hist_fit_e_b.Draw("bar") - hist_fit_e_s.Draw("barsame") - leg_rat = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg_rat.SetFillColor(0) - leg_rat.SetTextFont(42) - leg_rat.AddEntry(hist_fit_e_b, "B-only fit", "F") - leg_rat.AddEntry(hist_fit_e_s, "S+B fit", "F") - leg_rat.Draw() - line_one = ROOT.TLine(0, 1, hist_fit_e_s.GetXaxis().GetXmax(), 1) - line_one.SetLineColor(1) - line_one.SetLineStyle(2) - line_one.SetLineWidth(2) - line_one.Draw() - canvas_pferrs.RedrawAxis() - - fout.WriteTObject(canvas_pferrs) + for idx in range(len(hist_dnll)): + canvas_nuis = ROOT.TCanvas("nuisances_%s" % idx, "nuisances_%s" % idx, 900, 600) + hist_fit_e_s = hist_fit_s[idx].Clone("errors_s %s" % idx) + hist_fit_e_b = hist_fit_b[idx].Clone("errors_b %s" % idx) + # gr_fit_s = getGraph(hist_fit_s,-0.1) + # gr_fit_b = getGraph(hist_fit_b, 0.1) + gr_fit_s[idx].SetLineColor(ROOT.kRed) + gr_fit_s[idx].SetMarkerColor(ROOT.kRed) + gr_fit_b[idx].SetLineColor(ROOT.kBlue) + gr_fit_b[idx].SetMarkerColor(ROOT.kBlue) + gr_fit_b[idx].SetMarkerStyle(20) + gr_fit_s[idx].SetMarkerStyle(20) + gr_fit_b[idx].SetMarkerSize(1.0) + gr_fit_s[idx].SetMarkerSize(1.0) + gr_fit_b[idx].SetLineWidth(2) + gr_fit_s[idx].SetLineWidth(2) + hist_prefit[idx].SetLineWidth(2) + hist_prefit[idx].SetTitle("Nuisance Paramaeters") + hist_prefit[idx].SetLineColor(ROOT.kBlack) + hist_prefit[idx].SetFillColor(ROOT.kGray) + hist_prefit[idx].SetMaximum(3) + hist_prefit[idx].SetMinimum(-3) + hist_prefit[idx].Draw("E2") + hist_prefit[idx].Draw("histsame") + gr_fit_b[idx].Draw("EPsame") + gr_fit_s[idx].Draw("EPsame") + canvas_nuis.SetGridx() + canvas_nuis.RedrawAxis() + canvas_nuis.RedrawAxis("g") + leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg.SetFillColor(0) + leg.SetTextFont(42) + leg.AddEntry(hist_prefit[idx], "Prefit", "FL") + leg.AddEntry(gr_fit_b[idx], "B-only fit", "EPL") + leg.AddEntry(gr_fit_s[idx], "S+B fit", "EPL") + leg.Draw() + fout.WriteTObject(canvas_nuis) + + canvas_pferrs = ROOT.TCanvas("post_fit_errs_%s" % idx, "post_fit_errs_%s" % idx, 900, 600) + for b in range(1, hist_fit_e_s.GetNbinsX() + 1): + hist_fit_e_s.SetBinContent(b, hist_fit_s[idx].GetBinError(b) / hist_prefit[idx].GetBinError(b)) + hist_fit_e_b.SetBinContent(b, hist_fit_b[idx].GetBinError(b) / hist_prefit[idx].GetBinError(b)) + hist_fit_e_s.SetBinError(b, 0) + hist_fit_e_b.SetBinError(b, 0) + hist_fit_e_s.SetFillColor(ROOT.kRed) + hist_fit_e_b.SetFillColor(ROOT.kBlue) + hist_fit_e_s.SetBarWidth(0.4) + hist_fit_e_b.SetBarWidth(0.4) + hist_fit_e_b.SetBarOffset(0.45) + hist_fit_e_b.GetYaxis().SetTitle("#sigma_{#theta}/(#sigma_{#theta} prefit)") + hist_fit_e_b.SetTitle("Nuisance Parameter Uncertainty Reduction") + hist_fit_e_b.SetMaximum(1.5) + hist_fit_e_b.SetMinimum(0) + hist_fit_e_b.Draw("bar") + hist_fit_e_s.Draw("barsame") + leg_rat = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg_rat.SetFillColor(0) + leg_rat.SetTextFont(42) + leg_rat.AddEntry(hist_fit_e_b, "B-only fit", "F") + leg_rat.AddEntry(hist_fit_e_s, "S+B fit", "F") + leg_rat.Draw() + line_one = ROOT.TLine(0, 1, hist_fit_e_s.GetXaxis().GetXmax(), 1) + line_one.SetLineColor(1) + line_one.SetLineStyle(2) + line_one.SetLineWidth(2) + line_one.Draw() + canvas_pferrs.RedrawAxis() + + fout.WriteTObject(canvas_pferrs) + + + canvas_dnll = ROOT.TCanvas("dnlls_%s" % idx, "dnlls_%s" % idx, 900, 600) + hist_dnll[idx].SetStats(0) + hist_dnll[idx].SetLineWidth(2) + hist_dnll[idx].SetTitle("Nuisance Parameters") + hist_dnll[idx].SetLineColor(ROOT.kBlack) + hist_dnll[idx].SetFillColor(ROOT.kGray) + hist_dnll[idx].SetMaximum(2) + hist_dnll[idx].SetMinimum(-1) + hist_dnll[idx].GetYaxis().SetTitle("NLL(fit_{b}) - NLL(fit_{s+b})") + hist_dnll[idx].Draw() + + hist_cdnll[idx].SetLineWidth(2) + hist_cdnll[idx].SetLineColor(ROOT.kBlue) + hist_cdnll[idx].Draw("same") + canvas_dnll.SetGridx() + canvas_dnll.RedrawAxis() + canvas_dnll.RedrawAxis("g") + leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg.SetFillColor(0) + leg.SetTextFont(42) + leg.AddEntry(hist_dnll[idx], "individual nuisance", "FL") + leg.AddEntry(hist_cdnll[idx], "cumulative", "FL") + leg.Draw() + latex = ROOT.TLatex() + posx = 4 if hist_dnll[idx].GetNbinsX() > 10 else 0.5 + latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll) ) + fout.WriteTObject(canvas_dnll) From 7c91db93e6b4416d7f54b9eae08fe05ff4a0f319 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Tue, 21 Mar 2023 14:53:52 +0100 Subject: [PATCH 02/11] merge tutorial diffNuisances to test diffNuisances --- test/diffNuisances.py | 352 ++++++++++++++++++++++++++++-------------- 1 file changed, 235 insertions(+), 117 deletions(-) diff --git a/test/diffNuisances.py b/test/diffNuisances.py index b667706960d..0549ec14e25 100644 --- a/test/diffNuisances.py +++ b/test/diffNuisances.py @@ -143,6 +143,22 @@ type="string", help="Include only nuisance parameters that passes the following regex filter", ) +parser.add_option( + "-w", + "--workspace", + dest="workspace", + default="", + type="string", + help="Workspace to use for evaluating NLL differences." +) +parser.add_option( + "", + "--max-nuis", + dest="max_nuis", + default=65, + type="int", + help="Maximum nuisances for a single plot." +) (options, args) = parser.parse_args() if len(args) == 0: @@ -184,6 +200,11 @@ if prefit == None or prefit.ClassName() != "RooArgSet": raise RuntimeError("File %s does not contain the prefit nuisances 'nuisances_prefit'" % args[0]) +workspace = None +if options.workspace: + workspace_file = ROOT.TFile( options.workspace, "READ" ) + workspace = workspace_file.Get("w") + isFlagged = {} # maps from nuisance parameter name to the row to be printed in the table @@ -194,6 +215,7 @@ fpf_s = fit_s.floatParsFinal() pulls = [] +dnlls = [] nuis_p_i = 0 title = "pull" if options.pullDef else "#theta" @@ -219,18 +241,45 @@ def getGraph(hist,shift): for i in range(fpf_s.getSize()): nuis_s = fpf_s.at(i) name = nuis_s.GetName() + nuis_p = prefit.find(name) + if nuis_p == None: + if not options.absolute_values and not (options.pullDef == "unconstPullAsym"): + continue if bool(regex_NP_obj.match(name)): np_count += 1 -# Also make histograms for pull distributions: -hist_fit_b = ROOT.TH1F("fit_b", "B-only fit Nuisances;;%s " % title, np_count, 0, np_count) -hist_fit_s = ROOT.TH1F("fit_s", "S+B fit Nuisances ;;%s " % title, np_count, 0, np_count) -hist_prefit = ROOT.TH1F("prefit_nuisancs", "Prefit Nuisances ;;%s " % title, np_count, 0, np_count) -# Store also the *asymmetric* uncertainties -gr_fit_b = ROOT.TGraphAsymmErrors() -gr_fit_b.SetTitle("fit_b_g") -gr_fit_s = ROOT.TGraphAsymmErrors() -gr_fit_s.SetTitle("fit_b_s") + +n_hists = np_count // options.max_nuis +if np_count % options.max_nuis == 0: + n_hists -= 1 + +hist_fit_b = [] +hist_fit_s = [] +hist_prefit = [] +gr_fit_b = [] +gr_fit_s = [] +hist_dnll = [] +hist_cdnll = [] +for idx in range( ((np_count - 1) // options.max_nuis) + 1 ): + # Also make histograms for pull distributions: + nbins = min(np_count, options.max_nuis, np_count - (options.max_nuis * idx) ) + hist_fit_b += [ ROOT.TH1F("fit_b %s" % idx, "B-only fit Nuisances;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_fit_s += [ ROOT.TH1F("fit_s %s" % idx, "S+B fit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_prefit +=[ ROOT.TH1F( + "prefit_nuisancs %s" % idx, + "Prefit Nuisances ;;%s %s" % (title, idx), + nbins, + 0, + nbins, + ) ] + # Store also the *asymmetric* uncertainties + gr_fit_b += [ ROOT.TGraphAsymmErrors() ] + gr_fit_b[idx].SetTitle("fit_b_g % s" % idx) + gr_fit_s += [ ROOT.TGraphAsymmErrors() ] + gr_fit_s[idx].SetTitle("fit_b_s % s" % idx) + + hist_dnll += [ ROOT.TH1F("dnll %s" % idx, "delta log-likelihoods ;;%s %s" % (title, idx), nbins, 0, nbins) ] + hist_cdnll += [ ROOT.TH1F("cdnll %s" % idx, "cumulative delta log-likelihoods ;;%s %s" % (title, idx), nbins, 0, nbins) ] error_poi = fpf_s.find(options.poi).getError() @@ -309,46 +358,60 @@ def getGraph(hist,shift): if nuis_p != None: if options.plotfile: + if options.workspace: + pdf = workspace.pdf(name + "_Pdf") + var = workspace.var(name) + var_val = var.getVal() if fit_name == "b": nuis_p_i += 1 + hist_idx = (nuis_p_i - 1) // options.max_nuis + bin_idx = ((nuis_p_i - 1) % options.max_nuis) + 1 if options.pullDef and nuis_p != None: # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) - gr_fit_b.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 + 0.1, nx) - gr_fit_b.SetPointError(nuis_p_i - 1, 0, 0, ned, neu) + gr_fit_b[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 + 0.1, nx) + gr_fit_b[hist_idx].SetPointError(bin_idx - 1, 0, 0, ned, neu) else: - gr_fit_b.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 + 0.1, nuis_x.getVal()) - gr_fit_b.SetPointError( - nuis_p_i - 1, + gr_fit_b[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 + 0.1, nuis_x.getVal()) + gr_fit_b[hist_idx].SetPointError( + bin_idx - 1, 0, 0, abs(nuis_x.getErrorLo()), nuis_x.getErrorHi(), ) - hist_fit_b.SetBinContent(nuis_p_i, nuis_x.getVal()) - hist_fit_b.SetBinError(nuis_p_i, nuis_x.getError()) - hist_fit_b.GetXaxis().SetBinLabel(nuis_p_i, name) - gr_fit_b.GetXaxis().SetBinLabel(nuis_p_i, name) + hist_fit_b[hist_idx].SetBinContent(bin_idx, nuis_x.getVal()) + hist_fit_b[hist_idx].SetBinError(bin_idx, nuis_x.getError()) + hist_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + gr_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + if options.workspace: + var.setVal(nuis_x.getVal()) + bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) if fit_name == "s": if options.pullDef and nuis_p != None: # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) - gr_fit_s.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 - 0.1, nx) - gr_fit_s.SetPointError(nuis_p_i - 1, 0, 0, ned, neu) + gr_fit_s[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 - 0.1, nx) + gr_fit_s[hist_idx].SetPointError(bin_idx - 1, 0, 0, ned, neu) else: - gr_fit_s.SetPoint(nuis_p_i - 1, nuis_p_i - 0.5 - 0.1, nuis_x.getVal()) - gr_fit_s.SetPointError( - nuis_p_i - 1, + gr_fit_s[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 - 0.1, nuis_x.getVal()) + gr_fit_s[hist_idx].SetPointError( + bin_idx - 1, 0, 0, abs(nuis_x.getErrorLo()), nuis_x.getErrorHi(), ) - hist_fit_s.SetBinContent(nuis_p_i, nuis_x.getVal()) - hist_fit_s.SetBinError(nuis_p_i, nuis_x.getError()) - hist_fit_s.GetXaxis().SetBinLabel(nuis_p_i, name) - gr_fit_s.GetXaxis().SetBinLabel(nuis_p_i, name) - hist_prefit.SetBinContent(nuis_p_i, mean_p) - hist_prefit.SetBinError(nuis_p_i, sigma_p) - hist_prefit.GetXaxis().SetBinLabel(nuis_p_i, name) + hist_fit_s[hist_idx].SetBinContent(bin_idx, nuis_x.getVal()) + hist_fit_s[hist_idx].SetBinError(bin_idx, nuis_x.getError()) + hist_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + gr_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) + if options.workspace: + var.setVal(nuis_x.getVal()) + sfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) + hist_prefit[hist_idx].SetBinContent(bin_idx, mean_p) + hist_prefit[hist_idx].SetBinError(bin_idx, sigma_p) + hist_prefit[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) if sigma_p > 0: if options.pullDef: @@ -370,11 +433,11 @@ def getGraph(hist,shift): sigShift = nuis_x.getError() if options.pullDef: - row[-1] += "" + row += [""] elif options.absolute_values: - row[-1] += " (%+4.2fsig, %4.2f)" % (valShift, sigShift) + row += [" (%+4.2fsig, %4.2f)" % (valShift, sigShift)] else: - row[-1] = " %+4.2f, %4.2f" % (valShift, sigShift) + row += [" %+4.2f, %4.2f" % (valShift, sigShift)] if fit_name == "b": pulls.append(valShift) @@ -406,12 +469,33 @@ def getGraph(hist,shift): # end of loop over s and b row += ["%+4.2f" % fit_s.correlation(name, options.poi)] - row += ["%+4.3f" % (nuis_x.getError() * fit_s.correlation(name, options.poi) * error_poi)] + if options.workspace: + dnll = bfit_nll - sfit_nll + dnlls.append( (name, dnll) ) + row += ["%.4f" % dnll ] if flag or options.show_all_parameters: table[name] = row # end of loop over all fitted parameters +#fill dnll and cumulative dnll +len_dnll = len(dnlls) +hist_dnll = [] +hist_cdnll = [] +for hist_idx in range( ((len_dnll - 1) // options.max_nuis) + 1 ): + nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx ) + hist_dnll += [ ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] + hist_cdnll += [ ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] +dnlls.sort( key= lambda x: -x[1] ) +cdnll = 0 +for idx, (nm, val) in enumerate(dnlls): + hist_idx = idx // options.max_nuis + bin_idx = idx % options.max_nuis + hist_dnll[hist_idx].SetBinContent(bin_idx+1, val) + hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx+1, nm) + cdnll += val + hist_cdnll[hist_idx].SetBinContent( bin_idx+1, cdnll ) + # ---------- # print the results # ---------- @@ -431,12 +515,17 @@ def getGraph(hist,shift): print(" option '--skipFitB' set true. b-only Fit is just a copy of the s+b fit") if options.pullDef: fmtstring = "%-40s %30s %30s %10s %10s" - print(fmtstring % ("name", "b-only fit pull", "s+b fit pull", "rho", "approx impact")) + headers = ["name", "b-only fit pull", "s+b fit pull", "rho", "approx impact"] elif options.absolute_values: fmtstring = "%-40s %15s %30s %30s %10s %10s" - print(fmtstring % ("name", "pre fit", "b-only fit", "s+b fit", "rho", "approx impact")) + headers = ["name", "pre fit", "b-only fit", "s+b fit", "rho", "approx impact"] else: - print(fmtstring % ("name", "b-only fit", "s+b fit", "rho", "approx impact")) + headers = ["name", "b-only fit", "s+b fit", "rho", "approx impact"] + if options.workspace: + fmtstring = fmtstring + " %10s" + headers += ["dnll"] + print( fmtstring % tuple(headers) ) + elif options.format == "latex": pmsub = (r"(\S+) \+/- (\S+)", r"$\1 \\pm \2$") sigsub = ("sig", r"$\\sigma$") @@ -552,6 +641,8 @@ def getGraph(hist,shift): names = [n[1] for n in names] highlighters = {1: highlight, 2: morelight} +fl_b_idx = -4 if not options.workspace else -5 +fl_s_idx = -3 if not options.workspace else -4 for n in names: v = table[n] if pmsub != None: @@ -559,15 +650,16 @@ def getGraph(hist,shift): if sigsub != None: v = [re.sub(sigsub[0], sigsub[1], i) for i in v] if (n, "b") in isFlagged: - v[-3] = highlighters[isFlagged[(n, "b")]] % v[-3] + v[fl_b_idx] = highlighters[isFlagged[(n, "b")]] % v[fl_b_idx] if (n, "s") in isFlagged: - v[-2] = highlighters[isFlagged[(n, "s")]] % v[-2] + v[fl_s_idx] = highlighters[isFlagged[(n, "s")]] % v[fl_s_idx] if options.format == "latex": n = n.replace(r"_", r"\_") - if options.absolute_values: - print(fmtstring % (n, v[0], v[1], v[2], v[3], v[4])) - else: - print(fmtstring % (n, v[0], v[1], v[2], v[3])) + + j = [n] + v + if not options.absolute_values: + j = j[:1] + j[2:] + print(fmtstring % tuple(j)) if options.format == "latex": print(" \\hline\n\\end{tabular}") @@ -581,10 +673,11 @@ def getGraph(hist,shift): fout = ROOT.TFile(options.plotfile, "RECREATE") ROOT.gROOT.SetStyle("Plain") ROOT.gStyle.SetOptFit(1) + ROOT.gStyle.SetPadBottomMargin(0.3) histogram = ROOT.TH1F("pulls", "Pulls", 60, -3, 3) for pull in pulls: histogram.Fill(pull) - canvas = ROOT.TCanvas("asdf", "asdf", 800, 800) + canvas = ROOT.TCanvas("asdf","asdf", 800, 800) if options.pullDef: histogram.GetXaxis().SetTitle("pull") else: @@ -595,76 +688,101 @@ def getGraph(hist,shift): histogram.Draw("pe") fout.WriteTObject(canvas) - canvas_nuis = ROOT.TCanvas("nuisances", "nuisances", 900, 600) - hist_fit_e_s = hist_fit_s.Clone("errors_s") - hist_fit_e_b = hist_fit_b.Clone("errors_b") - # gr_fit_s = getGraph(hist_fit_s,-0.1) - # gr_fit_b = getGraph(hist_fit_b, 0.1) - gr_fit_s.SetLineColor(ROOT.kRed) - gr_fit_s.SetMarkerColor(ROOT.kRed) - gr_fit_b.SetLineColor(ROOT.kBlue) - gr_fit_b.SetMarkerColor(ROOT.kBlue) - gr_fit_b.SetMarkerStyle(20) - gr_fit_s.SetMarkerStyle(20) - gr_fit_b.SetMarkerSize(1.0) - gr_fit_s.SetMarkerSize(1.0) - gr_fit_b.SetLineWidth(2) - gr_fit_s.SetLineWidth(2) - hist_prefit.SetLineWidth(2) - hist_prefit.SetTitle("Nuisance Parameters") - hist_prefit.SetLineColor(ROOT.kBlack) - hist_prefit.SetFillColor(ROOT.kGray) - hist_prefit.SetMaximum(3) - hist_prefit.SetMinimum(-3) - hist_prefit.Draw("E2") - hist_prefit.Draw("histsame") - if not options.skipFitB: - gr_fit_b.Draw("EPsame") - if not options.skipFitS: - gr_fit_s.Draw("EPsame") - canvas_nuis.SetGridx() - canvas_nuis.RedrawAxis() - canvas_nuis.RedrawAxis("g") - leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg.SetFillColor(0) - leg.SetTextFont(42) - leg.AddEntry(hist_prefit, "Prefit", "FL") - if not options.skipFitB: - leg.AddEntry(gr_fit_b, "B-only fit", "EPL") - if not options.skipFitS: - leg.AddEntry(gr_fit_s, "S+B fit", "EPL") - leg.Draw() - fout.WriteTObject(canvas_nuis) - canvas_pferrs = ROOT.TCanvas("post_fit_errs", "post_fit_errs", 900, 600) - for b in range(1, hist_fit_e_s.GetNbinsX() + 1): - if hist_prefit.GetBinError(b) < 0.000001: - continue - hist_fit_e_s.SetBinContent(b, hist_fit_s.GetBinError(b) / hist_prefit.GetBinError(b)) - hist_fit_e_b.SetBinContent(b, hist_fit_b.GetBinError(b) / hist_prefit.GetBinError(b)) - hist_fit_e_s.SetBinError(b, 0) - hist_fit_e_b.SetBinError(b, 0) - hist_fit_e_s.SetFillColor(ROOT.kRed) - hist_fit_e_b.SetFillColor(ROOT.kBlue) - hist_fit_e_s.SetBarWidth(0.4) - hist_fit_e_b.SetBarWidth(0.4) - hist_fit_e_b.SetBarOffset(0.45) - hist_fit_e_b.GetYaxis().SetTitle("#sigma_{#theta}/(#sigma_{#theta} prefit)") - hist_fit_e_b.SetTitle("Nuisance Parameter Uncertainty Reduction") - hist_fit_e_b.SetMaximum(1.5) - hist_fit_e_b.SetMinimum(0) - hist_fit_e_b.Draw("bar") - hist_fit_e_s.Draw("barsame") - leg_rat = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg_rat.SetFillColor(0) - leg_rat.SetTextFont(42) - leg_rat.AddEntry(hist_fit_e_b, "B-only fit", "F") - leg_rat.AddEntry(hist_fit_e_s, "S+B fit", "F") - leg_rat.Draw() - line_one = ROOT.TLine(0, 1, hist_fit_e_s.GetXaxis().GetXmax(), 1) - line_one.SetLineColor(1) - line_one.SetLineStyle(2) - line_one.SetLineWidth(2) - line_one.Draw() - canvas_pferrs.RedrawAxis() - - fout.WriteTObject(canvas_pferrs) + for idx in range(len(hist_dnll)): + canvas_nuis = ROOT.TCanvas("nuisances_%s" % idx, "nuisances_%s" % idx, 900, 600) + hist_fit_e_s = hist_fit_s[idx].Clone("errors_s %s" % idx) + hist_fit_e_b = hist_fit_b[idx].Clone("errors_b %s" % idx) + # gr_fit_s = getGraph(hist_fit_s,-0.1) + # gr_fit_b = getGraph(hist_fit_b, 0.1) + gr_fit_s[idx].SetLineColor(ROOT.kRed) + gr_fit_s[idx].SetMarkerColor(ROOT.kRed) + gr_fit_b[idx].SetLineColor(ROOT.kBlue) + gr_fit_b[idx].SetMarkerColor(ROOT.kBlue) + gr_fit_b[idx].SetMarkerStyle(20) + gr_fit_s[idx].SetMarkerStyle(20) + gr_fit_b[idx].SetMarkerSize(1.0) + gr_fit_s[idx].SetMarkerSize(1.0) + gr_fit_b[idx].SetLineWidth(2) + gr_fit_s[idx].SetLineWidth(2) + hist_prefit[idx].SetLineWidth(2) + hist_prefit[idx].SetTitle("Nuisance Paramaeters") + hist_prefit[idx].SetLineColor(ROOT.kBlack) + hist_prefit[idx].SetFillColor(ROOT.kGray) + hist_prefit[idx].SetMaximum(3) + hist_prefit[idx].SetMinimum(-3) + hist_prefit[idx].Draw("E2") + hist_prefit[idx].Draw("histsame") + gr_fit_b[idx].Draw("EPsame") + gr_fit_s[idx].Draw("EPsame") + canvas_nuis.SetGridx() + canvas_nuis.RedrawAxis() + canvas_nuis.RedrawAxis("g") + leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg.SetFillColor(0) + leg.SetTextFont(42) + leg.AddEntry(hist_prefit[idx], "Prefit", "FL") + leg.AddEntry(gr_fit_b[idx], "B-only fit", "EPL") + leg.AddEntry(gr_fit_s[idx], "S+B fit", "EPL") + leg.Draw() + fout.WriteTObject(canvas_nuis) + + canvas_pferrs = ROOT.TCanvas("post_fit_errs_%s" % idx, "post_fit_errs_%s" % idx, 900, 600) + for b in range(1, hist_fit_e_s.GetNbinsX() + 1): + hist_fit_e_s.SetBinContent(b, hist_fit_s[idx].GetBinError(b) / hist_prefit[idx].GetBinError(b)) + hist_fit_e_b.SetBinContent(b, hist_fit_b[idx].GetBinError(b) / hist_prefit[idx].GetBinError(b)) + hist_fit_e_s.SetBinError(b, 0) + hist_fit_e_b.SetBinError(b, 0) + hist_fit_e_s.SetFillColor(ROOT.kRed) + hist_fit_e_b.SetFillColor(ROOT.kBlue) + hist_fit_e_s.SetBarWidth(0.4) + hist_fit_e_b.SetBarWidth(0.4) + hist_fit_e_b.SetBarOffset(0.45) + hist_fit_e_b.GetYaxis().SetTitle("#sigma_{#theta}/(#sigma_{#theta} prefit)") + hist_fit_e_b.SetTitle("Nuisance Parameter Uncertainty Reduction") + hist_fit_e_b.SetMaximum(1.5) + hist_fit_e_b.SetMinimum(0) + hist_fit_e_b.Draw("bar") + hist_fit_e_s.Draw("barsame") + leg_rat = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg_rat.SetFillColor(0) + leg_rat.SetTextFont(42) + leg_rat.AddEntry(hist_fit_e_b, "B-only fit", "F") + leg_rat.AddEntry(hist_fit_e_s, "S+B fit", "F") + leg_rat.Draw() + line_one = ROOT.TLine(0, 1, hist_fit_e_s.GetXaxis().GetXmax(), 1) + line_one.SetLineColor(1) + line_one.SetLineStyle(2) + line_one.SetLineWidth(2) + line_one.Draw() + canvas_pferrs.RedrawAxis() + + fout.WriteTObject(canvas_pferrs) + + + canvas_dnll = ROOT.TCanvas("dnlls_%s" % idx, "dnlls_%s" % idx, 900, 600) + hist_dnll[idx].SetStats(0) + hist_dnll[idx].SetLineWidth(2) + hist_dnll[idx].SetTitle("Nuisance Parameters") + hist_dnll[idx].SetLineColor(ROOT.kBlack) + hist_dnll[idx].SetFillColor(ROOT.kGray) + hist_dnll[idx].SetMaximum(2) + hist_dnll[idx].SetMinimum(-1) + hist_dnll[idx].GetYaxis().SetTitle("NLL(fit_{b}) - NLL(fit_{s+b})") + hist_dnll[idx].Draw() + + hist_cdnll[idx].SetLineWidth(2) + hist_cdnll[idx].SetLineColor(ROOT.kBlue) + hist_cdnll[idx].Draw("same") + canvas_dnll.SetGridx() + canvas_dnll.RedrawAxis() + canvas_dnll.RedrawAxis("g") + leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg.SetFillColor(0) + leg.SetTextFont(42) + leg.AddEntry(hist_dnll[idx], "individual nuisance", "FL") + leg.AddEntry(hist_cdnll[idx], "cumulative", "FL") + leg.Draw() + latex = ROOT.TLatex() + posx = 4 if hist_dnll[idx].GetNbinsX() > 10 else 0.5 + latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll) ) + fout.WriteTObject(canvas_dnll) From 7072001ae7ec05069c8afeea850d6f179d9b49f2 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Tue, 21 Mar 2023 15:51:06 +0100 Subject: [PATCH 03/11] move diffNuisances to scripts, make executable --- {test => scripts}/diffNuisances.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {test => scripts}/diffNuisances.py (100%) mode change 100644 => 100755 diff --git a/test/diffNuisances.py b/scripts/diffNuisances.py old mode 100644 new mode 100755 similarity index 100% rename from test/diffNuisances.py rename to scripts/diffNuisances.py From 0d191b7f9c22a3eb06810ea59b923d7c5050943b Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Tue, 21 Mar 2023 19:23:39 +0100 Subject: [PATCH 04/11] Update all formatting print options --- scripts/diffNuisances.py | 232 ++++++++++++++++++++------------------- 1 file changed, 118 insertions(+), 114 deletions(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index 0549ec14e25..4665fce2309 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -175,8 +175,10 @@ if options.pullDef: options.show_all_parameters = True -if options.sortBy not in ["correlation", "impact"]: - exit("choose one of [ %s ] for --sortBy" % (",".join()["correlation", "impact"])) +if options.sortBy not in ["correlation", "impact","dnll"]: + exit("choose one of [ %s ] for --sortBy" % (",".join()["correlation", "impact","dnll"])) +elif options.sortBy == "dnll" and (not options.workspace): + exit("can only sort by deltaNLL if a workspace is provided, otherwise no dnll can be calcualted") if options.regex != ".*": print("Including only nuisance parameters following this regex query:") @@ -356,12 +358,22 @@ def getGraph(hist,shift): else: row += ["%+.2f +%.2f %.2f" % (nx, neu, ned)] - if nuis_p != None: - if options.plotfile: - if options.workspace: + if options.workspace: + if fit_name == "b": pdf = workspace.pdf(name + "_Pdf") var = workspace.var(name) var_val = var.getVal() + var.setVal(nuis_x.getVal()) + bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) + + if options.workspace: + var.setVal(nuis_x.getVal()) + sfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) + + if nuis_p != None: + if options.plotfile: if fit_name == "b": nuis_p_i += 1 hist_idx = (nuis_p_i - 1) // options.max_nuis @@ -383,10 +395,6 @@ def getGraph(hist,shift): hist_fit_b[hist_idx].SetBinError(bin_idx, nuis_x.getError()) hist_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) gr_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - if options.workspace: - var.setVal(nuis_x.getVal()) - bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) - var.setVal(var_val) if fit_name == "s": if options.pullDef and nuis_p != None: # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) @@ -405,10 +413,6 @@ def getGraph(hist,shift): hist_fit_s[hist_idx].SetBinError(bin_idx, nuis_x.getError()) hist_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) gr_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - if options.workspace: - var.setVal(nuis_x.getVal()) - sfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) - var.setVal(var_val) hist_prefit[hist_idx].SetBinContent(bin_idx, mean_p) hist_prefit[hist_idx].SetBinError(bin_idx, sigma_p) hist_prefit[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) @@ -433,11 +437,11 @@ def getGraph(hist,shift): sigShift = nuis_x.getError() if options.pullDef: - row += [""] + row[-1] = "" elif options.absolute_values: - row += [" (%+4.2fsig, %4.2f)" % (valShift, sigShift)] + row[-1] = " (%+4.2fsig, %4.2f)" % (valShift, sigShift) else: - row += [" %+4.2f, %4.2f" % (valShift, sigShift)] + row[-1] = " %+4.2f, %4.2f" % (valShift, sigShift) if fit_name == "b": pulls.append(valShift) @@ -469,6 +473,7 @@ def getGraph(hist,shift): # end of loop over s and b row += ["%+4.2f" % fit_s.correlation(name, options.poi)] + row += ["%+4.3f" % (nuis_x.getError() * fit_s.correlation(name, options.poi) * error_poi)] if options.workspace: dnll = bfit_nll - sfit_nll dnlls.append( (name, dnll) ) @@ -478,23 +483,24 @@ def getGraph(hist,shift): # end of loop over all fitted parameters -#fill dnll and cumulative dnll -len_dnll = len(dnlls) -hist_dnll = [] -hist_cdnll = [] -for hist_idx in range( ((len_dnll - 1) // options.max_nuis) + 1 ): - nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx ) - hist_dnll += [ ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] - hist_cdnll += [ ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] -dnlls.sort( key= lambda x: -x[1] ) -cdnll = 0 -for idx, (nm, val) in enumerate(dnlls): - hist_idx = idx // options.max_nuis - bin_idx = idx % options.max_nuis - hist_dnll[hist_idx].SetBinContent(bin_idx+1, val) - hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx+1, nm) - cdnll += val - hist_cdnll[hist_idx].SetBinContent( bin_idx+1, cdnll ) +#fill dnll and cumulative dnll plots +if options.plotfile and options.workspace: + len_dnll = len(dnlls) + hist_dnll = [] + hist_cdnll = [] + for hist_idx in range( ((len_dnll - 1) // options.max_nuis) + 1 ): + nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx ) + hist_dnll += [ ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] + hist_cdnll += [ ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] + dnlls.sort( key= lambda x: -x[1] ) + cdnll = 0 + for idx, (nm, val) in enumerate(dnlls): + hist_idx = idx // options.max_nuis + bin_idx = idx % options.max_nuis + hist_dnll[hist_idx].SetBinContent(bin_idx+1, val) + hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx+1, nm) + cdnll += val + hist_cdnll[hist_idx].SetBinContent( bin_idx+1, cdnll ) # ---------- # print the results @@ -531,53 +537,37 @@ def getGraph(hist,shift): sigsub = ("sig", r"$\\sigma$") highlight = "\\textbf{%s}" morelight = "{{\\color{red}\\textbf{%s}}}" + headers_first = None if options.skipFitS: print(" option '--skipFitS' set true. $s+b$ Fit is just a copy of the $b$-only fit") if options.skipFitB: print(" option '--skipFitB' set true. $b$-only Fit is just a copy of the $s+b$ fit") if options.pullDef: - fmtstring = "%-40s & %30s & %30s & %6s & %6s \\\\" - print("\\begin{tabular}{|l|r|r|r|r|} \\hline ") - print( - ( - fmtstring - % ( - "name", - "$b$-only fit pull", - "$s+b$ fit pull", - r"$\rho(\theta, \mu)$", - r"I(\theta, \mu)", - ) - ), - " \\hline", - ) + fmtstring = "%-40s & %30s & %30s & %6s & %6s " + tab_header="\\begin{tabular}{|l|r|r|r|r|" + headers = ["name", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)" ] elif options.absolute_values: - fmtstring = "%-40s & %15s & %30s & %30s & %6s & %6s \\\\" - print("\\begin{tabular}{|l|r|r|r|r|r|} \\hline ") - print( - ( - fmtstring - % ( - "name", - "pre fit", - "$b$-only fit", - "$s+b$ fit", - r"$\rho(\theta, \mu)$", - r"I(\theta, \mu)", - ) - ), - " \\hline", - ) + fmtstring = "%-40s & %15s & %30s & %30s & %6s & %6s " + tab_header = "\\begin{tabular}{|l|r|r|r|r|r|" + headers = ["name", "pre fit", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)" ] else: - fmtstring = "%-40s & %15s & %15s & %6s & %6s \\\\" - print("\\begin{tabular}{|l|r|r|r|r|} \\hline ") + fmtstring = "%-40s & %15s & %15s & %6s & %6s " + tab_header = "\\begin{tabular}{|l|r|r|r|r|" # what = r"$(x_\text{out} - x_\text{in})/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$" what = r"\Delta x/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$" - print(fmtstring % ("", "$b$-only fit", "$s+b$ fit", "", "")) - print( - (fmtstring % ("name", what, what, r"$\rho(\theta, \mu)$", r"I(\theta, \mu)")), - " \\hline", - ) + headers_first = [ "", "$b$-only fit", "$s+b$ fit", "", "" ] + headers = ["name", what, what, r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] + if options.workspace: + tab_header += "r|" + fmtstring += " & %6s " + if headers_first is not None: + headers_first += [""] + headers += [r"\Delta NLL"] + fmtstring += " \\\\" + print( tab_header + "} \\hline " ) + if headers_first is not None: + print( fmtstring % tuple(headers_first) ) + print( fmtstring % tuple(headers), " \\hline") elif options.format == "twiki": pmsub = (r"(\S+) \+/- (\S+)", r"\1 ± \2") sigsub = ("sig", r"σ") @@ -589,13 +579,18 @@ def getGraph(hist,shift): print(" option '--skipFitB' set true. $b$-only Fit is just a copy of the $s+b$ fit") if options.pullDef: fmtstring = "| %-40s | %-30s | %-30s | %-15s | %-15s | %-15s |" - print("| *name* | *b-only fit pull* | *s+b fit pull* | *corr.* | *approx. impact* |") + header = "| *name* | *b-only fit pull* | *s+b fit pull* | *corr.* | *approx. impact* |" elif options.absolute_values: fmtstring = "| %-40s | %-15s | %-30s | %-30s | %-15s | %-15s |" - print("| *name* | *pre fit* | *b-only fit* | *s+b fit* | *corr.* | *approx. impact* |") + header = "| *name* | *pre fit* | *b-only fit* | *s+b fit* | *corr.* | *approx. impact* |" else: fmtstring = "| %-40s | %-15s | %-15s | %-15s | %-15s |" - print("| *name* | *b-only fit* | *s+b fit* | *corr.* | *approx. impact* |") + header = "| *name* | *b-only fit* | *s+b fit* | *corr.* | *approx. impact* |" + if options.workspace: + fmtstring += " %-15s |" + header += " *delta nll* |" + print(header) + elif options.format == "html": pmsub = (r"(\S+) \+/- (\S+)", r"\1 ± \2") sigsub = ("sig", r"σ") @@ -614,28 +609,38 @@ def getGraph(hist,shift): """ ) if options.pullDef: - print("nuisancebackground fit pull signal fit pullρ(μ, θ)I(μ, θ)") + header = "nuisancebackground fit pull signal fit pullρ(μ, θ)I(μ, θ)" fmtstring = "%-40s %-30s %-30s %-15s %-15s " elif options.absolute_values: - print("nuisancepre fitbackground fit signal fitcorrelation") - fmtstring = "%-40s %-15s %-30s %-30s %-15s %-15s " + header = "nuisancepre fitbackground fit signal fitcorrelation" + fmtstring = "%-40s %-15s %-30s %-30s %-15s %-15s " else: what = "Δx/σin, σout/σin" - print( - "nuisancebackground fit
%s signal fit
%sρ(μ, θ)I(μ, θ)" - % (what, what) - ) - fmtstring = "%-40s %-15s %-15s %-15s %-15s " + header = "nuisancebackground fit
%s signal fit
%sρ(μ, θ)I(μ, θ)" % (what, what) + fmtstring = "%-40s %-15s %-15s %-15s %-15s " + if options.workspace: + header += "ΔNLL" + fmtstring += " %-15s " + header += "" + print(header) + fmtstring += "" names = list(table.keys()) names.sort() +rho_idx = -2 if not options.workspace else -3 +imp_idx = -1 if not options.workspace else -2 if options.sortBy == "correlation": - names = [[abs(float(table[t][-2])), t] for t in table.keys()] + names = [[abs(float(table[t][rho_idx])), t] for t in table.keys()] names.sort() names.reverse() names = [n[1] for n in names] elif options.sortBy == "impact": - names = [[abs(float(table[t][-1])), t] for t in table.keys()] + names = [[abs(float(table[t][imp_idx])), t] for t in table.keys()] + names.sort() + names.reverse() + names = [n[1] for n in names] +elif options.sortBy == "dnll": + names = [[float(table[t][-1]), t] for t in table.keys()] names.sort() names.reverse() names = [n[1] for n in names] @@ -657,8 +662,6 @@ def getGraph(hist,shift): n = n.replace(r"_", r"\_") j = [n] + v - if not options.absolute_values: - j = j[:1] + j[2:] print(fmtstring % tuple(j)) if options.format == "latex": @@ -688,7 +691,7 @@ def getGraph(hist,shift): histogram.Draw("pe") fout.WriteTObject(canvas) - for idx in range(len(hist_dnll)): + for idx in range(len(hist_fit_s)): canvas_nuis = ROOT.TCanvas("nuisances_%s" % idx, "nuisances_%s" % idx, 900, 600) hist_fit_e_s = hist_fit_s[idx].Clone("errors_s %s" % idx) hist_fit_e_b = hist_fit_b[idx].Clone("errors_b %s" % idx) @@ -759,30 +762,31 @@ def getGraph(hist,shift): fout.WriteTObject(canvas_pferrs) - canvas_dnll = ROOT.TCanvas("dnlls_%s" % idx, "dnlls_%s" % idx, 900, 600) - hist_dnll[idx].SetStats(0) - hist_dnll[idx].SetLineWidth(2) - hist_dnll[idx].SetTitle("Nuisance Parameters") - hist_dnll[idx].SetLineColor(ROOT.kBlack) - hist_dnll[idx].SetFillColor(ROOT.kGray) - hist_dnll[idx].SetMaximum(2) - hist_dnll[idx].SetMinimum(-1) - hist_dnll[idx].GetYaxis().SetTitle("NLL(fit_{b}) - NLL(fit_{s+b})") - hist_dnll[idx].Draw() - - hist_cdnll[idx].SetLineWidth(2) - hist_cdnll[idx].SetLineColor(ROOT.kBlue) - hist_cdnll[idx].Draw("same") - canvas_dnll.SetGridx() - canvas_dnll.RedrawAxis() - canvas_dnll.RedrawAxis("g") - leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg.SetFillColor(0) - leg.SetTextFont(42) - leg.AddEntry(hist_dnll[idx], "individual nuisance", "FL") - leg.AddEntry(hist_cdnll[idx], "cumulative", "FL") - leg.Draw() - latex = ROOT.TLatex() - posx = 4 if hist_dnll[idx].GetNbinsX() > 10 else 0.5 - latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll) ) - fout.WriteTObject(canvas_dnll) + if options.workspace: + canvas_dnll = ROOT.TCanvas("dnlls_%s" % idx, "dnlls_%s" % idx, 900, 600) + hist_dnll[idx].SetStats(0) + hist_dnll[idx].SetLineWidth(2) + hist_dnll[idx].SetTitle("Nuisance Parameters") + hist_dnll[idx].SetLineColor(ROOT.kBlack) + hist_dnll[idx].SetFillColor(ROOT.kGray) + hist_dnll[idx].SetMaximum(2) + hist_dnll[idx].SetMinimum(-1) + hist_dnll[idx].GetYaxis().SetTitle("NLL(fit_{b}) - NLL(fit_{s+b})") + hist_dnll[idx].Draw() + + hist_cdnll[idx].SetLineWidth(2) + hist_cdnll[idx].SetLineColor(ROOT.kBlue) + hist_cdnll[idx].Draw("same") + canvas_dnll.SetGridx() + canvas_dnll.RedrawAxis() + canvas_dnll.RedrawAxis("g") + leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) + leg.SetFillColor(0) + leg.SetTextFont(42) + leg.AddEntry(hist_dnll[idx], "individual nuisance", "FL") + leg.AddEntry(hist_cdnll[idx], "cumulative", "FL") + leg.Draw() + latex = ROOT.TLatex() + posx = 4 if hist_dnll[idx].GetNbinsX() > 10 else 0.5 + latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll) ) + fout.WriteTObject(canvas_dnll) From d35b6e62854a0d60108bda3f27c23cc2d0493f34 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 08:57:13 +0100 Subject: [PATCH 05/11] Remove separate version of diffNusiances in tutorial, update doc --- data/tutorials/longexercise/diffNuisances.py | 680 ------------------- docs/part5/longexercise.md | 2 +- 2 files changed, 1 insertion(+), 681 deletions(-) delete mode 100644 data/tutorials/longexercise/diffNuisances.py diff --git a/data/tutorials/longexercise/diffNuisances.py b/data/tutorials/longexercise/diffNuisances.py deleted file mode 100644 index 86323da1793..00000000000 --- a/data/tutorials/longexercise/diffNuisances.py +++ /dev/null @@ -1,680 +0,0 @@ -#!/usr/bin/env python3 -from __future__ import absolute_import, print_function - -import datetime -import re -from optparse import OptionParser -from sys import argv, exit, stderr, stdout - -from six.moves import range - -import HiggsAnalysis.CombinedLimit.calculate_pulls as CP -import ROOT - -# tool to compare fitted nuisance parameters to prefit values. -# -# Also used to check for potential problems in RooFit workspaces to be used with combine -# (see https://twiki.cern.ch/twiki/bin/viewauth/CMS/HiggsWG/HiggsPAGPreapprovalChecks) - -# import ROOT with a fix to get batch mode (http://root.cern.ch/phpBB3/viewtopic.php?t=3198) -hasHelp = False -for X in ("-h", "-?", "--help"): - if X in argv: - hasHelp = True - argv.remove(X) -argv.append("-b-") - -ROOT.gROOT.SetBatch(True) -# ROOT.gSystem.Load("libHiggsAnalysisCombinedLimit") -argv.remove("-b-") -if hasHelp: - argv.append("-h") - -parser = OptionParser(usage="usage: %prog [options] in.root \nrun with --help to get list of options") -parser.add_option( - "--vtol", - "--val-tolerance", - dest="vtol", - default=0.30, - type="float", - help="Report nuisances whose value changes by more than this amount of sigmas", -) -parser.add_option( - "--stol", - "--sig-tolerance", - dest="stol", - default=0.10, - type="float", - help="Report nuisances whose sigma changes by more than this amount", -) -parser.add_option( - "--vtol2", - "--val-tolerance2", - dest="vtol2", - default=2.0, - type="float", - help="Report severely nuisances whose value changes by more than this amount of sigmas", -) -parser.add_option( - "--stol2", - "--sig-tolerance2", - dest="stol2", - default=0.50, - type="float", - help="Report severely nuisances whose sigma changes by more than this amount", -) -parser.add_option( - "-a", - "--all", - dest="show_all_parameters", - default=False, - action="store_true", - help="Print all nuisances, even the ones which are unchanged w.r.t. pre-fit values.", -) -parser.add_option( - "-A", - "--abs", - dest="absolute_values", - default=False, - action="store_true", - help="Report also absolute values of nuisance values and errors, not only the ones normalized to the input sigma", -) -parser.add_option( - "-p", - "--poi", - dest="poi", - default="r", - type="string", - help="Name of signal strength parameter (default is 'r' as per text2workspace.py)", -) -parser.add_option( - "-f", - "--format", - dest="format", - default="text", - type="string", - help="Output format ('text', 'latex', 'twiki'", -) -parser.add_option( - "-g", - "--histogram", - dest="plotfile", - default=None, - type="string", - help="If true, plot the pulls of the nuisances to the given file.", -) -parser.add_option( - "", - "--pullDef", - dest="pullDef", - default="", - type="string", - help="Choose the definition of the pull, see python/calculate_pulls.py for options", -) -parser.add_option( - "-w", - "--workspace", - dest="workspace", - default="", - type="string", - help="Workspace to use for evaluating NLL differences." -) -parser.add_option( - "", - "--max-nuis", - dest="max_nuis", - default=65, - type="int", - help="Maximum nuisances for a single plot" -) - -(options, args) = parser.parse_args() -if len(args) == 0: - parser.print_usage() - exit(1) - -if options.pullDef != "" and options.pullDef not in CP.allowed_methods(): - exit("Method %s not allowed, choose one of [%s]" % (options.pullDef, ",".join(CP.allowed_methods()))) - -if options.pullDef and options.absolute_values: - print("Pulls are always defined as absolute, will modify --absolute_values to False for you") - options.absolute_values = False - -if options.pullDef: - options.show_all_parameters = True - -setUpString = "diffNuisances run on %s, at %s with the following options ... " % ( - args[0], - datetime.datetime.utcnow(), -) + str(options) - -file = ROOT.TFile(args[0]) -workspace = None -if options.workspace: - workspace_file = ROOT.TFile( options.workspace, "READ" ) - workspace = workspace_file.Get("w") - -if file == None: - raise RuntimeError("Cannot open file %s" % args[0]) -fit_s = file.Get("fit_s") -fit_b = file.Get("fit_b") -prefit = file.Get("nuisances_prefit") -if fit_s == None or fit_s.ClassName() != "RooFitResult": - raise RuntimeError("File %s does not contain the output of the signal fit 'fit_s'" % args[0]) -if fit_b == None or fit_b.ClassName() != "RooFitResult": - raise RuntimeError("File %s does not contain the output of the background fit 'fit_b'" % args[0]) -if prefit == None or prefit.ClassName() != "RooArgSet": - raise RuntimeError("File %s does not contain the prefit nuisances 'nuisances_prefit'" % args[0]) - -isFlagged = {} - -# maps from nuisance parameter name to the row to be printed in the table -table = {} - -# get the fitted parameters -fpf_b = fit_b.floatParsFinal() -fpf_s = fit_s.floatParsFinal() - -pulls = [] -dnlls = [] - -nuis_p_i = 0 -title = "pull" if options.pullDef else "#theta" - -""" -def getGraph(hist,shift): - - gr = ROOT.TGraphAsymErrors() - gr.SetName(hist.GetName()) - for i in range(hist.GetNbinsX()): - x = hist.GetBinCenter(i+1)+shift - y = hist.GetBinContent(i+1) - e = hist.GetBinError(i+1) - gr.SetPoint(i,x,y) - gr.SetPointError(i,float(abs(shift))*0.8,e) - return gr -""" - -n_hists = prefit.getSize() // options.max_nuis -if prefit.getSize() % options.max_nuis == 0: - n_hists -= 1 - -hist_fit_b = [] -hist_fit_s = [] -hist_prefit = [] -gr_fit_b = [] -gr_fit_s = [] -for idx in range( ((prefit.getSize() - 1) // options.max_nuis) + 1 ): - # Also make histograms for pull distributions: - nbins = min(prefit.getSize(), options.max_nuis, prefit.getSize() - options.max_nuis * idx ) - hist_fit_b += [ ROOT.TH1F("fit_b %s" % idx, "B-only fit Nuisances;;%s %s" % (title, idx), nbins, 0, nbins)] - hist_fit_s += [ ROOT.TH1F("fit_s %s" % idx, "S+B fit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins)] - hist_prefit +=[ ROOT.TH1F( - "prefit_nuisancs %s" % idx, - "Prefit Nuisances ;;%s %s" % (title, idx), - nbins, - 0, - nbins, - ) ] - # Store also the *asymmetric* uncertainties - gr_fit_b += [ ROOT.TGraphAsymmErrors() ] - gr_fit_b[idx].SetTitle("fit_b_g % s" % idx) - gr_fit_s += [ ROOT.TGraphAsymmErrors() ] - gr_fit_s[idx].SetTitle("fit_b_s % s" % idx) - - -# loop over all fitted parameters -dnll = [] -for i in range(fpf_s.getSize()): - nuis_s = fpf_s.at(i) - name = nuis_s.GetName() - nuis_b = fpf_b.find(name) - nuis_p = prefit.find(name) - - # keeps information to be printed about the nuisance parameter - row = [] - - flag = False - mean_p, sigma_p, sigma_pu, sigma_pd = 0, 0, 0, 0 - - if nuis_p == None: - # nuisance parameter NOT present in the prefit result - if not options.absolute_values and not (options.pullDef == "unconstPullAsym"): - continue - row += ["[%.2f, %.2f]" % (nuis_s.getMin(), nuis_s.getMax())] - - else: - # get best-fit value and uncertainty at prefit for this - # nuisance parameter - if nuis_p.getErrorLo() == 0: - nuis_p.setError(nuis_p.getErrorHi()) - mean_p, sigma_p, sigma_pu, sigma_pd = ( - nuis_p.getVal(), - nuis_p.getError(), - nuis_p.getErrorHi(), - nuis_p.getErrorLo(), - ) - - if not sigma_p > 0: - sigma_p = (nuis_p.getMax() - nuis_p.getMin()) / 2 - nuisIsSymm = abs(abs(nuis_p.getErrorLo()) - abs(nuis_p.getErrorHi())) < 0.01 or nuis_p.getErrorLo() == 0 - if options.absolute_values: - if nuisIsSymm: - row += ["%.6f +/- %.6f" % (nuis_p.getVal(), nuis_p.getError())] - else: - row += ["%.6f +%.6f %.6f" % (nuis_p.getVal(), nuis_p.getErrorHi(), nuis_p.getErrorLo())] - - for fit_name, nuis_x in [("b", nuis_b), ("s", nuis_s)]: - if nuis_x == None: - row += [" n/a "] - else: - nuisIsSymm = abs(abs(nuis_x.getErrorLo()) - abs(nuis_x.getErrorHi())) < 0.01 or nuis_x.getErrorLo() == 0 - if nuisIsSymm: - nuis_x.setError(nuis_x.getErrorHi()) - nuiselo = abs(nuis_x.getErrorLo()) if nuis_x.getErrorLo() > 0 else nuis_x.getError() - nuisehi = nuis_x.getErrorHi() - if options.pullDef and nuis_p != None: - nx, ned, neu = CP.returnPullAsym( - options.pullDef, - nuis_x.getVal(), - mean_p, - nuisehi, - sigma_pu, - abs(nuiselo), - abs(sigma_pd), - ) - else: - nx, ned, neu = nuis_x.getVal(), nuiselo, nuisehi - - if nuisIsSymm: - row += ["%+.2f +/- %.2f" % (nx, (abs(ned) + abs(neu)) / 2)] - else: - row += ["%+.2f +%.2f %.2f" % (nx, neu, ned)] - - if nuis_p != None: - if options.plotfile: - pdf = workspace.pdf(name + "_Pdf") - var = workspace.var(name) - var_val = var.getVal() - if fit_name == "b": - nuis_p_i += 1 - hist_idx = (nuis_p_i - 1) // options.max_nuis - bin_idx = (nuis_p_i - 1) % options.max_nuis + 1 - print(hist_idx, bin_idx) - if options.pullDef and nuis_p != None: - # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) - gr_fit_b[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 + 0.1, nx) - gr_fit_b[hist_idx].SetPointError(bin_idx - 1, 0, 0, ned, neu) - else: - gr_fit_b[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 + 0.1, nuis_x.getVal()) - gr_fit_b[hist_idx].SetPointError( - bin_idx - 1, - 0, - 0, - abs(nuis_x.getErrorLo()), - nuis_x.getErrorHi(), - ) - hist_fit_b[hist_idx].SetBinContent(bin_idx, nuis_x.getVal()) - hist_fit_b[hist_idx].SetBinError(bin_idx, nuis_x.getError()) - hist_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - gr_fit_b[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - var.setVal(nuis_x.getVal()) - bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) - var.setVal(var_val) - if fit_name == "s": - if options.pullDef and nuis_p != None: - # nx,ned,neu = CP.returnPullAsym(options.pullDef,nuis_x.getVal(),mean_p,nuis_x.getErrorHi(),sigma_pu,abs(nuis_x.getErrorLo()),abs(sigma_pd)) - gr_fit_s[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 - 0.1, nx) - gr_fit_s[hist_idx].SetPointError(bin_idx - 1, 0, 0, ned, neu) - else: - gr_fit_s[hist_idx].SetPoint(bin_idx - 1, bin_idx - 0.5 - 0.1, nuis_x.getVal()) - gr_fit_s[hist_idx].SetPointError( - bin_idx - 1, - 0, - 0, - abs(nuis_x.getErrorLo()), - nuis_x.getErrorHi(), - ) - hist_fit_s[hist_idx].SetBinContent(bin_idx, nuis_x.getVal()) - hist_fit_s[hist_idx].SetBinError(bin_idx, nuis_x.getError()) - hist_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - gr_fit_s[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - var.setVal(nuis_x.getVal()) - sfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) - var.setVal(var_val) - hist_prefit[hist_idx].SetBinContent(bin_idx, mean_p) - hist_prefit[hist_idx].SetBinError(bin_idx, sigma_p) - hist_prefit[hist_idx].GetXaxis().SetBinLabel(bin_idx, name) - - if sigma_p > 0: - if options.pullDef: - valShift = nx - sigShift = 1 - else: - # calculate the difference of the nuisance parameter - # w.r.t to the prefit value in terms of the uncertainty - # on the prefit value - valShift = (nuis_x.getVal() - mean_p) / sigma_p - - # ratio of the nuisance parameter's uncertainty - # w.r.t the prefit uncertainty - sigShift = nuis_x.getError() / sigma_p - - else: - # print "No definition for prefit uncertainty %s. Printing absolute shifts"%(nuis_p.GetName()) - valShift = nuis_x.getVal() - mean_p - sigShift = nuis_x.getError() - - if options.pullDef: - row[-1] += "" - elif options.absolute_values: - row[-1] += " (%+4.2fsig, %4.2f)" % (valShift, sigShift) - else: - row[-1] = " %+4.2f, %4.2f" % (valShift, sigShift) - - if fit_name == "b": - pulls.append(valShift) - - if abs(valShift) > options.vtol2 or abs(sigShift - 1) > options.stol2: - # severely report this nuisance: - # - # the best fit moved by more than 2.0 sigma or the uncertainty (sigma) - # changed by more than 50% (default thresholds) w.r.t the prefit values - - isFlagged[(name, fit_name)] = 2 - - flag = True - - elif abs(valShift) > options.vtol or abs(sigShift - 1) > options.stol: - # report this nuisance: - # - # the best fit moved by more than 0.3 sigma or the uncertainty (sigma) - # changed by more than 10% (default thresholds) w.r.t the prefit values - - if options.show_all_parameters: - isFlagged[(name, fit_name)] = 1 - - flag = True - - elif options.show_all_parameters: - flag = True - - dnll.append( (name, bfit_nll - sfit_nll) ) - # end of loop over s and b - - row += ["%+4.2f" % fit_s.correlation(name, options.poi)] - row += ["%.4f" % (bfit_nll - sfit_nll) ] - if flag or options.show_all_parameters: - table[name] = row - -len_dnll = len(dnll) -hist_dnll = [] -hist_cdnll = [] -for hist_idx in range( ((len_dnll - 1) // options.max_nuis) + 1 ): - nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx ) - hist_dnll += [ ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] - hist_cdnll += [ ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] -dnll.sort( key= lambda x: -x[1] ) -cdnll = 0 -for idx, (nm, val) in enumerate(dnll): - hist_idx = idx // options.max_nuis - bin_idx = idx % options.max_nuis - hist_dnll[hist_idx].SetBinContent(bin_idx+1, val) - hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx+1, nm) - cdnll += val - hist_cdnll[hist_idx].SetBinContent( bin_idx+1, cdnll ) -# fmt_string = "%-40s %.3f" -# print(fmt_string % (nm, val ) ) -# end of loop over all fitted parameters - -# ---------- -# print the results -# ---------- - -# print details -print(setUpString) -print() - -fmtstring = "%-40s %15s %15s %10s %10s" -highlight = "*%s*" -morelight = "!%s!" -pmsub, sigsub = None, None -if options.format == "text": - if options.pullDef: - fmtstring = "%-40s %30s %30s %10s" - print(fmtstring % ("name", "b-only fit pull", "s+b fit pull", "rho")) - elif options.absolute_values: - fmtstring = "%-40s %15s %30s %30s %10s" - print(fmtstring % ("name", "pre fit", "b-only fit", "s+b fit", "rho")) - else: - print(fmtstring % ("name", "b-only fit", "s+b fit", "rho","dnll")) -elif options.format == "latex": - pmsub = (r"(\S+) \+/- (\S+)", r"$\1 \\pm \2$") - sigsub = ("sig", r"$\\sigma$") - highlight = "\\textbf{%s}" - morelight = "{{\\color{red}\\textbf{%s}}}" - if options.pullDef: - fmtstring = "%-40s & %30s & %30s & %6s \\\\" - print("\\begin{tabular}{|l|r|r|r|} \\hline ") - print( - ( - fmtstring - % ( - "name", - "$b$-only fit pull", - "$s+b$ fit pull", - r"$\rho(\theta, \mu)$", - ) - ), - " \\hline", - ) - elif options.absolute_values: - fmtstring = "%-40s & %15s & %30s & %30s & %6s \\\\" - print("\\begin{tabular}{|l|r|r|r|r|} \\hline ") - print( - ( - fmtstring - % ( - "name", - "pre fit", - "$b$-only fit", - "$s+b$ fit", - r"$\rho(\theta, \mu)$", - ) - ), - " \\hline", - ) - else: - fmtstring = "%-40s & %15s & %15s & %6s \\\\" - print("\\begin{tabular}{|l|r|r|r|} \\hline ") - # what = r"$(x_\text{out} - x_\text{in})/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$" - what = r"\Delta x/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$" - print(fmtstring % ("", "$b$-only fit", "$s+b$ fit", "")) - print((fmtstring % ("name", what, what, r"$\rho(\theta, \mu)$")), " \\hline") -elif options.format == "twiki": - pmsub = (r"(\S+) \+/- (\S+)", r"\1 ± \2") - sigsub = ("sig", r"σ") - highlight = "%s" - morelight = "%s" - if options.pullDef: - fmtstring = "| %-40s | %-30s | %-30s | %-15s |" - print("| *name* | *b-only fit pull* | *s+b fit pull* | ") - elif options.absolute_values: - fmtstring = "| %-40s | %-15s | %-30s | %-30s | %-15s |" - print("| *name* | *pre fit* | *b-only fit* | *s+b fit* | ") - else: - fmtstring = "| %-40s | %-15s | %-15s | %-15s |" - print("| *name* | *b-only fit* | *s+b fit* | *corr.* |") -elif options.format == "html": - pmsub = (r"(\S+) \+/- (\S+)", r"\1 ± \2") - sigsub = ("sig", r"σ") - highlight = "%s" - morelight = "%s" - print( - """ -Comparison of nuisances - -

Comparison of nuisances

- -""" - ) - if options.pullDef: - print("") - fmtstring = "" - elif options.absolute_values: - print("") - fmtstring = "" - else: - what = "Δx/σin, σout/σin" - print("" % (what, what)) - fmtstring = "" - -names = list(table.keys()) -names.sort() -highlighters = {1: highlight, 2: morelight} -for n in names: - v = table[n] - if options.format == "latex": - n = n.replace(r"_", r"\_") - if pmsub != None: - v = [re.sub(pmsub[0], pmsub[1], i) for i in v] - if sigsub != None: - v = [re.sub(sigsub[0], sigsub[1], i) for i in v] - if (n, "b") in isFlagged: - v[-4] = highlighters[isFlagged[(n, "b")]] % v[-4] - if (n, "s") in isFlagged: - v[-3] = highlighters[isFlagged[(n, "s")]] % v[-3] - if options.absolute_values: - print(fmtstring % (n, v[0], v[1], v[2], v[3])) - else: - print(fmtstring % (n, v[0], v[1], v[2], v[3])) - -if options.format == "latex": - print(" \\hline\n\\end{tabular}") -elif options.format == "html": - print("
nuisancebackground fit pull signal fit pullcorrelation
%-40s %-30s %-30s %-15s
nuisancepre fitbackground fit signal fitcorrelation
%-40s %-15s %-30s %-30s %-15s
nuisancebackground fit
%s 		signal fit
%s		ρ(μ, θ)
%-40s %-15s %-15s %-15s
") - - -if options.plotfile: - import ROOT - - fout = ROOT.TFile(options.plotfile, "RECREATE") - ROOT.gROOT.SetStyle("Plain") - ROOT.gStyle.SetOptFit(1) - ROOT.gStyle.SetPadBottomMargin(0.3) - histogram = ROOT.TH1F("pulls", "Pulls", 60, -3, 3) - for pull in pulls: - histogram.Fill(pull) - canvas = ROOT.TCanvas("asdf","asdf", 800, 800) - if options.pullDef: - histogram.GetXaxis().SetTitle("pull") - else: - histogram.GetXaxis().SetTitle("(#theta-#theta_{0})/#sigma_{pre-fit}") - histogram.SetTitle("Post-fit nuisance pull distribution") - histogram.SetMarkerStyle(20) - histogram.SetMarkerSize(2) - histogram.Draw("pe") - fout.WriteTObject(canvas) - - for idx in range(len(hist_dnll)): - canvas_nuis = ROOT.TCanvas("nuisances_%s" % idx, "nuisances_%s" % idx, 900, 600) - hist_fit_e_s = hist_fit_s[idx].Clone("errors_s %s" % idx) - hist_fit_e_b = hist_fit_b[idx].Clone("errors_b %s" % idx) - # gr_fit_s = getGraph(hist_fit_s,-0.1) - # gr_fit_b = getGraph(hist_fit_b, 0.1) - gr_fit_s[idx].SetLineColor(ROOT.kRed) - gr_fit_s[idx].SetMarkerColor(ROOT.kRed) - gr_fit_b[idx].SetLineColor(ROOT.kBlue) - gr_fit_b[idx].SetMarkerColor(ROOT.kBlue) - gr_fit_b[idx].SetMarkerStyle(20) - gr_fit_s[idx].SetMarkerStyle(20) - gr_fit_b[idx].SetMarkerSize(1.0) - gr_fit_s[idx].SetMarkerSize(1.0) - gr_fit_b[idx].SetLineWidth(2) - gr_fit_s[idx].SetLineWidth(2) - hist_prefit[idx].SetLineWidth(2) - hist_prefit[idx].SetTitle("Nuisance Paramaeters") - hist_prefit[idx].SetLineColor(ROOT.kBlack) - hist_prefit[idx].SetFillColor(ROOT.kGray) - hist_prefit[idx].SetMaximum(3) - hist_prefit[idx].SetMinimum(-3) - hist_prefit[idx].Draw("E2") - hist_prefit[idx].Draw("histsame") - gr_fit_b[idx].Draw("EPsame") - gr_fit_s[idx].Draw("EPsame") - canvas_nuis.SetGridx() - canvas_nuis.RedrawAxis() - canvas_nuis.RedrawAxis("g") - leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg.SetFillColor(0) - leg.SetTextFont(42) - leg.AddEntry(hist_prefit[idx], "Prefit", "FL") - leg.AddEntry(gr_fit_b[idx], "B-only fit", "EPL") - leg.AddEntry(gr_fit_s[idx], "S+B fit", "EPL") - leg.Draw() - fout.WriteTObject(canvas_nuis) - - canvas_pferrs = ROOT.TCanvas("post_fit_errs_%s" % idx, "post_fit_errs_%s" % idx, 900, 600) - for b in range(1, hist_fit_e_s.GetNbinsX() + 1): - hist_fit_e_s.SetBinContent(b, hist_fit_s[idx].GetBinError(b) / hist_prefit[idx].GetBinError(b)) - hist_fit_e_b.SetBinContent(b, hist_fit_b[idx].GetBinError(b) / hist_prefit[idx].GetBinError(b)) - hist_fit_e_s.SetBinError(b, 0) - hist_fit_e_b.SetBinError(b, 0) - hist_fit_e_s.SetFillColor(ROOT.kRed) - hist_fit_e_b.SetFillColor(ROOT.kBlue) - hist_fit_e_s.SetBarWidth(0.4) - hist_fit_e_b.SetBarWidth(0.4) - hist_fit_e_b.SetBarOffset(0.45) - hist_fit_e_b.GetYaxis().SetTitle("#sigma_{#theta}/(#sigma_{#theta} prefit)") - hist_fit_e_b.SetTitle("Nuisance Parameter Uncertainty Reduction") - hist_fit_e_b.SetMaximum(1.5) - hist_fit_e_b.SetMinimum(0) - hist_fit_e_b.Draw("bar") - hist_fit_e_s.Draw("barsame") - leg_rat = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg_rat.SetFillColor(0) - leg_rat.SetTextFont(42) - leg_rat.AddEntry(hist_fit_e_b, "B-only fit", "F") - leg_rat.AddEntry(hist_fit_e_s, "S+B fit", "F") - leg_rat.Draw() - line_one = ROOT.TLine(0, 1, hist_fit_e_s.GetXaxis().GetXmax(), 1) - line_one.SetLineColor(1) - line_one.SetLineStyle(2) - line_one.SetLineWidth(2) - line_one.Draw() - canvas_pferrs.RedrawAxis() - - fout.WriteTObject(canvas_pferrs) - - - canvas_dnll = ROOT.TCanvas("dnlls_%s" % idx, "dnlls_%s" % idx, 900, 600) - hist_dnll[idx].SetStats(0) - hist_dnll[idx].SetLineWidth(2) - hist_dnll[idx].SetTitle("Nuisance Parameters") - hist_dnll[idx].SetLineColor(ROOT.kBlack) - hist_dnll[idx].SetFillColor(ROOT.kGray) - hist_dnll[idx].SetMaximum(2) - hist_dnll[idx].SetMinimum(-1) - hist_dnll[idx].GetYaxis().SetTitle("NLL(fit_{b}) - NLL(fit_{s+b})") - hist_dnll[idx].Draw() - - hist_cdnll[idx].SetLineWidth(2) - hist_cdnll[idx].SetLineColor(ROOT.kBlue) - hist_cdnll[idx].Draw("same") - canvas_dnll.SetGridx() - canvas_dnll.RedrawAxis() - canvas_dnll.RedrawAxis("g") - leg = ROOT.TLegend(0.6, 0.7, 0.89, 0.89) - leg.SetFillColor(0) - leg.SetTextFont(42) - leg.AddEntry(hist_dnll[idx], "individual nuisance", "FL") - leg.AddEntry(hist_cdnll[idx], "cumulative", "FL") - leg.Draw() - latex = ROOT.TLatex() - posx = 4 if hist_dnll[idx].GetNbinsX() > 10 else 0.5 - latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll) ) - fout.WriteTObject(canvas_dnll) diff --git a/docs/part5/longexercise.md b/docs/part5/longexercise.md index 3549809b404..f3ab9415045 100644 --- a/docs/part5/longexercise.md +++ b/docs/part5/longexercise.md @@ -297,7 +297,7 @@ Underneath this the best-fit values ($\theta$) and symmetrised uncertainties for A more useful way of looking at this is to compare the pre- and post-fit values of the parameters, to see how much the fit to data has shifted and constrained these parameters with respect to the input uncertainty. The script `diffNuisances.py` can be used for this: ```shell -python diffNuisances.py fitDiagnosticsTest.root --all +diffNuisances.py fitDiagnosticsTest.root --all ```
Show output From 3dad6ac63634cc5238a39750d5dd57a373e32739 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 09:32:07 +0100 Subject: [PATCH 06/11] Fix linting problems --- scripts/diffNuisances.py | 88 ++++++++++++++++++++-------------------- 1 file changed, 44 insertions(+), 44 deletions(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index 4665fce2309..0b82707db40 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -175,8 +175,8 @@ if options.pullDef: options.show_all_parameters = True -if options.sortBy not in ["correlation", "impact","dnll"]: - exit("choose one of [ %s ] for --sortBy" % (",".join()["correlation", "impact","dnll"])) +if options.sortBy not in ["correlation", "impact", "dnll"]: + exit("choose one of [ %s ] for --sortBy" % (",".join()["correlation", "impact", "dnll"])) elif options.sortBy == "dnll" and (not options.workspace): exit("can only sort by deltaNLL if a workspace is provided, otherwise no dnll can be calcualted") @@ -204,7 +204,7 @@ workspace = None if options.workspace: - workspace_file = ROOT.TFile( options.workspace, "READ" ) + workspace_file = ROOT.TFile(options.workspace, "READ") workspace = workspace_file.Get("w") isFlagged = {} @@ -246,7 +246,7 @@ def getGraph(hist,shift): nuis_p = prefit.find(name) if nuis_p == None: if not options.absolute_values and not (options.pullDef == "unconstPullAsym"): - continue + continue if bool(regex_NP_obj.match(name)): np_count += 1 @@ -262,26 +262,26 @@ def getGraph(hist,shift): gr_fit_s = [] hist_dnll = [] hist_cdnll = [] -for idx in range( ((np_count - 1) // options.max_nuis) + 1 ): +for idx in range(((np_count - 1) // options.max_nuis) + 1): # Also make histograms for pull distributions: - nbins = min(np_count, options.max_nuis, np_count - (options.max_nuis * idx) ) - hist_fit_b += [ ROOT.TH1F("fit_b %s" % idx, "B-only fit Nuisances;;%s %s" % (title, idx), nbins, 0, nbins)] - hist_fit_s += [ ROOT.TH1F("fit_s %s" % idx, "S+B fit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins)] - hist_prefit +=[ ROOT.TH1F( + nbins = min(np_count, options.max_nuis, np_count - (options.max_nuis * idx)) + hist_fit_b += [ROOT.TH1F("fit_b %s" % idx, "B-only fit Nuisances;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_fit_s += [ROOT.TH1F("fit_s %s" % idx, "S+B fit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_prefit += [ROOT.TH1F( "prefit_nuisancs %s" % idx, "Prefit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins, - ) ] + )] # Store also the *asymmetric* uncertainties - gr_fit_b += [ ROOT.TGraphAsymmErrors() ] + gr_fit_b += [ROOT.TGraphAsymmErrors()] gr_fit_b[idx].SetTitle("fit_b_g % s" % idx) - gr_fit_s += [ ROOT.TGraphAsymmErrors() ] + gr_fit_s += [ROOT.TGraphAsymmErrors()] gr_fit_s[idx].SetTitle("fit_b_s % s" % idx) - hist_dnll += [ ROOT.TH1F("dnll %s" % idx, "delta log-likelihoods ;;%s %s" % (title, idx), nbins, 0, nbins) ] - hist_cdnll += [ ROOT.TH1F("cdnll %s" % idx, "cumulative delta log-likelihoods ;;%s %s" % (title, idx), nbins, 0, nbins) ] + hist_dnll += [ROOT.TH1F("dnll %s" % idx, "delta log-likelihoods ;;%s %s" % (title, idx), nbins, 0, nbins)] + hist_cdnll += [ROOT.TH1F("cdnll %s" % idx, "cumulative delta log-likelihoods ;;%s %s" % (title, idx), nbins, 0, nbins)] error_poi = fpf_s.find(options.poi).getError() @@ -360,12 +360,12 @@ def getGraph(hist,shift): if options.workspace: if fit_name == "b": - pdf = workspace.pdf(name + "_Pdf") - var = workspace.var(name) - var_val = var.getVal() - var.setVal(nuis_x.getVal()) - bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) - var.setVal(var_val) + pdf = workspace.pdf(name + "_Pdf") + var = workspace.var(name) + var_val = var.getVal() + var.setVal(nuis_x.getVal()) + bfit_nll = -pdf.getLogVal(ROOT.RooArgSet(var)) + var.setVal(var_val) if options.workspace: var.setVal(nuis_x.getVal()) @@ -476,8 +476,8 @@ def getGraph(hist,shift): row += ["%+4.3f" % (nuis_x.getError() * fit_s.correlation(name, options.poi) * error_poi)] if options.workspace: dnll = bfit_nll - sfit_nll - dnlls.append( (name, dnll) ) - row += ["%.4f" % dnll ] + dnlls.append((name, dnll)) + row += ["%.4f" % dnll] if flag or options.show_all_parameters: table[name] = row @@ -488,19 +488,19 @@ def getGraph(hist,shift): len_dnll = len(dnlls) hist_dnll = [] hist_cdnll = [] - for hist_idx in range( ((len_dnll - 1) // options.max_nuis) + 1 ): - nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx ) - hist_dnll += [ ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] - hist_cdnll += [ ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins) ] - dnlls.sort( key= lambda x: -x[1] ) + for hist_idx in range(((len_dnll - 1) // options.max_nuis) + 1): + nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx) + hist_dnll += [ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins)] + hist_cdnll += [ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins)] + dnlls.sort(key= lambda x: -x[1]) cdnll = 0 for idx, (nm, val) in enumerate(dnlls): - hist_idx = idx // options.max_nuis - bin_idx = idx % options.max_nuis - hist_dnll[hist_idx].SetBinContent(bin_idx+1, val) - hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx+1, nm) - cdnll += val - hist_cdnll[hist_idx].SetBinContent( bin_idx+1, cdnll ) + hist_idx = idx // options.max_nuis + bin_idx = idx % options.max_nuis + hist_dnll[hist_idx].SetBinContent(bin_idx + 1, val) + hist_dnll[hist_idx].GetXaxis().SetBinLabel(bin_idx + 1, nm) + cdnll += val + hist_cdnll[hist_idx].SetBinContent(bin_idx + 1, cdnll) # ---------- # print the results @@ -545,17 +545,17 @@ def getGraph(hist,shift): if options.pullDef: fmtstring = "%-40s & %30s & %30s & %6s & %6s " tab_header="\\begin{tabular}{|l|r|r|r|r|" - headers = ["name", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)" ] + headers = ["name", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] elif options.absolute_values: fmtstring = "%-40s & %15s & %30s & %30s & %6s & %6s " tab_header = "\\begin{tabular}{|l|r|r|r|r|r|" - headers = ["name", "pre fit", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)" ] + headers = ["name", "pre fit", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] else: fmtstring = "%-40s & %15s & %15s & %6s & %6s " tab_header = "\\begin{tabular}{|l|r|r|r|r|" # what = r"$(x_\text{out} - x_\text{in})/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$" what = r"\Delta x/\sigma_{\text{in}}$, $\sigma_{\text{out}}/\sigma_{\text{in}}$" - headers_first = [ "", "$b$-only fit", "$s+b$ fit", "", "" ] + headers_first = ["", "$b$-only fit", "$s+b$ fit", "", ""] headers = ["name", what, what, r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] if options.workspace: tab_header += "r|" @@ -564,10 +564,11 @@ def getGraph(hist,shift): headers_first += [""] headers += [r"\Delta NLL"] fmtstring += " \\\\" - print( tab_header + "} \\hline " ) + print(tab_header + "} \\hline ") if headers_first is not None: - print( fmtstring % tuple(headers_first) ) - print( fmtstring % tuple(headers), " \\hline") + print(fmtstring % tuple(headers_first)) + print(fmtstring % tuple(headers), " \\hline") + elif options.format == "twiki": pmsub = (r"(\S+) \+/- (\S+)", r"\1 ± \2") sigsub = ("sig", r"σ") @@ -579,7 +580,7 @@ def getGraph(hist,shift): print(" option '--skipFitB' set true. $b$-only Fit is just a copy of the $s+b$ fit") if options.pullDef: fmtstring = "| %-40s | %-30s | %-30s | %-15s | %-15s | %-15s |" - header = "| *name* | *b-only fit pull* | *s+b fit pull* | *corr.* | *approx. impact* |" + header = "| *name* | *b-only fit pull* | *s+b fit pull* | *corr.* | *approx. impact* |" elif options.absolute_values: fmtstring = "| %-40s | %-15s | %-30s | %-30s | %-15s | %-15s |" header = "| *name* | *pre fit* | *b-only fit* | *s+b fit* | *corr.* | *approx. impact* |" @@ -609,7 +610,7 @@ def getGraph(hist,shift): """ ) if options.pullDef: - header = "nuisancebackground fit pull signal fit pullρ(μ, θ)I(μ, θ)" + header = "nuisancebackground fit pull signal fit pullρ(μ, θ)I(μ, θ)" fmtstring = "%-40s %-30s %-30s %-15s %-15s " elif options.absolute_values: header = "nuisancepre fitbackground fit signal fitcorrelation" @@ -680,7 +681,7 @@ def getGraph(hist,shift): histogram = ROOT.TH1F("pulls", "Pulls", 60, -3, 3) for pull in pulls: histogram.Fill(pull) - canvas = ROOT.TCanvas("asdf","asdf", 800, 800) + canvas = ROOT.TCanvas("asdf", "asdf", 800, 800) if options.pullDef: histogram.GetXaxis().SetTitle("pull") else: @@ -761,7 +762,6 @@ def getGraph(hist,shift): fout.WriteTObject(canvas_pferrs) - if options.workspace: canvas_dnll = ROOT.TCanvas("dnlls_%s" % idx, "dnlls_%s" % idx, 900, 600) hist_dnll[idx].SetStats(0) @@ -788,5 +788,5 @@ def getGraph(hist,shift): leg.Draw() latex = ROOT.TLatex() posx = 4 if hist_dnll[idx].GetNbinsX() > 10 else 0.5 - latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll) ) + latex.DrawLatex(posx, 1.5, "Total #Delta NLL %.3f" % (cdnll)) fout.WriteTObject(canvas_dnll) From 4af3707776e7e1ac9f2d098c57c8863167f8304c Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 09:34:09 +0100 Subject: [PATCH 07/11] Missed some linting issues --- scripts/diffNuisances.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index 0b82707db40..7001ef0276b 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -492,7 +492,7 @@ def getGraph(hist,shift): nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx) hist_dnll += [ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins)] hist_cdnll += [ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins)] - dnlls.sort(key= lambda x: -x[1]) + dnlls.sort(key = lambda x: -x[1]) cdnll = 0 for idx, (nm, val) in enumerate(dnlls): hist_idx = idx // options.max_nuis @@ -530,7 +530,7 @@ def getGraph(hist,shift): if options.workspace: fmtstring = fmtstring + " %10s" headers += ["dnll"] - print( fmtstring % tuple(headers) ) + print(fmtstring % tuple(headers)) elif options.format == "latex": pmsub = (r"(\S+) \+/- (\S+)", r"$\1 \\pm \2$") @@ -544,12 +544,12 @@ def getGraph(hist,shift): print(" option '--skipFitB' set true. $b$-only Fit is just a copy of the $s+b$ fit") if options.pullDef: fmtstring = "%-40s & %30s & %30s & %6s & %6s " - tab_header="\\begin{tabular}{|l|r|r|r|r|" - headers = ["name", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] + tab_header = "\\begin{tabular}{|l|r|r|r|r|" + headers = ["name", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] elif options.absolute_values: fmtstring = "%-40s & %15s & %30s & %30s & %6s & %6s " tab_header = "\\begin{tabular}{|l|r|r|r|r|r|" - headers = ["name", "pre fit", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] + headers = ["name", "pre fit", "$b$-only fit pull", "$s+b$ fit pull", r"$\rho(\theta, \mu)$", r"I(\theta, \mu)"] else: fmtstring = "%-40s & %15s & %15s & %6s & %6s " tab_header = "\\begin{tabular}{|l|r|r|r|r|" From d38fa92dccf511458b040bb5db86fb7c84767582 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 09:36:52 +0100 Subject: [PATCH 08/11] really fix linting --- scripts/diffNuisances.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index 7001ef0276b..defdb7ae8a8 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -492,7 +492,7 @@ def getGraph(hist,shift): nbins = min(len_dnll, options.max_nuis, len_dnll - options.max_nuis * hist_idx) hist_dnll += [ROOT.TH1F("dnll %s" % hist_idx, "delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins)] hist_cdnll += [ROOT.TH1F("cdnll %s" % hist_idx, "cumulative delta log-likelihoods ;;%s %s" % (title, hist_idx), nbins, 0, nbins)] - dnlls.sort(key = lambda x: -x[1]) + dnlls.sort(key=lambda x: -x[1]) cdnll = 0 for idx, (nm, val) in enumerate(dnlls): hist_idx = idx // options.max_nuis From 4436e4bc9b8628959c803df5306906d1b9af9b3d Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 11:48:50 +0100 Subject: [PATCH 09/11] more verbose help, (also linting) --- scripts/diffNuisances.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index defdb7ae8a8..88d74f8cf65 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -149,7 +149,7 @@ dest="workspace", default="", type="string", - help="Workspace to use for evaluating NLL differences." + help="Workspace to use for evaluating NLL differences. If no workspace is passed, NLL differences won't be calculated.", ) parser.add_option( "", @@ -157,7 +157,7 @@ dest="max_nuis", default=65, type="int", - help="Maximum nuisances for a single plot." + help="Maximum nuisances for a single plot. If more than the specified number of nuisance parameters exist they will be split across multiple plots.", ) (options, args) = parser.parse_args() @@ -267,13 +267,15 @@ def getGraph(hist,shift): nbins = min(np_count, options.max_nuis, np_count - (options.max_nuis * idx)) hist_fit_b += [ROOT.TH1F("fit_b %s" % idx, "B-only fit Nuisances;;%s %s" % (title, idx), nbins, 0, nbins)] hist_fit_s += [ROOT.TH1F("fit_s %s" % idx, "S+B fit Nuisances ;;%s %s" % (title, idx), nbins, 0, nbins)] - hist_prefit += [ROOT.TH1F( - "prefit_nuisancs %s" % idx, - "Prefit Nuisances ;;%s %s" % (title, idx), - nbins, - 0, - nbins, - )] + hist_prefit += [ + ROOT.TH1F( + "prefit_nuisancs %s" % idx, + "Prefit Nuisances ;;%s %s" % (title, idx), + nbins, + 0, + nbins, + ) + ] # Store also the *asymmetric* uncertainties gr_fit_b += [ROOT.TGraphAsymmErrors()] gr_fit_b[idx].SetTitle("fit_b_g % s" % idx) @@ -617,7 +619,10 @@ def getGraph(hist,shift): fmtstring = "%-40s %-15s %-30s %-30s %-15s %-15s " else: what = "Δx/σin, σout/σin" - header = "nuisancebackground fit
%s signal fit
%sρ(μ, θ)I(μ, θ)" % (what, what) + header = "nuisancebackground fit
%s signal fit
%sρ(μ, θ)I(μ, θ)" % ( + what, + what, + ) fmtstring = "%-40s %-15s %-15s %-15s %-15s " if options.workspace: header += "ΔNLL" From c5c85b47da5dd50c9695419d28df78c2a4263dc9 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 11:50:37 +0100 Subject: [PATCH 10/11] remove whitespace --- scripts/diffNuisances.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index 88d74f8cf65..4f6e5f24c2b 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -620,7 +620,7 @@ def getGraph(hist,shift): else: what = "Δx/σin, σout/σin" header = "nuisancebackground fit
%s signal fit
%sρ(μ, θ)I(μ, θ)" % ( - what, + what, what, ) fmtstring = "%-40s %-15s %-15s %-15s %-15s " From 4466ce5e45d848abe912a0d93a498e06ed0af912 Mon Sep 17 00:00:00 2001 From: Kyle Cormier Date: Wed, 22 Mar 2023 11:53:48 +0100 Subject: [PATCH 11/11] :( okay, should be last linting issue Will improve my checks in future --- scripts/diffNuisances.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/diffNuisances.py b/scripts/diffNuisances.py index 4f6e5f24c2b..aa88db28b1e 100755 --- a/scripts/diffNuisances.py +++ b/scripts/diffNuisances.py @@ -485,7 +485,7 @@ def getGraph(hist,shift): # end of loop over all fitted parameters -#fill dnll and cumulative dnll plots +# fill dnll and cumulative dnll plots if options.plotfile and options.workspace: len_dnll = len(dnlls) hist_dnll = []