Limits: number counting experiment with uncertainty on both the background rate and signal efficiency.
The usage of a Confidence Interval Calculator to set a limit on the signal is illustrated
Author: Artem Busorgin, Kyle Cranmer (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, April 17, 2024 at 11:20 AM.
import ROOT
An example of setting a limit in a number counting experiment with uncertainty on background and signal
to time the macro
t = ROOT.TStopwatch()
t.Start()
wspace = ROOT.RooWorkspace()
wspace.factory(
"Poisson::countingModel(obs[150,0,300], " "sum(s[50,0,120]*ratioSigEff[1.,0,3.],b[100]*ratioBkgEff[1.,0.,3.]))"
) # counting model
wspace.factory("Gaussian::sigConstraint(gSigEff[1,0,3],ratioSigEff,0.05)") # 5% signal efficiency uncertainty
wspace.factory("Gaussian::bkgConstraint(gSigBkg[1,0,3],ratioBkgEff,0.2)") # 10% background efficiency uncertainty
wspace.factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)") # product of terms
wspace.Print()
modelWithConstraints = wspace["modelWithConstraints"] # get the model
obs = wspace["obs"] # get the observable
s = wspace["s"] # get the signal we care about
b = wspace["b"] # get the background and set it to a constant. Uncertainty included in ratioBkgEff
b.setConstant()
ratioSigEff = wspace["ratioSigEff"] # get uncertain parameter to constrain
ratioBkgEff = wspace["ratioBkgEff"] # get uncertain parameter to constrain
constrainedParams = {ratioSigEff, ratioBkgEff} # need to constrain these in the fit (should change default behavior)
gSigEff = wspace["gSigEff"] # global observables for signal efficiency
gSigBkg = wspace["gSigBkg"] # global obs for background efficiency
gSigEff.setConstant()
gSigBkg.setConstant()
RooWorkspace() contents variables --------- (b,gSigBkg,gSigEff,obs,ratioBkgEff,ratioSigEff,s) p.d.f.s ------- RooGaussian::bkgConstraint[ x=gSigBkg mean=ratioBkgEff sigma=0.2 ] = 1 RooPoisson::countingModel[ x=obs mean=countingModel_2 ] = 0.0325554 RooProdPdf::modelWithConstraints[ countingModel * sigConstraint * bkgConstraint ] = 0.0325554 RooGaussian::sigConstraint[ x=gSigEff mean=ratioSigEff sigma=0.05 ] = 1 functions -------- RooAddition::countingModel_2[ countingModel_2_[s_x_ratioSigEff] + countingModel_2_[b_x_ratioBkgEff] ] = 150 RooProduct::countingModel_2_[b_x_ratioBkgEff][ b * ratioBkgEff ] = 100 RooProduct::countingModel_2_[s_x_ratioSigEff][ s * ratioSigEff ] = 50
Create an example dataset with 160 observed events
obs.setVal(160.0)
dataOrig = ROOT.RooDataSet("exampleData", "exampleData", {obs})
dataOrig.add(obs)
not necessary
modelWithConstraints.fitTo(dataOrig, Constrain=constrainedParams, PrintLevel=-1)
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Minimization -- Including the following constraint terms in minimization: (sigConstraint,bkgConstraint) [#1] INFO:Minimization -- The global observables are not defined , normalize constraints with respect to the parameters (ratioBkgEff,ratioSigEff) [#1] INFO:Fitting -- RooAbsPdf::fitTo(modelWithConstraints) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2 [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_modelWithConstraints_exampleData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
Now let's make some confidence intervals for s, our parameter of interest
modelConfig = ROOT.RooStats.ModelConfig(wspace)
modelConfig.SetPdf(modelWithConstraints)
modelConfig.SetParametersOfInterest({s})
modelConfig.SetNuisanceParameters(constrainedParams)
modelConfig.SetObservables(obs)
modelConfig.SetGlobalObservables({gSigEff, gSigBkg})
modelConfig.SetName("ModelConfig")
wspace.Import(modelConfig)
wspace.Import(dataOrig)
wspace.SetName("w")
wspace.writeToFile("rs101_ws.root")
False
[#1] INFO:ObjectHandling -- RooWorkspace::import() importing dataset exampleData
Make sure we reference the data in the workspace from now on
data = wspace[dataOrig.GetName()]
First, let's use a Calculator based on the Profile Likelihood Ratio
plc = ROOT.RooStats.ProfileLikelihoodCalculator(data, modelConfig)
plc.SetTestSize(0.05)
lrinterval = plc.GetInterval()
[#1] INFO:InputArguments -- The deprecated RooFit::CloneData(1) option passed to createNLL() is ignored. [#1] INFO:Minimization -- Including the following constraint terms in minimization: (sigConstraint,bkgConstraint) [#1] INFO:Minimization -- The following global observables have been defined and their values are taken from the model: (gSigEff,gSigBkg) [#1] INFO:Fitting -- RooAbsPdf::fitTo(modelWithConstraints) fixing normalization set for coefficient determination to observables in data [#0] PROGRESS:Minimization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_modelWithConstraints_exampleData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#0] PROGRESS:Minimization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit2 / Migrad with strategy 1 [#1] INFO:Minimization -- RooFitResult: minimized FCN value: 0.689753, estimated distance to minimum: 2.26046e-16 covariance matrix quality: Full, accurate covariance matrix Status : MINIMIZE=0 Floating Parameter FinalValue +/- Error -------------------- -------------------------- ratioBkgEff 1.0000e+00 +/- 1.99e-01 ratioSigEff 1.0000e+00 +/- 5.00e-02 s 6.0000e+01 +/- 2.32e+01
Let's make a plot
dataCanvas = ROOT.TCanvas("dataCanvas")
dataCanvas.Divide(2, 1)
dataCanvas.cd(1)
plotInt = ROOT.RooStats.LikelihoodIntervalPlot(lrinterval)
plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S")
plotInt.Draw()
[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[s]) Creating instance of MINUIT [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_modelWithConstraints_exampleData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[s]) determining minimum likelihood for current configurations w.r.t all observable [#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[s]) minimum found at (s=60) . [#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[s]) Creating instance of MINUIT [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_modelWithConstraints_exampleData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[s]) determining minimum likelihood for current configurations w.r.t all observable [#0] ERROR:InputArguments -- RooArgSet::checkForDup: ERROR argument with name s is already in this set [#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[s]) minimum found at (s=60.0341) ..........................................................................................................................................................................................................
Second, use a Calculator based on the Feldman Cousins technique
fc = ROOT.RooStats.FeldmanCousins(data, modelConfig)
fc.UseAdaptiveSampling(True)
fc.FluctuateNumDataEntries(False) # number counting analysis: dataset always has 1 entry with N events observed
fc.SetNBins(100) # number of points to test per parameter
fc.SetTestSize(0.05)
fc.SaveBeltToFile(True) # optional
fcint = fc.GetInterval()
fit = modelWithConstraints.fitTo(data, Save=True, PrintLevel=-1)
=== Using the following for ModelConfig === Observables: RooArgSet:: = (obs) Parameters of Interest: RooArgSet:: = (s) Nuisance Parameters: RooArgSet:: = (ratioBkgEff,ratioSigEff) Global Observables: RooArgSet:: = (gSigEff,gSigBkg) PDF: RooProdPdf::modelWithConstraints[ countingModel * sigConstraint * bkgConstraint ] = 0.0036613 FeldmanCousins: ntoys per point: adaptive FeldmanCousins: nEvents per toy will not fluctuate, will always be 1 FeldmanCousins: Model has nuisance parameters, will do profile construction FeldmanCousins: # points to test = 100 NeymanConstruction: Prog: 1/100 total MC = 80 this test stat = 3.21178 s=0.6 ratioBkgEff=1.43644 ratioSigEff=0.999975 [-inf, 4.68588] in interval = 1 NeymanConstruction: Prog: 2/100 total MC = 80 this test stat = 3.08184 s=1.8 ratioBkgEff=1.42694 ratioSigEff=1.00018 [-inf, 3.7189] in interval = 1 NeymanConstruction: Prog: 3/100 total MC = 80 this test stat = 2.95515 s=3 ratioBkgEff=1.41769 ratioSigEff=1.00055 [-inf, 3.91249] in interval = 1 NeymanConstruction: Prog: 4/100 total MC = 80 this test stat = 2.83091 s=4.2 ratioBkgEff=1.40849 ratioSigEff=1.00106 [-inf, 3.65947] in interval = 1 NeymanConstruction: Prog: 5/100 total MC = 80 this test stat = 2.7092 s=5.4 ratioBkgEff=1.39955 ratioSigEff=1.00133 [-inf, 3.10323] in interval = 1 NeymanConstruction: Prog: 6/100 total MC = 80 this test stat = 2.59009 s=6.6 ratioBkgEff=1.3906 ratioSigEff=1.00159 [-inf, 2.97611] in interval = 1 NeymanConstruction: Prog: 7/100 total MC = 240 this test stat = 2.47431 s=7.8 ratioBkgEff=1.38084 ratioSigEff=1.00102 [-inf, 3.25381] in interval = 1 NeymanConstruction: Prog: 8/100 total MC = 80 this test stat = 2.36072 s=9 ratioBkgEff=1.37175 ratioSigEff=1.00134 [-inf, 2.92353] in interval = 1 NeymanConstruction: Prog: 9/100 total MC = 80 this test stat = 2.24982 s=10.2 ratioBkgEff=1.36272 ratioSigEff=1.00166 [-inf, 3.4069] in interval = 1 NeymanConstruction: Prog: 10/100 total MC = 80 this test stat = 2.14164 s=11.4 ratioBkgEff=1.35374 ratioSigEff=1.00196 [-inf, 3.07011] in interval = 1 NeymanConstruction: Prog: 11/100 total MC = 80 this test stat = 2.03623 s=12.6 ratioBkgEff=1.34405 ratioSigEff=1.00273 [-inf, 3.24614] in interval = 1 NeymanConstruction: Prog: 12/100 total MC = 80 this test stat = 1.93333 s=13.8 ratioBkgEff=1.33528 ratioSigEff=1.00294 [-inf, 2.4475] in interval = 1 NeymanConstruction: Prog: 13/100 total MC = 80 this test stat = 1.83331 s=15 ratioBkgEff=1.32684 ratioSigEff=1.0023 [-inf, 2.70017] in interval = 1 NeymanConstruction: Prog: 14/100 total MC = 80 this test stat = 1.73555 s=16.2 ratioBkgEff=1.31792 ratioSigEff=1.00257 [-inf, 2.77026] in interval = 1 NeymanConstruction: Prog: 15/100 total MC = 80 this test stat = 1.64101 s=17.4 ratioBkgEff=1.30904 ratioSigEff=1.00283 [-inf, 2.84104] in interval = 1 NeymanConstruction: Prog: 16/100 total MC = 80 this test stat = 1.54913 s=18.6 ratioBkgEff=1.2991 ratioSigEff=1.00364 [-inf, 2.44366] in interval = 1 NeymanConstruction: Prog: 17/100 total MC = 80 this test stat = 1.45967 s=19.8 ratioBkgEff=1.29046 ratioSigEff=1.0038 [-inf, 2.51046] in interval = 1 NeymanConstruction: Prog: 18/100 total MC = 80 this test stat = 1.3729 s=21 ratioBkgEff=1.28178 ratioSigEff=1.00395 [-inf, 2.3088] in interval = 1 NeymanConstruction: Prog: 19/100 total MC = 80 this test stat = 1.28897 s=22.2 ratioBkgEff=1.27309 ratioSigEff=1.00408 [-inf, 2.28612] in interval = 1 NeymanConstruction: Prog: 20/100 total MC = 80 this test stat = 1.20763 s=23.4 ratioBkgEff=1.26436 ratioSigEff=1.0042 [-inf, 2.0941] in interval = 1 NeymanConstruction: Prog: 21/100 total MC = 80 this test stat = 1.12898 s=24.6 ratioBkgEff=1.25562 ratioSigEff=1.00431 [-inf, 2.07266] in interval = 1 NeymanConstruction: Prog: 22/100 total MC = 80 this test stat = 1.05305 s=25.8 ratioBkgEff=1.24687 ratioSigEff=1.0044 [-inf, 2.4819] in interval = 1 NeymanConstruction: Prog: 23/100 total MC = 80 this test stat = 0.979804 s=27 ratioBkgEff=1.2381 ratioSigEff=1.00448 [-inf, 1.94868] in interval = 1 NeymanConstruction: Prog: 24/100 total MC = 80 this test stat = 0.909232 s=28.2 ratioBkgEff=1.22932 ratioSigEff=1.00454 [-inf, 2.09127] in interval = 1 NeymanConstruction: Prog: 25/100 total MC = 80 this test stat = 0.841264 s=29.4 ratioBkgEff=1.22097 ratioSigEff=1.00408 [-inf, 1.9082] in interval = 1 NeymanConstruction: Prog: 26/100 total MC = 80 this test stat = 0.776027 s=30.6 ratioBkgEff=1.21214 ratioSigEff=1.00408 [-inf, 2.04865] in interval = 1 NeymanConstruction: Prog: 27/100 total MC = 80 this test stat = 0.713453 s=31.8 ratioBkgEff=1.2033 ratioSigEff=1.00406 [-inf, 1.4258] in interval = 1 NeymanConstruction: Prog: 28/100 total MC = 80 this test stat = 0.653398 s=33 ratioBkgEff=1.19446 ratioSigEff=1.00403 [-inf, 1.92614] in interval = 1 NeymanConstruction: Prog: 29/100 total MC = 80 this test stat = 0.596234 s=34.2 ratioBkgEff=1.18562 ratioSigEff=1.00399 [-inf, 1.53064] in interval = 1 NeymanConstruction: Prog: 30/100 total MC = 80 this test stat = 0.541678 s=35.4 ratioBkgEff=1.17679 ratioSigEff=1.00393 [-inf, 1.65725] in interval = 1 NeymanConstruction: Prog: 31/100 total MC = 80 this test stat = 0.489749 s=36.6 ratioBkgEff=1.16795 ratioSigEff=1.00386 [-inf, 1.42567] in interval = 1 NeymanConstruction: Prog: 32/100 total MC = 80 this test stat = 0.440357 s=37.8 ratioBkgEff=1.15913 ratioSigEff=1.00378 [-inf, 1.40855] in interval = 1 NeymanConstruction: Prog: 33/100 total MC = 80 this test stat = 0.393811 s=39 ratioBkgEff=1.15032 ratioSigEff=1.00369 [-inf, 1.2593] in interval = 1 NeymanConstruction: Prog: 34/100 total MC = 80 this test stat = 0.349798 s=40.2 ratioBkgEff=1.14152 ratioSigEff=1.00358 [-inf, 1.30825] in interval = 1 NeymanConstruction: Prog: 35/100 total MC = 80 this test stat = 0.30842 s=41.4 ratioBkgEff=1.13274 ratioSigEff=1.00347 [-inf, 0.874179] in interval = 1 NeymanConstruction: Prog: 36/100 total MC = 80 this test stat = 0.269673 s=42.6 ratioBkgEff=1.12397 ratioSigEff=1.00334 [-inf, 1.14949] in interval = 1 NeymanConstruction: Prog: 37/100 total MC = 80 this test stat = 0.233551 s=43.8 ratioBkgEff=1.11574 ratioSigEff=1.00319 [-inf, 1.07403] in interval = 1 NeymanConstruction: Prog: 38/100 total MC = 80 this test stat = 0.200042 s=45 ratioBkgEff=1.10705 ratioSigEff=1.00303 [-inf, 1.30848] in interval = 1 NeymanConstruction: Prog: 39/100 total MC = 80 this test stat = 0.169154 s=46.2 ratioBkgEff=1.09838 ratioSigEff=1.00285 [-inf, 0.987229] in interval = 1 NeymanConstruction: Prog: 40/100 total MC = 80 this test stat = 0.140876 s=47.4 ratioBkgEff=1.08972 ratioSigEff=1.00267 [-inf, 0.917557] in interval = 1 NeymanConstruction: Prog: 41/100 total MC = 80 this test stat = 0.115202 s=48.6 ratioBkgEff=1.08109 ratioSigEff=1.00247 [-inf, 0.850481] in interval = 1 NeymanConstruction: Prog: 42/100 total MC = 80 this test stat = 0.0921253 s=49.8 ratioBkgEff=1.07247 ratioSigEff=1.00226 [-inf, 0.891152] in interval = 1 NeymanConstruction: Prog: 43/100 total MC = 80 this test stat = 0.0716406 s=51 ratioBkgEff=1.06387 ratioSigEff=1.00204 [-inf, 1.04631] in interval = 1 NeymanConstruction: Prog: 44/100 total MC = 80 this test stat = 0.0537413 s=52.2 ratioBkgEff=1.05529 ratioSigEff=1.00181 [-inf, 0.574131] in interval = 1 NeymanConstruction: Prog: 45/100 total MC = 80 this test stat = 0.0384208 s=53.4 ratioBkgEff=1.04672 ratioSigEff=1.00156 [-inf, 0.749926] in interval = 1 NeymanConstruction: Prog: 46/100 total MC = 80 this test stat = 0.0256726 s=54.6 ratioBkgEff=1.03818 ratioSigEff=1.00131 [-inf, 0.839739] in interval = 1 NeymanConstruction: Prog: 47/100 total MC = 80 this test stat = 0.0154896 s=55.8 ratioBkgEff=1.02966 ratioSigEff=1.00104 [-inf, 0.715031] in interval = 1 NeymanConstruction: Prog: 48/100 total MC = 80 this test stat = 0.00786498 s=57 ratioBkgEff=1.02116 ratioSigEff=1.00075 [-inf, 0.574092] in interval = 1 NeymanConstruction: Prog: 49/100 total MC = 80 this test stat = 0.00303984 s=58.2 ratioBkgEff=1.01164 ratioSigEff=1.00152 [-inf, 0.731975] in interval = 1 NeymanConstruction: Prog: 50/100 total MC = 80 this test stat = 0.000296144 s=59.4 ratioBkgEff=1.0038 ratioSigEff=1.00057 [-inf, 0.692058] in interval = 1 NeymanConstruction: Prog: 51/100 total MC = 80 this test stat = 0.000277433 s=60.6 ratioBkgEff=0.995997 ratioSigEff=0.999598 [-inf, 0.645915] in interval = 1 NeymanConstruction: Prog: 52/100 total MC = 80 this test stat = 0.00297545 s=61.8 ratioBkgEff=0.988241 ratioSigEff=0.998613 [-inf, 0.513395] in interval = 1 NeymanConstruction: Prog: 53/100 total MC = 80 this test stat = 0.00784733 s=63 ratioBkgEff=0.978981 ratioSigEff=0.999173 [-inf, 0.477501] in interval = 1 NeymanConstruction: Prog: 54/100 total MC = 80 this test stat = 0.0154172 s=64.2 ratioBkgEff=0.970614 ratioSigEff=0.998822 [-inf, 0.721422] in interval = 1 NeymanConstruction: Prog: 55/100 total MC = 80 this test stat = 0.0254769 s=65.4 ratioBkgEff=0.962271 ratioSigEff=0.998461 [-inf, 0.678264] in interval = 1 NeymanConstruction: Prog: 56/100 total MC = 80 this test stat = 0.0380399 s=66.6 ratioBkgEff=0.953952 ratioSigEff=0.998087 [-inf, 0.636603] in interval = 1 NeymanConstruction: Prog: 57/100 total MC = 80 this test stat = 0.0530874 s=67.8 ratioBkgEff=0.945657 ratioSigEff=0.997703 [-inf, 0.596425] in interval = 1 NeymanConstruction: Prog: 58/100 total MC = 80 this test stat = 0.0706058 s=69 ratioBkgEff=0.937388 ratioSigEff=0.997308 [-inf, 0.70275] in interval = 1 NeymanConstruction: Prog: 59/100 total MC = 80 this test stat = 0.0905894 s=70.2 ratioBkgEff=0.929144 ratioSigEff=0.996902 [-inf, 0.93255] in interval = 1 NeymanConstruction: Prog: 60/100 total MC = 80 this test stat = 0.11303 s=71.4 ratioBkgEff=0.920926 ratioSigEff=0.996485 [-inf, 1.06482] in interval = 1 NeymanConstruction: Prog: 61/100 total MC = 80 this test stat = 0.137918 s=72.6 ratioBkgEff=0.912734 ratioSigEff=0.996057 [-inf, 0.780337] in interval = 1 NeymanConstruction: Prog: 62/100 total MC = 80 this test stat = 0.165246 s=73.8 ratioBkgEff=0.904568 ratioSigEff=0.995619 [-inf, 0.684448] in interval = 1 NeymanConstruction: Prog: 63/100 total MC = 80 this test stat = 0.195001 s=75 ratioBkgEff=0.89643 ratioSigEff=0.995169 [-inf, 1.03027] in interval = 1 NeymanConstruction: Prog: 64/100 total MC = 80 this test stat = 0.227176 s=76.2 ratioBkgEff=0.888318 ratioSigEff=0.99471 [-inf, 0.978263] in interval = 1 NeymanConstruction: Prog: 65/100 total MC = 80 this test stat = 0.261776 s=77.4 ratioBkgEff=0.880234 ratioSigEff=0.99424 [-inf, 1.11256] in interval = 1 NeymanConstruction: Prog: 66/100 total MC = 80 this test stat = 0.298772 s=78.6 ratioBkgEff=0.872179 ratioSigEff=0.993759 [-inf, 1.18776] in interval = 1 NeymanConstruction: Prog: 67/100 total MC = 80 this test stat = 0.338144 s=79.8 ratioBkgEff=0.864151 ratioSigEff=0.993269 [-inf, 1.48101] in interval = 1 NeymanConstruction: Prog: 68/100 total MC = 80 this test stat = 0.379914 s=81 ratioBkgEff=0.856152 ratioSigEff=0.992768 [-inf, 1.56755] in interval = 1 NeymanConstruction: Prog: 69/100 total MC = 80 this test stat = 0.424071 s=82.2 ratioBkgEff=0.848182 ratioSigEff=0.992257 [-inf, 1.21708] in interval = 1 NeymanConstruction: Prog: 70/100 total MC = 80 this test stat = 0.470561 s=83.4 ratioBkgEff=0.840241 ratioSigEff=0.991736 [-inf, 1.43809] in interval = 1 NeymanConstruction: Prog: 71/100 total MC = 80 this test stat = 0.519388 s=84.6 ratioBkgEff=0.83233 ratioSigEff=0.991205 [-inf, 1.84094] in interval = 1 NeymanConstruction: Prog: 72/100 total MC = 80 this test stat = 0.570571 s=85.8 ratioBkgEff=0.824449 ratioSigEff=0.990664 [-inf, 1.31514] in interval = 1 NeymanConstruction: Prog: 73/100 total MC = 80 this test stat = 0.624119 s=87 ratioBkgEff=0.816598 ratioSigEff=0.990113 [-inf, 1.18918] in interval = 1 NeymanConstruction: Prog: 74/100 total MC = 80 this test stat = 0.679954 s=88.2 ratioBkgEff=0.808777 ratioSigEff=0.989553 [-inf, 1.33503] in interval = 1 NeymanConstruction: Prog: 75/100 total MC = 80 this test stat = 0.738105 s=89.4 ratioBkgEff=0.800988 ratioSigEff=0.988983 [-inf, 1.56456] in interval = 1 NeymanConstruction: Prog: 76/100 total MC = 80 this test stat = 0.798548 s=90.6 ratioBkgEff=0.79323 ratioSigEff=0.988403 [-inf, 1.49965] in interval = 1 NeymanConstruction: Prog: 77/100 total MC = 80 this test stat = 0.861306 s=91.8 ratioBkgEff=0.785503 ratioSigEff=0.987814 [-inf, 1.66275] in interval = 1 NeymanConstruction: Prog: 78/100 total MC = 80 this test stat = 0.926319 s=93 ratioBkgEff=0.777809 ratioSigEff=0.987216 [-inf, 1.44643] in interval = 1 NeymanConstruction: Prog: 79/100 total MC = 80 this test stat = 0.993583 s=94.2 ratioBkgEff=0.770146 ratioSigEff=0.986608 [-inf, 2.19051] in interval = 1 NeymanConstruction: Prog: 80/100 total MC = 80 this test stat = 1.06308 s=95.4 ratioBkgEff=0.762516 ratioSigEff=0.985991 [-inf, 1.54054] in interval = 1 NeymanConstruction: Prog: 81/100 total MC = 80 this test stat = 1.13458 s=96.6 ratioBkgEff=0.754919 ratioSigEff=0.985365 [-inf, 1.86615] in interval = 1 NeymanConstruction: Prog: 82/100 total MC = 80 this test stat = 1.20859 s=97.8 ratioBkgEff=0.747355 ratioSigEff=0.98473 [-inf, 2.40993] in interval = 1 NeymanConstruction: Prog: 83/100 total MC = 80 this test stat = 1.28501 s=99 ratioBkgEff=0.739824 ratioSigEff=0.984086 [-inf, 2.32724] in interval = 1 NeymanConstruction: Prog: 84/100 total MC = 80 this test stat = 1.36342 s=100.2 ratioBkgEff=0.732327 ratioSigEff=0.983433 [-inf, 2.43281] in interval = 1 NeymanConstruction: Prog: 85/100 total MC = 80 this test stat = 1.44405 s=101.4 ratioBkgEff=0.724864 ratioSigEff=0.982771 [-inf, 2.44403] in interval = 1 NeymanConstruction: Prog: 86/100 total MC = 80 this test stat = 1.52687 s=102.6 ratioBkgEff=0.717436 ratioSigEff=0.982101 [-inf, 2.09032] in interval = 1 NeymanConstruction: Prog: 87/100 total MC = 80 this test stat = 1.61182 s=103.8 ratioBkgEff=0.710041 ratioSigEff=0.981422 [-inf, 2.76135] in interval = 1 NeymanConstruction: Prog: 88/100 total MC = 80 this test stat = 1.69896 s=105 ratioBkgEff=0.702682 ratioSigEff=0.980734 [-inf, 2.38331] in interval = 1 NeymanConstruction: Prog: 89/100 total MC = 80 this test stat = 1.78823 s=106.2 ratioBkgEff=0.695358 ratioSigEff=0.980038 [-inf, 2.39444] in interval = 1 NeymanConstruction: Prog: 90/100 total MC = 80 this test stat = 1.87929 s=107.4 ratioBkgEff=0.688069 ratioSigEff=0.979333 [-inf, 2.89798] in interval = 1 NeymanConstruction: Prog: 91/100 total MC = 80 this test stat = 1.97313 s=108.6 ratioBkgEff=0.680815 ratioSigEff=0.97862 [-inf, 3.1204] in interval = 1 NeymanConstruction: Prog: 92/100 total MC = 80 this test stat = 2.06825 s=109.8 ratioBkgEff=0.67064 ratioSigEff=0.977383 [-inf, 3.13169] in interval = 1 NeymanConstruction: Prog: 93/100 total MC = 80 this test stat = 2.16592 s=111 ratioBkgEff=0.663285 ratioSigEff=0.976625 [-inf, 2.93195] in interval = 1 NeymanConstruction: Prog: 94/100 total MC = 80 this test stat = 2.26567 s=112.2 ratioBkgEff=0.655961 ratioSigEff=0.975857 [-inf, 2.54424] in interval = 1 NeymanConstruction: Prog: 95/100 total MC = 80 this test stat = 2.36747 s=113.4 ratioBkgEff=0.648488 ratioSigEff=0.975049 [-inf, 3.05911] in interval = 1 NeymanConstruction: Prog: 96/100 total MC = 80 this test stat = 2.47112 s=114.6 ratioBkgEff=0.641205 ratioSigEff=0.97426 [-inf, 2.96555] in interval = 1 NeymanConstruction: Prog: 97/100 total MC = 80 this test stat = 2.57725 s=115.8 ratioBkgEff=0.633952 ratioSigEff=0.973462 [-inf, 3.5205] in interval = 1 NeymanConstruction: Prog: 98/100 total MC = 240 this test stat = 2.68512 s=117 ratioBkgEff=0.626728 ratioSigEff=0.972655 [-inf, 3.19934] in interval = 1 NeymanConstruction: Prog: 99/100 total MC = 80 this test stat = 2.79507 s=118.2 ratioBkgEff=0.619534 ratioSigEff=0.971839 [-inf, 3.77433] in interval = 1 NeymanConstruction: Prog: 100/100 total MC = 80 this test stat = 2.90706 s=119.4 ratioBkgEff=0.61237 ratioSigEff=0.971014 [-inf, 3.66884] in interval = 1 [#1] INFO:Eval -- 100 points in interval [#1] INFO:Minimization -- Including the following constraint terms in minimization: (sigConstraint,bkgConstraint) [#1] INFO:Minimization -- The global observables are not defined , normalize constraints with respect to the parameters (b,ratioBkgEff,ratioSigEff,s) [#1] INFO:Fitting -- RooAbsPdf::fitTo(modelWithConstraints) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_modelWithConstraints_exampleData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
Third, use a Calculator based on Markov Chain monte carlo Before configuring the calculator, let's make a ProposalFunction that will achieve a high acceptance rate
ph = ROOT.RooStats.ProposalHelper()
ph.SetVariables(fit.floatParsFinal())
ph.SetCovMatrix(fit.covarianceMatrix())
ph.SetUpdateProposalParameters(True)
ph.SetCacheSize(100)
pdfProp = ph.GetProposalFunction()
mc = ROOT.RooStats.MCMCCalculator(data, modelConfig)
mc.SetNumIters(20000) # steps to propose in the chain
mc.SetTestSize(0.05) # 95% CL
mc.SetNumBurnInSteps(40) # ignore first N steps in chain as "burn in"
mc.SetProposalFunction(pdfProp)
mc.SetLeftSideTailFraction(0.5) # find a "central" interval
mcInt = mc.GetInterval()
[#1] INFO:Minimization -- Including the following constraint terms in minimization: (sigConstraint,bkgConstraint) [#1] INFO:Minimization -- The following global observables have been defined and their values are taken from the model: (gSigEff,gSigBkg) [#1] INFO:Fitting -- RooAbsPdf::fitTo(modelWithConstraints) fixing normalization set for coefficient determination to observables in data Metropolis-Hastings progress: .................................................................................................... [#1] INFO:Eval -- Proposal acceptance rate: 48.15% [#1] INFO:Eval -- Number of steps in chain: 9630
Get Lower and Upper limits from Profile Calculator
print("Profile lower limit on s = ", lrinterval.LowerLimit(s))
print("Profile upper limit on s = ", lrinterval.UpperLimit(s))
Profile lower limit on s = 13.949738060154317 Profile upper limit on s = 107.9503617757897
Get Lower and Upper limits from FeldmanCousins with profile construction
if fcint:
fcul = fcint.UpperLimit(s)
fcll = fcint.LowerLimit(s)
print("FC lower limit on s = ", fcll)
print("FC upper limit on s = ", fcul)
fcllLine = ROOT.TLine(fcll, 0, fcll, 1)
fculLine = ROOT.TLine(fcul, 0, fcul, 1)
fcllLine.SetLineColor(ROOT.kRed)
fculLine.SetLineColor(ROOT.kRed)
fcllLine.Draw("same")
fculLine.Draw("same")
dataCanvas.Update()
FC lower limit on s = 0.6 FC upper limit on s = 119.39999999999999
Plot MCMC interval and print some statistics
mcPlot = ROOT.RooStats.MCMCIntervalPlot(mcInt)
mcPlot.SetLineColor(ROOT.kMagenta)
mcPlot.SetLineWidth(2)
mcPlot.Draw("same")
mcul = mcInt.UpperLimit(s)
mcll = mcInt.LowerLimit(s)
print("MCMC lower limit on s = ", mcll)
print("MCMC upper limit on s = ", mcul)
print("MCMC Actual confidence level: ", mcInt.GetActualConfidenceLevel())
MCMC lower limit on s = 19.16837376625842 MCMC upper limit on s = 103.15897826304688 MCMC Actual confidence level: 0.9499146843320285
3-d plot of the parameter points
dataCanvas.cd(2)
<cppyy.gbl.TPad object at 0xacafbe0>
also plot the points in the markov chain
chainData = mcInt.GetChainAsDataSet()
print("plotting the chain data - nentries = ", chainData.numEntries())
chain = ROOT.RooStats.GetAsTTree("chainTreeData", "chainTreeData", chainData)
chain.SetMarkerStyle(6)
chain.SetMarkerColor(ROOT.kRed)
chain.Draw("s:ratioSigEff:ratioBkgEff", "nll_MarkovChain_local_", "box") # 3-d box proportional to posterior
plotting the chain data - nentries = 9630
9630
the points used in the profile construction
parScanData = fc.GetPointsToScan()
print("plotting the scanned points used in the frequentist construction - npoints = ", parScanData.numEntries())
gr = ROOT.TGraph2D(parScanData.numEntries())
for ievt in range(parScanData.numEntries()):
evt = parScanData.get(ievt)
x = evt.getRealValue("ratioBkgEff")
y = evt.getRealValue("ratioSigEff")
z = evt.getRealValue("s")
gr.SetPoint(ievt, x, y, z)
gr.SetMarkerStyle(24)
gr.Draw("P SAME")
plotting the scanned points used in the frequentist construction - npoints = 100
print timing info
t.Stop()
t.Print()
dataCanvas.SaveAs("rs101_limitexample.png")
Real time 0:00:20, CP time 17.630
Warning in <TAxis::TAxis::SetRangeUser>: ufirst < fXmin, fXmin is used Warning in <TAxis::TAxis::SetRangeUser>: ulast > fXmax, fXmax is used Warning in <TAxis::TAxis::SetRangeUser>: ufirst < fXmin, fXmin is used Warning in <TAxis::TAxis::SetRangeUser>: ulast > fXmax, fXmax is used Warning in <TAxis::TAxis::SetRangeUser>: ufirst < fXmin, fXmin is used Warning in <TAxis::TAxis::SetRangeUser>: ulast > fXmax, fXmax is used Warning in <TAxis::TAxis::SetRangeUser>: ufirst < fXmin, fXmin is used Warning in <TAxis::TAxis::SetRangeUser>: ulast > fXmax, fXmax is used Info in <TCanvas::Print>: png file rs101_limitexample.png has been created
Draw all canvases
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()
Warning in <TAxis::TAxis::SetRangeUser>: ufirst < fXmin, fXmin is used Warning in <TAxis::TAxis::SetRangeUser>: ulast > fXmax, fXmax is used Warning in <TAxis::TAxis::SetRangeUser>: ufirst < fXmin, fXmin is used Warning in <TAxis::TAxis::SetRangeUser>: ulast > fXmax, fXmax is used