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