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"
28 using namespace RooFit;
29 
30 void 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 }
RooFormulaVar.h
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
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooAbsData::SumW2
@ SumW2
Definition: RooAbsData.h:96
RooFit::DataError
RooCmdArg DataError(Int_t)
Definition: RooGlobalFunc.cxx:156
TNamed::GetTitle
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:54
RooFit::SumW2Error
RooCmdArg SumW2Error(Bool_t flag)
Definition: RooGlobalFunc.cxx:207
RooArgList
Definition: RooArgList.h:21
RooGaussian.h
RooChi2Var
Definition: RooChi2Var.h:25
RooAbsData::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:177
x
Double_t x[n]
Definition: legend1.C:17
TCanvas.h
RooGenericPdf.h
RooDataSet.h
RooFitResult
Definition: RooFitResult.h:40
RooChi2Var.h
RooDataSet::get
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Definition: RooDataSet.cxx:1038
RooPolynomial.h
RooDataHist
Definition: RooDataHist.h:39
RooFormulaVar
Definition: RooFormulaVar.h:30
RooFit
Definition: RooCFunction1Binding.h:29
RooPolynomial
Definition: RooPolynomial.h:28
RooDataHist.h
RooPlot.h
RooPlot::GetYaxis
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1256
RooPlot
Definition: RooPlot.h:44
RooRealVar.h
kRed
@ kRed
Definition: Rtypes.h:66
RooFitResult.h
RooConstVar.h
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:56
TCanvas
Definition: TCanvas.h:23
RooMinimizer.h
TAxis.h
RooGenericPdf
Definition: RooGenericPdf.h:25
RooFitResult::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
kDashed
@ kDashed
Definition: TAttLine.h:48
gPad
#define gPad
Definition: TVirtualPad.h:287
RooDataSet
Definition: RooDataSet.h:33
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition: RooGlobalFunc.cxx:57
TAttAxis::SetTitleOffset
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition: TAttAxis.cxx:293
RooRealVar
Definition: RooRealVar.h:35
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
RooMinimizer
Definition: RooMinimizer.h:40
RooFit::Title
RooCmdArg Title(const char *name)
Definition: RooGlobalFunc.cxx:173
RooDataSet::addColumn
virtual RooAbsArg * addColumn(RooAbsArg &var, Bool_t adjustRange=kTRUE)
Add a column with the values of the given (function) argument to this dataset.
Definition: RooDataSet.cxx:1383