Rf 4 0 3_Weightedevts

Data and categories: using weights in unbinned datasets

Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, January 19, 2022 at 10:21 AM.

In [ ]:
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooPolynomial.h"
#include "RooChi2Var.h"
#include "RooMinimizer.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooFitResult.h"
In [ ]:
%%cpp -d
// This is a workaround to make sure the namespace is used inside functions
using namespace RooFit;

Create observable and unweighted dataset

Declare observable

In [ ]:
RooRealVar x("x", "x", -10, 10);
x.setBins(40);

Construction a uniform pdf

In [ ]:
RooPolynomial p0("px", "px", x);

Sample 1000 events from pdf

In [ ]:
RooDataSet *data = p0.generate(x, 1000);

Calculate weight and make dataset weighted

Construct formula to calculate (fake) weight for events

In [ ]:
RooFormulaVar wFunc("w", "event weight", "(x*x+10)", x);

Add column with variable w to previously generated dataset

In [ ]:
RooRealVar *w = (RooRealVar *)data->addColumn(wFunc);

Dataset d is now a dataset with two observable (x,w) with 1000 entries

In [ ]:
data->Print();

Instruct dataset wdata in interpret w as event weight rather than as observable

In [ ]:
RooDataSet wdata(data->GetName(), data->GetTitle(), data, *data->get(), 0, w->GetName());

Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430k

In [ ]:
wdata.Print();

Unbinned ml fit to weighted data

Construction quadratic polynomial pdf for fitting

In [ ]:
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

In [ ]:
RooFitResult *r_ml_wgt = p2.fitTo(wdata, 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:

In [ ]:
RooFitResult *r_ml_wgt_corr = p2.fitTo(wdata, Save(), SumW2Error(kTRUE));

Plot weighed data and fit result

Construct plot frame

In [ ]:
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

In [ ]:
wdata.plotOn(frame, DataError(RooAbsData::SumW2));

Overlay result of 2nd order polynomial fit to weighted data

In [ ]:
p2.plotOn(frame);

Ml fit of pdf to equivalent unweighted dataset

Construct a pdf with the same shape as p0 after weighting

In [ ]:
RooGenericPdf genPdf("genPdf", "x*x+10", x);

Sample a dataset with the same number of events as data

In [ ]:
RooDataSet *data2 = genPdf.generate(x, 1000);

Sample a dataset with the same number of weights as data

In [ ]:
RooDataSet *data3 = genPdf.generate(x, 43000);

Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison

In [ ]:
RooFitResult *r_ml_unw10 = p2.fitTo(*data2, Save());
RooFitResult *r_ml_unw43 = p2.fitTo(*data3, Save());

Chi2 fit of pdf to binned weighted dataset

Construct binned clone of unbinned weighted dataset

In [ ]:
RooDataHist *binnedData = wdata.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

In [ ]:
RooChi2Var chi2("chi2", "chi2", p2, *binnedData, DataError(RooAbsData::SumW2));
RooMinimizer m(chi2);
m.migrad();
m.hesse();

Plot chi^2 fit result on frame as well

In [ ]:
RooFitResult *r_chi2_wgt = m.save();
p2.plotOn(frame, LineStyle(kDashed), LineColor(kRed));

Compare fit results of chi2,ml fits to (un)weighted data

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.

In [ ]:
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();

Draw all canvases

In [ ]:
%jsroot on
gROOT->GetListOfCanvases()->Draw()