#include "RooGlobalFunc.h"
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooFormulaVar.h"
#include "RooHist.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 "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include <iostream>
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "TFile.h"
#include "RooDataSet.h"
#include "TTree.h"
#include "myRooJohnsonSU.cpp"
#include "myRooPolBG.cpp"


using namespace RooFit ;
using namespace std;

void simFit_KK_final_2()
{

  //File and the tree containing the variable
  //TFile *f = new TFile("/home/sanjeeda/Desktop/thesis_new/analysis/D0_KK/complete_fit_strategy/new_sim_fit/for_2dfit_MC14.root");
  TFile *f = new TFile("for_2dfit_MC14.root");
  TTree *tree = (TTree*)f->Get("distributions");

  //Variables for 2Dfit
  RooRealVar mD0pi("mD0pi","#it{m}(#it{D}^{0}#it{#pi}^{+})", 2.005, 2.025, "GeV/#it{c}^{2}");
  RooRealVar mKK("mKK","#it{m}(#it{K}^{#minus}#it{K}^{+})", 1.75, 2.05, "GeV/#it{c}^{2}");

  //Variable used to separate the total sample into two (should be there in the tree)
  RooRealVar pisoft("pisoft","",-2,2,"");

  //Separated datasets
  RooDataSet *dataxy_d0 = new RooDataSet("data_d0","data",tree,RooArgSet(mKK,mD0pi,pisoft),"pisoft==1");
  RooDataSet *dataxy_d0bar = new RooDataSet("data_d0bar","data",tree,RooArgSet(mKK,mD0pi,pisoft),"pisoft==-1");  

  /*-----------------------S I G N A L----------------------*/
  
  //d0pi
  RooRealVar f_J_sig("f_{Jsig}", "rel_frac_j",0.8937);
  RooRealVar mu_J_sig("#mu_{Jsig}","mean_johnson",2.0102529);
  RooRealVar sigma_J_sig("#sigma_{Jsig}", "sigma_johnson",0.0001675);
  RooRealVar gamma_J_sig("#gamma_{Jsig}", "gamma", -0.0329);
  RooRealVar delta_J_sig("#delta_{Jsig}", "delta", 1.033);
  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.136);
  RooRealVar mu_G1_sig("#mu_{G1sig}", "mean_gauss1",2.010300);
  RooRealVar sigma_G1_sig("#sigma_{G1sig}", "sigma_gauss1",0.00385);
  RooRealVar sigma_G2_sig("#sigma_{G2sig", "sigma_gauss2",0.001224);
  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","sig_mD0pi", johnson_sig, doublegaussian_sig, f_J_sig);
  
  //kk
  RooRealVar f_J_sig_d0("f_{J_sig_d0}", "relative_frac_j",0.9624);
  RooRealVar mu_J_sig_d0("#mu_{J_sig_d0}","mean_johnson",1.864738);
  RooRealVar sigma_J_sig_d0("#sigma_{J_sig_d0}", "sigma_johnson", 0.00612);
  RooRealVar gamma_J_sig_d0("#gamma_{J_sig_d0}", "gamma",0.011);
  RooRealVar delta_J_sig_d0("#delta_{J_sig_d0}", "delta",1.689);
  myRooJohnsonSU johnson_sig_d0("johnson_sig_d0", "johnson_sig_d0", mKK, 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.451);
  RooRealVar mu_G1_sig_d0("#mu_{G1}_sig_d0", "mean_gauss1",1.8506);
  RooRealVar sigma_G1_sig_d0("#sigma_{G1}_sig_d0", "sigma_gauss1",0.0616);
  RooRealVar sigma_G2_sig_d0("#sigma_{G2}_sig_d0", "sigma_gauss2",0.0225);
  RooGaussian gauss1_sig_d0("gauss1_sig_d0", "1st gaussian PDF", mKK, mu_G1_sig_d0, sigma_G1_sig_d0);
  RooGaussian gauss2_sig_d0("gauss2_sig_d0", "2nd gaussian PDF", mKK, 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_mKK("sig_mKK","sig_mKK", johnson_sig_d0, doublegaussian_sig_d0, f_J_sig_d0);
  
  //signal peak
  RooProdPdf signal_peak("signal_peak", "signal_peak", sig_mKK, sig_mD0pi);
  
  //random pion (common for all)
  RooRealVar threshold("threshold", "no data beyond this", 2.00441);
  RooRealVar alpha("#alpha", "alpha", -16, -1e6, 0.);
  RooRealVar beta("#beta", "beta", 0.);//, -1e6, 1e6);
  myRooPolBG bkg_mD0pi("bkg_mD0pi","bkg_mD0pi", mD0pi, threshold, alpha, beta);
  
  //signal & random pion
  RooProdPdf signal_rnd("signal_rnd", "signal_rnd", sig_mKK, bkg_mD0pi);
  
  /*-----------------------M U L T I B O D Y-------------------------*/
  
  //d0pi
  RooRealVar delta_J_mult("#delta_{J}", "delta_mult",-1.2859);
  RooRealVar mu_J_mult("#mu_{J}_mult","mean_johnson",2.008198);
  RooRealVar sigma_J_mult("#sigma_{J}_mult", "sigma_johnson",0.00363);
  RooRealVar gamma_J_mult("#gamma_{J}_mult", "gamma", 1.);
  myRooJohnsonSU mult_mD0pi("mult_mD0pi","mult_mD0pi", mD0pi, mu_J_mult, sigma_J_mult, delta_J_mult, gamma_J_mult);
  
  //kk
  RooRealVar alpha1("a_{1}", "A1", -12.388);
  RooExponential mult_mKK("mult_mKK", "mult_mKK", mKK, alpha1);  

  RooProdPdf mult_peak("mult_peak", "mult_peak", mult_mKK, mult_mD0pi);
  
  //multibody & random pion
  RooProdPdf mult_rnd("mult_rnd","mult_rnd", mult_mKK, bkg_mD0pi);
  
  /*-----------------------K Pi-------------------------*/
  
  //kk
  RooRealVar f1_d0kpi("f_{G}d0kpi", "relative fraction gaussian",0.385);
  RooRealVar mu_G_d0kpi("#mu_{G}d0kpi","mean_gaus",1.94120);                                                          
  RooRealVar sigma_G_d0kpi("#sigma_{G}d0kpi", "sigma_gaus",0.00716);
  RooGaussian gauss_d0kpi("gauss_d0kpi", "gauss component", mKK, mu_G_d0kpi, sigma_G_d0kpi);
  
  RooRealVar gamma_J_d0kpi("#gamma_{J}d0kpi", "gamma",1.302);
  RooRealVar delta_J_d0kpi("#delta_{J}d0kpi", "delta",-0.3516);
  RooRealVar mu_J_d0kpi("#mu_{J}_d0kpi", "mean",1.95397);
  RooRealVar sigma_J_d0kpi("#sigma_{J}_d0kpi", "sigma",0.00760);
  
  myRooJohnsonSU johnson_d0kpi("johnson_d0kpi", "Johnson PDF", mKK, mu_J_d0kpi, sigma_J_d0kpi, delta_J_d0kpi, gamma_J_d0kpi);
  
  RooAddPdf kpi_mKK("kpi_mKK","kpi_mKK",gauss_d0kpi,johnson_d0kpi,f1_d0kpi);
  
  //kpi peaking background
  RooProdPdf kpi_peak("kpi_peak", "kpi_peak", kpi_mKK, sig_mD0pi);//as the shape of m(d0pi) is same as signal
  
  //kpi & random pion
  RooProdPdf kpi_rnd("kpi_rnd", "kpi_rnd", kpi_mKK, bkg_mD0pi);

  /*----------------------- Ds+ -> K+ K- pi+ -------------------------*/
  //d0pi
  RooRealVar slope_ds("slope_ds", "Slope1 of Polynomial",0.790);
  RooChebychev ds_mD0pi("ds_mD0pi","Chebshev Polynomial", mD0pi,RooArgList(slope_ds));

  //kk
  RooRealVar mu_0("#mu_{0}", "",0.000214269);
  RooRealVar constant("constant", "",3.83333);
  RooFormulaVar mu_ds("#mu_{ds}", "@0+@1-@2",RooArgList(mu_0,constant, mD0pi));
  RooRealVar sigma_ds("#sigma_{ds}", "",0.00684);
  RooRealVar gamma_ds("#gamma_{ds}", "",-0.1118);
  RooRealVar delta_ds("#delta_{ds}", "",1.125);      
  myRooJohnsonSU ds_mKK("ds_mKK","johnson",mKK, mu_ds,sigma_ds,delta_ds,gamma_ds);
 
  RooProdPdf ds_peak("ds_peak", "ds_peak", ds_mD0pi, Conditional(ds_mKK, mKK));
  
  /*------------------------C O M B I N A T O R I A L----------------------------*/
  //Chebychev
  RooRealVar slope1("slope1", "Slope1 of Polynomial",-0.2,-1,1);
  RooChebychev comb_mKK("chebpol","Chebshev Polynomial", mKK,RooArgList(slope1));
  RooProdPdf comb_bkg("comb_bkg", "comb_bkg", comb_mKK, 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",3.18353e+04,1e1,1e5);
  RooRealVar sig_randompi_yield("N_{sig_rnd}","signal + random-pi yield",1.31772e+03,1e1,1e5);
  RooRealVar mult_peak_yield("N_{mult_peak}","multibody yield",1.72079e+04,1e1,1e5);
  RooRealVar mult_randompi_yield("N_{mult_rnd}","multibody + random-pi yield",1.13956e+04,1e1,1e5);
  RooRealVar kpi_peak_yield("N_{kpi_peak}","kpi yield",2.23321e+04,1e1,1e5);
  RooRealVar kpi_randompi_yield("N_{kpi_rnd}","kpi + random-pi yield",1.14669e+03,1e1,1e5);
  RooRealVar ds_peak_yield("N_{ds_peak}","ds yield",65,50,1e5);
  RooRealVar comb_yield("N_{comb}","combinatorial yield",2.57367e+04,1e1,1e5);*/

  //Total yields
  RooRealVar sig_peak_yield("N_{sig_peak}","signal yield",36795,1e1,1e5);
  RooRealVar sig_randompi_yield("N_{sig_rnd}","signal + random-pi yield",2294,1e1,1e5);
  RooRealVar mult_peak_yield("N_{mult_peak}","multibody yield",18429,1e1,1e5);
  RooRealVar mult_randompi_yield("N_{mult_rnd}","multibody + random-pi yield",1202,1e1,1e5);
  RooRealVar kpi_peak_yield("N_{kpi_peak}","kpi yield",25380,1e1,1e5);
  RooRealVar kpi_randompi_yield("N_{kpi_rnd}","kpi + random-pi yield",1266,1e1,1e5);
  RooRealVar ds_peak_yield("N_{ds_peak}","ds yield",1097,50,1e5);
  RooRealVar comb_yield("N_{comb}","combinatorial yield",32922,1e1,1e5);

  //Asymmetries
  RooRealVar Araw_sig("A_{sig}","Araw",0.2,-1,1);
  RooRealVar Araw_sig_randompi("A_{sig_randompi}","Araw",0,-1,1);
  RooRealVar Araw_mult("A_{mult}","Araw",0,-1,1);
  RooRealVar Araw_mult_randompi("A_{mult_randompi}","Araw",0,-1,1);
  RooRealVar Araw_kpi("A_{kpi}","Araw",0,-1,1);
  RooRealVar Araw_kpi_randompi("A_{kpi_randompi}","Araw",0,-1,1);
  RooRealVar Araw_ds("A_{Ds}","Araw",0,-1,1);
  RooRealVar Araw_comb("A_{comb}","Araw",0,-1,1);

  //Yields for D0 sample: ND0 = N(1+A)/2
  RooFormulaVar sig_peak_yield_d0("ND0_{sig_peak}","(@0*(1 + @1))/2",RooArgList(sig_peak_yield,Araw_sig));
  RooFormulaVar sig_randompi_yield_d0("ND0_{sig_rnd}","(@0*(1 + @1))/2",RooArgList(sig_randompi_yield,Araw_sig_randompi));
  RooFormulaVar mult_peak_yield_d0("ND0_{mult_peak}","(@0*(1 + @1))/2",RooArgList(mult_peak_yield,Araw_mult));
  RooFormulaVar mult_randompi_yield_d0("ND0_{mult_rnd}","(@0*(1 + @1))/2",RooArgList(mult_randompi_yield,Araw_mult_randompi));
  RooFormulaVar kpi_peak_yield_d0("ND0_{kpi_peak}","(@0*(1 + @1))/2",RooArgList(kpi_peak_yield,Araw_kpi));
  RooFormulaVar kpi_randompi_yield_d0("ND0_{kpi_rnd}","(@0*(1 + @1))/2",RooArgList(kpi_randompi_yield,Araw_kpi_randompi));
  RooFormulaVar ds_peak_yield_d0("ND0_{ds_peak}","(@0*(1 + @1))/2",RooArgList(ds_peak_yield,Araw_ds));
  RooFormulaVar comb_yield_d0("ND0_{comb}","(@0*(1 + @1))/2",RooArgList(comb_yield,Araw_comb));

  //Yields for D0bar sample: ND0bar = N(1-A)/2
  RooFormulaVar sig_peak_yield_d0bar("ND0bar_{sig_peak}","(@0*(1 - @1))/2",RooArgList(sig_peak_yield,Araw_sig));
  RooFormulaVar sig_randompi_yield_d0bar("ND0bar_{sig_rnd}","(@0*(1 - @1))/2",RooArgList(sig_randompi_yield,Araw_sig_randompi));
  RooFormulaVar mult_peak_yield_d0bar("ND0bar_{mult_peak}","(@0*(1 - @1))/2",RooArgList(mult_peak_yield,Araw_mult));
  RooFormulaVar mult_randompi_yield_d0bar("ND0bar_{mult_rnd}","(@0*(1 - @1))/2",RooArgList(mult_randompi_yield,Araw_mult_randompi));
  RooFormulaVar kpi_peak_yield_d0bar("ND0bar_{kpi_peak}","(@0*(1 - @1))/2",RooArgList(kpi_peak_yield,Araw_kpi));
  RooFormulaVar kpi_randompi_yield_d0bar("ND0bar_{kpi_rnd}","(@0*(1 - @1))/2",RooArgList(kpi_randompi_yield,Araw_kpi_randompi));
  RooFormulaVar ds_peak_yield_d0bar("ND0bar_{ds_peak}","(@0*(1 - @1))/2",RooArgList(ds_peak_yield,Araw_ds));
  RooFormulaVar comb_yield_d0bar("ND0bar_{comb}","(@0*(1 - @1))/2",RooArgList(comb_yield,Araw_comb));

  /*------------------------C O M P L E T E  M O D E L----------------------------*/

  //Model for D0 sample
  RooAddPdf model_d0("model_d0","",RooArgList(signal_peak,signal_rnd, mult_peak,mult_rnd,kpi_peak, kpi_rnd,ds_peak, comb_bkg), RooArgList(sig_peak_yield_d0,sig_randompi_yield_d0,mult_peak_yield_d0, mult_randompi_yield_d0, kpi_peak_yield_d0,kpi_randompi_yield_d0,ds_peak_yield_d0,comb_yield_d0));

  //Model for D0bar sample
  RooAddPdf model_d0bar("model_d0bar","",RooArgList(signal_peak,signal_rnd, mult_peak,mult_rnd,kpi_peak, kpi_rnd,ds_peak, comb_bkg), RooArgList(sig_peak_yield_d0bar,sig_randompi_yield_d0bar,mult_peak_yield_d0bar, mult_randompi_yield_d0bar, kpi_peak_yield_d0bar,kpi_randompi_yield_d0bar,ds_peak_yield_d0bar,comb_yield_d0bar));

  //Define category to distinguish d0 and d0bar samples events
  RooCategory sample("sample","sample") ;
  sample.defineType("D0") ;
  sample.defineType("D0bar") ;

  //Construct combined dataset
  RooDataSet combData("combData","combined data",RooArgSet(mKK, mD0pi),Index(sample),Import("D0",*dataxy_d0),Import("D0bar",*dataxy_d0bar));
  
  //Construct a simultaneous pdf using category sample as index
  RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
  simPdf.addPdf(model_d0,"D0") ;
  simPdf.addPdf(model_d0bar,"D0bar");
  
  //Perform the simultaneous fit
  simPdf.fitTo(combData,Strategy(2));
  
  /*-----------------------------------PLOT THE RESULTS-------------------------------------*/

  // Plot "D0" slice of simultaneous pdf for m(KK). 
  //RooPlot *kk_d0 = mKK.frame(Title("D0"));
  RooPlot *kk_d0 = mKK.frame(Title(""));
  combData.plotOn(kk_d0,Cut("sample==sample::D0"));
  simPdf.plotOn(kk_d0,Slice(sample,"D0"),ProjWData(sample,combData));
  
  kk_d0->SetXTitle("#it{m}(#it{K}^{#minus}#it{K}^{+}) [GeV/#it{c}^{2}]");
  kk_d0->GetXaxis()->CenterTitle();
  kk_d0->GetXaxis()->SetTitleSize(0.047);
  kk_d0->GetXaxis()->SetLabelSize(0.045);
  
  kk_d0->SetYTitle("Candidates per 3 MeV/#it{c}^{2}");
  kk_d0->GetYaxis()->CenterTitle();
  kk_d0->GetXaxis()->SetTitleSize(0.047);
  kk_d0->GetXaxis()->SetLabelSize(0.045);
  kk_d0->GetYaxis()->SetTitleOffset(1.2);
 
  RooHist* hpullx_d0_mkk = kk_d0->pullHist();
  hpullx_d0_mkk->SetFillStyle(1001);
  hpullx_d0_mkk->SetFillColor(1);
  for(int i=0;i<hpullx_d0_mkk->GetN();++i) hpullx_d0_mkk->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullplotx_d0_mkk = mKK.frame(Title(" ")) ;
  pullplotx_d0_mkk->addPlotable(hpullx_d0_mkk,"B");
  pullplotx_d0_mkk->SetXTitle("#it{m}(#it{K}^{#minus}#it{K}^{+}) [GeV/#it{c}^{2}]");
  pullplotx_d0_mkk->GetXaxis()->CenterTitle();
  pullplotx_d0_mkk->SetMinimum(-8.);
  pullplotx_d0_mkk->SetMaximum(8.);
  pullplotx_d0_mkk->GetXaxis()->SetLabelSize(0.1);
  pullplotx_d0_mkk->GetXaxis()->SetTitleSize(0.1);
  pullplotx_d0_mkk->GetYaxis()->SetLabelSize(0.07);
  
  TCanvas* c1 = new TCanvas("c1", "c1",700, 700);
  Double_t xlow, ylow, xup, yup;
  c1->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
  c1->Divide(1,2);
  TVirtualPad *upPad1 = c1->GetPad(1);
  upPad1->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup);
  TVirtualPad *dwPad1 = c1->GetPad(2);
  dwPad1->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow));
  c1->Update();
  c1->cd(1);
  simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(signal_peak,signal_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(kpi_peak,kpi_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(mult_peak,mult_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(ds_peak)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
  simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(comb_bkg),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
  //simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(johnson_sig_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kViolet));
  //simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(gauss1_sig_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlue));
  //simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(gauss2_sig_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kOrange+5));
  
  kk_d0->Draw();
  upPad1->SetLogy();
  c1->cd(2);
  pullplotx_d0_mkk->Draw();
  c1->Update();
  c1->SaveAs("c1_mKK_d0.root");
    
  //Plot "D0bar" slice of simultaneous pdf for m(KK).
  RooPlot *kk_d0bar = mKK.frame(Title("D0bar"));
  combData.plotOn(kk_d0bar,Cut("sample==sample::D0bar"));
  simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),ProjWData(sample,combData));
  
  kk_d0bar->SetXTitle("#it{m}(#it{K}^{#minus}#it{K}^{+}) [GeV/#it{c}^{2}]");
  kk_d0bar->GetXaxis()->CenterTitle();
  kk_d0bar->GetXaxis()->SetTitleSize(0.047);
  kk_d0bar->GetXaxis()->SetLabelSize(0.045);
  
  kk_d0bar->SetYTitle("Candidates per 3 MeV/#it{c}^{2}");
  kk_d0bar->GetYaxis()->CenterTitle();
  kk_d0bar->GetXaxis()->SetTitleSize(0.047);
  kk_d0bar->GetXaxis()->SetLabelSize(0.045);
  kk_d0bar->GetYaxis()->SetTitleOffset(1.2);
 
  RooHist* hpullx_d0bar_mkk = kk_d0bar->pullHist();
  hpullx_d0bar_mkk->SetFillStyle(1001);
  hpullx_d0bar_mkk->SetFillColor(1);
  for(int i=0;i<hpullx_d0bar_mkk->GetN();++i) hpullx_d0bar_mkk->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullplotx_d0bar_mkk = mKK.frame(Title(" ")) ;
  pullplotx_d0bar_mkk->addPlotable(hpullx_d0bar_mkk,"B");
  pullplotx_d0bar_mkk->SetXTitle("#it{m}(#it{K}^{#minus}#it{K}^{+}) [GeV/#it{c}^{2}]");
  pullplotx_d0bar_mkk->SetMinimum(-8.);
  pullplotx_d0bar_mkk->SetMaximum(8.);
  pullplotx_d0bar_mkk->GetXaxis()->SetLabelSize(0.1);
  pullplotx_d0bar_mkk->GetXaxis()->SetTitleSize(0.1);
  pullplotx_d0bar_mkk->GetYaxis()->SetLabelSize(0.07);
  
  TCanvas* c2 = new TCanvas("c2", "c2",700, 700);
  c2->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
  c2->Divide(1,2);
  TVirtualPad *upPad2 = c2->GetPad(1);
  upPad2->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup);
  TVirtualPad *dwPad2 = c2->GetPad(2);
  dwPad2->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow));
  c2->Update();
  c2->cd(1);
  simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(signal_peak,signal_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(kpi_peak,kpi_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(mult_peak,mult_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(ds_peak)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
  simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(comb_bkg),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));

  //simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(johnson_sig_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kViolet));
  //simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(gauss1_sig_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlue));
  //simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(gauss2_sig_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kOrange+5));
  kk_d0bar->Draw();
  upPad2->SetLogy();
  c2->cd(2);
  pullplotx_d0bar_mkk->Draw();
  c2->Update();
  c2->SaveAs("c2_mKK_d0bar.root");
  
  // Plot "D0" slice of simultaneous pdf for m(D0pi).
  RooPlot *d0pi_d0 = mD0pi.frame(Title("D0"));
  combData.plotOn(d0pi_d0,Cut("sample==sample::D0"));
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),ProjWData(sample,combData));
  //simPdf.paramOn(d0pi_d0,dataxy_d0);
 
  d0pi_d0->SetXTitle("#it{m}(#it{D}^{0}#it{#pi}^{#plus}) [GeV/#it{c}^{2}]");
  d0pi_d0->GetXaxis()->CenterTitle();
  d0pi_d0->GetXaxis()->SetTitleSize(0.047);
  d0pi_d0->GetXaxis()->SetLabelSize(0.045);
  
  d0pi_d0->SetYTitle("Candidates per 0.2 MeV/#it{c}^{2}");
  d0pi_d0->GetYaxis()->CenterTitle();
  d0pi_d0->GetXaxis()->SetTitleSize(0.047);
  d0pi_d0->GetXaxis()->SetLabelSize(0.045);
  d0pi_d0->GetYaxis()->SetTitleOffset(1.2);

  RooHist* hpully_d0_md0pi = d0pi_d0->pullHist();
  hpully_d0_md0pi->SetFillStyle(1001);
  hpully_d0_md0pi->SetFillColor(1);
  for(int i=0;i<hpully_d0_md0pi->GetN();++i) hpully_d0_md0pi->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullploty_d0_md0pi = mD0pi.frame(Title(" ")) ;
  pullploty_d0_md0pi->addPlotable(hpully_d0_md0pi,"B");
  pullploty_d0_md0pi->SetXTitle("#it{m}(#it{D}^{0}#it{#pi}^{#plus}) [GeV/#it{c}^{2}]");
  pullploty_d0_md0pi->SetMinimum(-8.);
  pullploty_d0_md0pi->SetMaximum(8.);
  pullploty_d0_md0pi->GetXaxis()->SetLabelSize(0.1);
  pullploty_d0_md0pi->GetXaxis()->SetTitleSize(0.1);
  pullploty_d0_md0pi->GetYaxis()->SetLabelSize(0.07);
  
  TCanvas* c3 = new TCanvas("c3", "c3",700, 700);
  c3->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
  c3->Divide(1,2);
  TVirtualPad *upPad3 = c3->GetPad(1);
  upPad3->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup);
  TVirtualPad *dwPad3 = c3->GetPad(2);
  dwPad3->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow));
  c3->Update();
  c3->cd(1);
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(signal_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(kpi_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(mult_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(ds_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(comb_bkg),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
  simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(RooArgSet(signal_rnd,kpi_rnd,mult_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kOrange));
  d0pi_d0->Draw();
  upPad3->SetLogy();
  c3->cd(2);
  pullploty_d0_md0pi->Draw();
  c3->Update();
  c3->SaveAs("c3_md0pi_d0.root");
  
  //Plot "D0bar" slice of simultaneous pdf for m(D0pi). 
  RooPlot *d0pi_d0bar = mD0pi.frame(Title("D0bar"));
  combData.plotOn(d0pi_d0bar,Cut("sample==sample::D0bar"));
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),ProjWData(sample,combData));
  simPdf.paramOn(d0pi_d0bar,dataxy_d0bar);

  d0pi_d0bar->SetXTitle("#it{m}(#it{D}^{0}#it{#pi}^{#plus}) [GeV/#it{c}^{2}]");
  d0pi_d0bar->GetXaxis()->CenterTitle();
  d0pi_d0bar->GetXaxis()->SetTitleSize(0.047);
  d0pi_d0bar->GetXaxis()->SetLabelSize(0.045);
  
  d0pi_d0bar->SetYTitle("Candidates per 0.2 MeV/#it{c}^{2}");
  d0pi_d0bar->GetYaxis()->CenterTitle();
  d0pi_d0bar->GetXaxis()->SetTitleSize(0.047);
  d0pi_d0bar->GetXaxis()->SetLabelSize(0.045);
  d0pi_d0bar->GetYaxis()->SetTitleOffset(1.2);
 
  RooHist* hpully_d0bar_md0pi = d0pi_d0bar->pullHist();
  hpully_d0bar_md0pi->SetFillStyle(1001);
  hpully_d0bar_md0pi->SetFillColor(1);
  for(int i=0;i<hpully_d0bar_md0pi->GetN();++i) hpully_d0bar_md0pi->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullploty_d0bar_md0pi = mD0pi.frame(Title(" ")) ;
  pullploty_d0bar_md0pi->addPlotable(hpully_d0bar_md0pi,"B");
  pullploty_d0bar_md0pi->SetXTitle("#it{m}(#it{D}^{0}#it{#pi}^{#plus}) [GeV/#it{c}^{2}]");
  pullploty_d0bar_md0pi->SetMinimum(-8.);
  pullploty_d0bar_md0pi->SetMaximum(8.);
  pullploty_d0bar_md0pi->GetXaxis()->SetLabelSize(0.1);
  pullploty_d0bar_md0pi->GetXaxis()->SetTitleSize(0.1);
  pullploty_d0bar_md0pi->GetYaxis()->SetLabelSize(0.07);
  
  TCanvas* c4 = new TCanvas("c4", "c4",700, 700);
  c4->GetPad(0)->GetPadPar(xlow,ylow,xup,yup);
  c4->Divide(1,2);
  TVirtualPad *upPad4 = c4->GetPad(1);
  upPad4->SetPad(xlow,ylow+0.25*(yup-ylow),xup,yup);
  TVirtualPad *dwPad4 = c4->GetPad(2);
  dwPad4->SetPad(xlow,ylow,xup,ylow+0.25*(yup-ylow));
  c4->Update();
  c4->cd(1);
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(signal_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(kpi_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(mult_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(ds_peak),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(comb_bkg),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
  simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(signal_rnd,kpi_rnd,mult_rnd)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kOrange));
  d0pi_d0bar->Draw();
  upPad4->SetLogy();
  c4->cd(2);
  pullploty_d0bar_md0pi->Draw();
  c4->Update();
  c4->SaveAs("c4_md0pi_d0bar.root");
}