//#ifndef __CINT__ #include "RooGlobalFunc.h" //#endif #include "RooRealVar.h" #include "RooArgList.h" #include "RooFormulaVar.h" #include "RooDataSet.h" #include "RooGaussian.h" #include "RooChebychev.h" #include "RooExponential.h" #include "RooBifurGauss.h" #include "RooAddModel.h" #include "RooProdPdf.h" #include "TCanvas.h" #include "RooPlot.h" #include "RooHist.h" #include "RooCBShape.h" #include "RooPolynomial.h" #include "RooBinning.h" #include "RooPolyVar.h" #include "TH1.h" #include "TH2.h" #include "RooAddPdf.h" #include "RooProdPdf.h" #include "RooFitResult.h" #include "RooGenericPdf.h" #include "RooLandau.h" #include "TChain.h" #include "TCanvas.h" #include "TAxis.h" #include "RooPlot.h" #include "RooCategory.h" #include "RooSuperCategory.h" #include "RooSimultaneous.h" #include "RooNLLVar.h" #include "TFile.h" #include "TPaveLabel.h" #include "RooBinning.h" #include "RooDataHist.h" #include "RooJohnson.h" #include "RooHypatia2.h" #include "RooCrystalBall.h" #include "RooNovosibirsk.h" #include "/home/chanchal/Documents/chanchal phd/analysis_new/Dstatg_CFT_GenMC_feb_2023/fit/WS_final_fit/myRooJohnsonSU.cpp" #include "/home/chanchal/Documents/chanchal phd/analysis_new/Dstatg_CFT_GenMC_feb_2023/fit/WS_final_fit/myRooPolBG.cpp" using namespace RooFit ; using namespace std; #include "TLegend.h" #include "TLatex.h" #include "TROOT.h" #include "RooWorkspace.h" #include "TStyle.h" #include "TGaxis.h" const int belle2_font = 133; const double belle2_tsize = 26; void fit_2d_WS_check_unbinned() { TFile *file = new TFile("WS_sample_mc15ri_2dfit_variable_weight.root"); // // // // unbinned dataset /// /// /// // // TTree *tree = (TTree*)file->Get("DstD0PiKPiPi0WS"); cout << "Entries = " << tree->GetEntries() << endl; RooRealVar mD0("D0_M","", 1.8, 1.95, "GeV/#it{c}^{2}"); RooRealVar mD0pi("D0pi_M","", 2.0044, 2.022, "GeV/#it{c}^{2}"); RooRealVar weight("Weight","", 0, 5, ""); RooDataSet * dataxy = new RooDataSet("dataxy","dataxy", tree, RooArgSet(mD0, mD0pi, weight)); /*--------------------------------RS S I G N A L----------------------*/ ///////////////////////////////////////// ////////////////////////////////////////////dddddddddd // // // MD0pi: johnson + double gaussian RooRealVar f_J_sig("f_{Jsig}", "rel_frac_j", 0.985157); RooRealVar mu_J_sig("#mu_{Jsig}","mean_johnson", 2.01025); RooRealVar sigma_J_sig("#sigma_{Jsig}", "sigma_johnson", 0.000265024); RooRealVar gamma_J_sig("#gamma_{Jsig}", "gamma", -0.0599776); RooRealVar delta_J_sig("#delta_{Jsig}", "delta", 1.07406); myRooJohnsonSU johnson_sig("johnson_sig","johnson_sig", mD0pi, mu_J_sig, sigma_J_sig, delta_J_sig, gamma_J_sig); RooRealVar f_G1_sig("f_{G1sig}", "frac_gauss1", 0.238573); RooRealVar mu_G1_sig("#mu_{G1sig}", "mean_gauss1", 2.00928); RooRealVar sigma_G1_sig("#sigma_{G1sig}", "sigma_gauss1", 0.00437507); RooRealVar sigma_G2_sig("#sigma_{G2sig}", "sigma_gauss2", 0.00135064); RooGaussian gauss1_sig("gauss1_sig", "1st gaussian PDF", mD0pi, mu_G1_sig, sigma_G1_sig); RooGaussian gauss2_sig("gauss2_sig", "2nd gaussian PDF", mD0pi, mu_G1_sig, sigma_G2_sig); RooAddPdf doublegaussian_sig("doublegaussian_sig", "doublegaussian_sig", gauss1_sig, gauss2_sig, f_G1_sig); RooAddPdf sig_mD0pi("sig_mD0pi","johnson + double gaussian", johnson_sig, doublegaussian_sig, f_J_sig); // //////////////////////////////////////////// // // // /// // MD0: Novosibirsk + CB from binned method // // // // Novosibirsk + CB (Dz_isSignal==1) RooRealVar f_novo_sig_d0("f_{novo_sig_d0}", "", 0.204795); RooRealVar mean_N1_sig_d0("#mu_{N1_sig_d0}", "peak", 1.86569); RooRealVar sigma_N1_sig_d0("#sigma_{N1_sig_d0}", "width", 0.00653488); RooRealVar lambda_N1_sig_d0("#lambda_{N1_sig_d0}", "tail", -0.360276); RooNovosibirsk novo_sig_d0("novo_sig_d0", "1st Novo", mD0, mean_N1_sig_d0, sigma_N1_sig_d0, lambda_N1_sig_d0); RooRealVar x0_sig_d0("x0_sig_d0", "x0", 1.86310); RooRealVar sigmaL_sig_d0("sigmaL_sig_d0", "sigmaL", 0.0102188); RooRealVar sigmaR_sig_d0("sigmaR_sig_d0", "sigmaR", 0.00717711); RooRealVar alphaL_sig_d0("alphaL_sig_d0", "alphaL", 0.938547); RooRealVar alphaR_sig_d0("alphaR_sig_d0", "alphaR", 1.29445); RooRealVar nL_sig_d0("nL_sig_d0", "nL", 5.93616); RooRealVar nR_sig_d0("nR_sig_d0", "nR", 4.19775); RooCrystalBall CB_sig_d0("CB_sig_d0", "CrystalBall1 PDF", mD0, x0_sig_d0, sigmaL_sig_d0, sigmaR_sig_d0, alphaL_sig_d0, nL_sig_d0, alphaR_sig_d0, nR_sig_d0); RooAddPdf sig_mD0("sig_mD0", "Novosibirsk + CB", novo_sig_d0, CB_sig_d0, f_novo_sig_d0); // // // signal peak RooProdPdf signal_peak("signal_peak", "signal_peak", sig_mD0, sig_mD0pi); //random pion (common for all) RooRealVar threshold("threshold", "no data beyond this", 2.00441); RooRealVar alpha("#alpha", "alpha", -30.2406, -60, 60); //-16.0351, -50, 50 RooRealVar beta("#beta", "beta", 283.054, 0, 5500); // -22.5549, -50, 50 myRooPolBG bkg_mD0pi("bkg_mD0pi","bkg_mD0pi", mD0pi, threshold, alpha, beta); //signal & random pion RooProdPdf signal_rnd("signal_rnd", "signal_rnd", sig_mD0, bkg_mD0pi); /*-----------------------Singly MisID: pipipi0-------------------------*/ //mD0: Double Gaussian : binned swap RooRealVar f1_d0SMisid("f_{G}_d0SMisid", "", 0.142881); RooRealVar mu_G1_d0SMisid("#mu_{G1}_d0SMisid", "mean_gauss1", 1.96808); RooRealVar sigma_G1_d0SMisid("#sigma_{G1}_d0SMisid", "sigma_gauss1", 0.225933); RooRealVar sigma_G2_d0SMisid("#sigma_{G2}_d0SMisid", "sigma_gauss2", 0.0216328); RooGaussian gauss1_d0SMisid("gauss1_d0SMisid", "1st gaussian PDF", mD0, mu_G1_d0SMisid, sigma_G1_d0SMisid); RooGaussian gauss2_d0SMisid("gauss2_d0SMisid", "2nd gaussian PDF", mD0, mu_G1_d0SMisid, sigma_G2_d0SMisid); RooAddPdf SMisid_md0("SMisid_md0", "double_gaussian", gauss1_d0SMisid, gauss2_d0SMisid, f1_d0SMisid); //Singly MisID: pipipi0 peaking background RooProdPdf SMisid_peak("SMisid_peak", "SMisid_peak", SMisid_md0, sig_mD0pi); //as the shape of m(d0pi) is same as signa //Singly MisID: pipipi0 & random pion RooProdPdf SMisid_rnd("SMisid_rnd", "SMisid_rnd", SMisid_md0, bkg_mD0pi); /*-----------------------M U L T I B O D Y-------------------------*/ ///////////////////////////////// // // //MD0pi: Johnson: best RooRealVar mu_J("#mu_{J}","mean_johnson", 2.01368); RooRealVar sigma_J ("#sigma_{J}", "sigma_johnson", 0.00990412); RooRealVar gamma_J("#gamma_{J}", "gamma", 1.67362); RooRealVar delta_J("#delta_{J}", "delta", 4.26132); myRooJohnsonSU mult_mD0pi("mult_mD0pi", "Johnson PDF", mD0pi, mu_J, sigma_J, delta_J, gamma_J ); /////////////////////////////// // // //mD0: 2nd order Chebychev : Unbinned RooRealVar c1_d0mult("c1_d0mult", "c1", -0.244443, -1, 1); RooRealVar c2_d0mult("c2_d0mult", "c2", 0.0807256, -1, 1); RooChebychev mult_md0("mult_md0","Chebshev Polynomial", mD0, RooArgList(c1_d0mult, c2_d0mult)); RooProdPdf mult_peak("mult_peak", "mult_peak", mult_md0, mult_mD0pi); //multibody & random pion RooProdPdf mult_rnd("mult_rnd","mult_rnd", mult_md0, bkg_mD0pi); /*------------------------C O M B I N A T O R I A L----------------------------*/ RooProdPdf comb_bkg("comb_bkg", "comb_bkg", mult_md0, bkg_mD0pi); /*------------------------Y I E L D S F O R A L L C O M P O N E N T S----------------------------*/ ///////////////////////////////// // RooRealVar sig_peak_yield("N_{sig_peak}","signal yield", 1.0179e+04, 1e3, 1e5); RooRealVar sig_peak_yield("N_{sig_peak}","signal yield", 1.7768e+04, 1e3, 1e5); RooRealVar sig_randompi_yield("N_{sig_rnd}","signal + random-pi yield", 1.15849e+05, 1e3, 1e8); RooRealVar mult_peak_yield("N_{mult_peak}","multibody yield", 3.8124e+04, 1e2, 1e6); RooRealVar SMisid_peak_yield("N_{SMisid_peak}","SMisid yield", 6.393e+03, 1e2, 1e5); RooRealVar comb_yield("N_{comb}","combinatorial yield", 3.57367e+05, 1e3, 1e8); /*------------------------C O M P L E T E M O D E L----------------------------*/ RooAddPdf model("model","",RooArgList(signal_peak, signal_rnd, mult_peak, SMisid_peak, comb_bkg), RooArgList(sig_peak_yield, sig_randompi_yield, mult_peak_yield, SMisid_peak_yield, comb_yield)); /*----------------------------------- FIT THE DATA-------------------------------------*/ RooFitResult *fitresult = model.fitTo(*dataxy, Save(true), Strategy(2), Extended(true), RooFit::SumW2Error(true)); /*-----------------------------------PLOT THE RESULTS-------------------------------------*/ //Plot m(KK) projections RooPlot *xframe = mD0.frame(); // dataxy->plotOn(xframe); // dataxy->plotOn(xframe, DataError(RooAbsData::SumW2)); dataxy->plotOn(xframe, Name("data")); // xframe->SetXTitle("#it{m}(#it{K}^{+}#pi^{-}#pi^{0}) [GeV/#it{c}^{2}]"); xframe->GetXaxis()->SetLabelSize(0.05); xframe->SetYTitle("Candidates per 0.75 MeV/#it{c}^{2}"); xframe->GetYaxis()->CenterTitle(); model.plotOn(xframe); xframe->GetYaxis()->CenterTitle(); xframe->GetYaxis()->SetLabelSize(0.05); xframe->GetYaxis()->SetTitleSize(0.05); xframe->GetXaxis()->SetLabelOffset(999); RooHist* hpullx = xframe->pullHist(); hpullx->SetFillStyle(1001); hpullx->SetFillColor(1); for(int i=0;iGetN();++i) hpullx->SetPointError(i,0.0,0.0,0.0,0.0); RooPlot* pullplotx = mD0.frame(Title(" ")) ; pullplotx->addPlotable(hpullx,"B"); pullplotx->SetXTitle("#it{M}(#it{K}^{+}#pi^{-}#pi^{0}) [GeV/#it{c}^{2}]"); pullplotx->SetMinimum(-10.); pullplotx->SetMaximum(10.); pullplotx->GetXaxis()->SetLabelSize(0.11); pullplotx->GetXaxis()->CenterTitle(); pullplotx->GetXaxis()->SetTitleSize(0.13); pullplotx->GetYaxis()->SetLabelSize(0.08); TCanvas* canvas1 = new TCanvas("canvas1", "canvas1",700, 700); Double_t xlow, ylow, xup, yup; canvas1->GetPad(0)->GetPadPar(xlow,ylow,xup,yup); canvas1->Divide(1,2); TVirtualPad *upPad1 = canvas1->GetPad(1); upPad1->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup); TVirtualPad *dwPad1 = canvas1->GetPad(2); dwPad1->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow)); canvas1->Update(); canvas1->cd(1); model.plotOn(xframe,Components(RooArgSet(model)),Name("pdf")); model.plotOn(xframe,Components(RooArgSet(signal_peak, signal_rnd)),LineStyle(kDashed),LineColor(kRed), Name("sig")); model.plotOn(xframe,Components(RooArgSet(SMisid_peak)),LineStyle(kDashed),LineColor(kGreen), Name("SMisid")); model.plotOn(xframe,Components(RooArgSet(mult_peak)),LineStyle(kDashed),LineColor(kBlack), Name("multi")); model.plotOn(xframe,Components(comb_bkg),LineStyle(kDashed),LineColor(kMagenta), Name("comb")); xframe->Draw(); auto logo1 = new TLatex(); logo1->SetNDC(); logo1->DrawLatex(0.5133779,0.8148148,"#bf{#it{Belle II}} Simulation"); logo1->DrawLatex(0.5,0.7337963," #it{#scale[0.6]{#int} L dt} = 1 ab^{#minus1}"); // /// D0 Sample TLegend *legend_d0 = new TLegend(0.54,0.32,0.84,0.52); legend_d0->SetBorderSize(0); // no border legend_d0->AddEntry(xframe->findObject("data"),"Data","pe"); legend_d0->AddEntry(xframe->findObject("pdf"),"Total Fit","l"); legend_d0->AddEntry(xframe->findObject("sig"),"D^{0}#rightarrowK^{+}#pi^{-}#pi^{0}","l"); legend_d0->AddEntry(xframe->findObject("SMisid"),"D^{0}#rightarrow#pi^{+}#pi^{-}#pi^{0}","l"); legend_d0->AddEntry(xframe->findObject("multi"),"D^{0}#rightarrow Multibody","l"); legend_d0->AddEntry(xframe->findObject("comb"),"Combinatorial ","l"); legend_d0->Draw(); canvas1->cd(2); pullplotx->Draw(); canvas1->Update(); //Plot m(D0pi) projections RooPlot *yframe = mD0pi.frame(); // dataxy->plotOn(yframe); // dataxy->plotOn(yframe, DataError(RooAbsData::SumW2)); dataxy->plotOn(yframe, Name("data")); // yframe->SetXTitle("#it{m}(#it{D}^{0}#it{#pi}^{#plus}) [GeV/#it{c}^{2}]"); yframe->GetXaxis()->SetLabelSize(0.05); yframe->SetYTitle("Candidates per 0.09 MeV/#it{c}^{2}"); yframe->GetYaxis()->CenterTitle(); yframe->GetYaxis()->SetLabelSize(0.05); yframe->GetYaxis()->SetTitleSize(0.05); yframe->GetXaxis()->SetLabelOffset(999); model.plotOn(yframe); model.paramOn(yframe,dataxy); RooHist* hpully = yframe->pullHist(); hpully->SetFillStyle(1001); hpully->SetFillColor(1); for(int i=0;iGetN();++i) hpully->SetPointError(i,0.0,0.0,0.0,0.0); RooPlot* pullploty = mD0pi.frame(Title(" ")) ; pullploty->addPlotable(hpully,"B"); pullploty->SetXTitle("#it{M}(#it{D}^{0}#it{#pi}^{+}) [GeV/#it{c}^{2}]"); pullploty->SetMinimum(-10.); pullploty->SetMaximum(10.); pullploty->GetXaxis()->SetLabelSize(0.11); pullploty->GetXaxis()->CenterTitle(); pullploty->GetXaxis()->SetTitleSize(0.13); pullploty->GetYaxis()->SetLabelSize(0.08); TCanvas* canvas2 = new TCanvas("canvas2", "canvas2",700, 700); canvas2->GetPad(0)->GetPadPar(xlow,ylow,xup,yup); canvas2->Divide(1,2); TVirtualPad *upPad2 = canvas2->GetPad(1); upPad2->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup); TVirtualPad *dwPad2 = canvas2->GetPad(2); dwPad2->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow)); canvas2->Update(); canvas2->cd(1); model.plotOn(yframe,Components(RooArgSet(model)), Name("pdf")); model.plotOn(yframe,Components(signal_peak),LineStyle(kDashed),LineColor(kRed), Name("sig")); model.plotOn(yframe,Components(SMisid_peak),LineStyle(kDashed),LineColor(kGreen), Name("SMisid")); model.plotOn(yframe,Components(mult_peak),LineStyle(kDashed),LineColor(kBlack), Name("multi")); model.plotOn(yframe,Components(comb_bkg),LineStyle(kDashed),LineColor(kMagenta), Name("comb")); model.plotOn(yframe,Components(RooArgSet(signal_rnd)),LineStyle(kDashed),LineColor(kOrange), Name("rnd")); yframe->Draw(); auto logo2 = new TLatex(); logo2->SetNDC(); logo2->DrawLatex(0.5133779,0.8148148,"#bf{#it{Belle II}} Simulation"); logo2->DrawLatex(0.5,0.7337963," #it{#scale[0.6]{#int} L dt} = 1 ab^{#minus1}"); // // // D0 Sample TLegend *legend_d0pi = new TLegend(0.54,0.32,0.84,0.52); legend_d0pi->SetBorderSize(0); // no border legend_d0pi->AddEntry(yframe->findObject("data"),"Data","pe"); legend_d0pi->AddEntry(yframe->findObject("pdf"),"Total Fit","l"); legend_d0pi->AddEntry(yframe->findObject("sig"),"D^{*+}#rightarrowD^{0}(#rightarrowK^{+}#pi^{-}#pi^{0})#pi_{s}^{+}","l"); legend_d0pi->AddEntry(yframe->findObject("SMisid"),"D^{*+}#rightarrowD^{0}(#rightarrow#pi^{+}#pi^{-}#pi^{0})#pi_{s}^{+}","l"); legend_d0pi->AddEntry(yframe->findObject("multi"),"D^{*+}#rightarrowD^{0}(#rightarrowMultibody)#pi_{s}^{+}","l"); legend_d0pi->AddEntry(yframe->findObject("comb"),"Combinatorial","l"); legend_d0pi->AddEntry(yframe->findObject("rnd"),"Random pion","l"); legend_d0pi->Draw(); canvas2->cd(2); pullploty->Draw(); canvas2->Update(); }