#include <iostream>
#include <vector>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <sstream>

#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include <TROOT.h>
#include <TH1.h>
#include "TSystem.h"
#include "TString.h"
#include "TH2.h"

#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooAbsReal.h"
#include "RooAbsPdf.h"
#include "RooAddPdf.h"
#include "RooGaussian.h"
#include "RooUniform.h"
#include "RooExponential.h"
#include "RooKeysPdf.h"
#include "RooPolynomial.h"
#include "RooChebychev.h"
#include "RooUnblindUniform.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "RooWorkspace.h"
#include "RooMCStudy.h"

using namespace std;
using namespace RooFit;

RooFitResult* DataFit(RooDataSet* data_KPi_SS, 
                      RooDataSet* data_KPi_OS,
                      RooDataSet* data_KK,
                      RooDataSet* data_PiPi)
{
  gSystem->Load("libRooStats.so");

  //define some vars needed
  RooRealVar B_mass("Bd_REFITTED_M","Bd_REFITTED_M",5100,5600,"MeV/c^{2}");
  RooRealVar KstarK_ID("KstarK_ID","KstarK_ID",-322,322,"");
  RooRealVar D0K_ID("D0K_ID","D0K_ID",-322,322,"");
  B_mass.SetTitle("M(D^{0}K^{*0})");

  // GETTING VARIOUS BACKGROUND PDFS
  //--------------------------------
  //BdDRho_D02KPi
  TFile *_file0 = new TFile("~/BdDKstar_analysis/fitters/MCbackground/BdDRho_D02KPi/"
                            "BdDRho_D02KPi_fitresults.root");
  RooWorkspace *work_BdDRho = (RooWorkspace*)_file0->Get("BdDRho_D02KPi_BKG");
  RooArgSet BdDRho_pdfs = work_BdDRho->allPdfs();
  RooAbsPdf *BdDRho_D02KPi_bkg_shape = (RooAbsPdf*)BdDRho_pdfs.find("bkg_shape");  
  BdDRho_D02KPi_bkg_shape->SetName("BdDRho_D02KPi_bkg_shape");
  
  RooArgSet BdDRho_params = work_BdDRho->allVars();
  string word_BdDRho;
  stringstream stream_BdDRho(BdDRho_params.contentsString());
  while(getline(stream_BdDRho, word_BdDRho, ','))
    ((RooRealVar*)BdDRho_params.find(word_BdDRho.c_str()))->setConstant(kTRUE);  
  cout << "got Rho" << endl;
  
  //B2DstKst KEYS PDF
  //both gamma and Pi0 are combined in these pdfs, in the right proportions...
  TFile *_file3 = new TFile("~/BdDKstar_analysis/fitters/MCbackground/B2DstKst_ALL_keys/"
                            "B2DstKst_fitresults.root");

  RooWorkspace *work_B2DstKst = (RooWorkspace*)_file3->Get("B2DstKst_BKG");
  RooArgSet B2DstKst3_pdfs = work_B2DstKst->allPdfs();
  RooAbsPdf *BsDstKst_Dst2D0_bkg_shape001 = (RooAbsPdf*)B2DstKst3_pdfs.find("mass_pdf001Bs");
  BsDstKst_Dst2D0_bkg_shape001->SetName("BsDstKst_Dst2D0_bkg_shape001");
  RooAbsPdf *BsDstKst_Dst2D0_bkg_shape010 = (RooAbsPdf*)B2DstKst3_pdfs.find("mass_pdf010Bs");
  BsDstKst_Dst2D0_bkg_shape010->SetName("BsDstKst_Dst2D0_bkg_shape010");
  RooAbsPdf *BsDstKst_Dst2D0_bkg_shape100 = (RooAbsPdf*)B2DstKst3_pdfs.find("mass_pdf100Bs");
  BsDstKst_Dst2D0_bkg_shape100->SetName("BsDstKst_Dst2D0_bkg_shape100");
  RooAbsPdf *BdDstKst_Dst2D0_bkg_shape001 = (RooAbsPdf*)B2DstKst3_pdfs.find("mass_pdf001Bd");
  BdDstKst_Dst2D0_bkg_shape001->SetName("BdDstKst_Dst2D0_bkg_shape001");
  RooAbsPdf *BdDstKst_Dst2D0_bkg_shape010 = (RooAbsPdf*)B2DstKst3_pdfs.find("mass_pdf010Bd");
  BdDstKst_Dst2D0_bkg_shape010->SetName("BdDstKst_Dst2D0_bkg_shape010");
  RooAbsPdf *BdDstKst_Dst2D0_bkg_shape100 = (RooAbsPdf*)B2DstKst3_pdfs.find("mass_pdf100Bd");
  BdDstKst_Dst2D0_bkg_shape100->SetName("BdDstKst_Dst2D0_bkg_shape100");

  RooArgSet B2DstKst_params = work_B2DstKst->allVars();
  string word_B2DstKst;
  stringstream stream_B2DstKst(B2DstKst_params.contentsString());
  while(getline(stream_B2DstKst, word_B2DstKst, ','))
    ((RooRealVar*)B2DstKst_params.find(word_B2DstKst.c_str()))->setConstant(kTRUE);  
  
  //defining the pdfs...
  
  // DEFINING A MODEL FOR B->DKstar (SS)
  //-------------------------------------
    
  //signal
  RooRealVar Bd_mass_signal_mean("Bd_mass_signal_mean","Bd_mass_signal_mean",5250,5300,"MeV/c^{2}");
  RooRealVar signal_core_sigma_BDKstar("signal_core_sigma_BDKstar","signal_core_sigma_BDKstar",0,50,"MeV/c^{2}");
    
  RooRealVar width_ratio_BDKstar("width_ratio_BDKstar","width_ratio_BDKstar",1,4,"");
  RooFormulaVar signal_wide_sigma_BDKstar("signal_wide_sigma_BDKstar","signal_wide_sigma_BDKstar",
                                          "@0*@1",RooArgList(width_ratio_BDKstar,signal_core_sigma_BDKstar));
  
  RooAbsPdf *signal_core_BDKstar_SS = new RooGaussian("signal_core_BDKstar_SS", "signal_core_BDKstar_SS", 
                                                      B_mass, Bd_mass_signal_mean, signal_core_sigma_BDKstar);
  
  RooAbsPdf *signal_tail_BDKstar_SS = new RooGaussian("signal_tail_BDKstar_SS","signal_tail_BDKstar_SS", 
                                                      B_mass, Bd_mass_signal_mean, signal_wide_sigma_BDKstar);
  
  RooRealVar fCore_BDKstar("fCore_BDKstar","fCore_BDKstar",0,1,"");  

  RooAbsPdf *Bd_signal_BDKstar = new RooAddPdf("Bd_signal_BDKstar","Bd_signal_BDKstar",
                                               RooArgList(*signal_core_BDKstar_SS,*signal_tail_BDKstar_SS),
                                               RooArgList(fCore_BDKstar));  
    

  //combinatoric background
  RooRealVar combinatoric_slope("combinatoric_slope","combinatoric_slope",-1,0,"");
  RooAbsPdf *combkg_mass_BDKstar = new RooChebychev("comb_bkg_BDKstar", "comb_bkg_BDKstar",
                                                    B_mass,RooArgList(combinatoric_slope));
  
  // DEFINING A MODEL FOR B->DKstar (OS)
  //-------------------------------------
    
  //Bs signal peak
  //fixed to be the Bd mass plus the difference from the PDG
  //Bs signal peak is exactly the same shape as Bd just shifted in mass
  RooFormulaVar Bs_mass_signal_mean("Bs_mass_signal_mean","Bs_mass_signal_mean","@0+86.8",RooArgList(Bd_mass_signal_mean));

  //these have the same resolution parameters as the SS case
  RooAbsPdf *signal_core_BDKstar_OS = new RooGaussian("signal_core_BDKstar_OS", "signal_core_BDKstar_OS", 
                                                      B_mass, Bs_mass_signal_mean, signal_core_sigma_BDKstar);
    
  RooAbsPdf *signal_tail_BDKstar_OS = new RooGaussian("signal_tail_BDKstar_OS","signal_tail_BDKstar_OS", 
                                                      B_mass, Bs_mass_signal_mean, signal_wide_sigma_BDKstar);
  
  
  RooAbsPdf *Bs_signal_BDKstar = new RooAddPdf("Bs_signal_BDKstar","Bs_signal_BDKstar",
                                               RooArgList(*signal_core_BDKstar_OS,*signal_tail_BDKstar_OS),
                                               RooArgList(fCore_BDKstar));  
  
  //combinatoric background is the same as in the SS case
      

  // DEFINING A MODEL FOR B->D(KK)Kstar
  //-------------------------------------
  
  //defining the pdfs...
  //Bs signal is the same as that in the OS
  //Bd signal is the same as that in the SS
  //also reuse the low mass shapes for Bs(d) -> D*K*

  //combinatoric background is the same as in the SS case
    
  // DEFINING A MODEL FOR B->D(PiPi)Kstar
  //-------------------------------------
  
  //defining the pdfs...
  //Bs signal is the same as that in the OS
  //Bd signal is the same as that in the SS
  //also reuse the low mass shapes for Bs(d) -> D*K*

  //combinatoric background is the same as in the SS case    
  
  //CREATE INDEX CATEGORY AND JOIN SAMPLES
  //--------------------------------------

  RooCategory complete_sample("complete_sample","complete_sample");
  complete_sample.defineType("BDKstar_SS");
  complete_sample.defineType("BDKstar_OS");
  complete_sample.defineType("BDKstar_KK");
  complete_sample.defineType("BDKstar_PiPi");

  RooDataSet combinedData("combinedData","combinedData",B_mass,Index(complete_sample),Import("BDKstar_SS",*data_KPi_SS),
                          Import("BDKstar_OS",*data_KPi_OS),Import("BDKstar_KK",*data_KK),
                          Import("BDKstar_PiPi",*data_PiPi));
  
  //CONSTRUCT SIMULTANEOUS PDF IN B MASS
  //------------------------------------

  //making individual pdfs with yields
  //----------------------------------
  //---------------------------------------------------------------------------------------------  
  //B --> D(KPi)K* SS
  //---------------------------------------------------------------------------------------------  
  RooRealVar NsigBd_BDKstar_SS("NsigBd_BDKstar_SS","NsigBd_BDKstar_SS",0.0,data_KPi_SS->numEntries()*1.5,"");  
  RooRealVar Ncombkg_BDKstar_SS("Ncombkg_BDKstar_SS","Ncombkg_BDKstar_SS",0.0,data_KPi_SS->numEntries()*1.5,"");
  RooRealVar Nxfeed_BdDRho_BDKstar_SS("Nxfeed_BdDRho_BDKstar_SS","Nxfeed_BdDRho_BDKstar_SS",0.0,
                                      data_KPi_SS->numEntries()*1.5,"");
  
  //this is fixed by the branching ratio of D* -> D0Pi0 relative to D* -> D0gamma
  RooRealVar Dst_gammaPi0_ratio("Dst_gammaPi0_ratio","Dst_gammaPi0_ratio",0.616,"");
  
  RooRealVar Nbkg_BdDstKst_SS001("Nbkg_BdDstKst_SS001","Nbkg_BdDstKst_SS001",0.0,data_KPi_SS->numEntries()*1.5,"");
  RooRealVar Nbkg_BdDstKst_SS010("Nbkg_BdDstKst_SS010","Nbkg_BdDstKst_SS010",0.0,data_KPi_SS->numEntries()*1.5,"");
  RooRealVar Nbkg_BdDstKst_SS100("Nbkg_BdDstKst_SS100","Nbkg_BdDstKst_SS100",0.0,data_KPi_SS->numEntries()*1.5,"");
  

  RooFormulaVar BdDstKst_ratio001_010("BdDstKst_ratio001_010","BdDstKst_ratio001_010","@0/@1",
                                           RooArgList(Nbkg_BdDstKst_SS001,Nbkg_BdDstKst_SS010));  
  RooFormulaVar BdDstKst_ratio001_100("BdDstKst_ratio001_100","BdDstKst_ratio001_100","@0/@1",
                                           RooArgList(Nbkg_BdDstKst_SS001,Nbkg_BdDstKst_SS100));
  
  
  //the combined pdf
  RooAbsPdf *mass_pdf_BDKstar_SS = new RooAddPdf("mass_pdf_BDKstar_SS","mass_pdf_BDKstar_SS",
                                                 RooArgList(*Bd_signal_BDKstar,*combkg_mass_BDKstar,
                                                            *BdDRho_D02KPi_bkg_shape,
                                                            *BdDstKst_Dst2D0_bkg_shape001,
                                                            *BdDstKst_Dst2D0_bkg_shape010,
                                                            *BdDstKst_Dst2D0_bkg_shape100),
                                                 RooArgList(NsigBd_BDKstar_SS,Ncombkg_BDKstar_SS,
                                                            Nxfeed_BdDRho_BDKstar_SS,
                                                            Nbkg_BdDstKst_SS001,
                                                            Nbkg_BdDstKst_SS010,
                                                            Nbkg_BdDstKst_SS100));

  //---------------------------------------------------------------------------------------------    
  //B --> D(KPi)K* OS
  //---------------------------------------------------------------------------------------------  
  RooRealVar NsigBs_BDKstar_OS("NsigBs_BDKstar_OS","NsigBs_BDKstar_OS",0.0,data_KPi_OS->numEntries()*1.5,"");  
  //blind the suppressed ADS signal yield
  RooRealVar NsigBd_BDKstar_OS("NsigBd_BDKstar_OS","NsigBd_BDKstar_OS",0.0,data_KPi_OS->numEntries()*1.5,"");  
  RooUnblindUniform ub_NsigBd_BDKstar_OS("ub_NsigBd_BDKstar_OS","ub_NsigBd_BDKstar_OS","cat_is_minnie",100,NsigBd_BDKstar_OS);

  RooRealVar Ncombkg_BDKstar_OS("Ncombkg_BDKstar_OS","Ncombkg_BDKstar_OS",0.0,data_KPi_OS->numEntries()*1.5,"");
  RooRealVar Nxfeed_BdDRho_BDKstar_OS("Nxfeed_BdDRho_BDKstar_OS","Nxfeed_BdDRho_BDKstar_OS",0.0,
                                      data_KPi_OS->numEntries()*1.5,"");

  RooRealVar Nbkg_BsDstKst_OS001("Nbkg_BsDstKst_OS001","Nbkg_BsDstKst_OS001",0.0,data_KPi_OS->numEntries()*1.5,"");
  RooRealVar Nbkg_BsDstKst_OS010("Nbkg_BsDstKst_OS010","Nbkg_BsDstKst_OS010",0.0,data_KPi_OS->numEntries()*1.5,"");
  RooRealVar Nbkg_BsDstKst_OS100("Nbkg_BsDstKst_OS100","Nbkg_BsDstKst_OS100",0.0,data_KPi_OS->numEntries()*1.5,"");


  RooFormulaVar BsDstKst_ratio001_010("BsDstKst_ratio001_010","BsDstKst_ratio001_010","@0/@1",
                                           RooArgList(Nbkg_BsDstKst_OS001,Nbkg_BsDstKst_OS010));  
  RooFormulaVar BsDstKst_ratio001_100("BsDstKst_ratio001_100","BsDstKst_ratio001_100","@0/@1",
                                           RooArgList(Nbkg_BsDstKst_OS001,Nbkg_BsDstKst_OS100));  

  RooArgList yieldList_OS = RooArgList(NsigBs_BDKstar_OS,ub_NsigBd_BDKstar_OS,
                                       Ncombkg_BDKstar_OS,Nxfeed_BdDRho_BDKstar_OS,
                                       Nbkg_BsDstKst_OS001);
  yieldList_OS.add(RooArgList(Nbkg_BsDstKst_OS010,
                              Nbkg_BsDstKst_OS100));
  
  RooArgList pdfList_OS = RooArgList(*Bs_signal_BDKstar,*Bd_signal_BDKstar,
                                     *combkg_mass_BDKstar,*BdDRho_D02KPi_bkg_shape,
                                     *BsDstKst_Dst2D0_bkg_shape001);
  pdfList_OS.add(RooArgList(*BsDstKst_Dst2D0_bkg_shape010,
                            *BsDstKst_Dst2D0_bkg_shape100));
  
  //the combined pdf
  RooAbsPdf *mass_pdf_BDKstar_OS = new RooAddPdf("mass_pdf_BDKstar_OS","mass_pdf_BDKstar_OS",pdfList_OS,yieldList_OS);

  //---------------------------------------------------------------------------------------------  
  //B --> D(KK)K*
  //---------------------------------------------------------------------------------------------
  RooRealVar Ncombkg_BDKstar_KK("Ncombkg_BDKstar_KK","Ncombkg_BDKstar_KK",0.0,data_KK->numEntries()*1.5,"");
  RooRealVar Npartbkg_BDKstar_KK("Npartbkg_BDKstar_KK","Npartbkg_BDKstar_KK",0.0,data_KK->numEntries()*1.5,"");  
  RooRealVar NsigBs_BDKstar_KK("NsigBs_BDKstar_KK","NsigBs_BDKstar_KK",0.0,data_KK->numEntries()*1.5,"");  
  RooRealVar NsigBd_BDKstar_KK("NsigBd_BDKstar_KK","NsigBd_BDKstar_KK",0.0,data_KK->numEntries()*1.5,"");
  RooUnblindUniform ub_NsigBd_BDKstar_KK("ub_NsigBd_BDKstar_KK","ub_NsigBd_BDKstar_KK","cat_is_minnie",10,NsigBd_BDKstar_KK);
  
  RooRealVar Nbkg_BdDstKst_KK001("Nbkg_BdDstKst_KK001","Nbkg_BdDstKst_KK001",0.0,data_KK->numEntries()*1.5,"");
  RooFormulaVar Nbkg_BdDstKst_KK010("Nbkg_BdDstKst_KK010","Nbkg_BdDstKst_KK010","@0/@1",
                                         RooArgList(Nbkg_BdDstKst_KK001,BdDstKst_ratio001_010));
  RooFormulaVar Nbkg_BdDstKst_KK100("Nbkg_BdDstKst_KK100","Nbkg_BdDstKst_KK100","@0/@1",
                                         RooArgList(Nbkg_BdDstKst_KK001,BdDstKst_ratio001_100));

  RooRealVar Nbkg_BsDstKst_KK001("Nbkg_BsDstKst_KK001","Nbkg_BsDstKst_KK001",0.0,data_KK->numEntries()*1.5,"");
  RooFormulaVar Nbkg_BsDstKst_KK010("Nbkg_BsDstKst_KK010","Nbkg_BsDstKst_KK010","@0/@1",
                                         RooArgList(Nbkg_BsDstKst_KK001,BsDstKst_ratio001_010));
  RooFormulaVar Nbkg_BsDstKst_KK100("Nbkg_BsDstKst_KK100","Nbkg_BsDstKst_KK100","@0/@1",
                                         RooArgList(Nbkg_BsDstKst_KK001,BsDstKst_ratio001_100));


  //use this in PiPi to enforce the same relative contributions of Bs and Bd in the low mass backgrounds of KK and PiPi
  RooFormulaVar BstoBd_DstKst_ratio("BstoBd_DstKst_ratio","BstoBd_DstKst_ratio","@0/@1",
                                    RooArgList(Nbkg_BsDstKst_KK001,Nbkg_BdDstKst_KK001));
  
  RooArgList pdfList_KK = RooArgList(*Bs_signal_BDKstar,*Bd_signal_BDKstar,*combkg_mass_BDKstar,
                                     *BdDstKst_Dst2D0_bkg_shape001,
                                     *BdDstKst_Dst2D0_bkg_shape010,
                                     *BdDstKst_Dst2D0_bkg_shape100);
  pdfList_KK.add(RooArgList(*BsDstKst_Dst2D0_bkg_shape001,
                            *BsDstKst_Dst2D0_bkg_shape010,
                            *BsDstKst_Dst2D0_bkg_shape100));

  cout << "=================PDF LIST IN KK=====================" << endl;
  pdfList_KK.Print();
  
  RooArgList yieldList_KK = RooArgList(NsigBs_BDKstar_KK,ub_NsigBd_BDKstar_KK,Ncombkg_BDKstar_KK,
                                       Nbkg_BdDstKst_KK001,
                                       Nbkg_BdDstKst_KK010,
                                       Nbkg_BdDstKst_KK100);
  yieldList_KK.add(RooArgList(Nbkg_BsDstKst_KK001,
                              Nbkg_BsDstKst_KK010,
                              Nbkg_BsDstKst_KK100));
  //the combined pdf
  RooAbsPdf *mass_pdf_BDKstar_KK = new RooAddPdf("mass_pdf_BDKstar_KK","mass_pdf_BDKstar_KK",pdfList_KK,yieldList_KK);
  
  
  //-----------------------------------------------------------------------------------------------------------------------------

  //-----------------------------------------------------------------------------------------------------------------------------
  //B --> D(PiPi)K*
  //-----------------------------------------------------------------------------------------------------------------------------
  RooRealVar Ncombkg_BDKstar_PiPi("Ncombkg_BDKstar_PiPi","Ncombkg_BDKstar_PiPi",0.0,data_PiPi->numEntries()*1.5,"");
  RooRealVar Npartbkg_BDKstar_PiPi("Npartbkg_BDKstar_PiPi","Npartbkg_BDKstar_PiPi",0.0,data_PiPi->numEntries()*1.5,"");  
  RooRealVar NsigBs_BDKstar_PiPi("NsigBs_BDKstar_PiPi","NsigBs_BDKstar_PiPi",0.0,data_PiPi->numEntries()*1.5,"");  
  RooRealVar NsigBd_BDKstar_PiPi("NsigBd_BDKstar_PiPi","NsigBd_BDKstar_PiPi",0.0,data_PiPi->numEntries()*1.5,"");
  RooUnblindUniform ub_NsigBd_BDKstar_PiPi("ub_NsigBd_BDKstar_PiPi","ub_NsigBd_BDKstar_PiPi","cat_is_minnie",10,
                                           NsigBd_BDKstar_PiPi);
  
  RooRealVar Nbkg_BdDstKst_PiPi001("Nbkg_BdDstKst_PiPi001","Nbkg_BdDstKst_PiPi001",0.0,
                                        data_PiPi->numEntries()*1.5,"");
  RooFormulaVar Nbkg_BdDstKst_PiPi010("Nbkg_BdDstKst_PiPi010","Nbkg_BdDstKst_PiPi010","@0/@1",
                                           RooArgList(Nbkg_BdDstKst_PiPi001,BdDstKst_ratio001_010));
  RooFormulaVar Nbkg_BdDstKst_PiPi100("Nbkg_BdDstKst_PiPi100","Nbkg_BdDstKst_PiPi100","@0/@1",
                                           RooArgList(Nbkg_BdDstKst_PiPi001,BdDstKst_ratio001_100));  
  
  RooFormulaVar Nbkg_BsDstKst_PiPi001("Nbkg_BsDstKst_PiPi001","Nbkg_BsDstKst_PiPi001","@0*@1",
                                           RooArgList(Nbkg_BdDstKst_PiPi001,BstoBd_DstKst_ratio));
  RooFormulaVar Nbkg_BsDstKst_PiPi010("Nbkg_BsDstKst_PiPi010","Nbkg_BsDstKst_PiPi010","@0/@1",
                                           RooArgList(Nbkg_BsDstKst_PiPi001,BsDstKst_ratio001_010));
  RooFormulaVar Nbkg_BsDstKst_PiPi100("Nbkg_BsDstKst_PiPi100","Nbkg_BsDstKst_PiPi100","@0/@1",
                                           RooArgList(Nbkg_BsDstKst_PiPi001,BsDstKst_ratio001_100));
  
  RooArgList pdfList_PiPi = RooArgList(*Bs_signal_BDKstar,*Bd_signal_BDKstar,*combkg_mass_BDKstar,
                                       *BdDstKst_Dst2D0_bkg_shape001,
                                       *BdDstKst_Dst2D0_bkg_shape010,
                                       *BdDstKst_Dst2D0_bkg_shape100);
  pdfList_PiPi.add(RooArgList(*BsDstKst_Dst2D0_bkg_shape001,
                              *BsDstKst_Dst2D0_bkg_shape010,
                              *BsDstKst_Dst2D0_bkg_shape100));
  
  cout << "=================PDF LIST IN PiPi=====================" << endl;
  pdfList_PiPi.Print();
  
  RooArgList yieldList_PiPi = RooArgList(NsigBs_BDKstar_PiPi,ub_NsigBd_BDKstar_PiPi,Ncombkg_BDKstar_PiPi,
                                         Nbkg_BdDstKst_PiPi001,
                                         Nbkg_BdDstKst_PiPi010,
                                         Nbkg_BdDstKst_PiPi100);
  yieldList_PiPi.add(RooArgList(Nbkg_BsDstKst_PiPi001,
                                Nbkg_BsDstKst_PiPi010,
                                Nbkg_BsDstKst_PiPi100));
  //the combined pdf
  RooAbsPdf *mass_pdf_BDKstar_PiPi = new RooAddPdf("mass_pdf_BDKstar_PiPi","mass_pdf_BDKstar_PiPi",pdfList_PiPi,yieldList_PiPi);  
  
  //-----------------------------------------------------------------------------------------------------------------------

  // Construct a simultaneous pdf using category sample as index
  RooSimultaneous *simPdf = new RooSimultaneous("simPdf","simPdf",complete_sample);  

  // Associate pdfs with their types in the complete sample and add them to the simulataneous pdf
  simPdf->addPdf(*mass_pdf_BDKstar_SS,"BDKstar_SS");
  simPdf->addPdf(*mass_pdf_BDKstar_OS,"BDKstar_OS");
  simPdf->addPdf(*mass_pdf_BDKstar_KK,"BDKstar_KK");
  simPdf->addPdf(*mass_pdf_BDKstar_PiPi,"BDKstar_PiPi");
  cout << "made pdf" << endl;
  
  //extended likelihood fit: simultaneous fit to the data and save the fitresult:
  RooFitResult *mass_result = simPdf->fitTo(combinedData,Extended(),Save(),PrintLevel(3));
  //print the mass fit result to the screen:
  mass_result->Print();
  return mass_result;
}