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