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