ROOT logo
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #403
// 
// Using weights in unbinned datasets
//
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooPolynomial.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooFitResult.h"
using namespace RooFit ;


void rf403_weightedevts()
{
  // 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 
  // -------------------------------------------------------------------------------

  // Declare observable
  RooRealVar x("x","x",-10,10) ;
  x.setBins(40) ;

  // Construction a uniform pdf
  RooPolynomial p0("px","px",x) ;

  // Sample 1000 events from pdf
  RooDataSet* data = p0.generate(x,1000) ;

 

  // 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 
  // -----------------------------------------------------------------------------------

  // Construct formula to calculate (fake) weight for events
  RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;

  // Add column with variable w to previously generated dataset
  RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;

  // Dataset d is now a dataset with two observable (x,w) with 1000 entries
  data->Print() ;

  // Instruct dataset d in interpret w as event weight rather than as observable
  data->setWeightVar(*w) ;

  // Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430K
  data->Print() ;



  // 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 
  // ---------------------------------------------------------------

  // Construction quadratic polynomial pdf for fitting
  RooRealVar a0("a0","a0",1) ;
  RooRealVar a1("a1","a1",0,-1,1) ;
  RooRealVar a2("a2","a2",1,0,10) ;
  RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;

  // Fit quadratic polynomial to weighted data

  // NOTE: A plain Maximum likelihood fit to weighted data does in general 
  //       NOT result in correct error estimates, unless individual
  //       event weights represent Poisson statistics themselves.
  //       
  // Fit with 'wrong' errors
  RooFitResult* r_ml_wgt = p2.fitTo(*data,Save()) ;
  
  // A first order correction to estimated parameter errors in an 
  // (unbinned) ML fit can be obtained by calculating the
  // covariance matrix as
  //
  //    V' = V C-1 V
  //
  // where V is the covariance matrix calculated from a fit
  // to -logL = - sum [ w_i log f(x_i) ] and C is the covariance
  // matrix calculated from -logL' = -sum [ w_i^2 log f(x_i) ] 
  // (i.e. the weights are applied squared)
  //
  // A fit in this mode can be performed as follows:

  RooFitResult* r_ml_wgt_corr = p2.fitTo(*data,Save(),SumW2Error(kTRUE)) ;



  // 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 
  // ---------------------------------------------------------------

  // Construct plot frame
  RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;

  // Plot data using sum-of-weights-squared error rather than Poisson errors
  data->plotOn(frame,DataError(RooAbsData::SumW2)) ;

  // Overlay result of 2nd order polynomial fit to weighted data
  p2.plotOn(frame) ;



  // M L  F i t   o f   p d f   t o   e q u i v a l e n t  u n w e i g h t e d   d a t a s e t
  // -----------------------------------------------------------------------------------------
  
  // Construct a pdf with the same shape as p0 after weighting
  RooGenericPdf genPdf("genPdf","x*x+10",x) ;

  // Sample a dataset with the same number of events as data
  RooDataSet* data2 = genPdf.generate(x,1000) ;

  // Sample a dataset with the same number of weights as data
  RooDataSet* data3 = genPdf.generate(x,43000) ;

  // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
  RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
  RooFitResult* r_ml_unw43 = p2.fitTo(*data3,Save()) ;


  // 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
  // ------------------------------------------------------------------------------------

  // Construct binned clone of unbinned weighted dataset
  RooDataHist* binnedData = data->binnedClone() ;
  binnedData->Print("v") ;

  // Perform chi2 fit to binned weighted dataset using sum-of-weights errors
  // 
  // NB: Within the usual approximations of a chi2 fit, a chi2 fit to weighted
  // data using sum-of-weights-squared errors does give correct error
  // estimates
  RooChi2Var chi2("chi2","chi2",p2,*binnedData,DataError(RooAbsData::SumW2)) ;
  RooMinuit m(chi2) ;
  m.migrad() ;
  m.hesse() ;

  // Plot chi^2 fit result on frame as well
  RooFitResult* r_chi2_wgt = m.save() ;
  p2.plotOn(frame,LineStyle(kDashed),LineColor(kRed)) ;



  // 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 
  // ---------------------------------------------------------------------------------------------------------------

  // Note that ML fit on 1Kevt of weighted data is closer to result of ML fit on 43Kevt of unweighted data 
  // than to 1Kevt of unweighted data, whereas the reference chi^2 fit with SumW2 error gives a result closer to
  // that of an unbinned ML fit to 1Kevt of unweighted data. 

  cout << "==> ML Fit results on 1K unweighted events" << endl ;
  r_ml_unw10->Print() ;
  cout << "==> ML Fit results on 43K unweighted events" << endl ;
  r_ml_unw43->Print() ;
  cout << "==> ML Fit results on 1K weighted events with a summed weight of 43K" << endl ;
  r_ml_wgt->Print() ;
  cout << "==> Corrected ML Fit results on 1K weighted events with a summed weight of 43K" << endl ;
  r_ml_wgt_corr->Print() ;
  cout << "==> Chi2 Fit results on 1K weighted events with a summed weight of 43K" << endl ;
  r_chi2_wgt->Print() ;


  new TCanvas("rf403_weightedevts","rf403_weightedevts",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.8) ; frame->Draw() ;


}
 rf403_weightedevts.C:1
 rf403_weightedevts.C:2
 rf403_weightedevts.C:3
 rf403_weightedevts.C:4
 rf403_weightedevts.C:5
 rf403_weightedevts.C:6
 rf403_weightedevts.C:7
 rf403_weightedevts.C:8
 rf403_weightedevts.C:9
 rf403_weightedevts.C:10
 rf403_weightedevts.C:11
 rf403_weightedevts.C:12
 rf403_weightedevts.C:13
 rf403_weightedevts.C:14
 rf403_weightedevts.C:15
 rf403_weightedevts.C:16
 rf403_weightedevts.C:17
 rf403_weightedevts.C:18
 rf403_weightedevts.C:19
 rf403_weightedevts.C:20
 rf403_weightedevts.C:21
 rf403_weightedevts.C:22
 rf403_weightedevts.C:23
 rf403_weightedevts.C:24
 rf403_weightedevts.C:25
 rf403_weightedevts.C:26
 rf403_weightedevts.C:27
 rf403_weightedevts.C:28
 rf403_weightedevts.C:29
 rf403_weightedevts.C:30
 rf403_weightedevts.C:31
 rf403_weightedevts.C:32
 rf403_weightedevts.C:33
 rf403_weightedevts.C:34
 rf403_weightedevts.C:35
 rf403_weightedevts.C:36
 rf403_weightedevts.C:37
 rf403_weightedevts.C:38
 rf403_weightedevts.C:39
 rf403_weightedevts.C:40
 rf403_weightedevts.C:41
 rf403_weightedevts.C:42
 rf403_weightedevts.C:43
 rf403_weightedevts.C:44
 rf403_weightedevts.C:45
 rf403_weightedevts.C:46
 rf403_weightedevts.C:47
 rf403_weightedevts.C:48
 rf403_weightedevts.C:49
 rf403_weightedevts.C:50
 rf403_weightedevts.C:51
 rf403_weightedevts.C:52
 rf403_weightedevts.C:53
 rf403_weightedevts.C:54
 rf403_weightedevts.C:55
 rf403_weightedevts.C:56
 rf403_weightedevts.C:57
 rf403_weightedevts.C:58
 rf403_weightedevts.C:59
 rf403_weightedevts.C:60
 rf403_weightedevts.C:61
 rf403_weightedevts.C:62
 rf403_weightedevts.C:63
 rf403_weightedevts.C:64
 rf403_weightedevts.C:65
 rf403_weightedevts.C:66
 rf403_weightedevts.C:67
 rf403_weightedevts.C:68
 rf403_weightedevts.C:69
 rf403_weightedevts.C:70
 rf403_weightedevts.C:71
 rf403_weightedevts.C:72
 rf403_weightedevts.C:73
 rf403_weightedevts.C:74
 rf403_weightedevts.C:75
 rf403_weightedevts.C:76
 rf403_weightedevts.C:77
 rf403_weightedevts.C:78
 rf403_weightedevts.C:79
 rf403_weightedevts.C:80
 rf403_weightedevts.C:81
 rf403_weightedevts.C:82
 rf403_weightedevts.C:83
 rf403_weightedevts.C:84
 rf403_weightedevts.C:85
 rf403_weightedevts.C:86
 rf403_weightedevts.C:87
 rf403_weightedevts.C:88
 rf403_weightedevts.C:89
 rf403_weightedevts.C:90
 rf403_weightedevts.C:91
 rf403_weightedevts.C:92
 rf403_weightedevts.C:93
 rf403_weightedevts.C:94
 rf403_weightedevts.C:95
 rf403_weightedevts.C:96
 rf403_weightedevts.C:97
 rf403_weightedevts.C:98
 rf403_weightedevts.C:99
 rf403_weightedevts.C:100
 rf403_weightedevts.C:101
 rf403_weightedevts.C:102
 rf403_weightedevts.C:103
 rf403_weightedevts.C:104
 rf403_weightedevts.C:105
 rf403_weightedevts.C:106
 rf403_weightedevts.C:107
 rf403_weightedevts.C:108
 rf403_weightedevts.C:109
 rf403_weightedevts.C:110
 rf403_weightedevts.C:111
 rf403_weightedevts.C:112
 rf403_weightedevts.C:113
 rf403_weightedevts.C:114
 rf403_weightedevts.C:115
 rf403_weightedevts.C:116
 rf403_weightedevts.C:117
 rf403_weightedevts.C:118
 rf403_weightedevts.C:119
 rf403_weightedevts.C:120
 rf403_weightedevts.C:121
 rf403_weightedevts.C:122
 rf403_weightedevts.C:123
 rf403_weightedevts.C:124
 rf403_weightedevts.C:125
 rf403_weightedevts.C:126
 rf403_weightedevts.C:127
 rf403_weightedevts.C:128
 rf403_weightedevts.C:129
 rf403_weightedevts.C:130
 rf403_weightedevts.C:131
 rf403_weightedevts.C:132
 rf403_weightedevts.C:133
 rf403_weightedevts.C:134
 rf403_weightedevts.C:135
 rf403_weightedevts.C:136
 rf403_weightedevts.C:137
 rf403_weightedevts.C:138
 rf403_weightedevts.C:139
 rf403_weightedevts.C:140
 rf403_weightedevts.C:141
 rf403_weightedevts.C:142
 rf403_weightedevts.C:143
 rf403_weightedevts.C:144
 rf403_weightedevts.C:145
 rf403_weightedevts.C:146
 rf403_weightedevts.C:147
 rf403_weightedevts.C:148
 rf403_weightedevts.C:149
 rf403_weightedevts.C:150
 rf403_weightedevts.C:151
 rf403_weightedevts.C:152
 rf403_weightedevts.C:153
 rf403_weightedevts.C:154
 rf403_weightedevts.C:155
 rf403_weightedevts.C:156
 rf403_weightedevts.C:157
 rf403_weightedevts.C:158
 rf403_weightedevts.C:159
 rf403_weightedevts.C:160
 rf403_weightedevts.C:161
 rf403_weightedevts.C:162
 rf403_weightedevts.C:163
 rf403_weightedevts.C:164
 rf403_weightedevts.C:165
 rf403_weightedevts.C:166
 rf403_weightedevts.C:167
 rf403_weightedevts.C:168
 rf403_weightedevts.C:169
 rf403_weightedevts.C:170
 rf403_weightedevts.C:171
 rf403_weightedevts.C:172
 rf403_weightedevts.C:173
 rf403_weightedevts.C:174
 rf403_weightedevts.C:175
 rf403_weightedevts.C:176
 rf403_weightedevts.C:177
 rf403_weightedevts.C:178
 rf403_weightedevts.C:179
 rf403_weightedevts.C:180
 rf403_weightedevts.C:181
 rf403_weightedevts.C:182
 rf403_weightedevts.C:183