Example of the BernsteinCorrection utility in RooStats.
The idea is that one has a distribution coming either from data or Monte Carlo (called "reality" in the macro) and a nominal model that is not sufficiently flexible to take into account the real distribution. One wants to take into account the systematic associated with this imperfect modeling by augmenting the nominal model with some correction term (in this case a polynomial). The BernsteinCorrection utility will import into your workspace a corrected model given by nominal(x) * poly_N(x), where poly_N is an n-th order polynomial in the Bernstein basis. The degree N of the polynomial is chosen by specifying the tolerance one has in adding an extra term to the polynomial. The Bernstein basis is nice because it only has positive-definite terms and works well with PDFs. Finally, the macro makes a plot of:
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 sys
import ROOT
set range of observable
lowRange = -1
highRange = 5
make a RooRealVar for the observable
x = ROOT.RooRealVar("x", "x", lowRange, highRange)
true model
narrow = ROOT.RooGaussian("narrow", "", x, ROOT.RooFit.RooConst(0.0), ROOT.RooFit.RooConst(0.8))
wide = ROOT.RooGaussian("wide", "", x, ROOT.RooFit.RooConst(0.0), ROOT.RooFit.RooConst(2.0))
reality = ROOT.RooAddPdf("reality", "", [narrow, wide], ROOT.RooFit.RooConst(0.8))
data = reality.generate(x, 1000)
nominal model
sigma = ROOT.RooRealVar("sigma", "", 1.0, 0, 10)
nominal = ROOT.RooGaussian("nominal", "", x, ROOT.RooFit.RooConst(0.0), sigma)
wks = ROOT.RooWorkspace("myWorksspace")
wks.Import(data, Rename="data")
wks.Import(nominal)
if ROOT.TClass.GetClass("ROOT::Minuit2::Minuit2Minimizer"):
# use Minuit2 if ROOT was built with support for it:
ROOT.Math.MinimizerOptions.SetDefaultMinimizer("Minuit2")
[#0] WARNING:InputArguments -- The parameter 'sigma' with range [0, 10] of the RooGaussian 'nominal' exceeds the safe range of (0, inf). Advise to limit its range. [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing dataset realityData [#1] INFO:ObjectHandling -- RooWorkSpace::import(myWorksspace) changing name of dataset from realityData to data [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::x [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooGaussian::nominal [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooConstVar::0 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::sigma
The tolerance sets the probability to add an unnecessary term. lower tolerance will add fewer terms, while higher tolerance will add more terms and provide a more flexible function.
tolerance = 0.05
bernsteinCorrection = ROOT.RooStats.BernsteinCorrection(tolerance)
degree = bernsteinCorrection.ImportCorrectedPdf(wks, "nominal", "x", "data")
if degree < 0:
ROOT.Error("rs_bernsteinCorrection", "Bernstein correction failed !")
sys.exit()
print("Correction based on Bernstein Poly of degree ", degree)
frame = x.frame()
data.plotOn(frame)
Correction based on Bernstein Poly of degree 6
<cppyy.gbl.RooPlot object at 0xb972150>
BernsteinCorrection::ImportCorrectedPdf - Doing initial Fit with nominal model [#1] INFO:Fitting -- RooAbsPdf::fitTo(nominal_over_nominal_Int[x]) 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_nominal_over_nominal_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_data) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooEffProd::corrected [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooBernstein::poly [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_0 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_1 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_2 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_3 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_4 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_5 [#1] INFO:ObjectHandling -- RooWorkspace::import(myWorksspace) importing RooRealVar::c_6 ------ Begin Bernstein Correction Log -------- degree = 1 -log L(0) = 1216.78 -log L(1) = 1208.89 q = 15.7692 P(chi^2_1 > q) = 7.1557e-05 degree = 2 -log L(1) = 1208.89 -log L(2) = 1203.21 q = 11.3692 P(chi^2_1 > q) = 0.000746732 degree = 3 -log L(2) = 1203.21 -log L(3) = 1198.85 q = 8.72171 P(chi^2_1 > q) = 0.00314444 degree = 4 -log L(3) = 1198.85 -log L(4) = 1190.16 q = 17.3777 P(chi^2_1 > q) = 3.06393e-05 degree = 5 -log L(4) = 1190.16 -log L(5) = 1183.56 q = 13.1965 P(chi^2_1 > q) = 0.00028048 degree = 6 -log L(5) = 1183.56 -log L(6) = 1182.57 q = 1.98352 P(chi^2_1 > q) = 0.15902 ------ End Bernstein Correction Log --------
plot the best fit nominal model in blue
nominal.fitTo(data, PrintLevel=0)
nominal.plotOn(frame)
<cppyy.gbl.RooPlot object at 0xb972150>
[#1] INFO:Fitting -- RooAbsPdf::fitTo(nominal_over_nominal_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_nominal_over_nominal_Int[x]_realityData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization Minuit2Minimizer: Minimize with max-calls 500 convergence for edm < 1 strategy 1 Minuit2Minimizer : Valid minimum - status = 0 FVAL = 1216.77793416266286 Edm = 4.16187344562173991e-07 Nfcn = 19 sigma = 1.18138 +/- 0.0315451 (limited) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
plot the best fit corrected model in red
corrected = wks["corrected"]
if not corrected:
sys.exit()
fit corrected model
corrected.fitTo(data, PrintLevel=0)
corrected.plotOn(frame, LineColor="r")
<cppyy.gbl.RooPlot object at 0xb972150>
[#1] INFO:Fitting -- RooAbsPdf::fitTo(corrected_over_corrected_Int[x]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_corrected_over_corrected_Int[x]_realityData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization Minuit2Minimizer: Minimize with max-calls 3500 convergence for edm < 1 strategy 1 [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x) Minuit2Minimizer : Valid minimum - status = 0 FVAL = 1182.56856239518879 Edm = 0.000132754139401894118 Nfcn = 187 c_1 = 3.18337 +/- 0.786698 (limited) c_2 = 4.97861e-06 +/- 3.10992 (limited) c_3 = 1.17018e-06 +/- 1.51847 (limited) c_4 = 0.966521 +/- 2.16515 (limited) c_5 = 0.216513 +/- 67.8056 (limited) c_6 = 10.4402 +/- 20.6985 (limited) sigma = 1.2667 +/- 0.210846 (limited) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization [#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
plot the correction term (* norm constant) in dashed green should make norm constant just be 1, not depend on binning of data
poly = wks["poly"]
if poly:
poly.plotOn(frame, LineColor="g", LineStyle="--")
this is a switch to check the sampling distribution of -2 log LR for two comparisons: the first is for n-1 vs. n degree polynomial corrections the second is for n vs. n+1 degree polynomial corrections Here we choose n to be the one chosen by the tolerance criterion above, eg. n = "degree" in the code. Setting this to true is takes about 10 min.
checkSamplingDist = True
numToyMC = 20 # increase this value for sensible results
c1 = ROOT.TCanvas()
if checkSamplingDist:
c1.Divide(1, 2)
c1.cd(1)
frame.Draw()
ROOT.gPad.Update()
if checkSamplingDist:
# check sampling dist
ROOT.Math.MinimizerOptions.SetDefaultPrintLevel(-1)
samplingDist = ROOT.TH1F("samplingDist", "", 20, 0, 10)
samplingDistExtra = ROOT.TH1F("samplingDistExtra", "", 20, 0, 10)
bernsteinCorrection.CreateQSamplingDist(
wks, "nominal", "x", "data", samplingDist, samplingDistExtra, degree, numToyMC
)
c1.cd(2)
samplingDistExtra.SetLineColor(ROOT.kRed)
samplingDistExtra.Draw()
samplingDist.Draw("same")
c1.SaveAs("rs_bernsteinCorrection.png")
made pdfs, make toy generator on toy 0 on toy 1 on toy 2 on toy 3 on toy 4 on toy 5 on toy 6 on toy 7 on toy 8 on toy 9 on toy 10 on toy 11 on toy 12 on toy 13 on toy 14 on toy 15 on toy 16 on toy 17 on toy 18 on toy 19
Info in <TCanvas::Print>: png file rs_bernsteinCorrection.png has been created
Draw all canvases
%jsroot on
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()