doc/master/rf403__weightedevts_8C.html

#include "RooRealVar.h"

#include "RooDataSet.h"

#include "RooDataHist.h"

#include "RooGaussian.h"

#include "RooFormulaVar.h"

#include "RooGenericPdf.h"

#include "RooPolynomial.h"

#include "RooMinimizer.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

   std::unique_ptr<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 wdata in interpret w as event weight rather than as observable

   RooDataSet wdata(data->GetName(), data->GetTitle(), data.get(), *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

   wdata.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

   std::unique_ptr<RooFitResult> r_ml_wgt{p2.fitTo(wdata, Save(), PrintLevel(-1))};


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


   std::unique_ptr<RooFitResult> r_ml_wgt_corr{p2.fitTo(wdata, Save(), SumW2Error(true), PrintLevel(-1))};


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

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


   // Overlay result of 2nd order polynomial fit to weighted data

   p2.plotOn(frame);


   // ML Fit of pdf to equivalent unweighted dataset

   // -----------------------------------------------------------------------------------------


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

   std::unique_ptr<RooDataSet> data2{genPdf.generate(x, 1000)};


   // Sample a dataset with the same number of weights as data

   std::unique_ptr<RooDataSet> data3{genPdf.generate(x, 43000)};


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

   std::unique_ptr<RooFitResult> r_ml_unw10{p2.fitTo(*data2, Save(), PrintLevel(-1))};

   std::unique_ptr<RooFitResult> r_ml_unw43{p2.fitTo(*data3, Save(), PrintLevel(-1))};


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

   std::unique_ptr<RooAbsData> 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

   std::unique_ptr<RooAbsReal> chi2{

      p2.createChi2(static_cast<RooDataHist &>(*binnedData), DataError(RooAbsData::SumW2))};

   RooMinimizer m(*chi2);

   m.migrad();

   m.hesse();


   // Plot chi^2 fit result on frame as well

   std::unique_ptr<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();

}

RooDataHist.h

RooDataSet.h

RooFitResult.h

RooFormulaVar.h

RooGaussian.h

RooGenericPdf.h

RooMinimizer.h

RooPlot.h

RooPolynomial.h

RooRealVar.h

kRed
@ kRed
Definition Rtypes.h:66

kDashed
@ kDashed
Definition TAttLine.h:48

TAxis.h

TCanvas.h

w
winID w
Definition TGWin32VirtualGLProxy.cxx:39

data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Definition TGWin32VirtualXProxy.cxx:104

gPad
#define gPad
Definition TVirtualPad.h:305

RooAbsData::SumW2
@ SumW2
Definition RooAbsData.h:108

RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22

RooDataHist
Container class to hold N-dimensional binned data.
Definition RooDataHist.h:39

RooDataSet
Container class to hold unbinned data.
Definition RooDataSet.h:57

RooFormulaVar
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Definition RooFormulaVar.h:30

RooGenericPdf
Implementation of a probability density function that takes a RooArgList of servers and a C++ express...
Definition RooGenericPdf.h:25

RooMinimizer
Wrapper class around ROOT::Fit:Fitter that provides a seamless interface between the minimizer functi...
Definition RooMinimizer.h:43

RooPlot
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43

RooPlot::frame
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

RooPlot::GetYaxis
TAxis * GetYaxis() const
Definition RooPlot.cxx:1276

RooPlot::Draw
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:649

RooPolynomial
RooPolynomial implements a polynomial p.d.f of the form.
Definition RooPolynomial.h:25

RooRealVar
Variable that can be changed from the outside.
Definition RooRealVar.h:37

TAttAxis::SetTitleOffset
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298

TCanvas
The Canvas class.
Definition TCanvas.h:23

TMarker::Print
void Print(Option_t *option="") const override
Dump this marker with its attributes.
Definition TMarker.cxx:335

RooFit::Save
RooCmdArg Save(bool flag=true)
Definition RooGlobalFunc.cxx:566

RooFit::SumW2Error
RooCmdArg SumW2Error(bool flag)
Definition RooGlobalFunc.cxx:650

RooFit::DataError
RooCmdArg DataError(Int_t)
Definition RooGlobalFunc.cxx:387

RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition RooGlobalFunc.cxx:574

RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition RooGlobalFunc.cxx:173

RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition RooGlobalFunc.cxx:177

x
Double_t x[n]
Definition legend1.C:17

RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26

rf403_weightedevts
Definition rf403_weightedevts.py:1

xmlio::Title
const char * Title
Definition TXMLSetup.cxx:68

m
TMarker m
Definition textangle.C:8