Logo ROOT  
Reference Guide
rf610_visualerror.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 /// Likelihood and minimization: visualization of errors from a covariance matrix
5 ///
6 /// \macro_image
7 /// \macro_output
8 /// \macro_code
9 ///
10 /// \date April 2009
11 /// \author Wouter Verkerke
12 
13 #include "RooRealVar.h"
14 #include "RooDataHist.h"
15 #include "RooGaussian.h"
16 #include "RooConstVar.h"
17 #include "RooAddPdf.h"
18 #include "RooPlot.h"
19 #include "TCanvas.h"
20 #include "TAxis.h"
21 #include "TAxis.h"
22 using namespace RooFit;
23 
24 void rf610_visualerror()
25 {
26  // S e t u p e x a m p l e f i t
27  // ---------------------------------------
28 
29  // Create sum of two Gaussians pdf with factory
30  RooRealVar x("x", "x", -10, 10);
31 
32  RooRealVar m("m", "m", 0, -10, 10);
33  RooRealVar s("s", "s", 2, 1, 50);
34  RooGaussian sig("sig", "sig", x, m, s);
35 
36  RooRealVar m2("m2", "m2", -1, -10, 10);
37  RooRealVar s2("s2", "s2", 6, 1, 50);
38  RooGaussian bkg("bkg", "bkg", x, m2, s2);
39 
40  RooRealVar fsig("fsig", "fsig", 0.33, 0, 1);
41  RooAddPdf model("model", "model", RooArgList(sig, bkg), fsig);
42 
43  // Create binned dataset
44  x.setBins(25);
45  RooAbsData *d = model.generateBinned(x, 1000);
46 
47  // Perform fit and save fit result
48  RooFitResult *r = model.fitTo(*d, Save());
49 
50  // V i s u a l i z e f i t e r r o r
51  // -------------------------------------
52 
53  // Make plot frame
54  RooPlot *frame = x.frame(Bins(40), Title("P.d.f with visualized 1-sigma error band"));
55  d->plotOn(frame);
56 
57  // Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation
58  // This results in an error band that is by construction symmetric
59  //
60  // The linear error is calculated as
61  // error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
62  //
63  // where F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2,
64  //
65  // with f(x) = the plotted curve
66  // 'da' = error taken from the fit result
67  // Corr(a,a') = the correlation matrix from the fit result
68  // Z = requested significance 'Z sigma band'
69  //
70  // The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters),
71  // but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and
72  // Gaussian approximations made
73  //
74  model.plotOn(frame, VisualizeError(*r, 1), FillColor(kOrange));
75 
76  // Calculate error using sampling method and visualize as dashed red line.
77  //
78  // In this method a number of curves is calculated with variations of the parameter values, as sampled
79  // from a multi-variate Gaussian pdf that is constructed from the fit results covariance matrix.
80  // The error(x) is determined by calculating a central interval that capture N% of the variations
81  // for each value of x, where N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves
82  // is chosen to be such that at least 100 curves are expected to be outside the N% interval, and is minimally
83  // 100 (e.g. Z=1->Ncurve=356, Z=2->Ncurve=2156)) Intervals from the sampling method can be asymmetric,
84  // and may perform better in the presence of strong correlations, but may take (much) longer to calculate
85  model.plotOn(frame, VisualizeError(*r, 1, kFALSE), DrawOption("L"), LineWidth(2), LineColor(kRed));
86 
87  // Perform the same type of error visualization on the background component only.
88  // The VisualizeError() option can generally applied to _any_ kind of plot (components, asymmetries, efficiencies
89  // etc..)
90  model.plotOn(frame, VisualizeError(*r, 1), FillColor(kOrange), Components("bkg"));
91  model.plotOn(frame, VisualizeError(*r, 1, kFALSE), DrawOption("L"), LineWidth(2), LineColor(kRed), Components("bkg"),
93 
94  // Overlay central value
95  model.plotOn(frame);
96  model.plotOn(frame, Components("bkg"), LineStyle(kDashed));
97  d->plotOn(frame);
98  frame->SetMinimum(0);
99 
100  // V i s u a l i z e p a r t i a l f i t e r r o r
101  // ------------------------------------------------------
102 
103  // Make plot frame
104  RooPlot *frame2 = x.frame(Bins(40), Title("Visualization of 2-sigma partial error from (m,m2)"));
105 
106  // Visualize partial error. For partial error visualization the covariance matrix is first reduced as follows
107  // ___ -1
108  // Vred = V22 = V11 - V12 * V22 * V21
109  //
110  // Where V11,V12,V21,V22 represent a block decomposition of the covariance matrix into observables that
111  // are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), and V22bar
112  // is the Shur complement of V22, calculated as shown above
113  //
114  // (Note that Vred is _not_ a simple sub-matrix of V)
115 
116  // Propagate partial error due to shape parameters (m,m2) using linear and sampling method
117  model.plotOn(frame2, VisualizeError(*r, RooArgSet(m, m2), 2), FillColor(kCyan));
118  model.plotOn(frame2, Components("bkg"), VisualizeError(*r, RooArgSet(m, m2), 2), FillColor(kCyan));
119 
120  model.plotOn(frame2);
121  model.plotOn(frame2, Components("bkg"), LineStyle(kDashed));
122  frame2->SetMinimum(0);
123 
124  // Make plot frame
125  RooPlot *frame3 = x.frame(Bins(40), Title("Visualization of 2-sigma partial error from (s,s2)"));
126 
127  // Propagate partial error due to yield parameter using linear and sampling method
128  model.plotOn(frame3, VisualizeError(*r, RooArgSet(s, s2), 2), FillColor(kGreen));
129  model.plotOn(frame3, Components("bkg"), VisualizeError(*r, RooArgSet(s, s2), 2), FillColor(kGreen));
130 
131  model.plotOn(frame3);
132  model.plotOn(frame3, Components("bkg"), LineStyle(kDashed));
133  frame3->SetMinimum(0);
134 
135  // Make plot frame
136  RooPlot *frame4 = x.frame(Bins(40), Title("Visualization of 2-sigma partial error from fsig"));
137 
138  // Propagate partial error due to yield parameter using linear and sampling method
139  model.plotOn(frame4, VisualizeError(*r, RooArgSet(fsig), 2), FillColor(kMagenta));
140  model.plotOn(frame4, Components("bkg"), VisualizeError(*r, RooArgSet(fsig), 2), FillColor(kMagenta));
141 
142  model.plotOn(frame4);
143  model.plotOn(frame4, Components("bkg"), LineStyle(kDashed));
144  frame4->SetMinimum(0);
145 
146  TCanvas *c = new TCanvas("rf610_visualerror", "rf610_visualerror", 800, 800);
147  c->Divide(2, 2);
148  c->cd(1);
149  gPad->SetLeftMargin(0.15);
150  frame->GetYaxis()->SetTitleOffset(1.4);
151  frame->Draw();
152  c->cd(2);
153  gPad->SetLeftMargin(0.15);
154  frame2->GetYaxis()->SetTitleOffset(1.6);
155  frame2->Draw();
156  c->cd(3);
157  gPad->SetLeftMargin(0.15);
158  frame3->GetYaxis()->SetTitleOffset(1.6);
159  frame3->Draw();
160  c->cd(4);
161  gPad->SetLeftMargin(0.15);
162  frame4->GetYaxis()->SetTitleOffset(1.6);
163  frame4->Draw();
164 }
c
#define c(i)
Definition: RSha256.hxx:101
m
auto * m
Definition: textangle.C:8
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
RooAddPdf
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
Definition: RooAddPdf.h:32
RooFit::Bins
RooCmdArg Bins(Int_t nbin)
Definition: RooGlobalFunc.cxx:177
RooAbsData
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:47
kGreen
@ kGreen
Definition: Rtypes.h:66
RooFit::VisualizeError
RooCmdArg VisualizeError(const RooDataSet &paramData, Double_t Z=1)
Definition: RooGlobalFunc.cxx:71
r
ROOT::R::TRInterface & r
Definition: Object.C:4
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooGaussian.h
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
kOrange
@ kOrange
Definition: Rtypes.h:67
x
Double_t x[n]
Definition: legend1.C:17
RooGaussian
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
RooAddPdf.h
TCanvas.h
RooFit::FillColor
RooCmdArg FillColor(Color_t color)
Definition: RooGlobalFunc.cxx:60
rf610_visualerror
Definition: rf610_visualerror.py:1
RooFitResult
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
TGeant4Unit::m2
static constexpr double m2
Definition: TGeant4SystemOfUnits.h:123
kCyan
@ kCyan
Definition: Rtypes.h:66
RooFit::DrawOption
RooCmdArg DrawOption(const char *opt)
Definition: RooGlobalFunc.cxx:39
kMagenta
@ kMagenta
Definition: Rtypes.h:66
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooDataHist.h
RooPlot.h
RooPlot::GetYaxis
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1258
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
RooRealVar.h
kRed
@ kRed
Definition: Rtypes.h:66
RooConstVar.h
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:57
TCanvas
The Canvas class.
Definition: TCanvas.h:23
TAxis.h
RooPlot::SetMinimum
virtual void SetMinimum(Double_t minimum=-1111)
Set minimum value of Y axis.
Definition: RooPlot.cxx:1091
d
#define d(i)
Definition: RSha256.hxx:102
kDashed
@ kDashed
Definition: TAttLine.h:48
gPad
#define gPad
Definition: TVirtualPad.h:287
make_cnn_model.model
model
Definition: make_cnn_model.py:6
RooFit::LineWidth
RooCmdArg LineWidth(Width_t width)
Definition: RooGlobalFunc.cxx:59
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition: RooGlobalFunc.cxx:58
TAttAxis::SetTitleOffset
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition: TAttAxis.cxx:293
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:37
RooFit::Components
RooCmdArg Components(const RooArgSet &compSet)
Definition: RooGlobalFunc.cxx:75
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:190
RooFit::Title
RooCmdArg Title(const char *name)
Definition: RooGlobalFunc.cxx:176
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29