///////////////////////////////////////////////////////////////////////// // // 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #305 // // Multi-dimensional p.d.f.s with conditional p.d.fs in product // // pdf = gauss(x,f(y),sx | y ) * gauss(y,ms,sx) with f(y) = a0 + a1*y // // // 07/2008 - Wouter Verkerke // ///////////////////////////////////////////////////////////////////////// #ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooRealVar.h" #include "RooDataSet.h" #include "RooGaussian.h" #include "RooConstVar.h" #include "RooPolyVar.h" #include "RooProdPdf.h" #include "RooPlot.h" #include "TCanvas.h" #include "TAxis.h" #include "TH1.h" #include "TH1F.h" using namespace RooFit ; void prodPdf_eg() { RooRealVar x_1("x1","x1",3,0,3) ; x_1.setBins(3); TH1F sig_h_1 ("sig_histo", "signal_histogram", 3,0,3); TH1F bkg_h_1("bkg_histo", "background_histogram", 3,0,3); RooRealVar x_2("x2","x2",4,0,4) ; TH1F sig_h_2 ("sig2_histo", "signal_histogram", 4,0,4); TH1F bkg_h_2("bkg2_histo", "background_histogram", 4,0,4); sig_h_1.SetBinContent(1,3); sig_h_1.SetBinContent(2,8); sig_h_1.SetBinContent(3,1); bkg_h_1.SetBinContent(1,3); bkg_h_1.SetBinContent(2,3); bkg_h_1.SetBinContent(3,4); sig_h_2.SetBinContent(1,2); sig_h_2.SetBinContent(2,9); sig_h_2.SetBinContent(3,10); sig_h_2.SetBinContent(4,2); bkg_h_2.SetBinContent(1,9); bkg_h_2.SetBinContent(2,12); bkg_h_2.SetBinContent(3,12); bkg_h_2.SetBinContent(4,11); //build model1 RooRealVar nsig1("nsig1","fitted number of signal events",12, 0, 250) ; RooRealVar nbkg1("nbkg1","fitted number of bkg events",8,0,100); nsig1.setBins(3); nbkg1.setBins(3); RooDataHist sig1_hist("sig_hist","signal histogram ",RooArgList(x_1),&sig_h_1); RooDataHist bkg1_hist("sig_hist","signal histogram ",RooArgList(x_1),&bkg_h_1); RooHistPdf sig1_pdf("sig_pdf", "" , x_1, sig1_hist); RooHistPdf bkg1_pdf("bkg_pdf", "" , x_1, bkg1_hist); RooAddPdf model1("model1","model1",RooArgList(sig1_pdf,bkg1_pdf),RooArgList(nsig1,nbkg1)) ; //build model2 RooConstVar const_a("ca","ca",2.); RooConstVar const_b("cb","cb",5.); RooFormulaVar nsig2("nsig2","@0*@1",RooArgList(nsig1,const_a)); RooFormulaVar nbkg2("nbkg2","@0*@1",RooArgList(nbkg1,const_b)); RooDataHist sig2_hist("sig_hist","signal histogram ",RooArgList(x_2),&sig_h_2); RooDataHist bkg2_hist("sig_hist","signal histogram ",RooArgList(x_2),&bkg_h_2); RooHistPdf sig2_pdf("sig_pdf", "" , x_2, sig2_hist); RooHistPdf bkg2_pdf("bkg_pdf", "" , x_2, bkg2_hist); RooAddPdf model2("model2","model2",RooArgList(sig2_pdf,bkg2_pdf),RooArgList(nsig2,nbkg2)) ; RooProdPdf model("model","model",model1,Conditional(model2,x_2)); // double nExpGen = model.expectedEvents(RooArgSet(x_1,x_2)); // double nEvt = RooRandom::randomGenerator()->Poisson(nExpGen) ; RooDataHist *data = model.generateBinned(RooArgSet(x_1,x_2),500); //RooDataSet *data = model.generate(RooArgSet(x_1,x_2),500); model.fitTo(*data,PrintLevel(-1)); std::cout << "nsig1: " << nsig1.getVal() << std::endl; std::cout << "nsig2: " << nsig2.getVal() << std::endl; std::cout << "nbkg1: " << nbkg1.getVal() << std::endl; std::cout << "nbkg2: " << nbkg2.getVal() << std::endl; // Plot x distribution of data and projection of model on x = Int(dy) model(x,y) RooPlot* x1frame = x_1.frame() ; data->plotOn(x1frame) ; model1.plotOn(x1frame) ; //model.plotOn(x1frame) ; // Plot x distribution of data and projection of model on y = Int(dx) model(x,y) RooPlot* x2frame = x_2.frame() ; data->plotOn(x2frame) ; model2.plotOn(x2frame) ; // model.plotOn(x2frame) ; // Make two-dimensional plot in x vs y TH1* hh_model = model.createHistogram("hh_model",x_1,Binning(50),YVar(x_2,Binning(50))) ; hh_model->SetLineColor(kBlue) ; // Make canvas and draw RooPlots TCanvas *c = new TCanvas("rf305_condcorrprod","rf05_condcorrprod",1200, 400); c->Divide(3); c->cd(1) ; gPad->SetLeftMargin(0.15) ; x1frame->GetYaxis()->SetTitleOffset(1.6) ; x1frame->Draw() ; c->cd(2) ; gPad->SetLeftMargin(0.15) ; x2frame->GetYaxis()->SetTitleOffset(1.6) ; x2frame->Draw() ; c->cd(3) ; gPad->SetLeftMargin(0.20) ; hh_model->GetZaxis()->SetTitleOffset(2.5) ; hh_model->Draw("surf") ; // // C r e a t e c o n d i t i o n a l p d f g x ( x | y ) // // ----------------------------------------------------------- // // Create observables // RooRealVar x("x","x",-5,5) ; // RooRealVar y("y","y",-5,5) ; // // Create function f(y) = a0 + a1*y // RooRealVar a0("a0","a0",-0.5,-5,5) ; // RooRealVar a1("a1","a1",-0.5,-1,1) ; // // RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ; // RooFormulaVar fy("fy","@0/@1",RooArgList(a0,a1)) ; // // Create gaussx(x,f(y),sx) // RooRealVar sigmax("sigma","width of gaussian",0.5) ; // RooGaussian gaussx("gaussx","Gaussian in x with shifting mean in y",x,fy,sigmax) ; // // C r e a t e p d f g y ( y ) // // ----------------------------------------------------------- // // Create gaussy(y,0,5) // RooGaussian gaussy("gaussy","Gaussian in y",y,RooConst(0),RooConst(3)) ; // // C r e a t e p r o d u c t g x ( x | y ) * g y ( y ) // // ------------------------------------------------------- // // Create gaussx(x,sx|y) * gaussy(y) // RooProdPdf model("model","gaussx(x|y)*gaussy(y)",gaussy,Conditional(gaussx,x)) ; // // S a m p l e , f i t a n d p l o t p r o d u c t p d f // // --------------------------------------------------------------- // // Generate 1000 events in x and y from model // RooDataSet *data = model.generate(RooArgSet(x,y),10000) ; // model.fitTo(*data); // // Plot x distribution of data and projection of model on x = Int(dy) model(x,y) // RooPlot* xframe = x.frame() ; // data->plotOn(xframe) ; // model.plotOn(xframe) ; // // Plot x distribution of data and projection of model on y = Int(dx) model(x,y) // RooPlot* yframe = y.frame() ; // data->plotOn(yframe) ; // model.plotOn(yframe) ; // // Make two-dimensional plot in x vs y // TH1* hh_model = model.createHistogram("hh_model",x,Binning(50),YVar(y,Binning(50))) ; // hh_model->SetLineColor(kBlue) ; // // Make canvas and draw RooPlots // TCanvas *c = new TCanvas("rf305_condcorrprod","rf05_condcorrprod",1200, 400); // c->Divide(3); // c->cd(1) ; gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.6) ; xframe->Draw() ; // c->cd(2) ; gPad->SetLeftMargin(0.15) ; yframe->GetYaxis()->SetTitleOffset(1.6) ; yframe->Draw() ; // c->cd(3) ; gPad->SetLeftMargin(0.20) ; hh_model->GetZaxis()->SetTitleOffset(2.5) ; hh_model->Draw("surf") ; }