Skip to content

Commit ca5b243

Browse files
guitargeekguitargeek
committed
[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
1 parent 9f5dd4b commit ca5b243

File tree

2 files changed

+27
-68
lines changed

2 files changed

+27
-68
lines changed

‎tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.C

+15-37
Original file line numberDiff line numberDiff line change
@@ -176,14 +176,14 @@ void StandardHistFactoryPlotsWithCategories(const char *infile = "", const char
176176
RooPlot *frame = obs->frame();
177177
frame->SetYTitle(var->GetName());
178178
data->plotOn(frame, MarkerSize(1));
179-
var->setVal(0);
179+
const double value = var->getVal();
180180
mc->GetPdf()->plotOn(frame, LineWidth(1.));
181-
var->setVal(1);
181+
var->setVal(value + var->getError());
182182
mc->GetPdf()->plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1));
183-
var->setVal(-1);
183+
var->setVal(value - var->getError());
184184
mc->GetPdf()->plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(1));
185185
frameList.push_back(frame);
186-
var->setVal(0);
186+
var->setVal(value);
187187
}
188188

189189
} else {
@@ -223,57 +223,35 @@ void StandardHistFactoryPlotsWithCategories(const char *infile = "", const char
223223
Double_t normCount =
224224
data->sumEntries(Form("%s==%s::%s", channelCat->GetName(), channelCat->GetName(), catName.c_str()));
225225

226-
if (strcmp(var->GetName(), "Lumi") == 0) {
227-
cout << "working on lumi" << endl;
228-
var->setVal(w->var("nominalLumi")->getVal());
229-
var->Print();
230-
} else {
231-
var->setVal(0);
232-
}
226+
// remember the nominal value
227+
const double value = var->getVal();
228+
233229
// w->allVars().Print("v");
234230
// mc->GetNuisanceParameters()->Print("v");
235231
// pdftmp->plotOn(frame,LineWidth(2.));
236232
// mc->GetPdf()->plotOn(frame,LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
237233
// pdftmp->plotOn(frame,LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
238234
normCount = pdftmp->expectedEvents(*obs);
239-
pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent));
240-
241-
if (strcmp(var->GetName(), "Lumi") == 0) {
242-
cout << "working on lumi" << endl;
243-
var->setVal(w->var("nominalLumi")->getVal() + 0.05);
244-
var->Print();
245-
} else {
246-
var->setVal(nSigmaToVary);
247-
}
235+
pdftmp->plotOn(frame, LineWidth(2.), Normalization(normCount, RooAbsReal::NumEvent)); // nominal
236+
237+
var->setVal(value + nSigmaToVary * var->getError());
248238
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2));
249239
// mc->GetPdf()->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
250240
// pdftmp->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(*channelCat,catName.c_str()),ProjWData(*data));
251241
normCount = pdftmp->expectedEvents(*obs);
252242
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kRed), LineStyle(kDashed),
253-
Normalization(normCount, RooAbsReal::NumEvent));
254-
255-
if (strcmp(var->GetName(), "Lumi") == 0) {
256-
cout << "working on lumi" << endl;
257-
var->setVal(w->var("nominalLumi")->getVal() - 0.05);
258-
var->Print();
259-
} else {
260-
var->setVal(-nSigmaToVary);
261-
}
243+
Normalization(normCount, RooAbsReal::NumEvent)); // +n sigma
244+
245+
var->setVal(value - nSigmaToVary * var->getError());
262246
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2));
263247
// mc->GetPdf()->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,catName.c_str()),ProjWData(*data));
264248
// pdftmp->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(*channelCat,catName.c_str()),ProjWData(*data));
265249
normCount = pdftmp->expectedEvents(*obs);
266250
pdftmp->plotOn(frame, LineWidth(2.), LineColor(kGreen), LineStyle(kDashed),
267-
Normalization(normCount, RooAbsReal::NumEvent));
251+
Normalization(normCount, RooAbsReal::NumEvent)); // -n sigma
268252

269253
// set them back to normal
270-
if (strcmp(var->GetName(), "Lumi") == 0) {
271-
cout << "working on lumi" << endl;
272-
var->setVal(w->var("nominalLumi")->getVal());
273-
var->Print();
274-
} else {
275-
var->setVal(0);
276-
}
254+
var->setVal(value);
277255

278256
frameList.push_back(frame);
279257

‎tutorials/roofit/roostats/StandardHistFactoryPlotsWithCategories.py

+12-31
Original file line numberDiff line numberDiff line change
@@ -154,14 +154,14 @@ def StandardHistFactoryPlotsWithCategories(
154154
frame = obs.frame()
155155
frame.SetYTitle(var.GetName())
156156
data.plotOn(frame, MarkerSize(1))
157-
var.setVal(0)
157+
value = var.getVal()
158158
mc.GetPdf().plotOn(frame, LineWidth(1))
159-
var.setVal(1)
159+
var.setVal(value + var.getError());
160160
mc.GetPdf().plotOn(frame, LineColor(kRed), LineStyle(kDashed), LineWidth(1))
161-
var.setVal(-1)
161+
var.setVal(value - var.getError());
162162
mc.GetPdf().plotOn(frame, LineColor(kGreen), LineStyle(kDashed), LineWidth(1))
163163
frameList.append(frame)
164-
var.setVal(0)
164+
var.setVal(value)
165165

166166
else:
167167
channelCat = simPdf.indexCat()
@@ -204,12 +204,8 @@ def StandardHistFactoryPlotsWithCategories(
204204

205205
normCount = data.sumEntries(cut)
206206

207-
if str(var.GetName()) == "Lumi":
208-
print(f"working on lumi")
209-
var.setVal(w.var("nominalLumi").getVal())
210-
var.Print()
211-
else:
212-
var.setVal(0)
207+
# remember the nominal value
208+
value = var.getVal()
213209

214210
# w.allVars().Print("v")
215211
# mc.GetNuisanceParameters().Print("v")
@@ -231,14 +227,9 @@ def StandardHistFactoryPlotsWithCategories(
231227
# instead, we have to use obstmp
232228
# normCount = pdftmp.expectedEvents(RooArgSet(obstmp)) #doesn´t work properly
233229
normCount = pdftmp.expectedEvents(obstmp)
234-
pdftmp.plotOn(frame, ROOT.RooFit.Normalization(normCount, ROOT.RooAbsReal.NumEvent), LineWidth=2)
230+
pdftmp.plotOn(frame, ROOT.RooFit.Normalization(normCount, ROOT.RooAbsReal.NumEvent), LineWidth=2) # nominal
235231

236-
if str(var.GetName()) == "Lumi":
237-
print(f"working on lumi")
238-
var.setVal(w.var("nominalLumi").getVal() + 0.05)
239-
var.Print()
240-
else:
241-
var.setVal(nSigmaToVary)
232+
var.setVal(value + nSigmaToVary * var.getError())
242233

243234
# pdftmp.plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2))
244235
# mc.GetPdf().plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(2.),Slice(channelCat,catName.c_str()),ProjWData(data))
@@ -250,14 +241,9 @@ def StandardHistFactoryPlotsWithCategories(
250241
LineWidth=2,
251242
LineColor="r",
252243
LineStyle="--",
253-
)
244+
) # +n sigma
254245

255-
if str(var.GetName()) == "Lumi":
256-
print(f"working on lumi")
257-
var.setVal(w.var("nominalLumi").getVal() - 0.05)
258-
var.Print()
259-
else:
260-
var.setVal(-nSigmaToVary)
246+
var.setVal(value - nSigmaToVary * var.getError())
261247

262248
# pdftmp.plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2))
263249
# mc.GetPdf().plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(2),Slice(channelCat,catName.c_str()),ProjWData(data))
@@ -269,15 +255,10 @@ def StandardHistFactoryPlotsWithCategories(
269255
LineWidth=2,
270256
LineColor="g",
271257
LineStyle="--",
272-
)
258+
) # -n sigma
273259

274260
# set them back to normal
275-
if str(var.GetName()) == "Lumi":
276-
print(f"working on lumi")
277-
var.setVal(w.var("nominalLumi").getVal())
278-
var.Print()
279-
else:
280-
var.setVal(0)
261+
var.setVal(value)
281262

282263
frameList.append(frame)
283264

0 commit comments

Comments
 (0)