#ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooRealVar.h" #include "RooDataSet.h" #include "RooGaussian.h" #include "RooPlot.h" #include "RooFitResult.h" #include "RooConstVar.h" #include "RooDataHist.h" #include "RooHistPdf.h" #include "RooFitResult.h" #include "RooArgSet.h" #include "RooAddPdf.h" #include "RooCategory.h" #include "RooSimultaneous.h" #include "RooChebychev.h" #include "RooAbsReal.h" #include "RooMinuit.h" #include "RooProdPdf.h" #include "RooGaussian.h" #include "RooLinearVar.h" #include #include #include #include #include #include #include #include #include #include #include #include using namespace RooFit; using namespace std; void MLeventrate(); // modification of tutorial rf_202 void test1(); //snippet of tutorial rf_202; for comparison and reference int main(){ MLeventrate(); //test1(); return 0; } void MLeventrate(){ // Define observable ----------------- RooRealVar x("x","x",0,10); // Data set 1 --------------------------- RooRealVar mean("mean", "mean of Gaussians", 5); RooRealVar sigma1("sigma1", "sigma of Gaussian 1", 0.5); RooRealVar sigma2("sigma2", "sigma of Gaussian 2", 1.0); RooGaussian gaus1 ("gaus1", "gaus1", x, mean, sigma1); RooGaussian gaus2 ("gaus2", "gaus2", x, mean, sigma2); RooRealVar a0("a0","a0",0.5,0.,1.) ; RooRealVar a1("a1","a1",-0.2,0.,1.) ; RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ; RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ; RooAddPdf sig("sig","Signal",RooArgList(gaus1,gaus1),sig1frac) ; // Data set 2 --------------------------- RooRealVar mean_ds2("mean_ds2", "mean of Gaussians", 5); RooRealVar sigma1_ds2("sigma1_ds2", "sigma of Gaussian 1", 0.5); RooRealVar sigma2_ds2("sigma2_ds2", "sigma of Gaussian 2", 1.0); RooGaussian gaus1_ds2 ("gaus1_ds2", "gaus1", x, mean_ds2, sigma1_ds2); RooGaussian gaus2_ds2 ("gaus2_ds2", "gaus2", x, mean_ds2, sigma2_ds2); RooRealVar a0_ds2("a0_ds2","a0_ds2",0.5,0.,1.) ; RooRealVar a1_ds2("a1_ds2","a1_ds2",-0.2,0.,1.) ; RooChebychev bkg_ds2("bkg_ds2","Background",x,RooArgSet(a0_ds2,a1_ds2)) ; RooRealVar sig1frac_ds2("sig1frac_ds2","fraction of component 1 in signal",0.8,0.,1.) ; RooAddPdf sig_ds2("sig_ds2","Signal",RooArgList(gaus1_ds2,gaus1_ds2),sig1frac_ds2) ; // Define fit parameters of interest --------------- RooRealVar nsig("nsig","number of signal events",500,0.,10000) ; RooRealVar nbkg("nbkg","number of background events",500,0,10000) ; RooRealVar bkgrate("bkgrate","bkgrate", 100, 0, 200); RooRealVar sigrate("sigrate","sigrate", 250, 0, 400); RooRealVar eff1("eff1", "eff1", 0.2); RooRealVar eff2("eff2", "eff2", 0.5); RooLinearVar events_bkg("events_bkg","events_bkg", bkgrate, eff1, RooConst(0.0)); RooLinearVar events_sig("events_sig","events_sig", sigrate, eff2, RooConst(0.0)); RooRealVar eff1_ds2("eff1_ds2", "eff1_ds2", 0.1); RooRealVar eff2_ds2("eff2_ds2", "eff2_ds2", 0.6); RooLinearVar events_bkg_ds2("events_bkg_ds2","events_bkg_ds2", bkgrate, eff1_ds2, RooConst(0.0)); RooLinearVar events_sig_ds2("events_sig_ds2","events_sig_ds2", sigrate, eff2_ds2, RooConst(0.0)); // Create models to fit ------------- RooAddPdf model("model","(g1+g2)+a",RooArgList(bkg,sig),RooArgList(events_bkg,events_sig)) ; RooAddPdf model_ds2("model_ds2","(g1+g2)+a",RooArgList(bkg_ds2,sig_ds2),RooArgList(events_bkg_ds2,events_sig_ds2)) ; sig1frac.setConstant(kTRUE); sig1frac_ds2.setConstant(kTRUE); // Define the two datasets and corresponding models------ RooCategory sample("sample", "sample"); sample.defineType("ds1"); sample.defineType("ds2"); RooSimultaneous simpdf("simpdf","simultaneous sample", sample); simpdf.addPdf(model,"ds1"); simpdf.addPdf(model_ds2, "ds2"); // Generate fake data sets and associate them with corresponding labels---- RooRealVar weight("weight", "weight", 1.); RooArgSet datasetVars(x, sample, weight); RooDataSet combData("combData","combData",datasetVars, WeightVar(weight)); RooDataSet *data1 = model.generate(x); TH1* hdata1 = data1->createHistogram("hdata1",x); sample.setLabel("ds1"); for (int i=1; i<=hdata1->GetNbinsX()+1;++i){ x.setVal(hdata1->GetBinCenter(i)); combData.add(datasetVars,hdata1->GetBinContent(i)); } RooDataSet *data2 = model_ds2.generate(x); TH1* hdata2 = data2->createHistogram("hdata2",x); sample.setLabel("ds2"); for (int i=1; i<=hdata2->GetNbinsX()+1;++i){ x.setVal(hdata2->GetBinCenter(i)); combData.add(datasetVars,hdata2->GetBinContent(i)); } // Define constraints RooGaussian constraint_bkg ("constraint_bkg", "constraint on background", bkgrate, RooConst(80), RooConst(80*0.1)); // Fit! RooFitResult *r = simpdf.fitTo(combData, ExternalConstraints(constraint_bkg), Extended(), Save()); r->Print("V"); // Plot! RooPlot* frame1 = x.frame(); combData.plotOn(frame1, Cut("sample==sample::ds1")); simpdf.plotOn(frame1, Slice(sample,"ds1"), ProjWData(sample, combData)); simpdf.plotOn(frame1, Slice(sample,"ds1"), ProjWData(sample, combData), Components(bkg), LineStyle(kDashed), LineColor(kGreen+2)); simpdf.plotOn(frame1, Slice(sample,"ds1"), ProjWData(sample, combData), Components(sig), LineStyle(kDashed), LineColor(kRed)); RooPlot* frame2 = x.frame(); combData.plotOn(frame2, Cut("sample==sample::ds2")); simpdf.plotOn(frame2, Slice(sample,"ds2"), ProjWData(sample, combData)); simpdf.plotOn(frame2, Slice(sample,"ds2"), ProjWData(sample, combData), Components(bkg_ds2), LineStyle(kDashed), LineColor(kGreen+2)); simpdf.plotOn(frame2, Slice(sample,"ds2"), ProjWData(sample, combData), Components(sig_ds2), LineStyle(kDashed), LineColor(kRed)); TCanvas *c1 = new TCanvas("c1","plot canvas", 800,800); c1->Divide(2,1); c1->cd(1); frame1->Draw(); c1->cd(2); frame2->Draw(); c1->SaveAs("test.png"); } void test1(){ // Declare observable x RooRealVar x("x","x",0,10) ; // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters RooRealVar mean("mean","mean of gaussians",5) ; RooRealVar sigma1("sigma1","width of gaussians",0.5) ; RooRealVar sigma2("sigma2","width of gaussians",1) ; RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ; RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ; // Build Chebychev polynomial p.d.f. RooRealVar a0("a0","a0",0.5,0.,1.) ; RooRealVar a1("a1","a1",-0.2,0.,1.) ; RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ; // Sum the signal components into a composite signal p.d.f. RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ; RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ; ///////////////////// // M E T H O D 1 // ///////////////////// // C o n s t r u c t e x t e n d e d c o m p o s i t e m o d e l // ------------------------------------------------------------------- // Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg RooRealVar nsig("nsig","number of signal events",500,0.,10000) ; RooRealVar nbkg("nbkg","number of background events",500,0,10000) ; RooAddPdf model("model","(g1+g2)+a",RooArgList(bkg,sig),RooArgList(nbkg,nsig)) ; // S a m p l e , f i t a n d p l o t e x t e n d e d m o d e l // --------------------------------------------------------------------- // Generate a data sample of expected number events in x from model // = model.expectedEvents() = nsig+nbkg RooDataSet *data = model.generate(x) ; // Fit model to data, extended ML term automatically included model.fitTo(*data) ; // Plot data and PDF overlaid, use expected number of events for p.d.f projection normalization // rather than observed number of events (==data->numEntries()) RooPlot* xframe = x.frame(Title("extended ML fit example")) ; data->plotOn(xframe) ; model.plotOn(xframe,Normalization(1.0,RooAbsReal::RelativeExpected)) ; // Overlay the background component of model with a dashed line model.plotOn(xframe,Components(bkg),LineStyle(kDashed),Normalization(1.0,RooAbsReal::RelativeExpected)) ; // Overlay the background+sig2 components of model with a dotted line model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted),Normalization(1.0,RooAbsReal::RelativeExpected)) ; // Print structure of composite p.d.f. model.Print("t") ; TCanvas *c2 = new TCanvas(); c2->cd(); xframe->Draw(); c2->SaveAs("test1.png"); }