From f2b020379f50e9bc3b0800b5751eef0794292265 Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Fri, 14 Mar 2025 19:39:39 +0100
Subject: [PATCH] [RF][RS] Fix StandardHistFactoryPlotsWithCategories tutorial
The StandardHistFactoryPlotsWithCategories tutorial assumed that the
nominal and plus-minus one sigma values of all nuisance parameters was
zero, plus one, and minus one respectively. This is not true, especially
for `gamma` parameters, where the nominal value is one.
The actual values and uncertainties should be taken from the parameters
with `getVal()` and `getError()`. This also fixes the output plot in the
tutorial for the gamma parameters:
https://root.cern/doc/master/StandardHistFactoryPlotsWithCategories_8C.html
---
.../StandardHistFactoryPlotsWithCategories.C | 52 ++++++-------------
.../StandardHistFactoryPlotsWithCategories.py | 43 +++++----------
2 files changed, 27 insertions(+), 68 deletions(-)
diff --git a/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.C b/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.C
index ea7983100f0a6..13e248619b03f 100644
--- a/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.C
+++ b/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.C
@@ -176,14 +176,14 @@ void StandardHistFactoryPlotsWithCategories(const char *infile = "", const char
RooPlot *frame = obs->frame();
frame->SetYTitle(var->GetName());
data->plotOn(frame, MarkerSize(1));
- var->setVal(0);
+ const double value = var->getVal();
mc->GetPdf()->plotOn(frame, LineWidth(1.));
- var->setVal(1);
+ var->setVal(value + var->getError());
mc->GetPdf()->plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1));
- var->setVal(-1);
+ var->setVal(value - var->getError());
mc->GetPdf()->plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(1));
frameList.push_back(frame);
- var->setVal(0);
+ var->setVal(value);
}
} else {
@@ -223,57 +223,35 @@ void StandardHistFactoryPlotsWithCategories(const char *infile = "", const char
Double_t normCount =
data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str()));
- if (strcmp(var->GetName(), "Lumi") == 0) {
- cout << "working on lumi" << endl;
- var->setVal(w->var("nominalLumi")->getVal());
- var->Print();
- } else {
- var->setVal(0);
- }
+ // remember the nominal value
+ const double value = var->getVal();
+
// w->allVars().Print("v");
// mc->GetNuisanceParameters()->Print("v");
// pdftmp->plotOn(frame,LineWidth(2.));
// mc->GetPdf()->plotOn(frame,LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
// pdftmp->plotOn(frame,LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
- pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent));
-
- if (strcmp(var->GetName(), "Lumi") == 0) {
- cout << "working on lumi" << endl;
- var->setVal(w->var("nominalLumi")->getVal() + 0.05);
- var->Print();
- } else {
- var->setVal(nSigmaToVary);
- }
+ pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent)); // nominal
+
+ var->setVal(value + nSigmaToVary * var->getError());
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2));
// mc->GetPdf()->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kRed), LineStyle(kDashed),
- Normalization(normCount, RooAbsReal::NumEvent));
-
- if (strcmp(var->GetName(), "Lumi") == 0) {
- cout << "working on lumi" << endl;
- var->setVal(w->var("nominalLumi")->getVal() - 0.05);
- var->Print();
- } else {
- var->setVal(-nSigmaToVary);
- }
+ Normalization(normCount, RooAbsReal::NumEvent)); // +n sigma
+
+ var->setVal(value - nSigmaToVary * var->getError());
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2));
// mc->GetPdf()->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,catName.c_str()),ProjWData(*data));
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,catName.c_str()),ProjWData(*data));
normCount = pdftmp->expectedEvents(*obs);
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kGreen), LineStyle(kDashed),
- Normalization(normCount, RooAbsReal::NumEvent));
+ Normalization(normCount, RooAbsReal::NumEvent)); // -n sigma
// set them back to normal
- if (strcmp(var->GetName(), "Lumi") == 0) {
- cout << "working on lumi" << endl;
- var->setVal(w->var("nominalLumi")->getVal());
- var->Print();
- } else {
- var->setVal(0);
- }
+ var->setVal(value);
frameList.push_back(frame);
diff --git a/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.py b/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.py
index 61d86eef82560..4983e5117d76a 100644
--- a/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.py
+++ b/tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.py
@@ -154,14 +154,14 @@ def StandardHistFactoryPlotsWithCategories(
frame = obs.frame()
frame.SetYTitle(var.GetName())
data.plotOn(frame, MarkerSize(1))
- var.setVal(0)
+ value = var.getVal()
mc.GetPdf().plotOn(frame, LineWidth(1))
- var.setVal(1)
+ var.setVal(value + var.getError());
mc.GetPdf().plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1))
- var.setVal(-1)
+ var.setVal(value - var.getError());
mc.GetPdf().plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(1))
frameList.append(frame)
- var.setVal(0)
+ var.setVal(value)
else:
channelCat = simPdf.indexCat()
@@ -204,12 +204,8 @@ def StandardHistFactoryPlotsWithCategories(
normCount = data.sumEntries(cut)
- if str(var.GetName()) == "Lumi":
- print(f"working on lumi")
- var.setVal(w.var("nominalLumi").getVal())
- var.Print()
- else:
- var.setVal(0)
+ # remember the nominal value
+ value = var.getVal()
# w.allVars().Print("v")
# mc.GetNuisanceParameters().Print("v")
@@ -231,14 +227,9 @@ def StandardHistFactoryPlotsWithCategories(
# instead, we have to use obstmp
# normCount = pdftmp.expectedEvents(RooArgSet(obstmp)) #doesn´t work properly
normCount = pdftmp.expectedEvents(obstmp)
- pdftmp.plotOn(frame, ROOT.RooFit.Normalization(normCount, ROOT.RooAbsReal.NumEvent), LineWidth=2)
+ pdftmp.plotOn(frame, ROOT.RooFit.Normalization(normCount, ROOT.RooAbsReal.NumEvent), LineWidth=2) # nominal
- if str(var.GetName()) == "Lumi":
- print(f"working on lumi")
- var.setVal(w.var("nominalLumi").getVal() + 0.05)
- var.Print()
- else:
- var.setVal(nSigmaToVary)
+ var.setVal(value + nSigmaToVary * var.getError())
# pdftmp.plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2))
# mc.GetPdf().plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(channelCat,catName.c_str()),ProjWData(data))
@@ -250,14 +241,9 @@ def StandardHistFactoryPlotsWithCategories(
LineWidth=2,
LineColor="r",
LineStyle="--",
- )
+ ) # +n sigma
- if str(var.GetName()) == "Lumi":
- print(f"working on lumi")
- var.setVal(w.var("nominalLumi").getVal() - 0.05)
- var.Print()
- else:
- var.setVal(-nSigmaToVary)
+ var.setVal(value - nSigmaToVary * var.getError())
# pdftmp.plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2))
# mc.GetPdf().plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(channelCat,catName.c_str()),ProjWData(data))
@@ -269,15 +255,10 @@ def StandardHistFactoryPlotsWithCategories(
LineWidth=2,
LineColor="g",
LineStyle="--",
- )
+ ) # -n sigma
# set them back to normal
- if str(var.GetName()) == "Lumi":
- print(f"working on lumi")
- var.setVal(w.var("nominalLumi").getVal())
- var.Print()
- else:
- var.setVal(0)
+ var.setVal(value)
frameList.append(frame)