ROOT   Reference Guide
rs_bernsteinCorrection.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// Example of the BernsteinCorrection utility in RooStats.
5///
6/// The idea is that one has a distribution coming either from data or Monte Carlo
7/// (called "reality" in the macro) and a nominal model that is not sufficiently
8/// flexible to take into account the real distribution. One wants to take into
9/// account the systematic associated with this imperfect modeling by augmenting
10/// the nominal model with some correction term (in this case a polynomial).
11/// The BernsteinCorrection utility will import into your workspace a corrected model
12/// given by nominal(x) * poly_N(x), where poly_N is an n-th order polynomial in
13/// the Bernstein basis. The degree N of the polynomial is chosen by specifying the tolerance
14/// one has in adding an extra term to the polynomial.
15/// The Bernstein basis is nice because it only has positive-definite terms
16/// and works well with PDFs.
17/// Finally, the macro makes a plot of:
18/// - the data (drawn from 'reality'),
19/// - the best fit of the nominal model (blue)
20/// - and the best fit corrected model.
21///
22/// \macro_image
23/// \macro_output
24/// \macro_code
25///
26/// \author Kyle Cranmer
27
28#include "RooDataSet.h"
29#include "RooRealVar.h"
30#include "RooConstVar.h"
31#include "RooBernstein.h"
32#include "TCanvas.h"
33#include "RooAbsPdf.h"
34#include "RooFit.h"
35#include "RooFitResult.h"
36#include "RooPlot.h"
37#include <string>
38#include <vector>
39#include <stdio.h>
40#include <sstream>
41#include <iostream>
42
43#include "RooProdPdf.h"
45#include "RooGaussian.h"
46#include "RooNLLVar.h"
47#include "RooMinuit.h"
48#include "RooProfileLL.h"
49#include "RooWorkspace.h"
50
52
54using namespace RooFit;
55using namespace RooStats;
56
57//____________________________________
58void rs_bernsteinCorrection()
59{
60
61 // set range of observable
62 Double_t lowRange = -1, highRange = 5;
63
64 // make a RooRealVar for the observable
65 RooRealVar x("x", "x", lowRange, highRange);
66
67 // true model
68 RooGaussian narrow("narrow", "", x, RooConst(0.), RooConst(.8));
69 RooGaussian wide("wide", "", x, RooConst(0.), RooConst(2.));
70 RooAddPdf reality("reality", "", RooArgList(narrow, wide), RooConst(0.8));
71
72 RooDataSet *data = reality.generate(x, 1000);
73
74 // nominal model
75 RooRealVar sigma("sigma", "", 1., 0, 10);
76 RooGaussian nominal("nominal", "", x, RooConst(0.), sigma);
77
78 RooWorkspace *wks = new RooWorkspace("myWorksspace");
79
80 wks->import(*data, Rename("data"));
81 wks->import(nominal);
82
83 if (TClass::GetClass("ROOT::Minuit2::Minuit2Minimizer")) {
84 // use Minuit2 if ROOT was built with support for it:
86 }
87
88 // The tolerance sets the probability to add an unnecessary term.
89 // lower tolerance will add fewer terms, while higher tolerance
90 // will add more terms and provide a more flexible function.
91 Double_t tolerance = 0.05;
92 BernsteinCorrection bernsteinCorrection(tolerance);
93 Int_t degree = bernsteinCorrection.ImportCorrectedPdf(wks, "nominal", "x", "data");
94
95 if (degree < 0) {
96 Error("rs_bernsteinCorrection", "Bernstein correction failed ! ");
97 return;
98 }
99
100 cout << " Correction based on Bernstein Poly of degree " << degree << endl;
101
102 RooPlot *frame = x.frame();
103 data->plotOn(frame);
104 // plot the best fit nominal model in blue
106 nominal.fitTo(*data, PrintLevel(0), Minimizer(minimType));
107 nominal.plotOn(frame);
108
109 // plot the best fit corrected model in red
110 RooAbsPdf *corrected = wks->pdf("corrected");
111 if (!corrected)
112 return;
113
114 // fit corrected model
115 corrected->fitTo(*data, PrintLevel(0), Minimizer(minimType));
116 corrected->plotOn(frame, LineColor(kRed));
117
118 // plot the correction term (* norm constant) in dashed green
119 // should make norm constant just be 1, not depend on binning of data
120 RooAbsPdf *poly = wks->pdf("poly");
121 if (poly)
122 poly->plotOn(frame, LineColor(kGreen), LineStyle(kDashed));
123
124 // this is a switch to check the sampling distribution
125 // of -2 log LR for two comparisons:
126 // the first is for n-1 vs. n degree polynomial corrections
127 // the second is for n vs. n+1 degree polynomial corrections
128 // Here we choose n to be the one chosen by the tolerance
129 // criterion above, eg. n = "degree" in the code.
130 // Setting this to true is takes about 10 min.
131 bool checkSamplingDist = true;
132 int numToyMC = 20; // increase this value for sensible results
133
134 TCanvas *c1 = new TCanvas();
135 if (checkSamplingDist) {
136 c1->Divide(1, 2);
137 c1->cd(1);
138 }
139 frame->Draw();
141
142 if (checkSamplingDist) {
143 // check sampling dist
145 TH1F *samplingDist = new TH1F("samplingDist", "", 20, 0, 10);
146 TH1F *samplingDistExtra = new TH1F("samplingDistExtra", "", 20, 0, 10);
147 bernsteinCorrection.CreateQSamplingDist(wks, "nominal", "x", "data", samplingDist, samplingDistExtra, degree,
148 numToyMC);
149
150 c1->cd(2);
151 samplingDistExtra->SetLineColor(kRed);
152 samplingDistExtra->Draw();
153 samplingDist->Draw("same");
154 }
155}
@ kRed
Definition: Rtypes.h:64
@ kGreen
Definition: Rtypes.h:64
@ kDashed
Definition: TAttLine.h:48
void Error(const char *location, const char *msgfmt,...)
static void SetDefaultMinimizer(const char *type, const char *algo=0)
static const std::string & DefaultMinimizerType()
static void SetDefaultPrintLevel(int level)
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Calls RooPlot* plotOn(RooPlot* frame, const RooLinkedList& cmdList) const ;.
Definition: RooAbsData.cxx:546
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1254
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Definition: RooAbsPdf.h:118
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
BernsteinCorrection is a utility in RooStats to augment a nominal PDF with a polynomial correction te...
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
The Canvas class.
Definition: TCanvas.h:27
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:2948
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
Basic string class.
Definition: TString.h:131
RooCmdArg Rename(const char *suffix)
RooConstVar & RooConst(Double_t val)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Minimizer(const char *type, const char *alg=0)
RooCmdArg LineColor(Color_t color)
RooCmdArg LineStyle(Style_t style)
const Double_t sigma
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19
static constexpr double degree