//#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 "RooMCStudy.h" #include "/home/chanchal/Documents/chanchal phd/analysis_new/bkg_study_mc15rib/WS_sample/myRooJohnsonSU.cpp" #include "/home/chanchal/Documents/chanchal phd/analysis_new/bkg_study_mc15rib/WS_sample/myRooPolBG.cpp" using namespace RooFit ; using namespace std; void toy_data_new() { // gROOT->Reset(); RooRealVar mD0pi("mD0pi","", 2.0044, 2.022,""); RooRealVar mD0("mD0","", 1.8, 1.95,""); /*-----------------------S I G N A L----------------------*/ // // // RS Signal //MD0pi: johnson + double gaussian RooRealVar f_J_sig("f_{Jsig}", "rel_frac_j", 0.986614); RooRealVar mu_J_sig("#mu_{Jsig}","mean_johnson", 2.01025); RooRealVar sigma_J_sig("#sigma_{Jsig}", "sigma_johnson", 0.000270142); RooRealVar gamma_J_sig("#gamma_{Jsig}", "gamma", -0.0648578); RooRealVar delta_J_sig("#delta_{Jsig}", "delta", 1.06486); 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.272020); RooRealVar mu_G1_sig("#mu_{G1sig}", "mean_gauss1", 2.00924); RooRealVar sigma_G1_sig("#sigma_{G1sig}", "sigma_gauss1", 0.00460804); RooRealVar sigma_G2_sig("#sigma_{G2sig", "sigma_gauss2", 0.00135548); 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: johnson + double gaussian RooRealVar f_J_sig_d0("f_{J_sig_d0}", "relative_frac_j", 0.404385); RooRealVar mu_J_sig_d0("#mu_{J_sig_d0}","mean_johnson", 1.86006); RooRealVar sigma_J_sig_d0("#sigma_{J_sig_d0}", "sigma_johnson", 0.0174644); RooRealVar gamma_J_sig_d0("#gamma_{J_sig_d0}", "gamma", 0.370282); RooRealVar delta_J_sig_d0("#delta_{J_sig_d0}", "delta", 1.17478); myRooJohnsonSU johnson_sig_d0("johnson_sig_d0", "johnson_sig_d0", mD0, mu_J_sig_d0, sigma_J_sig_d0, delta_J_sig_d0, gamma_J_sig_d0 ); RooRealVar f_G1_sig_d0("f_{G1}_sig_d0", "relative_frac_g", 0.522563); RooRealVar mu_G1_sig_d0("#mu_{G1}_sig_d0", "mean_gauss1", 1.86461); RooRealVar sigma_G1_sig_d0("#sigma_{G1}_sig_d0", "sigma_gauss1", 0.0115506); RooRealVar sigma_G2_sig_d0("#sigma_{G2}_sig_d0", "sigma_gauss2", 0.00561775); // check RooGaussian gauss1_sig_d0("gauss1_sig_d0", "1st gaussian PDF", mD0, mu_G1_sig_d0, sigma_G1_sig_d0); RooGaussian gauss2_sig_d0("gauss2_sig_d0", "2nd gaussian PDF", mD0, mu_G1_sig_d0, sigma_G2_sig_d0); RooAddPdf doublegaussian_sig_d0("doublegaussian_sig_d0", "doublegaussian_sig_d0", gauss1_sig_d0, gauss2_sig_d0, f_G1_sig_d0); RooAddPdf sig_mD0("sig_mD0","johnson + double gaussian", johnson_sig_d0, doublegaussian_sig_d0, f_J_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); RooRealVar beta("#beta", "beta", 283.054, 0, 5500); 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 RooRealVar f1_d0SMisid("f_{G}_d0SMisid", "", 0.143521); RooRealVar mu_G1_d0SMisid("#mu_{G1}_d0SMisid", "mean_gauss1", 1.97570); RooRealVar sigma_G1_d0SMisid("#sigma_{G1}_d0SMisid", "sigma_gauss1", 0.147941); RooRealVar sigma_G2_d0SMisid("#sigma_{G2}_d0SMisid", "sigma_gauss2", 0.0245542); 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 RooProdPdf SMisid_peak("SMisid_peak", "SMisid_peak", SMisid_md0, sig_mD0pi); //as the shape of m(d0pi) is same as signal /*-----------------------Doubly MisID: piKpi0-------------------------*/ RooRealVar c1_d0DMisid("c1_DMisid", "c1", 0.0789491); RooChebychev DMisid_md0("DMisid_md0","Chebshev Polynomial", mD0, RooArgList(c1_d0DMisid)); //Doubly MisID RooProdPdf DMisid_peak("DMisid_peak", "DMisid_peak", DMisid_md0, sig_mD0pi); //as the shape of m(d0pi) is same as signal /*-----------------------M U L T I B O D Y-------------------------*/ //mD0pi: johnson + gauss RooRealVar mu_J_mult("#mu_{J}_mult","mean_johnson", 2.00941); RooRealVar sigma_J_mult("#sigma_{J}_mult", "sigma_johnson", 0.00186565); RooRealVar gamma_J_mult("#gamma_{J}_mult", "gamma", 1.27576); RooRealVar delta_J_mult("#delta_{J}_mult", "delta_mult", -0.884165); myRooJohnsonSU johnson("johnson","johnson", mD0pi, mu_J_mult, sigma_J_mult, delta_J_mult, gamma_J_mult); RooRealVar f_J_mult("f_{J}_mult", "", 0.662795); RooRealVar mu_G_mult("#mu_{G}_mult", "mean_gauss1", 2.00835); RooRealVar sigma_G_mult("#sigma_{G}_mult", "sigma_gauss1", 0.00204865); RooGaussian gauss("gauss", "1st gaussian PDF", mD0pi, mu_G_mult, sigma_G_mult); RooAddPdf mult_mD0pi("mult_mD0pi", "johnson + gaussian", johnson, gauss, f_J_mult); //mD0: 2nd order Chebychev RooRealVar c1_d0mult("c1_d0mult", "c1", -0.227318); RooRealVar c2_d0mult("c2_d0mult", "c2", 0.0765964); RooChebychev mult_md0("mult_md0","Chebshev Polynomial", mD0, RooArgList(c1_d0mult, c2_d0mult)); RooProdPdf mult_peak("mult_peak", "mult_peak", mult_md0, mult_mD0pi); // /*------------------------C O M B I N A T O R I A L----------------------------*/ // //mD0: Chebychev RooRealVar slope1("slope1", "Slope2 of Polynomial", -0.2, -1, 1); RooChebychev comb_mD0("comb_mD0","Chebshev Polynomial", mD0, RooArgList(slope1)); //mD0pi RooProdPdf comb_bkg("comb_bkg", "comb_bkg", comb_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.1116e+04, 1e1, 1e5); RooRealVar sig_randompi_yield("N_{sig_rnd}","signal + random-pi yield", 1.15849e+05, 1e1, 1e8); RooRealVar mult_peak_yield("N_{mult_peak}","multibody yield", 3.8124e+04, 1e1, 1e8); RooRealVar SMisid_peak_yield("N_{SMisid_peak}","SMisid yield", 6.393e+03, 1e1, 1e5); RooRealVar DMisid_peak_yield("N_{DMisid_peak}","DMisid yield", 2.3048e+04, 1e1, 1e5); RooRealVar comb_yield("N_{comb}","combinatorial yield", 2.57367e+04, 1e1, 1e9); /*------------------------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, DMisid_peak, comb_bkg), RooArgList(sig_peak_yield, sig_randompi_yield, mult_peak_yield, SMisid_peak_yield, DMisid_peak_yield, comb_yield)); //set initial parameters of the model //random pion (common for all) alpha.setVal(-42.3827); beta.setVal(756.662); // C O M B I N A T O R I A L //mD0: Chebychev slope1.setVal(-0.243772); // Y I E L D S F O R A L L C O M P O N E N T S sig_peak_yield.setVal(11116); sig_randompi_yield.setVal(22737); mult_peak_yield.setVal(50302); SMisid_peak_yield.setVal(6662); DMisid_peak_yield.setVal(15162); comb_yield.setVal(231395); /////////////////////////// RooMCStudy ////////////////////////////////////// RooMCStudy *mcs = new RooMCStudy(model, RooArgSet(mD0, mD0pi), Binned(true), Extended(), Silence(), FitOptions(Save(kTRUE), PrintEvalErrors(0))); mcs->generateAndFit(500, 145072, true); mcs->fitParDataSet().Print("v"); RooDataSet* toySample = (RooDataSet*)mcs->genData(0); model.fitTo(*toySample); // Plot m(D0) projections RooPlot *xframe = mD0.frame(); toySample->plotOn(xframe); xframe->SetXTitle("#it{m}(#it{K}^{+}#pi^{-}#pi^{0}) [GeV/#it{c}^{2}]"); xframe->SetYTitle("Candidates per 1.5 MeV/#it{c}^{2}"); model.plotOn(xframe); 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(-4.); pullplotx->SetMaximum(4.); pullplotx->GetXaxis()->SetLabelSize(0.1); pullplotx->GetXaxis()->SetTitleSize(0.1); pullplotx->GetYaxis()->SetLabelSize(0.07); 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(signal_peak, signal_rnd)),LineStyle(kDashed),LineColor(kRed)); model.plotOn(xframe,Components(RooArgSet(SMisid_peak)),LineStyle(kDashed),LineColor(kGreen)); model.plotOn(xframe,Components(RooArgSet(DMisid_peak)),LineStyle(kDashed),LineColor(kCyan)); model.plotOn(xframe,Components(RooArgSet(mult_peak)),LineStyle(kDashed),LineColor(kBlack)); model.plotOn(xframe,Components(comb_bkg),LineStyle(kDashed),LineColor(kMagenta)); xframe->Draw(); canvas1->cd(2); pullplotx->Draw(); canvas1->Update(); //Plot m(D0pi) projections RooPlot *yframe = mD0pi.frame(); toySample->plotOn(yframe); yframe->SetXTitle("#it{m}(#it{D}^{0}#it{#pi}^{#plus}) [GeV/#it{c}^{2}]"); yframe->SetYTitle("Candidates per 19 MeV/#it{c}^{2}"); model.plotOn(yframe); model.paramOn(yframe,toySample); 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(-4.); pullploty->SetMaximum(4.); pullploty->GetXaxis()->SetLabelSize(0.1); pullploty->GetXaxis()->SetTitleSize(0.1); pullploty->GetYaxis()->SetLabelSize(0.07); 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(signal_peak),LineStyle(kDashed),LineColor(kRed)); model.plotOn(yframe,Components(SMisid_peak),LineStyle(kDashed),LineColor(kGreen)); model.plotOn(yframe,Components(DMisid_peak),LineStyle(kDashed),LineColor(kCyan)); model.plotOn(yframe,Components(mult_peak),LineStyle(kDashed),LineColor(kBlack)); model.plotOn(yframe,Components(comb_bkg),LineStyle(kDashed),LineColor(kMagenta)); model.plotOn(yframe,Components(RooArgSet(signal_rnd)),LineStyle(kDashed),LineColor(kOrange)); yframe->Draw(); canvas2->cd(2); pullploty->Draw(); canvas2->Update(); /////////////////////////////// // Make plot of distribution of generated value of N_{sig_peak} parameter TH1* h_f_gen = mcs->fitParDataSet().createHistogram("N_{sig_peak}", -40) ; RooPlot* frame1 = mcs->plotParam(sig_peak_yield, Bins(100)) ; frame1->SetTitle("Distribution of fitted signal yield") ; RooPlot* frame2 = mcs->plotError(sig_peak_yield, Bins(100),FitGauss(false)) ; frame2->SetTitle("Distribution of fit uncertainty on signal yield") ; RooPlot* frame3 = mcs->plotPull(sig_peak_yield, Bins(100),FitGauss(false)) ; frame3->SetTitle("Distribution of signal yield pull") ; TCanvas* c = new TCanvas("TestPulls","TestPulls",1200,400); c->Divide(3, 2) ; c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ; c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ; c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.4) ; frame3->Draw() ; c->cd(4) ; gPad->SetLeftMargin(0.15) ; h_f_gen->GetYaxis()->SetTitleOffset(1.4) ; h_f_gen->Draw() ; // Make RooMCStudy object available on command line after macro finishes gDirectory->Add(mcs); }