107 infile="", workspaceName="combined", modelConfigName="ModelConfig", dataName="obsData"
110 confidenceLevel = 0.95
119 filename =
"results/example_combined_GaussExample_model.root"
120 fileExist =
not ROOT.gSystem.AccessPathName(filename)
124 print(f
"will run standard hist2workspace example")
125 ROOT.gROOT.ProcessLine(
".not prepareHistFactory .")
126 ROOT.gROOT.ProcessLine(
".not hist2workspace config/example.xml")
127 print(f
"\n\n---------------------")
128 print(f
"Done creating example input")
129 print(f
"---------------------\n\n")
135 file = ROOT.TFile.Open(filename)
139 print(f
"StandardRooStatsDemoMacro: Input file {filename} is not found")
146 w = file.Get(workspaceName)
153 print(f
"workspace not found")
157 mc = w.obj(modelConfigName)
160 data = w.data(dataName)
163 if not data
or not mc:
165 print(f
"data or ModelConfig was not found")
172 firstPOI = mc.GetParametersOfInterest().first()
183 fc = ROOT.RooStats.FeldmanCousins(data, mc)
184 fc.SetConfidenceLevel(confidenceLevel)
185 fc.AdditionalNToysFactor(
189 fc.SetNBins(nPointsToScan)
190 fc.CreateConfBelt(
True)
201 toymcsampler = fc.GetTestStatSampler()
202 testStat = toymcsampler.GetTestStatistic()
203 testStat.SetOneSided(
True)
212 if not mc.GetPdf().canBeExtended():
213 if data.numEntries() == 1:
214 fc.FluctuateNumDataEntries(
False)
216 print(f
"Not sure what to do about this model")
218 if mc.GetGlobalObservables():
219 print(f
"will use global observables for unconditional ensemble")
220 mc.GetGlobalObservables().
Print()
221 toymcsampler.SetGlobalObservables(mc.GetGlobalObservables())
224 interval = fc.GetInterval()
225 belt = fc.GetConfidenceBelt()
229 f
"\n95% interval on {firstPOI.GetName()} is : [{interval.LowerLimit(firstPOI)}, {interval.UpperLimit(firstPOI)} ]"
233 tmpPOI = ROOT.RooArgSet(firstPOI)
234 observedUL = interval.UpperLimit(firstPOI)
235 firstPOI.setVal(observedUL)
236 obsTSatObsUL = fc.GetTestStatSampler().EvaluateTestStatistic(data, tmpPOI)
239 parameterScan = fc.GetPointsToScan()
240 tmpPoint = ROOT.RooArgSet()
243 histOfThresholds = ROOT.TH1F(
244 "histOfThresholds",
"", parameterScan.numEntries(), firstPOI.getMin(), firstPOI.getMax()
246 histOfThresholds.GetXaxis().SetTitle(firstPOI.GetName())
247 histOfThresholds.GetYaxis().SetTitle(
"Threshold")
252 for i
in range(parameterScan.numEntries()):
253 tmpPoint = parameterScan.get(i).clone(
"temp")
255 arMax = belt.GetAcceptanceRegionMax(tmpPoint)
256 poiVal = tmpPoint.getRealValue(firstPOI.GetName())
257 histOfThresholds.Fill(poiVal, arMax)
262 histOfThresholds.SetMinimum(0)
263 histOfThresholds.Draw()
272 nll = mc.GetPdf().createNLL(data)
273 profile = nll.createProfile(mc.GetParametersOfInterest())
276 poiAndNuisance = ROOT.RooArgSet()
277 if mc.GetNuisanceParameters():
278 poiAndNuisance.add(mc.GetNuisanceParameters())
279 poiAndNuisance.add(mc.GetParametersOfInterest())
280 w.saveSnapshot(
"paramsToGenerateData", poiAndNuisance)
281 paramsToGenerateData = poiAndNuisance.snapshot()
282 print(f
"\nWill use these parameter points to generate pseudo data for bkg only")
283 paramsToGenerateData.Print(
"v")
285 unconditionalObs = ROOT.RooArgSet()
286 unconditionalObs.add(mc.GetObservables())
287 unconditionalObs.add(mc.GetGlobalObservables())
293 histOfUL = ROOT.TH1F(
"histOfUL",
"", 100, 0, firstPOI.getMax())
294 histOfUL.GetXaxis().SetTitle(
"Upper Limit (background only)")
295 histOfUL.GetYaxis().SetTitle(
"Entries")
296 for imc
in range(nToyMC):
300 w.loadSnapshot(
"paramsToGenerateData")
303 toyData = ROOT.RooDataSet()
309 if not mc.GetPdf().canBeExtended():
310 if data.numEntries() == 1:
311 toyData = mc.GetPdf().generate(mc.GetObservables(), 1)
313 print(f
"Not sure what to do about this model")
316 toyData = mc.GetPdf().generate(mc.GetObservables(), Extended=
True)
324 one = mc.GetPdf().generate(mc.GetGlobalObservables(), 1)
326 allVars = mc.GetPdf().getVariables()
327 allVars.assign(values)
333 for tt
in simPdf.indexCat():
339 pdftmp = simPdf.getPdf(str(catName))
342 globtmp = pdftmp.getObservables(mc.GetGlobalObservables())
343 tmp = pdftmp.generate(globtmp, 1)
346 globtmp.assign(tmp.get(0))
356 firstPOI.setVal(observedUL)
357 toyTSatObsUL = fc.GetTestStatSampler().EvaluateTestStatistic(toyData, tmpPOI)
360 if obsTSatObsUL < toyTSatObsUL:
361 CLb += (1.0) / nToyMC
362 if obsTSatObsUL <= toyTSatObsUL:
363 CLbinclusive += (1.0) / nToyMC
366 thisUL = ROOT.Double_t(0)
367 for i
in range(parameterScan.numEntries()):
368 tmpPoint = parameterScan.get(i).clone(
"temp")
369 arMax = belt.GetAcceptanceRegionMax(tmpPoint)
370 firstPOI.setVal(tmpPoint.getRealValue(firstPOI.GetName()))
372 thisTS = fc.GetTestStatSampler().EvaluateTestStatistic(toyData, tmpPOI)
378 thisUL = firstPOI.getVal()
384 # loop over points in belt to find upper limit for this toy data
386 for i in range(histOfThresholds.GetNbinsX() ++i)
387 tmpPoint = (RooArgSet) parameterScan.get(i).clone("temp")
388 print("---------------- ", i)
390 print(f"from hist ", histOfThresholds.GetBinCenter(i+1) )
391 arMax = histOfThresholds.GetBinContent(i+1)
392 # cout << " threshold from Hist = aMax " << arMax<<endl;
393 # double arMax2 = belt->GetAcceptanceRegionMax(*tmpPoint);
394 # cout << "from scan arMax2 = "<< arMax2 << endl; # not the same due to TH1F not TH1D
395 # cout << "scan - hist" << arMax2-arMax << endl;
396 firstPOI.setVal( histOfThresholds.GetBinCenter(i+1))
397 # double thisTS = profile->getVal();
398 thisTS = fc.GetTestStatSampler().EvaluateTestStatistic(toyData,tmpPOI)
400 # cout << "poi = " << firstPOI->getVal()
401 # = ROOT.Double_t() << " max is " << arMax << " this profile = " << thisTS << endl;
402 # cout << "thisTS = " << thisTS<<endl;
404 # NOTE: need to add a small epsilon term for single precision vs. double precision
405# if(thisTS<=arMax + 1e-7){
406# thisUL = firstPOI->getVal();
415 histOfUL.Fill(thisUL)
425 c1.SaveAs(
"OneSidedFrequentistUpperLimitWithBands.png")
430 SamplingDistPlot sampPlot
432 tmpPoint = (RooArgSet) parameterScan.get(indexInScan).clone("temp")
433 firstPOI.setVal( tmpPoint.getRealValue(firstPOI.GetName()) )
434 toymcsampler.SetParametersForTestStat(tmpPOI)
435 samp = toymcsampler.GetSamplingDistribution(tmpPoint)
436 sampPlot.AddSamplingDistribution(samp)
441 bins = histOfUL.GetIntegral()
442 cumulative = histOfUL.Clone(
"cumulative")
443 cumulative.SetContent(bins)
444 band2sigDown = band1sigDown = bandMedian = band1sigUp = band2sigUp = ROOT.Double_t()
445 for i
in range(cumulative.GetNbinsX()):
446 if bins[i] < ROOT.RooStats.SignificanceToPValue(2):
447 band2sigDown = cumulative.GetBinCenter(i)
448 if bins[i] < ROOT.RooStats.SignificanceToPValue(1):
449 band1sigDown = cumulative.GetBinCenter(i)
451 bandMedian = cumulative.GetBinCenter(i)
452 if bins[i] < ROOT.RooStats.SignificanceToPValue(-1):
453 band1sigUp = cumulative.GetBinCenter(i)
454 if bins[i] < ROOT.RooStats.SignificanceToPValue(-2):
455 band2sigUp = cumulative.GetBinCenter(i)
457 print(f
"-2 sigma band ", band2sigDown)
458 print(f
"-1 sigma band {band1sigDown} [Power Constraint)]")
459 print(f
"median of band ", bandMedian)
460 print(f
"+1 sigma band ", band1sigUp)
461 print(f
"+2 sigma band ", band2sigUp)
464 print(f
"\nObserved 95% upper-limit ", interval.UpperLimit(firstPOI))
465 print(f
"CLb strict [P(toy>obs|0)] for observed 95% upper-limit ", CLb)
466 print(
"inclusive [P(toy>=obs|0)] for observed 95% upper-limit ", CLbinclusive)
470 infile=
"", workspaceName=
"combined", modelConfigName=
"ModelConfig", dataName=
"obsData"
void Print(GNN_Data &d, std::string txt="")