See StandardHypoTestInvDemo.C in RooStats examples for more information.
Standard tutorial macro for performing an inverted hypothesis test.
This macro will perform a scan of the p-values for computing the limit
import ROOT
from ROOT import RooStats, RooFit, TFile, RooArgSet
import rootnotes
import rootprint
import utils
plotHypoTestResult = True
useProof = True
optimize = False
writeResult = False
nworkers = 2
def RunInverter(w, modelSBName, modelBName, \
dataName, calculatorType, testStatType, \
npoints, poimin, poimax, \
ntoys, useCls ):
print "Running HypoTestInverter on the workspace", w.GetName()
w.Print()
data = w.data(dataName)
if (not data):
print "StandardHypoTestDemo","Not existing data",dataName
return 0
else:
print "Using data set", dataName
# get models from WS
# get the modelConfig out of the file
bModel = w.obj(modelBName)
sbModel = w.obj(modelSBName)
if(not sbModel.GetSnapshot()):
print "StandardHypoTestInvDemo","Model",modelSBName,"has no snapshot - make one using model poi"
sbModel.SetSnapshot( sbModel.GetParametersOfInterest() )
# if(not bModel or bModel and sbModel):
print "StandardHypoTestInvDemo","The background model",modelBName,"does not exist"
print "StandardHypoTestInvDemo","Copy it from ModelConfig",modelSBName,"and set POI to zero"
bModel = sbModel.Clone();
bModel.SetName(modelSBName+"_with_poi_0")
var = bModel.GetParametersOfInterest().first()
if(not var): return 0
oldval = var.getVal()
var.setVal(0)
bModel.SetSnapshot( RooArgSet(var) )
var.setVal(oldval)
# Define profile likelihood and likelihood ratio
# ----------------------------------------------
slrts = RooStats.SimpleLikelihoodRatioTestStat(sbModel.GetPdf(),bModel.GetPdf())
if(sbModel.GetSnapshot()): slrts.SetNullParameters(sbModel.GetSnapshot())
if(bModel.GetSnapshot()): slrts.SetAltParameters(bModel.GetSnapshot())
# ratio of profile likelihood - need to pass snapshot for the alt
ropl = RooStats.RatioOfProfiledLikelihoodsTestStat(sbModel.GetPdf(), bModel.GetPdf(), bModel.GetSnapshot())
ropl.SetSubtractMLE(False)
profll = RooStats.ProfileLikelihoodTestStat(sbModel.GetPdf())
if(testStatType == 3): profll.SetOneSided(1)
if(optimize): profll.SetReuseNLL(True)
testStat = slrts
if(testStatType == 1): testStat = ropl
if(testStatType == 2 or testStatType == 3): testStat = profll
if(calculatorType == 0): hc = RooStats.FrequentistCalculator(data, bModel, sbModel)
else: hc = RooStats.HybridCalculator(data, bModel, sbModel)
# Generate toy MC
# ---------------
toymcs = hc.GetTestStatSampler()
# toymcs.SetNEventsPerToy(1)
toymcs.SetTestStatistic(testStat)
if(optimize): toymcs.SetUseMultiGen(True)
if(calculatorType == 1):
hhc = hc
hhc.SetToys(ntoys,ntoys)
# check for nuisance prior pdf
if(bModel.GetPriorPdf() and sbModel.GetPriorPdf()):
hhc.ForcePriorNuisanceAlt(bModel.GetPriorPdf())
hhc.ForcePriorNuisanceNull(sbModel.GetPriorPdf())
else:
if(bModel.GetNuisanceParameters() or sbModel.GetNuisanceParameters() ):
print "StandardHypoTestInvDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified"
return 0
else:
hc.SetToys(ntoys,ntoys)
# Get the result
# RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
# TStopwatch tw; tw.Start();
poiSet = sbModel.GetParametersOfInterest()
poi = poiSet.first()
# fit the data first
sbModel.GetPdf().fitTo(data)
poihat = poi.getVal()
calc = RooStats.HypoTestInverter(hc)
calc.SetConfidenceLevel(0.95)
calc.UseCLs(useCls)
calc.SetVerbose(True)
# can speed up using proof-lite
if(useProof and nworkers > 1):
pc = RooStats.ProofConfig(w, nworkers, "", ROOT.kFALSE);
toymcs.SetProofConfig(pc) # enable proof
if(npoints > 0):
if(poimin >= poimax):
# if no min/max given scan between MLE and +4 sigma
poimin = int(poihat)
poimax = int(poihat + 4 * poi.getError())
print "Doing a fixed scan in interval :",poimin,",",poimax
calc.SetFixedScan(npoints,poimin,poimax)
else:
# poi->setMax(10*int( (poihat+ 10 *poi->getError() )/10 ) );
print "Doing an automatic scan in interval :", poi.getMin(), ",", poi.getMax()
r = calc.GetInterval()
return r
%%rootprint
wsName = "combined"
modelSBName = "ModelConfig"
modelBName = ""
dataName = "obsData"
calculatorType = 0
testStatType = 3
useCls = True
npoints = 5
poimin = 0
poimax = 5
ntoys=1000
# Use standard file generated with HistFactory
inputFile = TFile("/afs/cern.ch/user/d/demattia/public/HATS/example_combined_GaussExample_model.root")
w = inputFile.Get(wsName)
print w, "\t"
r = RunInverter(w, modelSBName, modelBName, dataName, calculatorType, testStatType, npoints, poimin, poimax, ntoys, useCls)
upperLimit = r.UpperLimit()
ulError = r.UpperLimitEstimatedError()
print "The computed upper limit is:", upperLimit, "+/-", ulError
<ROOT.RooWorkspace object ("combined") at 0x10496230> Running HypoTestInverter on the workspace combined Using data set obsData StandardHypoTestInvDemo Model ModelConfig has no snapshot - make one using model poi StandardHypoTestInvDemo The background model does not exist StandardHypoTestInvDemo Copy it from ModelConfig ModelConfig and set POI to zero Doing a fixed scan in interval : 0 , 5
Plotting
nEntries = r.ArraySize()
c1 = rootnotes.canvas("roostats1", (800, 800))
typeName = ("Frequentist" if (calculatorType == 0) else "Hybrid")
resultName = (w.GetName() if w else r.GetName())
plotTitle = typeName+" CL Scan for workspace "+resultName
plot = RooStats.HypoTestInverterPlot("HTI_Result_Plot",plotTitle,r)
plot.Draw("CLb 2CL") # plot all and Clb
c1.SetLogy(False)
c2 = rootnotes.canvas("roostats2", (800, 800))
if(plotHypoTestResult):
# c2.Divide( 2, ROOT.TMath.Ceil(nEntries/2))
c2.Divide( 2, nEntries/2+1)
for i in range(nEntries):
c2.cd(i+1)
pl = plot.MakeTestStatPlot(i)
pl.SetLogYaxis(True)
pl.Draw()
print " expected limit (median)", r.GetExpectedUpperLimit(0)
print " expected limit (-1 sig)", r.GetExpectedUpperLimit(-1)
print " expected limit (+1 sig)", r.GetExpectedUpperLimit(1)
if (w != ROOT.NULL and writeResult):
# write to a file the results
calcType = ("Freq" if (calculatorType == 0) else "Hybr")
limitType = ("CLs" if (useCls) else "Cls+b")
scanType = ("auto" if (npoints < 0) else "grid")
resultFileName = calcType+"_"+limitType+"_"+scanType+"_"+testStatType
resultFileName += fileName
fileOut = TFile(resultFileName,"RECREATE")
r.Write()
fileOut.Close()
TCanvas::Constructor:0: RuntimeWarning: Deleting canvas with same name: roostats1 TROOT::Append:0: RuntimeWarning: Replacing existing TH1: CLs_observed (Potential memory leak). TCanvas::Constructor:0: RuntimeWarning: Deleting canvas with same name: roostats2 TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_with_poi_0_SigXsecOverSM_0.00 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_SigXsecOverSM_0.00 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_with_poi_0_SigXsecOverSM_0.75 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_SigXsecOverSM_0.75 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_with_poi_0_SigXsecOverSM_1.50 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_SigXsecOverSM_1.50 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_with_poi_0_SigXsecOverSM_2.25 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_SigXsecOverSM_2.25 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_with_poi_0_SigXsecOverSM_3.00 (Potential memory leak). TROOT::Append:0: RuntimeWarning: Replacing existing TH1: ModelConfig_SigXsecOverSM_3.00 (Potential memory leak).
c1
c2