Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
rs_bernsteinCorrection.C File Reference
Tutorials » RooStats Tutorials

Detailed Description

View in nbviewer Open in SWAN
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:

  • the data (drawn from 'reality'),
  • the best fit of the nominal model (blue)
  • and the best fit corrected model.
[#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
BernsteinCorrection::ImportCorrectedPdf - Doing initial Fit with nominal model
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#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 --------
Correction based on Bernstein Poly of degree 6
[#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
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
Minuit2Minimizer: Minimize with max-calls 3500 convergence for edm < 1 strategy 1
Minuit2Minimizer : Valid minimum - status = 0
FVAL = 1182.56856195034334
Edm = 0.000132254574581617092
Nfcn = 185
c_1 = 3.18345 +/- 0.786682 (limited)
c_2 = 4.73727e-06 +/- 3.10512 (limited)
c_3 = 1.11806e-06 +/- 1.51813 (limited)
c_4 = 0.966535 +/- 2.16507 (limited)
c_5 = 0.216439 +/- 67.8113 (limited)
c_6 = 10.4409 +/- 20.6997 (limited)
sigma = 1.26669 +/- 0.210849 (limited)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(corrected_Int[x]) using numeric integrator RooRombergIntegrator to calculate Int(x)
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
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooBernstein.h"
#include "TCanvas.h"
#include "RooAbsPdf.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include <string>
#include <vector>
#include <stdio.h>
#include <sstream>
#include <iostream>
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooGaussian.h"
#include "RooProfileLL.h"
#include "RooWorkspace.h"
#include "RooStats/BernsteinCorrection.h"
// use this order for safety on library loading
using namespace RooFit;
using namespace RooStats;
//____________________________________
void rs_bernsteinCorrection()
{
// set range of observable
Double_t lowRange = -1, highRange = 5;
// make a RooRealVar for the observable
RooRealVar x("x", "x", lowRange, highRange);
// true model
RooGaussian narrow("narrow", "", x, RooConst(0.), RooConst(.8));
RooGaussian wide("wide", "", x, RooConst(0.), RooConst(2.));
RooAddPdf reality("reality", "", RooArgList(narrow, wide), RooConst(0.8));
std::unique_ptr<RooDataSet> data{reality.generate(x, 1000)};
// nominal model
RooRealVar sigma("sigma", "", 1., 0, 10);
RooGaussian nominal("nominal", "", x, RooConst(0.), sigma);
RooWorkspace *wks = new RooWorkspace("myWorksspace");
wks->import(*data, Rename("data"));
wks->import(nominal);
if (TClass::GetClass("ROOT::Minuit2::Minuit2Minimizer")) {
// use Minuit2 if ROOT was built with support for it:
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
}
// 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.
Double_t tolerance = 0.05;
BernsteinCorrection bernsteinCorrection(tolerance);
Int_t degree = bernsteinCorrection.ImportCorrectedPdf(wks, "nominal", "x", "data");
if (degree < 0) {
Error("rs_bernsteinCorrection", "Bernstein correction failed ! ");
return;
}
cout << " Correction based on Bernstein Poly of degree " << degree << endl;
RooPlot *frame = x.frame();
data->plotOn(frame);
// plot the best fit nominal model in blue
nominal.fitTo(*data, PrintLevel(0));
nominal.plotOn(frame);
// plot the best fit corrected model in red
RooAbsPdf *corrected = wks->pdf("corrected");
if (!corrected)
return;
// fit corrected model
corrected->fitTo(*data, PrintLevel(0));
corrected->plotOn(frame, LineColor(kRed));
// plot the correction term (* norm constant) in dashed green
// should make norm constant just be 1, not depend on binning of data
RooAbsPdf *poly = wks->pdf("poly");
if (poly)
poly->plotOn(frame, LineColor(kGreen), LineStyle(kDashed));
// 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.
bool checkSamplingDist = true;
int numToyMC = 20; // increase this value for sensible results
TCanvas *c1 = new TCanvas();
if (checkSamplingDist) {
c1->Divide(1, 2);
c1->cd(1);
}
frame->Draw();
gPad->Update();
if (checkSamplingDist) {
// check sampling dist
ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(-1);
TH1F *samplingDist = new TH1F("samplingDist", "", 20, 0, 10);
TH1F *samplingDistExtra = new TH1F("samplingDistExtra", "", 20, 0, 10);
bernsteinCorrection.CreateQSamplingDist(wks, "nominal", "x", "data", samplingDist, samplingDistExtra, degree,
numToyMC);
c1->cd(2);
samplingDistExtra->SetLineColor(kRed);
samplingDistExtra->Draw();
samplingDist->Draw("same");
}
}
Int_t
int Int_t
Definition RtypesCore.h:45
Double_t
double Double_t
Definition RtypesCore.h:59
kRed
@ kRed
Definition Rtypes.h:66
kGreen
@ kGreen
Definition Rtypes.h:66
kDashed
@ kDashed
Definition TAttLine.h:48
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Definition TGWin32VirtualXProxy.cxx:104
gPad
#define gPad
Definition TVirtualPad.h:305
ROOT::Math::MinimizerOptions::SetDefaultMinimizer
static void SetDefaultMinimizer(const char *type, const char *algo=nullptr)
Set the default Minimizer type and corresponding algorithms.
Definition MinimizerOptions.cxx:43
ROOT::Math::MinimizerOptions::SetDefaultPrintLevel
static void SetDefaultPrintLevel(int level)
Set the default Print Level.
Definition MinimizerOptions.cxx:75
RooAbsPdf
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooAbsPdf::plotOn
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:123
RooAbsPdf::fitTo
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:156
RooAddPdf
Efficient implementation of a sum of PDFs of the form.
Definition RooAddPdf.h:33
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooGaussian
Plain Gaussian p.d.f.
Definition RooGaussian.h:24
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
RooPlot::frame
static RooPlot * frame(const RooAbsRealLValue &var, double xmin, double xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:239
RooPlot::Draw
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:652
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:37
RooStats::BernsteinCorrection
BernsteinCorrection is a utility in RooStats to augment a nominal PDF with a polynomial correction te...
Definition BernsteinCorrection.h:22
RooWorkspace
Persistable container for RooFit projects.
Definition RooWorkspace.h:43
RooWorkspace::pdf
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Definition RooWorkspace.cxx:1212
RooWorkspace::import
bool import(const RooAbsArg &arg, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={})
Import a RooAbsArg object, e.g.
Definition RooWorkspace.cxx:368
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
TCanvas
The Canvas class.
Definition TCanvas.h:23
TClass::GetClass
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition TClass.cxx:2968
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:577
TH1::Draw
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3067
RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition RooGlobalFunc.cxx:572
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition RooGlobalFunc.cxx:173
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition RooGlobalFunc.cxx:177
sigma
const Double_t sigma
Definition h1analysisProxy.h:11
c1
return c1
Definition legend1.C:41
x
Double_t x[n]
Definition legend1.C:17
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
RooStats
Namespace for the RooStats classes.
Definition Asimov.h:19
rs_bernsteinCorrection
Definition rs_bernsteinCorrection.py:1
Author
Kyle Cranmer

Definition in file rs_bernsteinCorrection.C.