////////////////////////////////////////////////////////////////////////// // // '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 "RooRealIntegral.h" #include "RooFormulaVar.h" #include "RooGenericPdf.h" #include "RooPolynomial.h" #include "RooChi2Var.h" #include "RooMinuit.h" #include "TCanvas.h" #include "RooPlot.h" #include "RooFitResult.h" #include "RooRandom.h" using namespace RooFit ; void pullsgaussian() { // -- load libraries gSystem->Load("libPhysics.so"); gSystem->Load("libRooFit.so"); gStyle->SetOptFit(1); gStyle->SetOptStat(0); gStyle->SetOptTitle(0); RooRandom::randomGenerator()->SetSeed(0); // 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",-15,15) ; x.setBins(40) ; // Construction quadratic polynomial pdf for fitting RooRealVar a0("a0","a0",1,0,100) ; RooRealVar a1("a1","a1",0,-1,1) ; RooGaussian g("g","g",x,a1,a0); // Construct formula to calculate (fake) weight for events RooRealVar wtest("wtest","wtest",-2,4); RooRealVar w2test("w2test","w2test",0,9); TH1D* hist0 = new TH1D("hist0","hist0",100,-6,6); TH1D* hist1 = new TH1D("hist1","hist1",100,-6,6); TH1D* hist2 = new TH1D("hist2","hist2",100,-6,6); TH1D* hist3 = new TH1D("hist3","hist3",100,-6,6); // ------------------------------------------------------------------------------- for (int k=0; k<1000; ++k) { a0.setVal(1.); a1.setVal(0.); a0.setConstant(false); a1.setConstant(false); // Sample 1000 events from pdf RooDataSet* data = g.generate(x,1000) ; RooDataSet* ds = new RooDataSet("ds","ds",RooArgSet(wtest,w2test)); // Fill dataset by hand for (Int_t i=0 ; inumEntries() ; i++) { double myrand = 6.*RooRandom::uniform()-2.; wtest.setVal(myrand); w2test.setVal(myrand*myrand); ds->add(RooArgSet(wtest,w2test)) ; } data->merge(ds); // Instruct dataset d in interpret w as event weight rather than as observable data->setWeightVar(wtest) ; data->Print("v"); // NOTE: Maximum likelihood fit to weighted data does in general // NOT result in correct error estimates, unless individual // event weights represent Poisson statistics themselves. // In general, parameter error reflect precision of SumOfWeights // events rather than NumEvents events. See comparisons below g.fitTo(*data,SumW2Error(kFALSE),Save(kTRUE),Minos(kFALSE)) ; hist0->Fill( (a0.getVal()-1.)/a0.getError() ); hist1->Fill( (a1.getVal())/a1.getError() ); g.fitTo(*data,SumW2Error(kTRUE),Save(kTRUE),Minos(kFALSE)) ; hist2->Fill( (a0.getVal()-1.)/a0.getError() ); hist3->Fill( (a1.getVal())/a1.getError() ); delete data; delete ds; } // ------------------------------------------------------------------------------- TCanvas* c1 = new TCanvas("weightfix","weightfix",1200,800) ; c1->Divide(2,2); c1->cd(1); hist0->Fit("gaus"); c1->cd(2); hist1->Fit("gaus"); c1->cd(3); hist2->Fit("gaus"); c1->cd(4); hist3->Fit("gaus"); c1->Draw(); c1->Print("gaussfit.eps"); }