#include "RooRealVar.h" #include "RooDataSet.h" #include "RooDataHist.h" #include "RooGaussian.h" #include "TCanvas.h" #include "RooPlot.h" #include "TTree.h" #include "TH1D.h" #include "TRandom.h" using namespace RooFit; using namespace std; void BmFit_DataFile_exp_doublgauss() { // --------------------------------------------------- // I m p o r t i n g R O O T h i s t o g r a m s // =================================================== // I m p o r t T H 1 i n t o a R o o D a t a H i s t // --------------------------------------------------------- TChain* chain_data = new TChain("tree"); chain_data->Add("/Users/sadhanaverma/Documents/PhD/Roofit/Roofit1/sel_Combine_2016_Mini_Presel_data_cut_bdt-Correction-s04.root"); Int_t nevt = (int)chain_data->GetEntries(); std::cout<<"Number of total events: "<reduce(bm && Q_2 && triggers && AR && B_dt && M_eta && pt && K_eta && trk && dr); RooPlot* frame = Bmass->frame(/*Title(""),*/Bins(100)); data->plotOn(frame, RooFit::Name("t_data")); // Fit with Exponential and double_Gaussian p.d.f to the data Bmass->setRange("signalRange",5.15,5.55); RooRealVar* mean = new RooRealVar("mean","mean",5.36,5.15,5.55); RooRealVar* sigma1 = new RooRealVar("sigma1","sigma1",0.0248,0.,.1); RooRealVar* sigma2 = new RooRealVar("sigma2","sigma2",0.0543,0.,.1); RooRealVar* excon = new RooRealVar("excon","excon",-5.59,3.5); RooRealVar* ns = new RooRealVar("ns","number signal events",45880,0.,1000000000); RooRealVar* nb = new RooRealVar("nb","number background events",2113,0.,10000000); RooGaussian* gauss1 = new RooGaussian("gauss1","gauss1",*Bmass,*mean,*sigma1); RooGaussian* gauss2 = new RooGaussian("gauss2","gauss2",*Bmass,*mean,*sigma2); RooExponential* expo = new RooExponential("expo","expo",*Bmass,*excon); RooRealVar* frac1 = new RooRealVar("frac1","frac1",0.52985,0.0,1.0); RooRealVar* frac2 = new RooRealVar("frac2","frac2",0.3,0.0,1.0); RooAddPdf* siggaugau = new RooAddPdf("siggaugau","siggaugau",RooArgList(*gauss1,*gauss2),RooArgList(*frac1)); RooExtendPdf* esig = new RooExtendPdf("esig","extended signal p.d.f",*siggaugau,*ns,"signalRange") ; //https://root.cern.ch/doc/master/classRooExtendPdf.html RooExtendPdf* ebkg = new RooExtendPdf("ebkg","extended background p.d.f",*expo,*nb,"signalRange") ; //https://root.cern.ch/doc/master/classRooExtendPdf.html RooAddPdf* gaugauexpo = new RooAddPdf("gaugauexpo","gaugauexpo",RooArgList(*esig,*ebkg)); RooFitResult* fitRes = gaugauexpo->fitTo(*data,Save(),Extended(1)); // print the results of the fit fitRes->Print("v"); gaugauexpo->plotOn(frame, RooFit::Name("t_gaugauexpo")) ; gaugauexpo->paramOn(frame,Layout(1,0.64,0.94)); gaugauexpo->plotOn(frame, LineColor(kBlue), LineWidth(2)); // You have to get the chi-square of the fit from MassFrame - the 7 tells RooFit to take into account the 7 fit parameters when calculating the number of degrees of freedom //cout<<" Fit chi square/dof = "<< frame->chiSquare(9)<chiSquare(7); cout<<"Fit chi_square/dof = "<< chisquare << endl; gaugauexpo->plotOn(frame, RooFit::Components(*gauss1), RooFit::LineColor(3), /*RooFit::Name("signal-1")*/ LineWidth(4), LineStyle(2)); gaugauexpo->plotOn(frame, RooFit::Components(*gauss2), RooFit::LineColor(2), /*RooFit::Name("signal-2")*/ LineWidth(4), LineStyle(2)); gaugauexpo->plotOn(frame, RooFit::Components(*expo), RooFit::LineColor(6), RooFit::Name("t_Background"), LineWidth(4), LineStyle(2)); TLegend *leg = new TLegend(0.7,0.7,0.9,0.9); //leg->AddEntry(frame->findObject("signal-1"),"B^{+}#rightarrow J/#psi K^{+}","l"); //leg->AddEntry(frame->findObject("signal-2"),"Combinatorial","l"); leg->AddEntry(frame->findObject("t_data"),"Data","EP"); leg->AddEntry(frame->findObject("t_gaugauexpo"),"Fit","L"); leg->AddEntry(frame->findObject("t_Background"),"Background","L"); // Put the Chi2 label on the plot TPaveLabel *t1=new TPaveLabel(0.38,0.82,0.63,0.7, Form("#chi^{2}/ndf = %f",chisquare),"brNDC"); t1->SetFillColor(0); t1->SetTextSize(0.4); frame->addObject(t1); //----- //TLegend *leg = new TLegend(0.7,0.7,0.9,0.9); //leg->AddEntry(frame->findObject("gauss1"),"B^{+}#rightarrow J/#psi K^{+}","l"); //leg->AddEntry(frame->findObject("gauss2"),"Combinatorial","l"); //leg->AddEntry(frame->findObject("expo"),"B^{+}#rightarrow J/#psi X","l"); //legend->AddEntry(cbgauexpo, "first gaussian pdf", "l"); //legend->AddEntry(frame->FindObject("gauss2"), "second gaussian pdf", "l"); //legend->AddEntry(frame->FindObject("expo"), "exponential pdf", "l"); //legend->Draw(); //frame->Draw(); //RooPlot* pullframe = Bmass->frame(RooFit::Title("Mass pull")); //RooHist* hpull1 = frame->pullHist(); //pullframe->addPlotable(hpull1,"P0") ; //pullframe->SetMinimum(-3) ; //pullframe->SetMaximum(+3) ; //pullframe->SetYTitle("Pull"); //pullframe->SetMarkerStyle(20); //pullframe->SetNdivisions(10); //pullframe->GetXaxis()->CenterTitle(); //pullframe->GetXaxis()->SetTitleOffset(0.95); //pullframe->GetXaxis()->SetTitleSize(0.12); //pullframe->GetXaxis()->SetLabelSize(0.12); //pullframe->GetYaxis()->CenterTitle(); //pullframe->GetYaxis()->SetTitleOffset(0.24); //pullframe->GetYaxis()->SetTitleSize(0.18); //pullframe->GetYaxis()->SetLabelSize(0.13); //TCanvas *c = new TCanvas("c", "c",0,0,600,600); //TPad *pad1 = new TPad("pad1","pad1",0,0.33,1,1); //TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.33); //pad1->SetBottomMargin(0.11); //pad1->SetBorderMode(0); //pad2->SetTopMargin(0.15); //pad2->SetBottomMargin(0.25); //pad2->SetBorderMode(0); //pad1->Draw(); //pad2->Draw(); //pad1->cd(); //gStyle->SetOptTitle(0); //c->SetFillColor(0); //c->SetBorderSize(2); //c->SetLeftMargin(0.1422222); //c->SetRightMargin(0.04444445); frame->SetStats(0); frame->Draw(); leg->Draw("same"); //TLatex* cms1 = new TLatex(5.0, 49050, " "); //cms1->SetNDC(true); //cms1->SetTextColor(12); //cms1->SetTextFont(42); //cms1->SetTextSize(0.055); //cms1-> DrawLatex(0.15,0.95,"Fit with Exponential and double_Gaussian pdf"); //pad2->cd(); //pullframe->SetStats(0); //pullframe->Draw(); //c->cd(); }