#include "RooRealVar.h" #include "RooGaussian.h" #include "RooUniform.h" #include "RooDataSet.h" #include "RooDataHist.h" #include "RooHistFunc.h" #include "RooRealSumPdf.h" #include "RooParamHistFunc.h" #include "RooHistConstraint.h" #include "RooProdPdf.h" #include "RooPlot.h" #include "TCanvas.h" #include "TPaveText.h" #include #include using namespace RooFit; void rf709_BarlowBeeston() { // First, construct a likelihood model with a Gaussian signal on top of a uniform background RooRealVar x("x", "x", -20, 20); x.setBins(25); RooRealVar meanG("meanG", "meanG", 1, -10, 10); RooRealVar sigG("sigG", "sigG", 1.5, -10, 10); RooGaussian g("g", "Gauss", x, meanG, sigG); RooUniform u("u", "Uniform", x); // Generate the data to be fitted std::unique_ptr sigData(g.generate(x, 50)); std::unique_ptr bkgData(u.generate(x, 1000)); RooDataSet sumData("sumData", "Gauss + Uniform", x); sumData.append(*sigData); sumData.append(*bkgData); // Make histogram templates for signal and background. // Let's take a signal distribution with low statistics and a more accurate // background distribution. // Normally, these come from Monte Carlo simulations, but we will just generate them. std::unique_ptr dh_sig( g.generateBinned(x, 50) ); std::unique_ptr dh_bkg( u.generateBinned(x, 10000) ); // ***** Case 0 - 'Rigid templates' ***** // Construct histogram shapes for signal and background RooHistFunc p_h_sig("p_h_sig","p_h_sig",x,*dh_sig); RooHistFunc p_h_bkg("p_h_bkg","p_h_bkg",x,*dh_bkg); // ***** Case 1 - 'Barlow Beeston' ***** // Construct parameterized histogram shapes for signal and background RooParamHistFunc p_ph_sig1("p_ph_sig","p_ph_sig",*dh_sig); RooParamHistFunc p_ph_bkg1("p_ph_bkg","p_ph_bkg",*dh_bkg); RooRealVar Asig1("Asig","Asig",1,0.01,5000); RooRealVar Abkg1("Abkg","Abkg",1,0.01,5000); // Construct the sum of these RooRealSumPdf model_tmp("sp_ph", "sp_ph", RooArgList(p_ph_sig1,p_ph_bkg1), RooArgList(Asig1,Abkg1), true); // Construct the subsidiary poisson measurements constraining the histogram parameters // These ensure that the bin contents of the histograms are only allowed to vary within // the statistical uncertainty of the Monte Carlo. RooHistConstraint hc_sig("hc_sig","hc_sig",p_ph_sig1); RooHistConstraint hc_bkg("hc_bkg","hc_bkg",p_ph_bkg1); // Construct the joint model with template PDFs and constraints RooProdPdf model1("model1","model1",RooArgSet(hc_sig,hc_bkg),Conditional(model_tmp,x)); auto result1 = model1.fitTo(sumData, PrintLevel(0), Save()); TCanvas* can = new TCanvas("can", "", 1500, 600); can->cd(); auto frame = x.frame(Title("")); sumData.plotOn(frame); model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1)); // Plot data again to show it on top of error bands: sumData.plotOn(frame); model1.plotOn(frame, LineColor(kBlue)); model1.plotOn(frame, Components(p_ph_sig1), LineColor(kAzure)); model1.plotOn(frame, Components(p_ph_bkg1), LineColor(kRed)); model1.paramOn(frame, Parameters(RooArgSet(Asig1, Abkg1))); sigData->plotOn(frame, MarkerColor(kBlue)); frame->Draw(); can->SaveAs("test.png"); x.setRange("window",-20,20); RooAbsReal * integral_sig1 = p_ph_sig1.createIntegral(RooArgSet(x), RooArgSet(x), "window"); Double_t Asig1val = Asig1.getVal(); Double_t Asig1val_err = Asig1.getError(); Double_t integval_sig = integral_sig1->getVal(); Double_t dintegval_sig = integral_sig1->getPropagatedError(*result1, RooArgSet(x)); cout << "Asig1 : " << Asig1val << " +/- " << Asig1val_err << endl; cout << "Signal Integral in window: " << integval_sig << " +/- " << dintegval_sig << endl; x.setRange("window1",-5,0); RooAbsReal * integral_sig2 = p_ph_sig1.createIntegral(RooArgSet(x), RooArgSet(x), "window1"); integval_sig = integral_sig2->getVal(); dintegval_sig = integral_sig2->getPropagatedError(*result1, RooArgSet(x)); cout << "Signal Integral in window1: " << integval_sig << " +/- " << dintegval_sig << endl; }