#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 <fstream>
#include <string>
#include "TMath.h"
#include <algorithm>
#include <functional>
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "TFile.h"
#include "RooDataSet.h"
#include "RooArgSet.h"
#include "myRooJohnsonSU.cpp"
#include "myRooPolBG.cpp"
#include "TTree.h"
#include "TROOT.h"
#include "TSystemDirectory.h"

using namespace RooFit ;
using namespace std;


void asymmetry_all_params_arg(const char* fname)
{

  //File and the tree containing the variable
  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_dst("f_{J_sig_dst}", "rel_frac_j",0.8937);
  RooRealVar A_fJ_sig_dst("A_{#fJ_sig_dst}","Asymmetry in rel frac for d0pi",0);
  RooFormulaVar fJ_sig_dst_d0("fJ_sig_dst_d0","rel frac for d0pi","(@0*(1 + @1))",RooArgList(f_J_sig_dst,A_fJ_sig_dst));
  RooFormulaVar fJ_sig_dst_d0bar("fJ_sig_dst_d0bar","rel frac for d0pi","(@0*(1- @1))",RooArgList(f_J_sig_dst,A_fJ_sig_dst));

  RooRealVar mu_J_sig_dst("#mu_{J_sig_dst}","mean_johnson",2.0102529);
  RooRealVar A_muJ_sig_dst("A_{#muJ_sig_dst}","asymm. in mean",0);
  RooFormulaVar mu_J_sig_dst_d0("mu_J_sig_dst_d0","mu_J_sig_dst_d0","(@0*(1 + @1))",RooArgList(mu_J_sig_dst, A_muJ_sig_dst));
  RooFormulaVar mu_J_sig_dst_d0bar("mu_J_sig_dst_d0bar","mu_J_sig_dst_d0bar","(@0*(1- @1))",RooArgList(mu_J_sig_dst, A_muJ_sig_dst));
  
  RooRealVar sigma_J_sig_dst("#sigma_{J_sig_dst}", "sigma_johnson",0.0001675);
  RooRealVar A_sigmaJ_sig_dst("A_{#sigma_Jsig_dst}","asymm. in sigma",0);
  RooFormulaVar sigma_J_sig_dst_d0("sigma_J_sig_dst_d0","sigma_J_sig_dst_d0","(@0*(1 + @1))",RooArgList(sigma_J_sig_dst, A_sigmaJ_sig_dst));
  RooFormulaVar sigma_J_sig_dst_d0bar("sigma_J_sig_dst_d0bar","sigma_J_sig_dst_d0bar","(@0*(1- @1))",RooArgList(sigma_J_sig_dst, A_sigmaJ_sig_dst));
  
  RooRealVar gamma_J_sig_dst("#gamma_{J_sig_dst}", "gamma", -0.0329);
  RooRealVar A_gammaJ_sig_dst("A_{#gamma_Jsig_dst}","asymm. in sigma",0);
  RooFormulaVar gamma_J_sig_dst_d0("gamma_J_sig_dst_d0","gamma_J_sig_dst_d0","(@0*(1 + @1))",RooArgList(gamma_J_sig_dst, A_gammaJ_sig_dst));
  RooFormulaVar gamma_J_sig_dst_d0bar("gamma_J_sig_dst_d0bar","gamma_J_sig_dst_d0bar","(@0*(1- @1))",RooArgList(gamma_J_sig_dst, A_gammaJ_sig_dst));

  RooRealVar delta_J_sig_dst("#delta_{J_sig_dst}", "delta", 1.033);
  RooRealVar A_deltaJ_sig_dst("A_{#delta_Jsig_dst}","asymm. in delta",0);
  RooFormulaVar delta_J_sig_dst_d0("delta_J_sig_dst_d0","delta_J_sig_dst_d0","(@0*(1 + @1))",RooArgList(delta_J_sig_dst, A_deltaJ_sig_dst));
  RooFormulaVar delta_J_sig_dst_d0bar("delta_J_sig_dst_d0bar","delta_J_sig_dst_d0bar","(@0*(1- @1))",RooArgList(delta_J_sig_dst, A_deltaJ_sig_dst));
  
  myRooJohnsonSU johnson_sig_dst_d0("johnson_sig_dst_d0","johnson_sig_d0", mD0pi, mu_J_sig_dst_d0, sigma_J_sig_dst_d0,  delta_J_sig_dst_d0,  gamma_J_sig_dst_d0);
  myRooJohnsonSU johnson_sig_dst_d0bar("johnson_sig_dst_d0bar","johnson_sig_d0bar", mD0pi, mu_J_sig_dst_d0bar, sigma_J_sig_dst_d0bar,  delta_J_sig_dst_d0bar,  gamma_J_sig_dst_d0bar);

  RooRealVar f_G1_sig_dst("f_{G1_sig_dst}", "frac_gauss1",0.136);
  RooRealVar A_fG1_sig_dst("A_{#fG1_sig_dst}","Asymmetry in rel frac (gaus) d0pi",0);
  RooFormulaVar fG1_sig_dst_d0("fG1_sig_dst_d0","rel frac for d0pi","(@0*(1 + @1))",RooArgList(f_G1_sig_dst, A_fG1_sig_dst));
  RooFormulaVar fG1_sig_dst_d0bar("fG1_sig_dst_d0bar","rel frac for d0pi","(@0*(1- @1))",RooArgList(f_G1_sig_dst, A_fG1_sig_dst));
  
  RooRealVar mu_G1_sig_dst("#mu_{G1_sig_dst}", "mean_gauss1",2.010300);
  RooRealVar A_muG1_sig_dst("A_{#muG1_sig_dst}","asymm. in mean",0);
  RooFormulaVar muG1_sig_dst_d0("muG1_sig_dst_d0","muG1_sig_dst_d0","(@0*(1 + @1))",RooArgList(mu_G1_sig_dst, A_muG1_sig_dst));
  RooFormulaVar muG1_sig_dst_d0bar("muG1_sig_dst_d0bar","muG1_sig_dst_d0bar","(@0*(1- @1))",RooArgList(mu_G1_sig_dst, A_muG1_sig_dst));
  
  RooRealVar sigma_G1_sig_dst("#sigma_{G1_sig_dst}", "sigma_gauss1",0.00385);
  RooRealVar A_sigmaG1_sig_dst("A_{#sigma_G1sig_dst}","asymm. in mean",0);
  RooFormulaVar sigmaG1_sig_dst_d0("sigma_G1sig_dst_d0","sigma_G1sig_dst_d0","(@0*(1 + @1))",RooArgList(sigma_G1_sig_dst, A_sigmaG1_sig_dst));
  RooFormulaVar sigmaG1_sig_dst_d0bar("sigma_G1sig_dst_d0bar","sigma_G1sig_dst_d0bar","(@0*(1- @1))",RooArgList(sigma_G1_sig_dst, A_sigmaG1_sig_dst));
  
  RooRealVar sigma_G2_sig_dst("#sigma_{G2_sig_dst}", "sigma_gauss2",0.001224);
  RooRealVar A_sigmaG2_sig_dst("A_{#sigma_G2sig_dst}","asymm. in sigma",0);
  RooFormulaVar sigmaG2_sig_dst_d0("sigma_G2sig_dst_d0","sigma_G2sig_dst_d0","(@0*(1 + @1))",RooArgList(sigma_G2_sig_dst, A_sigmaG2_sig_dst));
  RooFormulaVar sigmaG2_sig_dst_d0bar("sigma_G2sig_dst_d0bar","sigma_G2sig_dst_d0bar","(@0*(1- @1))",RooArgList(sigma_G2_sig_dst, A_sigmaG2_sig_dst));
  
  RooGaussian gauss1_sig_dst_d0("gauss1_sig_d0", "1st gaussian PDF_dst_ds0", mD0pi, muG1_sig_dst_d0, sigmaG1_sig_dst_d0);
  RooGaussian gauss1_sig_dst_d0bar("gauss1_sig_d0bar", "1st gaussian PDF_dst_d0", mD0pi, muG1_sig_dst_d0bar, sigmaG1_sig_dst_d0bar);

  RooGaussian gauss2_sig_dst_d0("gauss2_sig_dst_d0", "2nd gaussian PDF_dst_d0", mD0pi, muG1_sig_dst_d0, sigmaG2_sig_dst_d0);
  RooGaussian gauss2_sig_dst_d0bar("gauss2_sig_dst_d0bar", "2nd gaussian PDF_dst_d0bar", mD0pi, muG1_sig_dst_d0bar, sigmaG2_sig_dst_d0bar);
  
  RooAddPdf doublegaussian_sig_dst_d0("doublegaussian_sig_dst_d0", "doublegaussian_sig_dst_d0", gauss1_sig_dst_d0, gauss2_sig_dst_d0, fG1_sig_dst_d0);
  RooAddPdf doublegaussian_sig_dst_d0bar("doublegaussian_sig_dst_d0bar", "doublegaussian_sig_dst_d0bar", gauss1_sig_dst_d0bar, gauss2_sig_dst_d0bar,  fG1_sig_dst_d0bar);
  
  RooAddPdf sig_mD0pi_d0("sig_mD0pi_d0","sig_mD0pi_d0", johnson_sig_dst_d0, doublegaussian_sig_dst_d0, fJ_sig_dst_d0);
  RooAddPdf sig_mD0pi_d0bar("sig_mD0pi_d0bar","sig_mD0pi_d0bar", johnson_sig_dst_d0bar, doublegaussian_sig_dst_d0bar, fJ_sig_dst_d0bar);
  
  //kk
  RooRealVar f_J_sig_dz("f_{J_sig_dz}", "relative_frac_j",0.9624);
  RooRealVar A_fJ_sig_dz("A_{#fJ_sig_dz}","Asymmetry in rel frac for dz",0);
  RooFormulaVar fJ_sig_dz_d0("fJ_sig_dz_d0","fJ_sig_dz_d0","(@0*(1 + @1))",RooArgList(f_J_sig_dz,A_fJ_sig_dz));
  RooFormulaVar fJ_sig_dz_d0bar("fJ_sig_dz_d0bar","fJ_sig_dz_d0bar","(@0*(1- @1))",RooArgList(f_J_sig_dz,A_fJ_sig_dz));

  RooRealVar mu_J_sig_dz("#mu_{J_sig_dz}","mean_johnson",1.864738);
  RooRealVar A_muJ_sig_dz("A_{#muJ_sig_dz}","asymm. in mean",0);
  RooFormulaVar mu_J_sig_dz_d0("mu_J_sig_dz_d0","mu_J_sig_dz_d0","(@0*(1 + @1))",RooArgList(mu_J_sig_dz, A_muJ_sig_dz));
  RooFormulaVar mu_J_sig_dz_d0bar("mu_J_sig_dz_d0bar","mu_J_sig_dz_d0bar","(@0*(1- @1))",RooArgList(mu_J_sig_dz, A_muJ_sig_dz));
  
  RooRealVar sigma_J_sig_dz("#sigma_{J_sig_dz}", "sigma_johnson", 0.00612);
  RooRealVar A_sigmaJ_sig_dz("A_{#sigma_Jsig_dz}","asymm. in sigma",0);
  RooFormulaVar sigma_J_sig_dz_d0("sigma_J_sig_dz_d0","sigma_J_sig_dz_d0","(@0*(1 + @1))",RooArgList(sigma_J_sig_dz, A_sigmaJ_sig_dz));
  RooFormulaVar sigma_J_sig_dz_d0bar("sigma_J_sig_dz_d0bar","sigma_J_sig_dz_d0bar","(@0*(1- @1))",RooArgList(sigma_J_sig_dz, A_sigmaJ_sig_dz));
  RooRealVar gamma_J_sig_dz("#gamma_{J_sig_dz}", "gamma",0.011);
  RooRealVar A_gammaJ_sig_dz("A_{#gamma_Jsig_dz}","asymm. in gamma",0);
  RooFormulaVar gamma_J_sig_dz_d0("gamma_J_sig_dz_d0","gamma_J_sig_dz_d0","(@0*(1 + @1))",RooArgList(gamma_J_sig_dz, A_gammaJ_sig_dz));
  RooFormulaVar gamma_J_sig_dz_d0bar("gamma_J_sig_dz_d0bar","gamma_J_sig_dz_d0bar","(@0*(1- @1))",RooArgList(gamma_J_sig_dz, A_gammaJ_sig_dz));
  
  RooRealVar delta_J_sig_dz("#delta_{J_sig_dz}", "delta",1.689);
  RooRealVar A_deltaJ_sig_dz("A_{#delta_Jsig_dz}","asymm. in delta",0);
  RooFormulaVar delta_J_sig_dz_d0("delta_J_sig_dz_d0","delta_J_sig_dz_d0","(@0*(1 + @1))",RooArgList(delta_J_sig_dz, A_deltaJ_sig_dz));
  RooFormulaVar delta_J_sig_dz_d0bar("delta_J_sig_dz_d0bar","delta_J_sig_dz_d0bar","(@0*(1- @1))",RooArgList(delta_J_sig_dz, A_deltaJ_sig_dz));
  
  myRooJohnsonSU johnson_sig_d0("johnson_sig_d0", "johnson_sig_d0", mKK, mu_J_sig_dz_d0,sigma_J_sig_dz_d0, delta_J_sig_dz_d0, gamma_J_sig_dz_d0);
  myRooJohnsonSU johnson_sig_d0bar("johnson_sig_d0bar", "johnson_sig_d0bar", mKK, mu_J_sig_dz_d0bar, sigma_J_sig_dz_d0bar, delta_J_sig_dz_d0bar, gamma_J_sig_dz_d0bar); 
  RooRealVar f_G1_sig_dz("f_{G1}_sig_dz", "relative_frac_g",0.451);
  RooRealVar A_fG1_sig_dz("A_{#fG1_sig_dz}","asymm. in fG",0);
  RooFormulaVar fG1_sig_dz_d0("fG1_sig_dz_d0","fG1_sig_dz_d0","(@0*(1 + @1))",RooArgList(f_G1_sig_dz, A_fG1_sig_dz));
  RooFormulaVar fG1_sig_dz_d0bar("fG1_sig_dz_d0bar","fG1_sig_dz_d0bar","(@0*(1- @1))",RooArgList(f_G1_sig_dz, A_fG1_sig_dz));
  
  RooRealVar mu_G1_sig_dz("#mu_{G1}_sig_dz", "mean_gauss1",1.8506);
  RooRealVar A_muG1_sig_dz("A_{#muG1_sig_dz}","asymm. in fG",0);
  RooFormulaVar muG1_sig_dz_d0("muG1_sig_dz_d0","muG1_sig_dz_d0","(@0*(1 + @1))",RooArgList(mu_G1_sig_dz, A_muG1_sig_dz));
  RooFormulaVar muG1_sig_dz_d0bar("muG1_sig_dz_d0bar","muG1_sig_dz_d0bar","(@0*(1- @1))",RooArgList(mu_G1_sig_dz, A_muG1_sig_dz));
  
  RooRealVar sigma_G1_sig_dz("#sigma_{G1}_sig_dz", "sigma_gauss1",0.0616);
  RooRealVar A_sigmaG1_sig_dz("A_{#sigma_G1sig_dz}","asymm. in sigma",0);
  RooFormulaVar sigmaG1_sig_dz_d0("sigma_G1sig_dz_d0","sigma_G1sig_dz_d0","(@0*(1 + @1))",RooArgList(sigma_G1_sig_dz, A_sigmaG1_sig_dz));
  RooFormulaVar sigmaG1_sig_dz_d0bar("sigma_G1sig_dz_d0bar","sigma_G1sig_dz_d0bar","(@0*(1- @1))",RooArgList(sigma_G1_sig_dz, A_sigmaG1_sig_dz));
  
  RooRealVar sigma_G2_sig_dz("#sigma_{G2}_sig_dz", "sigma_gauss2",0.0225);
  RooRealVar A_sigmaG2_sig_dz("A_{#sigma_G2sig_dz}","asymm. in sigma",0);
  RooFormulaVar sigmaG2_sig_dz_d0("sigma_G2sig_dz_d0","sigma_G2sig_dz_d0","(@0*(1 + @1))",RooArgList(sigma_G2_sig_dz, A_sigmaG2_sig_dz));
  RooFormulaVar sigmaG2_sig_dz_d0bar("sigma_G2sig_dz_d0bar","sigma_G2sig_dz_d0","(@0*(1- @1))",RooArgList(sigma_G2_sig_dz, A_sigmaG2_sig_dz));
  
  RooGaussian gauss1_sig_dz_d0("gauss1_sig_dz_d0", "1st gaussian PDF D0", mKK, muG1_sig_dz_d0, sigmaG1_sig_dz_d0);
  RooGaussian gauss1_sig_dz_d0bar("gauss1_sig_dz_d0bar", "1st gaussian PDF D0bar", mKK, muG1_sig_dz_d0, sigmaG1_sig_dz_d0);
  
  RooGaussian gauss2_sig_dz_d0("gauss2_sig_dz_d0", "2nd gaussian PDF D0", mKK, muG1_sig_dz_d0, sigmaG2_sig_dz_d0);
  RooGaussian gauss2_sig_dz_d0bar("gauss2_sig_dz_d0bar", "2nd gaussian PDF D0bar", mKK, muG1_sig_dz_d0bar, sigmaG2_sig_dz_d0bar);
  
  RooAddPdf doublegaussian_sig_dz_d0("doublegaussian_sig_dz_d0", "doublegaussian_sig_dz_d0", gauss1_sig_dz_d0, gauss2_sig_dz_d0, fG1_sig_dz_d0);
  RooAddPdf doublegaussian_sig_dz_d0bar("doublegaussian_sig_dz_d0bar", "doublegaussian_sig_dz_d0bar", gauss1_sig_dz_d0, gauss2_sig_dz_d0bar, fG1_sig_dz_d0bar);
  RooAddPdf sig_mKK_d0("sig_mKK_d0","sig_mKK_d0", johnson_sig_d0, doublegaussian_sig_dz_d0,fJ_sig_dz_d0);
  RooAddPdf sig_mKK_d0bar("sig_mKK_d0bar","sig_mKK_d0bar", johnson_sig_d0bar, doublegaussian_sig_dz_d0bar, fJ_sig_dz_d0bar);
  
  //signal peak
  RooProdPdf signal_peak_d0("signal_peak_d0", "signal_peak_d0", sig_mKK_d0, sig_mD0pi_d0);
  RooProdPdf signal_peak_d0bar("signal_peak_d0bar", "signal_peak_d0bar", sig_mKK_d0bar, sig_mD0pi_d0bar);
  
  //random pion (common for all)
  RooRealVar threshold("threshold", "no data beyond this", 2.00441);
  RooRealVar alpha("#alpha", "alpha", -16, -1e6, 0.);
  RooRealVar A_alpha("A_{#alpha}","Asymmetry in random pion",0);
  RooFormulaVar alpha_d0("alpha_d0","alpha for Ds peak (D0)","(@0*(1 + @1))",RooArgList(alpha, A_alpha));
  RooFormulaVar alpha_d0bar("alpha_d0bar","alpha for Ds peak (D0bar)","(@0*(1- @1))",RooArgList(alpha, A_alpha));
  RooRealVar beta("#beta", "beta", 0.);//, -1e6, 1e6);
  myRooPolBG bkg_mD0pi_d0("bkg_mD0pi_d0","bkg_mD0pi", mD0pi, threshold, alpha_d0, beta);
  myRooPolBG bkg_mD0pi_d0bar("bkg_mD0pi_d0bar","bkg_mD0pi", mD0pi, threshold, alpha_d0bar, beta);
  
  //signal & random pion
  RooProdPdf signal_rnd_d0("signal_rnd_d0", "signal_rnd_d0", sig_mKK_d0, bkg_mD0pi_d0);
  RooProdPdf signal_rnd_d0bar("signal_rnd_d0bar", "signal_rnd_d0bar", sig_mKK_d0bar, bkg_mD0pi_d0bar);
  
  
  /*-----------------------M U L T I B O D Y-------------------------*/
  
  //d0pi
  RooRealVar delta_J_mult_dst("#delta_{J}_mult", "delta_mult",-1.2859);
  RooRealVar A_deltaJ_mult_dst("A_{#delta_J_mult_dst}","asym. delta",0);
  RooFormulaVar delta_J_mult_dst_d0("delta_J_mult_dst_d0","delta_J_mult_dst_d0","(@0*(1 + @1))",RooArgList(delta_J_mult_dst,A_deltaJ_mult_dst));
  RooFormulaVar delta_J_mult_dst_d0bar("delta_J_mult_dst_d0bar","delta_J_mult_dst_d0bar","(@0*(1- @1))",RooArgList(delta_J_mult_dst,A_deltaJ_mult_dst));
  
  RooRealVar mu_J_mult_dst("#mu_{J}_mult","mean_johnson",2.008198);
  RooRealVar A_muJ_mult_dst("A_{#muJ_mult_dst}","asym. mu",0);
  RooFormulaVar mu_J_mult_dst_d0("mu_J_mult_dst_d0","mu_J_mult_dst_d0","(@0*(1 + @1))",RooArgList(mu_J_mult_dst, A_muJ_mult_dst));
  RooFormulaVar mu_J_mult_dst_d0bar("mu_J_mult_dst_d0bar","mu_J_mult_dst_d0bar","(@0*(1- @1))",RooArgList(mu_J_mult_dst, A_muJ_mult_dst));
  
  RooRealVar sigma_J_mult_dst("#sigma_{J}_mult", "sigma_johnson",0.00363);
  RooRealVar A_sigmaJ_mult_dst("A_{#sigma_J_mult_dst}","asym. sigma",0);
  RooFormulaVar sigma_J_mult_dst_d0("sigma_J_mult_dst_d0","sigma_J_mult_dst_d0","(@0*(1 + @1))",RooArgList(sigma_J_mult_dst, A_sigmaJ_mult_dst));
  RooFormulaVar sigma_J_mult_dst_d0bar("sigma_J_mult_dst_d0bar","sigma_J_mult_dst_d0bar","(@0*(1- @1))",RooArgList(sigma_J_mult_dst, A_sigmaJ_mult_dst));

  RooRealVar gamma_J_mult_dst("#gamma_{J}_mult", "gamma", 1.);
  RooRealVar A_gammaJ_mult_dst("A_{#gamma_J_mult_dst}","asym. gamma",0);
  RooFormulaVar gamma_J_mult_dst_d0("gamma_J_mult_dst_d0","gamma_J_mult_dst_d0","(@0*(1 + @1))",RooArgList(gamma_J_mult_dst,  A_gammaJ_mult_dst));
  RooFormulaVar gamma_J_mult_dst_d0bar("gamma_J_mult_dst_d0bar","gamma_J_mult_dst_d0bar","(@0*(1- @1))",RooArgList(gamma_J_mult_dst, A_gammaJ_mult_dst));
  
  myRooJohnsonSU mult_mD0pi_dst_d0("mult_mD0pi_dst_d0","mult_mD0pi_dst_d0", mD0pi, mu_J_mult_dst_d0,sigma_J_mult_dst_d0, delta_J_mult_dst_d0,  gamma_J_mult_dst_d0);
  myRooJohnsonSU mult_mD0pi_dst_d0bar("mult_mD0pi_dst_d0bar","mult_mD0pi_dst_d0bar", mD0pi,mu_J_mult_dst_d0bar,sigma_J_mult_dst_d0bar, delta_J_mult_dst_d0bar, gamma_J_mult_dst_d0bar);
  
  //kk
  RooRealVar a1_dz("a_{mult}", "A1", -12.388);
  RooRealVar A_a1_mult_dz("A_{#a1_mult_dz}","asym. a",0);
  RooFormulaVar a1_mult_dz_d0("a1_mult_dz_d0","a1_mult_dz_d0","(@0*(1 + @1))",RooArgList(a1_dz,  A_a1_mult_dz));
  RooFormulaVar a1_mult_dz_d0bar("a1_mult_dz_d0bar","a1_mult_dz_d0bar","(@0*(1- @1))",RooArgList(a1_dz,A_a1_mult_dz));
  
  RooExponential mult_mKK_d0("mult_mKK_d0", "mult_mKK_d0", mKK, a1_mult_dz_d0);
  RooExponential mult_mKK_d0bar("mult_mKK_d0bar", "mult_mKK_d0bar", mKK, a1_mult_dz_d0);  

  RooProdPdf mult_peak_d0("mult_peak_dz_d0", "mult_peak", mult_mKK_d0, mult_mD0pi_dst_d0);
  RooProdPdf mult_peak_d0bar("mult_peak_dz_d0bar", "mult_peak", mult_mKK_d0bar, mult_mD0pi_dst_d0bar);
  
  //multibody & random pion	
  RooProdPdf mult_rnd_d0("mult_rnd_d0","mult_rnd_d0", mult_mKK_d0, bkg_mD0pi_d0);
  RooProdPdf mult_rnd_d0bar("mult_rnd_d0bar","mult_rnd_d0bar", mult_mKK_d0bar, bkg_mD0pi_d0bar);
  
  /*-----------------------K Pi-------------------------*/
  
  //kk
  RooRealVar fG_d0kpi_dz("f_{G}d0kpi_dz", "relative fraction gaussian",0.385);
  RooRealVar A_fG_d0kpi_dz("A_{fG_d0kpi_dz}","asym mean",0);
  RooFormulaVar fG_d0kpi_dz_d0("fG_d0kpi_dz_d0","fG_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(fG_d0kpi_dz, A_fG_d0kpi_dz));
  RooFormulaVar fG_d0kpi_dz_d0bar("fG_d0kpi_dz_d0bar","fG_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(fG_d0kpi_dz, A_fG_d0kpi_dz));
  
  RooRealVar mu_G_d0kpi("#mu_{G}d0kpi","mean_gaus",1.94120);
  RooRealVar A_muG_d0kpi_dz("A_{#muG1_d0kpi_dz}","asym mean",0);
  RooFormulaVar muG_d0kpi_dz_d0("muG_d0kpi_dz_d0","muG_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(mu_G_d0kpi, A_muG_d0kpi_dz));
  RooFormulaVar muG_d0kpi_dz_d0bar("muG_d0kpi_dz_d0bar","muG_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(mu_G_d0kpi, A_muG_d0kpi_dz));
  
  RooRealVar sigma_G_d0kpi("#sigma_{G}d0kpi", "sigma_gaus",0.00716);
  RooRealVar A_sigmaG_d0kpi_dz("A_{#sigmaG_d0kpi_dz}","asym sigma",0);
  RooFormulaVar sigmaG_d0kpi_dz_d0("sigmaG_d0kpi_dz_d0","sigmaG_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(sigma_G_d0kpi, A_sigmaG_d0kpi_dz));
  RooFormulaVar sigmaG_d0kpi_dz_d0bar("sigmaG_d0kpi_dz_d0bar","sigmaG_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(sigma_G_d0kpi, A_sigmaG_d0kpi_dz));
  
  RooGaussian gauss_d0kpi_dz_d0("gauss_d0kpi_dz_d0", "gauss component d0", mKK,  muG_d0kpi_dz_d0, sigmaG_d0kpi_dz_d0);
  RooGaussian gauss_d0kpi_dz_d0bar("gauss_d0kpi_dz_d0bar", "gauss component d0bar", mKK,  muG_d0kpi_dz_d0bar, sigmaG_d0kpi_dz_d0bar);
  
  RooRealVar gamma_J_d0kpi_dz("#gamma_{J}d0kpi", "gamma",1.302);
  RooRealVar A_gammaJ_d0kpi_dz("A_{#gamma_J_d0kpi_dz}","asym gamma",0);
  RooFormulaVar gamma_J_d0kpi_dz_d0("gamma_J_d0kpi_dz_d0","gamma_J_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(gamma_J_d0kpi_dz, A_gammaJ_d0kpi_dz));
  RooFormulaVar gamma_J_d0kpi_dz_d0bar("gamma_J_d0kpi_dz_d0bar","gamma_J_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(gamma_J_d0kpi_dz, A_gammaJ_d0kpi_dz));
  
  RooRealVar delta_J_d0kpi_dz("#delta_{J}d0kpi", "delta",-0.3516);
  RooRealVar A_deltaJ_d0kpi_dz("A_{#delta_J_d0kpi_dz}","asym delta",0);
  RooFormulaVar delta_J_d0kpi_dz_d0("delta_J_d0kpi_dz_d0","delta_J_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(delta_J_d0kpi_dz,A_deltaJ_d0kpi_dz));
  RooFormulaVar delta_J_d0kpi_dz_d0bar("delta_J_d0kpi_dz_d0bar","delta_J_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(delta_J_d0kpi_dz,A_deltaJ_d0kpi_dz));
  
  RooRealVar mu_J_d0kpi_dz("#mu_{J}_d0kpi", "mean",1.95397);
  RooRealVar A_muJ_d0kpi_dz("A_{#muJ_d0kpi_dz}","asym mean",0);
  RooFormulaVar mu_J_d0kpi_dz_d0("mu_J_d0kpi_dz_d0","mu_J_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(mu_J_d0kpi_dz, A_muJ_d0kpi_dz));
  RooFormulaVar mu_J_d0kpi_dz_d0bar("mu_J_d0kpi_dz_d0bar","mu_J_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(mu_J_d0kpi_dz, A_muJ_d0kpi_dz));

  RooRealVar sigma_J_d0kpi_dz("#sigma_{J}_d0kpi", "sigma",0.00760);
  RooRealVar A_sigmaJ_d0kpi_dz("A_{#sigma_Jd0kpi_dz}","asym sigma",0);
  RooFormulaVar sigma_J_d0kpi_dz_d0("sigma_J_d0kpi_dz_d0","sigma_J_d0kpi_dz_d0","(@0*(1 + @1))",RooArgList(sigma_J_d0kpi_dz, A_sigmaJ_d0kpi_dz));
  RooFormulaVar sigma_J_d0kpi_dz_d0bar("sigma_J_d0kpi_dz_d0bar","sigma_J_d0kpi_dz_d0bar","(@0*(1- @1))",RooArgList(sigma_J_d0kpi_dz, A_sigmaJ_d0kpi_dz));
  
  myRooJohnsonSU johnson_d0kpi_dz_d0("johnson_d0kpi_dz_d0bar", "Johnson PDF", mKK,mu_J_d0kpi_dz_d0,  sigma_J_d0kpi_dz_d0,delta_J_d0kpi_dz_d0,gamma_J_d0kpi_dz_d0);
  myRooJohnsonSU johnson_d0kpi_dz_d0bar("johnson_d0kpi_dz_d0bar", "Johnson PDF", mKK,mu_J_d0kpi_dz_d0bar, sigma_J_d0kpi_dz_d0bar, delta_J_d0kpi_dz_d0bar, gamma_J_d0kpi_dz_d0bar);  
  RooAddPdf kpi_mKK_dz_d0("kpi_mKK_dz_d0","kpi_mKK",gauss_d0kpi_dz_d0,johnson_d0kpi_dz_d0,fG_d0kpi_dz_d0);
  RooAddPdf kpi_mKK_dz_d0bar("kpi_mKK_dz_d0bar","kpi_mKK",gauss_d0kpi_dz_d0bar,johnson_d0kpi_dz_d0bar,fG_d0kpi_dz_d0bar);
  
  //kpi peaking background
  RooProdPdf kpi_peak_d0("kpi_peak_d0", "kpi_peak_d0", kpi_mKK_dz_d0, sig_mD0pi_d0);//as the shape of m(d0pi) is same as signal
  RooProdPdf kpi_peak_d0bar("kpi_peak_d0bar", "kpi_peakbar_d0",  kpi_mKK_dz_d0bar, sig_mD0pi_d0bar);

  //kpi & random pion
  RooProdPdf kpi_rnd_d0("kpi_rnd_d0", "kpi_rnd_d0",  kpi_mKK_dz_d0, bkg_mD0pi_d0);
  RooProdPdf kpi_rnd_d0bar("kpi_rnd_d0bar", "kpi_rnd_d0", kpi_mKK_dz_d0bar, bkg_mD0pi_d0bar);

  /*----------------------- Ds+ -> K+ K- pi+ -------------------------*/
  //d0pi
  RooRealVar d_dst("d_{s}", "Slope of Polynomial for Ds",0.790);
  RooRealVar A_d_ds_dst("A_d_{ds}"," asym in slope of Polynomial for Ds",0);
  RooFormulaVar d_ds_dst_d0("d_ds_dst_d0","d_ds_dst_d0","(@0*(1 + @1))",RooArgList(d_dst, A_d_ds_dst));
  RooFormulaVar d_ds_dst_d0bar("d_ds_dst_d0bar","d_ds_dst_d0bar","(@0*(1- @1))",RooArgList(d_dst, A_d_ds_dst));
  RooChebychev ds_mD0pi_d0("ds_mD0pi_d0","Chebshev Polynomial", mD0pi,RooArgList(d_ds_dst_d0));
  RooChebychev ds_mD0pi_d0bar("ds_mD0pi_d0bar","Chebshev Polynomial", mD0pi,RooArgList(d_ds_dst_d0bar));

  //kk
  RooRealVar mu_0_ds("#mu_{0}_ds", "",0.000214269);
  RooRealVar A_mu0_ds_dz("A_{mu0}_ds","asym mu0",0);
  RooFormulaVar mu0_ds_dz_d0("mu0_ds_dz_d0","mu0_ds_dz_d0","(@0*(1 + @1))",RooArgList(mu_0_ds,A_mu0_ds_dz));
  RooFormulaVar mu0_ds_dz_d0bar("mu0_ds_dz_d0bar","mu0_ds_dz_d0bar","(@0*(1- @1))",RooArgList(mu_0_ds, A_mu0_ds_dz));
    
  RooRealVar constant_ds("constant_ds", "slope",3.83333);
  RooFormulaVar mu_ds_dz_d0("mu_ds_dz_d0", "@0+@1-@2",RooArgList(mu0_ds_dz_d0,constant_ds, mD0pi));
  RooFormulaVar mu_ds_dz_d0bar("mu_ds_dz_d0", "@0+@1-@2",RooArgList(mu0_ds_dz_d0bar,constant_ds, mD0pi));
  
  RooRealVar sigma_ds("#sigma_{ds}", "",0.00684);
  RooRealVar A_sigmads("A_{#sigma_ds}","Asymmetry in sigma for Ds peak",0);
  RooFormulaVar sigmads_d0("sigmads_d0","sigma for Ds peak (D0)","(@0*(1 + @1))",RooArgList(sigma_ds,A_sigmads));
  RooFormulaVar sigmads_d0bar("sigmads_d0bar","sigma for Ds peak (D0bar)","(@0*(1 - @1))",RooArgList(sigma_ds,A_sigmads));
  
  RooRealVar gamma_ds("#gamma_{ds}", "",-0.1118);
  RooRealVar A_gammads("A_{#gamma_ds}","Asymmetry in gamma for Ds peak",0);
  RooFormulaVar gammads_d0("gammads_d0","gamma for Ds peak (D0)","(@0*(1 + @1))",RooArgList(gamma_ds,A_gammads));
  RooFormulaVar gammads_d0bar("gammads_d0bar","gammma for Ds peak (D0bar)","(@0*(1 - @1))",RooArgList(gamma_ds,A_gammads));

  RooRealVar delta_ds("#delta_{ds}", "",1.125);
  RooRealVar A_deltads("A_{#delta_ds}","Asymmetry in delta for Ds peak",0);
  RooFormulaVar deltads_d0("deltads_d0","delta for Ds peak (D0)","(@0*(1 + @1))",RooArgList(delta_ds,A_deltads));
  RooFormulaVar deltads_d0bar("deltads_d0bar","delta for Ds peak (D0bar)","(@0*(1- @1))",RooArgList(delta_ds, A_deltads));
  
  myRooJohnsonSU ds_mKK_d0("ds_mKK_d0","johnson",mKK, mu_ds_dz_d0,sigmads_d0,deltads_d0,gammads_d0);
  myRooJohnsonSU ds_mKK_d0bar("ds_mKK_d0bar","johnson",mKK, mu_ds_dz_d0bar,sigmads_d0bar,deltads_d0bar,gammads_d0bar);
  
  RooProdPdf ds_peak_d0("ds_peak_d0", "ds_peak", ds_mD0pi_d0, Conditional(ds_mKK_d0, mKK));
  RooProdPdf ds_peak_d0bar("ds_peak_d0bar", "ds_peak", ds_mD0pi_d0bar, Conditional(ds_mKK_d0bar, mKK)); 
  
  
  /*------------------------C O M B I N A T O R I A L----------------------------*/
  //Chebychev
  RooRealVar c("c", "Slope1 of Polynomial",-0.2,-1,1);
  RooRealVar A_c("A_{#c}","Asymmetry in combinatorial",0);
  RooFormulaVar c_d0("c_d0","c for comb (D0)","(@0*(1 + @1))",RooArgList(c, A_c));
  RooFormulaVar c_d0bar("c_d0bar","c for comb (D0bar)","(@0*(1- @1))",RooArgList(c, A_c));
  RooChebychev comb_mKK_d0("chebpol_d0","Chebshev Polynomial d0", mKK,RooArgList(c_d0));
  RooChebychev comb_mKK_d0bar("chebpol_d0bar","Chebshev Polynomial d0bar", mKK,RooArgList(c_d0bar));

  RooProdPdf comb_bkg_d0("comb_bkg_d0", "comb_bkg_d0", comb_mKK_d0, bkg_mD0pi_d0);
  RooProdPdf comb_bkg_d0bar("comb_bkg_d0bar", "comb_bkg_d0", comb_mKK_d0, bkg_mD0pi_d0bar);
  
  /*------------------------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_d0,signal_rnd_d0, mult_peak_d0,mult_rnd_d0,kpi_peak_d0, kpi_rnd_d0,ds_peak_d0, comb_bkg_d0), 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_d0bar,signal_rnd_d0bar, mult_peak_d0bar,mult_rnd_d0bar,kpi_peak_d0bar, kpi_rnd_d0bar,ds_peak_d0bar, comb_bkg_d0bar), 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));

  //cout << "---------------------------------------------------------------------------------------------------------------------" << endl;
  //cout << "------------------------------------------- First fit ends here  ----------------------------------------------------" << endl;
  //cout << "---------------------------------------------------------------------------------------------------------------------" << endl;


  //Obtain set of parameters
  RooArgSet* params = model_d0.getParameters(combData);
  //Write parameters to file
  //params->writeToFile("params.txt");
   
  TString dir1 = gROOT->GetTutorialsDir();
  dir1.Append(fname);
  params->readFromFile(fname);
    
  //params->Print("v");
  RooFitResult *res = simPdf.fitTo(combData,Save(true),Strategy(2));
  if (res->status() == 0&&res->covQual() == 3){
    TString kk_d01, kk_d0bar1, d0pi_d01, d0pi_d0bar1;
    char *ext1=".pdf";
    char *ext2=".log";
    //TString parameters_after_fit =  "new_params_" + fname +ext2;  
    std::string parameters_after_fit  = std::string("new_params_") + fname + ext2;
    
    RooArgSet* params_new = model_d0.getParameters(combData);
    params_new->writeToFile(parameters_after_fit.c_str());
    
    std::string output_kk_d0  = std::string("kk_d01_") + fname + ext1;
    std::string output_kk_d0bar  = std::string("kk_d0bar1_") + fname + ext1;
    std::string output_d0pi_d0  = std::string("d0pi_d01_") + fname + ext1;
    std::string output_d0pi_d0bar  = std::string("d0pi_d0bar1_") + fname + ext1;
    
    /*-----------------------------------PLOT THE RESULTS-------------------------------------*/
      
      // Plot "D0" slice of simultaneous pdf for m(KK). 
      RooPlot *kk_d0 = mKK.frame(Title("D0"));
      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->SetYTitle("Candidates per 3 MeV/#it{c}^{2}");
      kk_d0->GetYaxis()->SetTitleOffset(1.5);
      
      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->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",500, 500);
      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_d0,signal_rnd_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
      simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(kpi_peak_d0,kpi_rnd_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
      simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(mult_peak_d0,mult_rnd_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
      simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(RooArgSet(ds_peak_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
      simPdf.plotOn(kk_d0,Slice(sample,"D0"),Components(comb_bkg_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
      
      kk_d0->Draw();
      upPad1->SetLogy(1);
      c1->cd(2);
      pullplotx_d0_mkk->Draw();
      c1->Update();
      c1->SaveAs(output_kk_d0.c_str());
      
      //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->SetYTitle("Candidates per 3 MeV/#it{c}^{2}");
      kk_d0bar->GetYaxis()->SetTitleOffset(1.5);
      
      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",500, 500);
      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_d0bar,signal_rnd_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
      simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(kpi_peak_d0bar,kpi_rnd_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
      simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(mult_peak_d0bar,mult_rnd_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
      simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(ds_peak_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
      simPdf.plotOn(kk_d0bar,Slice(sample,"D0bar"),Components(comb_bkg_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
      
      kk_d0bar->Draw();
      upPad2->SetLogy();
      c2->cd(2);
      pullplotx_d0bar_mkk->Draw();
      c2->Update();
      c2->SaveAs(output_kk_d0bar.c_str());
      
      // 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->SetYTitle("Candidates per 0.2 MeV/#it{c}^{2}");
      d0pi_d0->GetYaxis()->SetTitleOffset(1.5);
      
      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",500, 500);
      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_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
      simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(kpi_peak_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
      simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(mult_peak_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
      simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(ds_peak_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
      simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(comb_bkg_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
      simPdf.plotOn(d0pi_d0,Slice(sample,"D0"),Components(RooArgSet(signal_rnd_d0,kpi_rnd_d0,mult_rnd_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kOrange));
      d0pi_d0->Draw();
      upPad3->SetLogy();
      c3->cd(2);
      pullploty_d0_md0pi->Draw();
      c3->Update();
      c3->SaveAs(output_d0pi_d0.c_str());
      
      //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->SetYTitle("Candidates per 0.2 MeV/#it{c}^{2}");
      d0pi_d0bar->GetYaxis()->SetTitleOffset(1.5);
      
      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",500, 500);
      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_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
      simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(kpi_peak_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
      simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(mult_peak_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
      simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(ds_peak_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kCyan));
      simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(comb_bkg_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kMagenta));
      simPdf.plotOn(d0pi_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(signal_rnd_d0bar,kpi_rnd_d0bar,mult_rnd_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kOrange));
      d0pi_d0bar->Draw();
      upPad4->SetLogy();
      c4->cd(2);
      pullploty_d0bar_md0pi->Draw();
      c4->Update();
      c4->SaveAs(output_d0pi_d0bar.c_str());
      
      cout << "-----------------------------------------------------------------------------------------------------------------" << endl;
      cout << "-----------------------------------------------------------------------------------------------------------------" << endl;
    }
  }
}
}
//}
void asymmetry_all_params_batch()
{
  asymmetry_all_params_arg();
  return 0;
}