Logo ROOT  
Reference Guide
rf403_weightedevts.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roofit
3/// \notebook -js
4///
5/// Data and categories: using weights in unbinned datasets
6///
7/// \macro_image
8/// \macro_output
9/// \macro_code
10///
11/// \date 07/2008
12/// \author Wouter Verkerke
13
14#include "RooRealVar.h"
15#include "RooDataSet.h"
16#include "RooDataHist.h"
17#include "RooGaussian.h"
18#include "RooConstVar.h"
19#include "RooFormulaVar.h"
20#include "RooGenericPdf.h"
21#include "RooPolynomial.h"
22#include "RooChi2Var.h"
23#include "RooMinimizer.h"
24#include "TCanvas.h"
25#include "TAxis.h"
26#include "RooPlot.h"
27#include "RooFitResult.h"
28using namespace RooFit;
29
30void rf403_weightedevts()
31{
32 // C r e a t e o b s e r v a b l e a n d u n w e i g h t e d d a t a s e t
33 // -------------------------------------------------------------------------------
34
35 // Declare observable
36 RooRealVar x("x", "x", -10, 10);
37 x.setBins(40);
38
39 // Construction a uniform pdf
40 RooPolynomial p0("px", "px", x);
41
42 // Sample 1000 events from pdf
43 RooDataSet *data = p0.generate(x, 1000);
44
45 // C a l c u l a t e w e i g h t a n d m a k e d a t a s e t w e i g h t e d
46 // -----------------------------------------------------------------------------------
47
48 // Construct formula to calculate (fake) weight for events
49 RooFormulaVar wFunc("w", "event weight", "(x*x+10)", x);
50
51 // Add column with variable w to previously generated dataset
52 RooRealVar *w = (RooRealVar *)data->addColumn(wFunc);
53
54 // Dataset d is now a dataset with two observable (x,w) with 1000 entries
55 data->Print();
56
57 // Instruct dataset wdata in interpret w as event weight rather than as observable
58 RooDataSet wdata(data->GetName(), data->GetTitle(), data, *data->get(), 0, w->GetName());
59
60 // Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430K
61 wdata.Print();
62
63 // U n b i n n e d M L f i t t o w e i g h t e d d a t a
64 // ---------------------------------------------------------------
65
66 // Construction quadratic polynomial pdf for fitting
67 RooRealVar a0("a0", "a0", 1);
68 RooRealVar a1("a1", "a1", 0, -1, 1);
69 RooRealVar a2("a2", "a2", 1, 0, 10);
70 RooPolynomial p2("p2", "p2", x, RooArgList(a0, a1, a2), 0);
71
72 // Fit quadratic polynomial to weighted data
73
74 // NOTE: A plain Maximum likelihood fit to weighted data does in general
75 // NOT result in correct error estimates, unless individual
76 // event weights represent Poisson statistics themselves.
77 //
78 // Fit with 'wrong' errors
79 RooFitResult *r_ml_wgt = p2.fitTo(wdata, Save());
80
81 // A first order correction to estimated parameter errors in an
82 // (unbinned) ML fit can be obtained by calculating the
83 // covariance matrix as
84 //
85 // V' = V C-1 V
86 //
87 // where V is the covariance matrix calculated from a fit
88 // to -logL = - sum [ w_i log f(x_i) ] and C is the covariance
89 // matrix calculated from -logL' = -sum [ w_i^2 log f(x_i) ]
90 // (i.e. the weights are applied squared)
91 //
92 // A fit in this mode can be performed as follows:
93
94 RooFitResult *r_ml_wgt_corr = p2.fitTo(wdata, Save(), SumW2Error(kTRUE));
95
96 // P l o t w e i g h e d d a t a a n d f i t r e s u l t
97 // ---------------------------------------------------------------
98
99 // Construct plot frame
100 RooPlot *frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data"));
101
102 // Plot data using sum-of-weights-squared error rather than Poisson errors
103 wdata.plotOn(frame, DataError(RooAbsData::SumW2));
104
105 // Overlay result of 2nd order polynomial fit to weighted data
106 p2.plotOn(frame);
107
108 // ML Fit of pdf to equivalent unweighted dataset
109 // -----------------------------------------------------------------------------------------
110
111 // Construct a pdf with the same shape as p0 after weighting
112 RooGenericPdf genPdf("genPdf", "x*x+10", x);
113
114 // Sample a dataset with the same number of events as data
115 RooDataSet *data2 = genPdf.generate(x, 1000);
116
117 // Sample a dataset with the same number of weights as data
118 RooDataSet *data3 = genPdf.generate(x, 43000);
119
120 // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
121 RooFitResult *r_ml_unw10 = p2.fitTo(*data2, Save());
122 RooFitResult *r_ml_unw43 = p2.fitTo(*data3, Save());
123
124 // C h i 2 f i t o f p d f t o b i n n e d w e i g h t e d d a t a s e t
125 // ------------------------------------------------------------------------------------
126
127 // Construct binned clone of unbinned weighted dataset
128 RooDataHist *binnedData = wdata.binnedClone();
129 binnedData->Print("v");
130
131 // Perform chi2 fit to binned weighted dataset using sum-of-weights errors
132 //
133 // NB: Within the usual approximations of a chi2 fit, a chi2 fit to weighted
134 // data using sum-of-weights-squared errors does give correct error
135 // estimates
136 RooChi2Var chi2("chi2", "chi2", p2, *binnedData, DataError(RooAbsData::SumW2));
137 RooMinimizer m(chi2);
138 m.migrad();
139 m.hesse();
140
141 // Plot chi^2 fit result on frame as well
142 RooFitResult *r_chi2_wgt = m.save();
143 p2.plotOn(frame, LineStyle(kDashed), LineColor(kRed));
144
145 // C o m p a r e f i t r e s u l t s o f c h i 2 , M L f i t s t o ( u n ) w e i g h t e d d a t a
146 // ---------------------------------------------------------------------------------------------------------------
147
148 // Note that ML fit on 1Kevt of weighted data is closer to result of ML fit on 43Kevt of unweighted data
149 // than to 1Kevt of unweighted data, whereas the reference chi^2 fit with SumW2 error gives a result closer to
150 // that of an unbinned ML fit to 1Kevt of unweighted data.
151
152 cout << "==> ML Fit results on 1K unweighted events" << endl;
153 r_ml_unw10->Print();
154 cout << "==> ML Fit results on 43K unweighted events" << endl;
155 r_ml_unw43->Print();
156 cout << "==> ML Fit results on 1K weighted events with a summed weight of 43K" << endl;
157 r_ml_wgt->Print();
158 cout << "==> Corrected ML Fit results on 1K weighted events with a summed weight of 43K" << endl;
159 r_ml_wgt_corr->Print();
160 cout << "==> Chi2 Fit results on 1K weighted events with a summed weight of 43K" << endl;
161 r_chi2_wgt->Print();
162
163 new TCanvas("rf403_weightedevts", "rf403_weightedevts", 600, 600);
164 gPad->SetLeftMargin(0.15);
165 frame->GetYaxis()->SetTitleOffset(1.8);
166 frame->Draw();
167}
const Bool_t kTRUE
Definition: RtypesCore.h:89
@ kRed
Definition: Rtypes.h:64
@ kDashed
Definition: TAttLine.h:48
#define gPad
Definition: TVirtualPad.h:287
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:175
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Class RooChi2Var implements a simple chi^2 calculation from a binned dataset and a PDF.
Definition: RooChi2Var.h:25
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
virtual RooAbsArg * addColumn(RooAbsArg &var, Bool_t adjustRange=kTRUE)
Add a column with the values of the given (function) argument to this dataset.
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Definition: RooFormulaVar.h:29
RooGenericPdf is a concrete implementation of a probability density function, which takes a RooArgLis...
Definition: RooGenericPdf.h:25
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:38
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
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
RooPolynomial implements a polynomial p.d.f of the form.
Definition: RooPolynomial.h:28
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
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooCmdArg SumW2Error(Bool_t flag)
RooCmdArg Save(Bool_t flag=kTRUE)
RooCmdArg DataError(Int_t)
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...
const char * Title
Definition: TXMLSetup.cxx:67
auto * m
Definition: textangle.C:8