//#ifndef __CINT__
#include "RooGlobalFunc.h"
//#endif
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooFormulaVar.h"
#include "RooDataSet.h"
#include "RooGaussian.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 "TLegend.h"
#include "RooPlot.h"
#include "RooCategory.h"
#include "RooSuperCategory.h"
#include "RooSimultaneous.h"
#include "RooNLLVar.h"
#include "TFile.h"
#include "RooJohnson.h"
#include "RooKeysPdf.h"
#include "TRandom.h"
#include "myRooPolBG.cpp"

using namespace RooFit ;
using namespace std;

void mc15_asymmetry_fit_check()
{
  TFile *f = new TFile("component_wise_trees_for_fit.root");
  TTree *tree = (TTree*)f->Get("complete");
  RooRealVar log_min_significance("log_min_significance","", -6, 6.7, "");
  RooRealVar deltaM("deltaM","", 0.139, 0.159, "");

  //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_d0",tree,RooArgSet(log_min_significance,deltaM,pisoft),"pisoft==1");
  RooDataSet *dataxy_d0bar = new RooDataSet("data_d0bar","data_d0bar",tree,RooArgSet(log_min_significance,deltaM,pisoft),"pisoft==-1");  
  /*-----------------------S I G N A L----------------------*/
  //deltaM
  RooRealVar mean_sig("#mu_{sig}", "mean",0.1454197);
  RooRealVar A_mean_sig("A_mean_sig","Asymmetry in mean",0);
  RooFormulaVar mean_sig_d0("mean_sig_d0","mean","(@0*(1 + @1))",RooArgList(mean_sig, A_mean_sig));
  RooFormulaVar mean_sig_d0bar("mean_sig_d0bar","mean","(@0*(1- @1))",RooArgList(mean_sig, A_mean_sig));

  RooRealVar lambda_sig("#lambda_{sig}", "sigma", 0.0001581);
  RooRealVar A_lambda_sig("A_lambda_sig","Asymmetry in lambda",0);
  RooFormulaVar lambda_sig_d0("lambda_sig_d0","lambda","(@0*(1 + @1))",RooArgList(lambda_sig, A_lambda_sig));
  RooFormulaVar lambda_sig_d0bar("lambda_sig_d0bar","lambda","(@0*(1- @1))",RooArgList(lambda_sig, A_lambda_sig));

  RooRealVar delta_sig("#delta_{sig}", "delta", -0.0559);
  RooRealVar A_delta_sig("A_delta_sig","Asymmetry in delta",0);
  RooFormulaVar delta_sig_d0("delta_sig_d0","delta","(@0*(1 + @1))",RooArgList(delta_sig, A_delta_sig));
  RooFormulaVar delta_sig_d0bar("delta_sig_d0bar","delta","(@0*(1- @1))",RooArgList(delta_sig, A_delta_sig));

  RooRealVar gamma_sig("#gamma_{sig}", "gamma", 0.776);
  RooRealVar A_gamma_sig("A_gamma_sig","Asymmetry in gamma",0);
  RooFormulaVar gamma_sig_d0("gamma_sig_d0","gamma","(@0*(1 + @1))",RooArgList(gamma_sig, A_gamma_sig));
  RooFormulaVar gamma_sig_d0bar("gamma_sig_d0bar","gamma","(@0*(1- @1))",RooArgList(gamma_sig, A_gamma_sig));

  RooJohnson deltam_sig_d0("deltam_sig_d0", "signal component_d0", deltaM, mean_sig_d0, lambda_sig_d0, delta_sig_d0, gamma_sig_d0);
  RooJohnson deltam_sig_d0bar("deltam_sig_d0bar", "signal component_d0bar", deltaM, mean_sig_d0bar, lambda_sig_d0bar, delta_sig_d0bar, gamma_sig_d0bar);

  //log_min_significance
  RooRealVar gamma_mean_sig("#mu_{sig_gamma}", "mean",7.021);
  RooRealVar A_gamma_mean_sig("A_gamma_mean_sig","Asymmetry in gamma_mean_sig",0);
  RooFormulaVar gamma_mean_sig_d0("gamma_mean_sig_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_mean_sig, A_gamma_mean_sig));
  RooFormulaVar gamma_mean_sig_d0bar("gamma_mean_sig_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_mean_sig, A_gamma_mean_sig));

  RooRealVar gamma_sigma_sig("#lambda_{sig_gamma}", "sigma",0.301);
  RooRealVar A_gamma_sigma_sig("A_gamma_sigma_sig","Asymmetry in gamma_sigma_sig",0);
  RooFormulaVar gamma_sigma_sig_d0("gamma_sigma_sig_d0","gamma_sigma","(@0*(1 + @1))",RooArgList(gamma_sigma_sig, A_gamma_sigma_sig));
  RooFormulaVar gamma_sigma_sig_d0bar("gamma_sigma_sig_d0bar","gamma_sigma","(@0*(1- @1))",RooArgList(gamma_sigma_sig, A_gamma_sigma_sig));

  RooRealVar gamma_delta_sig("#delta_{sig_gamma}", "delta", 5);
  RooRealVar A_gamma_delta_sig("A_gamma_delta_sig","Asymmetry in gamma_delta_sig",0);
  RooFormulaVar gamma_delta_sig_d0("gamma_delta_sig_d0","gamma_delta","(@0*(1 + @1))",RooArgList(gamma_delta_sig, A_gamma_delta_sig));
  RooFormulaVar gamma_delta_sig_d0bar("gamma_delta_sig_d0bar","gamma_delta","(@0*(1- @1))",RooArgList(gamma_delta_sig, A_gamma_delta_sig));

  RooRealVar gamma_gamma_sig("#gamma_{sig_gamma}", "gamma", 1.954);
  RooRealVar A_gamma_gamma_sig("A_gamma_gamma_sig","Asymmetry in gamma_gamma_sig",0);
  RooFormulaVar gamma_gamma_sig_d0("gamma_gamma_sig_d0","gamma_gamma","(@0*(1 + @1))",RooArgList(gamma_gamma_sig, A_gamma_gamma_sig));
  RooFormulaVar gamma_gamma_sig_d0bar("gamma_gamma_sig_d0bar","gamma_gamma","(@0*(1- @1))",RooArgList(gamma_gamma_sig, A_gamma_gamma_sig));

  RooJohnson gammapdf_sig_d0("gammapdf_sig_d0", "signal component_d0", log_min_significance, gamma_mean_sig_d0, gamma_sigma_sig_d0, gamma_delta_sig, gamma_gamma_sig_d0);
  RooJohnson gammapdf_sig_d0bar("gammapdf_sig_d0bar", "signal component_d0bar", log_min_significance, gamma_mean_sig_d0bar, gamma_sigma_sig_d0bar, gamma_delta_sig_d0bar, gamma_gamma_sig_d0bar);
  
  //total signal component
  RooProdPdf signal_peak_d0("signal_peak_d0","signal_peak_d0", RooArgSet(deltam_sig_d0,gammapdf_sig_d0));
  RooProdPdf signal_peak_d0bar("signal_peak_d0bar","signal_peak_d0bar", RooArgSet(deltam_sig_d0bar,gammapdf_sig_d0bar));
  
  // /*-----------------------P E A K I N G     B A C K G R O U N D     (KS Pi Pi)-------------------------*/
  
  //deltaM
  RooRealVar f_G_kspipi("f_{G_kspipi}", "",0.234);
  RooRealVar A_f_G_kspipi("A_f_G_kspipi","Asymmetry in f_G_kspipi",0);
  RooFormulaVar f_G_kspipi_d0("f_G_kspipi_d0","f_G_kspipi","(@0*(1 + @1))",RooArgList(f_G_kspipi, A_f_G_kspipi));
  RooFormulaVar f_G_kspipi_d0bar("f_G_kspipi_d0bar","f_G_kspipi","(@0*(1- @1))",RooArgList(f_G_kspipi, A_f_G_kspipi));

  RooRealVar mean_G_kspipi("#mu_{G_kspipi}", "mean",0.1454229);
  RooRealVar A_mean_G_kspipi("A_mean_G_kspipi","Asymmetry in mean_G_kspipi",0);
  RooFormulaVar mean_G_kspipi_d0("mean_G_kspipi_d0","mean_G_kspipi","(@0*(1 + @1))",RooArgList(mean_G_kspipi, A_mean_G_kspipi));
  RooFormulaVar mean_G_kspipi_d0bar("mean_G_kspipi_d0bar","mean_G_kspipi","(@0*(1- @1))",RooArgList(mean_G_kspipi, A_mean_G_kspipi));

  RooRealVar sigma_G_kspipi("#sigma_{G_kspipi}", "sigma", 0.000185);
  RooRealVar A_sigma_G_kspipi("A_sigma_G_kspipi","Asymmetry in sigma_G_kspipi",0);
  RooFormulaVar sigma_G_kspipi_d0("sigma_G_kspipi_d0","sigma_G_kspipi","(@0*(1 + @1))",RooArgList(sigma_G_kspipi, A_sigma_G_kspipi));
  RooFormulaVar sigma_G_kspipi_d0bar("sigma_G_kspipi_d0bar","sigma_G_kspipi","(@0*(1- @1))",RooArgList(sigma_G_kspipi, A_sigma_G_kspipi));
  
  RooGaussian gauss_kspipi_d0("gauss_kspipi_d0", "gauss component_d0", deltaM, mean_G_kspipi_d0, sigma_G_kspipi_d0);
  RooGaussian gauss_kspipi_d0bar("gauss_kspipi_d0bar", "gauss component_d0bar", deltaM, mean_G_kspipi_d0bar, sigma_G_kspipi_d0bar);
 
  RooRealVar mean_J_kspipi("#mu_{J_kspipi}", "mean",0.1454229);
  RooRealVar A_mean_J_kspipi("A_mean_J_kspipi","Asymmetry in mean_J_kspipi",0);
  RooFormulaVar mean_J_kspipi_d0("mean_J_kspipi_d0","mean_J_kspipi","(@0*(1 + @1))",RooArgList(mean_J_kspipi, A_mean_J_kspipi));
  RooFormulaVar mean_J_kspipi_d0bar("mean_J_kspipi_d0bar","mean_J_kspipi","(@0*(1- @1))",RooArgList(mean_J_kspipi, A_mean_J_kspipi));

  RooRealVar lambda_J_kspipi("#lambda_{J_kspipi}", "lambda", 0.0001408);
  RooRealVar A_lambda_J_kspipi("A_lambda_J_kspipi","Asymmetry in lambda_J_kspipi",0);
  RooFormulaVar lambda_J_kspipi_d0("lambda_J_kspipi_d0","lambda_J_kspipi","(@0*(1 + @1))",RooArgList(lambda_J_kspipi, A_lambda_J_kspipi));
  RooFormulaVar lambda_J_kspipi_d0bar("lambda_J_kspipi_d0bar","lambda_J_kspipi","(@0*(1- @1))",RooArgList(lambda_J_kspipi, A_lambda_J_kspipi));

  RooRealVar delta_J_kspipi("#delta_{J_kspipi}", "delta", -0.0404);
  RooRealVar A_delta_J_kspipi("A_delta_J_kspipi","Asymmetry in delta_J_kspipi",0);
  RooFormulaVar delta_J_kspipi_d0("delta_J_kspipi_d0","delta_J_kspipi","(@0*(1 + @1))",RooArgList(delta_J_kspipi, A_delta_J_kspipi));
  RooFormulaVar delta_J_kspipi_d0bar("delta_J_kspipi_d0bar","delta_J_kspipi","(@0*(1- @1))",RooArgList(delta_J_kspipi, A_delta_J_kspipi));

  RooRealVar gamma_J_kspipi("#gamma_{J_kspipi}", "gamma", 0.684);
  RooRealVar A_gamma_J_kspipi("A_gamma_J_kspipi","Asymmetry in gamma_J_kspipi",0);
  RooFormulaVar gamma_J_kspipi_d0("gamma_J_kspipi_d0","gamma_J_kspipi","(@0*(1 + @1))",RooArgList(gamma_J_kspipi, A_gamma_J_kspipi));
  RooFormulaVar gamma_J_kspipi_d0bar("gamma_J_kspipi_d0bar","gamma_J_kspipi","(@0*(1- @1))",RooArgList(gamma_J_kspipi, A_gamma_J_kspipi));

  RooJohnson johnson_kspipi_d0("johnson_kspipi_d0", "johnson component_d0", deltaM, mean_G_kspipi_d0, lambda_J_kspipi_d0, delta_J_kspipi_d0, gamma_J_kspipi_d0);
  RooJohnson johnson_kspipi_d0bar("johnson_kspipi_d0bar", "johnson component_d0bar", deltaM, mean_G_kspipi_d0bar, lambda_J_kspipi_d0bar, delta_J_kspipi_d0bar, gamma_J_kspipi_d0bar);

  RooAddPdf deltam_kspipi_d0("deltam_kspipi_d0" ,"gauss_kspipi_d0 + johnson_kspipi_d0",gauss_kspipi_d0, johnson_kspipi_d0,f_G_kspipi_d0);
  RooAddPdf deltam_kspipi_d0bar("deltam_kspipi_d0bar" ,"gauss_kspipi_d0bar + johnson_kspipi_d0bar",gauss_kspipi_d0bar, johnson_kspipi_d0bar,f_G_kspipi_d0bar);
  
  //log_min_significance
  RooRealVar gamma_mean_kspipi("#mu_{kspipi_gamma}", "mean",2.242);
  RooRealVar A_gamma_mean_kspipi("A_gamma_mean_kspipi","Asymmetry in gamma_mean_kspipi",0);
  RooFormulaVar gamma_mean_kspipi_d0("gamma_mean_kspipi_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_mean_kspipi, A_gamma_mean_kspipi));
  RooFormulaVar gamma_mean_kspipi_d0bar("gamma_mean_kspipi_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_mean_kspipi, A_gamma_mean_kspipi));

  RooRealVar gamma_sigma_kspipi("#lambda_{kspipi_gamma}", "sigma",1.390);
  RooRealVar A_gamma_sigma_kspipi("A_gamma_sigma_kspipi","Asymmetry in gamma_sigma_kspipi",0);
  RooFormulaVar gamma_sigma_kspipi_d0("gamma_sigma_kspipi_d0","gamma_sigma","(@0*(1 + @1))",RooArgList(gamma_sigma_kspipi, A_gamma_sigma_kspipi));
  RooFormulaVar gamma_sigma_kspipi_d0bar("gamma_sigma_kspipi_d0bar","gamma_sigma","(@0*(1- @1))",RooArgList(gamma_sigma_kspipi, A_gamma_sigma_kspipi));

  RooRealVar gamma_delta_kspipi("#delta_{kspipi_gamma}", "delta", 1.790);
  RooRealVar A_gamma_delta_kspipi("A_gamma_delta_kspipi","Asymmetry in gamma_delta_kspipi",0);
  RooFormulaVar gamma_delta_kspipi_d0("gamma_delta_kspipi_d0","gamma_delta","(@0*(1 + @1))",RooArgList(gamma_delta_kspipi, A_gamma_delta_kspipi));
  RooFormulaVar gamma_delta_kspipi_d0bar("gamma_delta_kspipi_d0bar","gamma_delta","(@0*(1- @1))",RooArgList(gamma_delta_kspipi, A_gamma_delta_kspipi));

  RooRealVar gamma_gamma_kspipi("#gamma_{kspipi_gamma}", "gamma", 1.926);
  RooRealVar A_gamma_gamma_kspipi("A_gamma_gamma_kspipi","Asymmetry in gamma_gamma_kspipi",0);
  RooFormulaVar gamma_gamma_kspipi_d0("gamma_gamma_kspipi_d0","gamma_gamma","(@0*(1 + @1))",RooArgList(gamma_gamma_kspipi, A_gamma_gamma_kspipi));
  RooFormulaVar gamma_gamma_kspipi_d0bar("gamma_gamma_kspipi_d0bar","gamma_gamma","(@0*(1- @1))",RooArgList(gamma_gamma_kspipi, A_gamma_gamma_kspipi));

  RooJohnson gamma_kspipi_d0("gamma_kspipi_d0", "signal component_d0", log_min_significance, gamma_mean_kspipi_d0, gamma_sigma_kspipi_d0, gamma_delta_kspipi, gamma_gamma_kspipi_d0);
  RooJohnson gamma_kspipi_d0bar("gamma_kspipi_d0bar", "signal component_d0bar", log_min_significance, gamma_mean_kspipi_d0bar, gamma_sigma_kspipi_d0bar, gamma_delta_kspipi_d0bar, gamma_gamma_kspipi_d0bar);   
  
  //kspipi component (peaking background)
  RooProdPdf kspipi_peak_d0("kspipi_peak_d0", "kspipi_peak_d0", deltam_kspipi_d0, gamma_kspipi_d0);
  RooProdPdf kspipi_peak_d0bar("kspipi_peak_d0bar", "kspipi_peak_d0bar", deltam_kspipi_d0bar, gamma_kspipi_d0bar);
  
  
  /*----------------------- N O N     P E A K I N G     B A C K G R O U N D-------------------------*/
  
  //deltaM
  RooRealVar mpi_non_peaking("mpi_non_peaking", "nominal pion mass", 0.13957039);
  RooRealVar alpha_non_peaking("#alpha_{non_peaking}", "alpha", -152, -1e3, 1e3);
  RooRealVar A_alpha_non_peaking("A_alpha_non_peaking","Asymmetry in alpha_non_peaking",0);
  RooFormulaVar alpha_non_peaking_d0("alpha_non_peaking_d0","alpha_non_peaking_d0","(@0*(1 + @1))",RooArgList(alpha_non_peaking, A_alpha_non_peaking));
  RooFormulaVar alpha_non_peaking_d0bar("alpha_non_peaking_d0bar","alpha_non_peaking_d0bar","(@0*(1- @1))",RooArgList(alpha_non_peaking, A_alpha_non_peaking));
  RooRealVar beta_non_peaking("#beta_{non_peaking}", "beta", 0);
  
  myRooPolBG deltam_non_peaking_d0("deltam_non_peaking_d0","deltam_non_peaking_d0", deltaM, mpi_non_peaking, alpha_non_peaking_d0, beta_non_peaking);
  myRooPolBG deltam_non_peaking_d0bar("deltam_non_peaking_d0bar","deltam_non_peaking_d0bar", deltaM, mpi_non_peaking, alpha_non_peaking_d0bar, beta_non_peaking);

  //log_min_significance
  RooRealVar f_J_non_peaking("fJ_{non_peaking}", "frac",0.691);
  //RooRealVar f_J_non_peaking("fJ_{non_peaking}", "frac",0.3282);
  RooRealVar A_fJ_non_peaking(" A_fJ_non_peaking","Asymmetry in gamma_mean_non_peaking",0);
  RooFormulaVar f_J_non_peaking_d0("f_J_non_peaking_d0","gamma_mean","(@0*(1 + @1))",RooArgList(f_J_non_peaking, A_fJ_non_peaking));
  RooFormulaVar f_J_non_peaking_d0bar("f_J_non_peaking_d0","gamma_mean","(@0*(1- @1))",RooArgList(f_J_non_peaking, A_fJ_non_peaking));

  RooRealVar gamma_meanJ_non_peaking("#muJ_{non_peaking_gamma}", "mean",0.86);
  RooRealVar A_gamma_meanJ_non_peaking("A_gamma_meanJ_non_peaking","Asymmetry in gamma_mean_non_peaking",0);
  RooFormulaVar gamma_meanJ_non_peaking_d0("gamma_meanJ_non_peaking_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_meanJ_non_peaking, A_gamma_meanJ_non_peaking));
  RooFormulaVar gamma_meanJ_non_peaking_d0bar("gamma_meanJ_non_peaking_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_meanJ_non_peaking, A_gamma_meanJ_non_peaking));

  RooRealVar gamma_sigmaJ_non_peaking("#lambdaJ_{non_peaking_gamma}", "sigma", 1.81);
  RooRealVar A_gamma_sigmaJ_non_peaking("A_gamma_sigmaJ_non_peaking","Asymmetry in gamma_sigma_non_peaking",0);
  RooFormulaVar gamma_sigmaJ_non_peaking_d0("gamma_sigmaJ_non_peaking_d0","gamma_sigma","(@0*(1 + @1))",RooArgList(gamma_sigmaJ_non_peaking, A_gamma_sigmaJ_non_peaking));
  RooFormulaVar gamma_sigmaJ_non_peaking_d0bar("gamma_sigmaJ_non_peaking_d0bar","gamma_sigma","(@0*(1- @1))",RooArgList(gamma_sigmaJ_non_peaking, A_gamma_sigmaJ_non_peaking));

  RooRealVar gamma_deltaJ_non_peaking("#deltaJ_{sig_gamma}", "delta", 0.50);
  RooRealVar A_gamma_deltaJ_non_peaking("A_gamma_deltaJ_non_peaking","Asymmetry in gamma_delta_non_peaking",0);
  RooFormulaVar gamma_deltaJ_non_peaking_d0("gamma_deltaJ_non_peaking_d0","gamma_delta","(@0*(1 + @1))",RooArgList(gamma_deltaJ_non_peaking, A_gamma_deltaJ_non_peaking));
  RooFormulaVar gamma_deltaJ_non_peaking_d0bar("gamma_deltaJ_non_peaking_d0bar","gamma_delta","(@0*(1- @1))",RooArgList(gamma_deltaJ_non_peaking, A_gamma_deltaJ_non_peaking));

  RooRealVar gamma_gammaJ_non_peaking("#gammaJ_{non_peaking_gamma}", "gamma", 1.68);
  RooRealVar A_gamma_gammaJ_non_peaking("A_gamma_gammaJ_non_peaking","Asymmetry in gamma_gamma_non_peaking",0);
  RooFormulaVar gamma_gammaJ_non_peaking_d0("gamma_gammaJ_non_peaking_d0","gamma_gamma","(@0*(1 + @1))",RooArgList(gamma_gammaJ_non_peaking, A_gamma_gammaJ_non_peaking));
  RooFormulaVar gamma_gammaJ_non_peaking_d0bar("gamma_gammaJ_non_peaking_d0bar","gamma_gamma","(@0*(1- @1))",RooArgList(gamma_gammaJ_non_peaking, A_gamma_gammaJ_non_peaking));

  RooJohnson gammaJ_non_peaking_d0("gammaJ_non_peaking_d0", "signal component_d0", log_min_significance, gamma_meanJ_non_peaking_d0, gamma_sigmaJ_non_peaking_d0, gamma_deltaJ_non_peaking, gamma_gammaJ_non_peaking_d0);
  RooJohnson gammaJ_non_peaking_d0bar("gammaJ_non_peaking_d0bar", "signal component_d0bar", log_min_significance, gamma_meanJ_non_peaking_d0bar, gamma_sigmaJ_non_peaking_d0bar, gamma_deltaJ_non_peaking_d0bar, gamma_gammaJ_non_peaking_d0bar);   
  
  ///////////////////////// 
  //RooRealVar gamma_meanG_non_peaking_shift("#muG_{non_peaking_gamma_shift}", "mean",6.632);
  RooRealVar gamma_meanG_non_peaking_shift("#muG_{non_peaking_gamma_shift}", "mean",6.69);
  //RooFormulaVar gamma_meanG_non_peaking("#mu_{sig_gamma}","(@0 + @1)",RooArgList(gamma_meanG_non_peaking_shift,gamma_meanJ_non_peaking));
  RooRealVar A_gamma_meanG_non_peaking("A_gamma_meanG_non_peaking","Asymmetry in gamma_mean_non_peaking",0);
  RooFormulaVar gamma_meanG_non_peaking_shift_d0("gamma_meanG_non_peaking_shift_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_meanG_non_peaking_shift, A_gamma_meanG_non_peaking));
  RooFormulaVar gamma_meanG_non_peaking_shift_d0bar("gamma_meanG_non_peaking_shift_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_meanG_non_peaking_shift, A_gamma_meanG_non_peaking));

  //RooRealVar gamma_sigmaG_non_peaking_shift("lambdaG_{non_peaking_gamma_shift}", "mean",0.016);
  RooRealVar gamma_sigmaG_non_peaking_shift("lambdaG_{non_peaking_gamma_shift}", "mean",0.87);
  //RooFormulaVar gamma_sigmaG_non_peaking("#lambda_{sig_gamma}","(@0 + @1)",RooArgList(gamma_sigmaG_non_peaking_shift,gamma_sigmaJ_non_peaking));
  RooRealVar A_gamma_sigmaG_non_peaking("A_gamma_sigmaG_non_peaking","Asymmetry in gamma_mean_non_peaking",0);
  RooFormulaVar gamma_sigmaG_non_peaking_shift_d0("gamma_sigmaG_non_peaking_shift_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_sigmaG_non_peaking_shift, A_gamma_sigmaG_non_peaking));
  RooFormulaVar gamma_sigmaG_non_peaking_shift_d0bar("gamma_sigmaG_non_peaking_shift_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_sigmaG_non_peaking_shift, A_gamma_sigmaG_non_peaking));

  //RooRealVar gamma_deltaG_non_peaking_shift("deltaG_{non_peaking_gamma_shift}", "mean",7.5);
  RooRealVar gamma_deltaG_non_peaking_shift("deltaG_{non_peaking_gamma_shift}", "mean",2.64);
  //RooFormulaVar gamma_deltaG_non_peaking("#delta_{sig_gamma}","(@0 + @1)",RooArgList(gamma_deltaG_non_peaking_shift,gamma_deltaJ_non_peaking));
  RooRealVar A_gamma_deltaG_non_peaking("A_delta_sigmaG_non_peaking","Asymmetry in gamma_mean_non_peaking",0);
  RooFormulaVar gamma_deltaG_non_peaking_shift_d0("gamma_deltaG_non_peaking_shift_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_deltaG_non_peaking_shift, A_gamma_deltaG_non_peaking));
  RooFormulaVar gamma_deltaG_non_peaking_shift_d0bar("gamma_deltaG_non_peaking_shift_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_deltaG_non_peaking_shift, A_gamma_deltaG_non_peaking));

  RooRealVar gamma_gammaG_non_peaking_shift("gammaG_{non_peaking_gamma_shift}", "mean",1.61);
  //RooRealVar gamma_gammaG_non_peaking_shift("gammaG_{non_peaking_gamma_shift}", "mean",1.329);
  //RooFormulaVar gamma_gammaG_non_peaking("#gamma_{sig_gamma}","(@0 + @1)",RooArgList(gamma_gammaG_non_peaking_shift,gamma_gammaJ_non_peaking));
  RooRealVar A_gamma_gammaG_non_peaking("A_gamma_gammaG_non_peaking","Asymmetry in gamma_mean_non_peaking",0);
  RooFormulaVar gamma_gammaG_non_peaking_shift_d0("gamma_gammaG_non_peaking_shift_d0","gamma_mean","(@0*(1 + @1))",RooArgList(gamma_gammaG_non_peaking_shift, A_gamma_gammaG_non_peaking));
  RooFormulaVar gamma_gammaG_non_peaking_shift_d0bar("gamma_gammaG_non_peaking_shift_d0bar","gamma_mean","(@0*(1- @1))",RooArgList(gamma_gammaG_non_peaking_shift, A_gamma_gammaG_non_peaking));
 
  RooJohnson gammaG_non_peaking_d0("gammaG_non_peaking_d0", "signal component_d0", log_min_significance, gamma_meanG_non_peaking_shift_d0, gamma_sigmaG_non_peaking_shift_d0, gamma_deltaG_non_peaking_shift_d0, gamma_gammaG_non_peaking_shift_d0);
  RooJohnson gammaG_non_peaking_d0bar("gammaG_non_peaking_d0bar", "signal component_d0bar", log_min_significance, gamma_meanG_non_peaking_shift_d0bar, gamma_sigmaG_non_peaking_shift_d0bar, gamma_deltaG_non_peaking_shift_d0,gamma_gammaG_non_peaking_shift_d0);   
  
  RooAddPdf gamma_non_peaking_d0("gamma_non_peaking_d0", "gammaJ_non_peaking_d0 + gammaG_non_peaking_d0",gammaJ_non_peaking_d0,gammaG_non_peaking_d0,f_J_non_peaking_d0);
  RooAddPdf gamma_non_peaking_d0bar("gamma_non_peaking_d0bar", "gammaJ_non_peaking_d0bar + gammaG_non_peaking_d0bar",gammaJ_non_peaking_d0,gammaG_non_peaking_d0,f_J_non_peaking_d0bar); 

  //non_peaking component
  RooProdPdf non_peaking_d0("non_peaking_d0", "non_peaking_d0", deltam_non_peaking_d0, gamma_non_peaking_d0);
  RooProdPdf non_peaking_d0bar("non_peaking_d0bar", "non_peaking_d0bar", deltam_non_peaking_d0bar, gamma_non_peaking_d0bar);
  
  // /*------------------------Y I E L D S  F O R  A L L  C O M P O N E N T S----------------------------*/
  RooRealVar signal_yield("N_{signal}","signal yield",2000,1e0,1e4);
  RooRealVar kspipi_yield("N_{kspipi}","kspipi yield",5000,1e0,1e5);
  RooRealVar nonpeaking_yield("N_{nonpeaking}","nonpeaking yield",5000,1e0,1e5);

 //For Blinding signal asymmetry
  TRandom3 *r=new TRandom3(1994);
  Double_t shift = r->Uniform(-0.1,0.1); 
  cout << "random = "<< shift << endl; 
  RooRealVar rnd("rnd","",shift,"");

  RooRealVar Araw_sig("A_{sig}","Araw",0.01,-1,1);
  //RooFormulaVar Araw_sig_blind("A_{sig_blind}","(@0 + @1)",RooArgList(rnd,Araw_sig));
  RooRealVar Araw_kspipi("A_{kspipi}","Araw",0,-1,1);
  RooRealVar Araw_nonpeaking("A_{non_peaking}","Araw",0,-1,1);

  //Yields for D0 sample: ND0 = N(1+A)/2
  RooFormulaVar signal_yield_d0("ND0_{signal}","(@0*(1 + @1))/2",RooArgList(signal_yield,Araw_sig));
  RooFormulaVar kspipi_yield_d0("ND0_{kspipi}","(@0*(1 + @1))/2",RooArgList(kspipi_yield,Araw_kspipi));
  RooFormulaVar nonpeaking_yield_d0("ND0_{nonpeaking}","(@0*(1 + @1))/2",RooArgList(nonpeaking_yield,Araw_nonpeaking));

  //Yields for D0 sample: ND0 = N(1-A)/2
  RooFormulaVar signal_yield_d0bar("ND0bar_{signal}","(@0*(1 - @1))/2",RooArgList(signal_yield,Araw_sig));
  RooFormulaVar kspipi_yield_d0bar("ND0bar_{kspipi}","(@0*(1 - @1))/2",RooArgList(kspipi_yield,Araw_kspipi));
  RooFormulaVar nonpeaking_yield_d0bar("ND0bar_{nonpeaking}","(@0*(1 - @1))/2",RooArgList(nonpeaking_yield,Araw_nonpeaking));

  
  /*------------------------C O M P L E T E  M O D E L----------------------------*/
  RooAddPdf model_d0("model_d0","",RooArgList(signal_peak_d0,kspipi_peak_d0,non_peaking_d0), RooArgList(signal_yield_d0,kspipi_yield_d0,nonpeaking_yield_d0));
  RooAddPdf model_d0bar("model_d0bar","",RooArgList(signal_peak_d0bar,kspipi_peak_d0bar,non_peaking_d0bar), RooArgList(signal_yield_d0bar,kspipi_yield_d0bar,nonpeaking_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(log_min_significance,deltaM),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 deltaM 
  RooPlot *deltaM_d0 = deltaM.frame(Title("D0"));
  combData.plotOn(deltaM_d0,Cut("sample==sample::D0"));
  simPdf.plotOn(deltaM_d0,Slice(sample,"D0"),ProjWData(sample,combData));
  simPdf.paramOn(deltaM_d0,dataxy_d0);
  deltaM_d0->SetXTitle("#it{#Deltam} [GeV/c^{2}]");
  deltaM_d0->SetYTitle("Candidates per 0.2 MeV/#it{c}^{2}");
  deltaM_d0->GetYaxis()->SetTitleOffset(1.5);
 
  RooHist* hpullx_d0_deltaM = deltaM_d0->pullHist();
  hpullx_d0_deltaM->SetFillStyle(1001);
  hpullx_d0_deltaM->SetFillColor(1);
  for(int i=0;i<hpullx_d0_deltaM->GetN();++i) hpullx_d0_deltaM->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullplotx_d0_deltaM = deltaM.frame(Title(" ")) ;
  pullplotx_d0_deltaM->addPlotable(hpullx_d0_deltaM,"B");
  pullplotx_d0_deltaM->SetXTitle("#it{#Deltam} [GeV/c^{2}]");
  pullplotx_d0_deltaM->SetMinimum(-4.);
  pullplotx_d0_deltaM->SetMaximum(4.);
  pullplotx_d0_deltaM->GetXaxis()->SetLabelSize(0.1);
  pullplotx_d0_deltaM->GetXaxis()->SetTitleSize(0.1);
  pullplotx_d0_deltaM->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(deltaM_d0,Slice(sample,"D0"),Components(RooArgSet(signal_peak_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(deltaM_d0,Slice(sample,"D0"),Components(RooArgSet(kspipi_peak_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(deltaM_d0,Slice(sample,"D0"),Components(RooArgSet(non_peaking_d0)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  
  deltaM_d0->Draw();

  TLegend *leg_deltaM_d0 = new TLegend(0.54,0.32,0.84,0.52);
  leg_deltaM_d0->SetBorderSize(0);  // no border
  leg_deltaM_d0->AddEntry(deltaM_d0->findObject("deltaM"),"Data","pe");
  leg_deltaM_d0->AddEntry(deltaM_d0->findObject("simpdf"),"Total Fit","l");
  leg_deltaM_d0->AddEntry(deltaM_d0->findObject("signal_peak_d0"),"Signal","l");
  leg_deltaM_d0->AddEntry(deltaM_d0->findObject("kspipi_peak_d0"),"K_{s}#pi^{+}#pi^{-} background","l");
  leg_deltaM_d0->AddEntry(deltaM_d0->findObject("non_peaking_d0"),"Other backgrounds","l");
  leg_deltaM_d0->Draw();
  c1->cd(2);
  pullplotx_d0_deltaM->Draw();
  c1->Update();

  // Plot "D0bar" slice of simultaneous pdf for deltaM 
  RooPlot *deltaM_d0bar = deltaM.frame(Title("D0bar"));
  combData.plotOn(deltaM_d0bar,Cut("sample==sample::D0bar"));
  simPdf.plotOn(deltaM_d0bar,Slice(sample,"D0bar"),ProjWData(sample,combData));
  deltaM_d0bar->SetXTitle("#it{#Deltam} [GeV/c^{2}]");
  deltaM_d0bar->SetYTitle("Candidates per 0.2 MeV/#it{c}^{2}");
  deltaM_d0bar->GetYaxis()->SetTitleOffset(1.5);
 
  RooHist* hpullx_d0bar_deltaM = deltaM_d0bar->pullHist();
  hpullx_d0bar_deltaM->SetFillStyle(1001);
  hpullx_d0bar_deltaM->SetFillColor(1);
  for(int i=0;i<hpullx_d0bar_deltaM->GetN();++i) hpullx_d0bar_deltaM->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullplotx_d0bar_deltaM = deltaM.frame(Title(" ")) ;
  pullplotx_d0bar_deltaM->addPlotable(hpullx_d0bar_deltaM,"B");
  pullplotx_d0bar_deltaM->SetXTitle("#it{#Deltam} [GeV/c^{2}]");
  pullplotx_d0bar_deltaM->SetMinimum(-4.);
  pullplotx_d0bar_deltaM->SetMaximum(4.);
  pullplotx_d0bar_deltaM->GetXaxis()->SetLabelSize(0.1);
  pullplotx_d0bar_deltaM->GetXaxis()->SetTitleSize(0.1);
  pullplotx_d0bar_deltaM->GetYaxis()->SetLabelSize(0.07);
  
  TCanvas* c2 = new TCanvas("c2", "c2",700, 700);
  //Double_t xlow, ylow, xup, yup;
  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(deltaM_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(signal_peak_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(deltaM_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(kspipi_peak_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(deltaM_d0bar,Slice(sample,"D0bar"),Components(RooArgSet(non_peaking_d0bar)),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  
  deltaM_d0bar->Draw();
  c2->cd(2);
  pullplotx_d0bar_deltaM->Draw();
  c2->Update();
//////////////////////////////////////////////////////////////////////////////////////////////

// Plot "D0" slice of simultaneous pdf for max_significance.
  RooPlot *log_min_significance_d0 = log_min_significance.frame(Title("D0"));
  combData.plotOn(log_min_significance_d0,Cut("sample==sample::D0"));
  simPdf.plotOn(log_min_significance_d0,Slice(sample,"D0"),ProjWData(sample,combData));
  //simPdf.paramOn(max_significance_d0,dataxy_d0);
  log_min_significance_d0->SetXTitle("#gamma");
  log_min_significance_d0->SetYTitle("Candidates per 0.127");
  log_min_significance_d0->GetYaxis()->SetTitleOffset(1.5);
 
  RooHist* hpully_d0_log_min_significance = log_min_significance_d0->pullHist();
  hpully_d0_log_min_significance->SetFillStyle(1001);
  hpully_d0_log_min_significance->SetFillColor(1);
  for(int i=0;i<hpully_d0_log_min_significance->GetN();++i) hpully_d0_log_min_significance->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullploty_d0_log_min_significance = log_min_significance.frame(Title(" ")) ;
  pullploty_d0_log_min_significance->addPlotable(hpully_d0_log_min_significance,"B");
  pullploty_d0_log_min_significance->SetXTitle("#gamma");
  pullploty_d0_log_min_significance->SetMinimum(-4.);
  pullploty_d0_log_min_significance->SetMaximum(4.);
  pullploty_d0_log_min_significance->GetXaxis()->SetLabelSize(0.1);
  pullploty_d0_log_min_significance->GetXaxis()->SetTitleSize(0.1);
  pullploty_d0_log_min_significance->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(log_min_significance_d0,Slice(sample,"D0"),Components(signal_peak_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(log_min_significance_d0,Slice(sample,"D0"),Components(kspipi_peak_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(log_min_significance_d0,Slice(sample,"D0"),Components(non_peaking_d0),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  
  log_min_significance_d0->Draw();
  c3->cd(2);
  pullploty_d0_log_min_significance->Draw();
  c3->Update();

  // Plot "D0bar" slice of simultaneous pdf for max_significance.
  RooPlot *log_min_significance_d0bar = log_min_significance.frame(Title("D0bar"));
  combData.plotOn(log_min_significance_d0bar,Cut("sample==sample::D0bar"));
  simPdf.plotOn(log_min_significance_d0bar,Slice(sample,"D0bar"),ProjWData(sample,combData));
  simPdf.paramOn(log_min_significance_d0bar,dataxy_d0bar);
  log_min_significance_d0bar->SetXTitle("#gamma");
  log_min_significance_d0bar->SetYTitle("Candidates per 0.127");
  log_min_significance_d0bar->GetYaxis()->SetTitleOffset(1.5);
 
  RooHist* hpully_d0bar_log_min_significance = log_min_significance_d0bar->pullHist();
  hpully_d0bar_log_min_significance->SetFillStyle(1001);
  hpully_d0bar_log_min_significance->SetFillColor(1);
  for(int i=0;i<hpully_d0bar_log_min_significance->GetN();++i) hpully_d0bar_log_min_significance->SetPointError(i,0.0,0.0,0.0,0.0);
  RooPlot* pullploty_d0bar_log_min_significance = log_min_significance.frame(Title(" ")) ;
  pullploty_d0bar_log_min_significance->addPlotable(hpully_d0bar_log_min_significance,"B");
  pullploty_d0bar_log_min_significance->SetXTitle("#gamma");
  pullploty_d0bar_log_min_significance->SetMinimum(-4.);
  pullploty_d0bar_log_min_significance->SetMaximum(4.);
  pullploty_d0bar_log_min_significance->GetXaxis()->SetLabelSize(0.1);
  pullploty_d0bar_log_min_significance->GetXaxis()->SetTitleSize(0.1);
  pullploty_d0bar_log_min_significance->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(log_min_significance_d0bar,Slice(sample,"D0bar"),Components(signal_peak_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kRed)) ;
  simPdf.plotOn(log_min_significance_d0bar,Slice(sample,"D0bar"),Components(kspipi_peak_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kGreen));
  simPdf.plotOn(log_min_significance_d0bar,Slice(sample,"D0bar"),Components(non_peaking_d0bar),ProjWData(sample,combData),LineStyle(kDashed),LineColor(kBlack));
  
  log_min_significance_d0bar->Draw();
  c4->cd(2);
  pullploty_d0bar_log_min_significance->Draw();
  c4->Update();

}