ROOT   6.21/01 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
54 using namespace RooFit;
55 using namespace RooStats;
56
57 //____________________________________
58 void 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 }
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
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:549
RooCmdArg LineColor(Color_t color)
return c1
Definition: legend1.C:41
Definition: Rtypes.h:64
RooCmdArg PrintLevel(Int_t code)
Basic string class.
Definition: TString.h:131
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
int Int_t
Definition: RtypesCore.h:41
Definition: Rtypes.h:64
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
static void SetDefaultPrintLevel(int level)
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
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
static constexpr double degree
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:2981
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:31
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
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
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:2906
RooConstVar & RooConst(Double_t val)
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
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.