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