## From \$ROOTSYS/tutorials/roofit/rf610_visualerror.C

```/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #610
//
// Visualization of errors from a covariance matrix
//
//
//
// 04/2009 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TAxis.h"
using namespace RooFit ;

void rf610_visualerror()
{
// S e t u p   e x a m p l e   f i t
// ---------------------------------------

// Create sum of two Gaussians p.d.f. with factory
RooRealVar x("x","x",-10,10) ;

RooRealVar m("m","m",0,-10,10) ;
RooRealVar s("s","s",2,1,50) ;
RooGaussian sig("sig","sig",x,m,s) ;

RooRealVar m2("m2","m2",-1,-10,10) ;
RooRealVar s2("s2","s2",6,1,50) ;
RooGaussian bkg("bkg","bkg",x,m2,s2) ;

RooRealVar fsig("fsig","fsig",0.33,0,1) ;

// Create binned dataset
x.setBins(25) ;
RooAbsData* d = model.generateBinned(x,1000) ;

// Perform fit and save fit result
RooFitResult* r = model.fitTo(*d,Save()) ;

// V i s u a l i z e   f i t   e r r o r
// -------------------------------------

// Make plot frame
RooPlot* frame = x.frame(Bins(40),Title("P.d.f with visualized 1-sigma error band")) ;
d->plotOn(frame) ;

// Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation
// This results in an error band that is by construction symmetric
//
// The linear error is calculated as
// error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
//
// where     F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2,
//
//         with f(x) = the plotted curve
//              'da' = error taken from the fit result
//        Corr(a,a') = the correlation matrix from the fit result
//                Z = requested significance 'Z sigma band'
//
// The linear method is fast (required 2*N evaluations of the curve, where N is the number of parameters),
// but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and
// Gaussian approximations made
//
model.plotOn(frame,VisualizeError(*r,1),FillColor(kOrange)) ;

// Calculate error using sampling method and visualize as dashed red line.
//
// In this method a number of curves is calculated with variations of the parameter values, as sampled
// from a multi-variate Gaussian p.d.f. that is constructed from the fit results covariance matrix.
// The error(x) is determined by calculating a central interval that capture N% of the variations
// for each valye of x, where N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves
// is chosen to be such that at least 100 curves are expected to be outside the N% interval, and is minimally
// 100 (e.g. Z=1->Ncurve=356, Z=2->Ncurve=2156)) Intervals from the sampling method can be asymmetric,
// and may perform better in the presence of strong correlations, but may take (much) longer to calculate
model.plotOn(frame,VisualizeError(*r,1,kFALSE),DrawOption("L"),LineWidth(2),LineColor(kRed)) ;

// Perform the same type of error visualization on the background component only.
// The VisualizeError() option can generally applied to _any_ kind of plot (components, asymmetries, efficiencies etc..)
model.plotOn(frame,VisualizeError(*r,1),FillColor(kOrange),Components("bkg")) ;
model.plotOn(frame,VisualizeError(*r,1,kFALSE),DrawOption("L"),LineWidth(2),LineColor(kRed),Components("bkg"),LineStyle(kDashed)) ;

// Overlay central value
model.plotOn(frame) ;
model.plotOn(frame,Components("bkg"),LineStyle(kDashed)) ;
d->plotOn(frame) ;
frame->SetMinimum(0) ;

// V i s u a l i z e   p a r t i a l   f i t   e r r o r
// ------------------------------------------------------

// Make plot frame
RooPlot* frame2 = x.frame(Bins(40),Title("Visualization of 2-sigma partial error from (m,m2)")) ;

// Visualize partial error. For partial error visualization the covariance matrix is first reduced as follows
//        ___                   -1
// Vred = V22  = V11 - V12 * V22   * V21
//
// Where V11,V12,V21,V22 represent a block decomposition of the covariance matrix into observables that
// are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), and V22bar
// is the Shur complement of V22, calculated as shown above
//
// (Note that Vred is _not_ a simple sub-matrix of V)

// Propagate partial error due to shape parameters (m,m2) using linear and sampling method
model.plotOn(frame2,VisualizeError(*r,RooArgSet(m,m2),2),FillColor(kCyan)) ;
model.plotOn(frame2,Components("bkg"),VisualizeError(*r,RooArgSet(m,m2),2),FillColor(kCyan)) ;

model.plotOn(frame2) ;
model.plotOn(frame2,Components("bkg"),LineStyle(kDashed)) ;
frame2->SetMinimum(0) ;

// Make plot frame
RooPlot* frame3 = x.frame(Bins(40),Title("Visualization of 2-sigma partial error from (s,s2)")) ;

// Propagate partial error due to yield parameter using linear and sampling method
model.plotOn(frame3,VisualizeError(*r,RooArgSet(s,s2),2),FillColor(kGreen)) ;
model.plotOn(frame3,Components("bkg"),VisualizeError(*r,RooArgSet(s,s2),2),FillColor(kGreen)) ;

model.plotOn(frame3) ;
model.plotOn(frame3,Components("bkg"),LineStyle(kDashed)) ;
frame3->SetMinimum(0) ;

// Make plot frame
RooPlot* frame4 = x.frame(Bins(40),Title("Visualization of 2-sigma partial error from fsig")) ;

// Propagate partial error due to yield parameter using linear and sampling method
model.plotOn(frame4,VisualizeError(*r,RooArgSet(fsig),2),FillColor(kMagenta)) ;
model.plotOn(frame4,Components("bkg"),VisualizeError(*r,RooArgSet(fsig),2),FillColor(kMagenta)) ;

model.plotOn(frame4) ;
model.plotOn(frame4,Components("bkg"),LineStyle(kDashed)) ;
frame4->SetMinimum(0) ;

TCanvas* c = new TCanvas("rf610_visualerror","rf610_visualerror",800,800) ;
c->Divide(2,2) ;
c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4)  ; frame->Draw() ;
c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; frame2->Draw() ;
c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; frame3->Draw() ;
c->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.6) ; frame4->Draw() ;
}
```
rf610_visualerror.C:1
rf610_visualerror.C:2
rf610_visualerror.C:3
rf610_visualerror.C:4
rf610_visualerror.C:5
rf610_visualerror.C:6
rf610_visualerror.C:7
rf610_visualerror.C:8
rf610_visualerror.C:9
rf610_visualerror.C:10
rf610_visualerror.C:11
rf610_visualerror.C:12
rf610_visualerror.C:13
rf610_visualerror.C:14
rf610_visualerror.C:15
rf610_visualerror.C:16
rf610_visualerror.C:17
rf610_visualerror.C:18
rf610_visualerror.C:19
rf610_visualerror.C:20
rf610_visualerror.C:21
rf610_visualerror.C:22
rf610_visualerror.C:23
rf610_visualerror.C:24
rf610_visualerror.C:25
rf610_visualerror.C:26
rf610_visualerror.C:27
rf610_visualerror.C:28
rf610_visualerror.C:29
rf610_visualerror.C:30
rf610_visualerror.C:31
rf610_visualerror.C:32
rf610_visualerror.C:33
rf610_visualerror.C:34
rf610_visualerror.C:35
rf610_visualerror.C:36
rf610_visualerror.C:37
rf610_visualerror.C:38
rf610_visualerror.C:39
rf610_visualerror.C:40
rf610_visualerror.C:41
rf610_visualerror.C:42
rf610_visualerror.C:43
rf610_visualerror.C:44
rf610_visualerror.C:45
rf610_visualerror.C:46
rf610_visualerror.C:47
rf610_visualerror.C:48
rf610_visualerror.C:49
rf610_visualerror.C:50
rf610_visualerror.C:51
rf610_visualerror.C:52
rf610_visualerror.C:53
rf610_visualerror.C:54
rf610_visualerror.C:55
rf610_visualerror.C:56
rf610_visualerror.C:57
rf610_visualerror.C:58
rf610_visualerror.C:59
rf610_visualerror.C:60
rf610_visualerror.C:61
rf610_visualerror.C:62
rf610_visualerror.C:63
rf610_visualerror.C:64
rf610_visualerror.C:65
rf610_visualerror.C:66
rf610_visualerror.C:67
rf610_visualerror.C:68
rf610_visualerror.C:69
rf610_visualerror.C:70
rf610_visualerror.C:71
rf610_visualerror.C:72
rf610_visualerror.C:73
rf610_visualerror.C:74
rf610_visualerror.C:75
rf610_visualerror.C:76
rf610_visualerror.C:77
rf610_visualerror.C:78
rf610_visualerror.C:79
rf610_visualerror.C:80
rf610_visualerror.C:81
rf610_visualerror.C:82
rf610_visualerror.C:83
rf610_visualerror.C:84
rf610_visualerror.C:85
rf610_visualerror.C:86
rf610_visualerror.C:87
rf610_visualerror.C:88
rf610_visualerror.C:89
rf610_visualerror.C:90
rf610_visualerror.C:91
rf610_visualerror.C:92
rf610_visualerror.C:93
rf610_visualerror.C:94
rf610_visualerror.C:95
rf610_visualerror.C:96
rf610_visualerror.C:97
rf610_visualerror.C:98
rf610_visualerror.C:99
rf610_visualerror.C:100
rf610_visualerror.C:101
rf610_visualerror.C:102
rf610_visualerror.C:103
rf610_visualerror.C:104
rf610_visualerror.C:105
rf610_visualerror.C:106
rf610_visualerror.C:107
rf610_visualerror.C:108
rf610_visualerror.C:109
rf610_visualerror.C:110
rf610_visualerror.C:111
rf610_visualerror.C:112
rf610_visualerror.C:113
rf610_visualerror.C:114
rf610_visualerror.C:115
rf610_visualerror.C:116
rf610_visualerror.C:117
rf610_visualerror.C:118
rf610_visualerror.C:119
rf610_visualerror.C:120
rf610_visualerror.C:121
rf610_visualerror.C:122
rf610_visualerror.C:123
rf610_visualerror.C:124
rf610_visualerror.C:125
rf610_visualerror.C:126
rf610_visualerror.C:127
rf610_visualerror.C:128
rf610_visualerror.C:129
rf610_visualerror.C:130
rf610_visualerror.C:131
rf610_visualerror.C:132
rf610_visualerror.C:133
rf610_visualerror.C:134
rf610_visualerror.C:135
rf610_visualerror.C:136
rf610_visualerror.C:137
rf610_visualerror.C:138
rf610_visualerror.C:139
rf610_visualerror.C:140
rf610_visualerror.C:141
rf610_visualerror.C:142
rf610_visualerror.C:143
rf610_visualerror.C:144
rf610_visualerror.C:145
rf610_visualerror.C:146
rf610_visualerror.C:147
rf610_visualerror.C:148
rf610_visualerror.C:149
rf610_visualerror.C:150
rf610_visualerror.C:151
rf610_visualerror.C:152
rf610_visualerror.C:153
rf610_visualerror.C:154
rf610_visualerror.C:155
rf610_visualerror.C:156
rf610_visualerror.C:157
rf610_visualerror.C:158
rf610_visualerror.C:159
rf610_visualerror.C:160
rf610_visualerror.C:161
rf610_visualerror.C:162