void Bmass_Fit_Bin3(int year=2016, const char* bin="Bin3", const char* w_s="bin1C", double peak_frac=0.064)
{

  gROOT->Reset();
  gROOT->SetBatch(1);
  gStyle->SetOptFit(kTRUE);
  Color_t white=10;
  gStyle->SetCanvasColor(white);gSystem->Load("libRooFit");
  using namespace RooFit;
  using namespace RooStats;


  //gSystem->Load("$LIB");

  char data_ws[200];
  TFile *F_data = new TFile("CMC_2016_Bin3_062023_1.root","UPDATED");
  strcpy(data_ws, "data_cmc");

  RooDataSet *Data = (RooDataSet *)F_data->Get(data_ws);

           TFile *F = new TFile("wspace_2016_bin1C.root");
           TFile *F_sigMDCB = new TFile("wspace_updated_sigMDCBFitter_2016_bin1C.root");
           TFile *F_MDCB_Kstar = new TFile("wspace_updated_bkgM_KStarFitter_2016_bin1C.root");
           TFile *F_bkgCombM = new TFile("wspace_updated_massSB_2016_bin1C.root");
           TFile *F_sig2D = new TFile("wspace_updated_sig2DFitter_2016_bin1C.root");

          char WS[200], F_ws[200];

          char F_sigMDCB_ws[200], F_MDCB_Kstar_ws[200], F_sig2D_ws[200], F_bkgCombM_ws[200];

          strcpy(WS, Form(".%s",w_s));
          strcpy(F_ws, Form("wspace.%d", year));
          strcpy(F_sigMDCB_ws, Form("wspace_updated.sigMDCBFitter_%d", year));
          strcpy(F_MDCB_Kstar_ws, Form("wspace_updated.bkgM_KStarFitter_%d", year));
          strcpy(F_sig2D_ws, Form("wspace_updated.sig2DFitter_%d", year));
          strcpy(F_bkgCombM_ws, Form("wspace_updated.massSB_%d",year));

          strcat (F_ws, WS);
          strcat (F_sigMDCB_ws, WS);
          strcat (F_MDCB_Kstar_ws, WS);
          strcat (F_sig2D_ws, WS);
          strcat (F_bkgCombM_ws, WS);

          printf("%s",F_ws);
          printf("\n");

          RooWorkspace *w = (RooWorkspace *)F->Get(F_ws);

          RooWorkspace *w_sigMDCB = (RooWorkspace *)F_sigMDCB->Get(F_sigMDCB_ws);

          RooWorkspace *w_MDCB_Kstar = (RooWorkspace *)F_MDCB_Kstar->Get(F_MDCB_Kstar_ws);

          RooWorkspace *w_sig2D = (RooWorkspace *)F_sig2D->Get(F_sig2D_ws);

          RooWorkspace *w_bkgCombM = (RooWorkspace *)F_bkgCombM->Get(F_bkgCombM_ws);
         
          //w->Print("v");

   	  RooRealVar *Bmass = w_sigMDCB->var("Bmass");

          Bmass->setMin(5.25); Bmass->setMax(5.55);

          RooRealVar *CosThetaK = w_sig2D->var("CosThetaK");

          RooRealVar *CosThetaL = w_sig2D->var("CosThetaL");

          RooArgSet variables(*Bmass, *CosThetaK, *CosThetaL);


          RooArgSet p_sM = w_sigMDCB->allVars();

          RooArgSet p_pM = w_MDCB_Kstar->allVars();

          RooArgSet p_eff = w_sig2D->allVars();

          cout<<"#Parameters in Signal Mass Model = "<<p_sM.getSize()<<endl;

          cout<<"#Parameters in Peak. Mass Model = "<<p_pM.getSize()<<endl;

          cout<<"#Parameters in Efficiency Model = "<<p_eff.getSize()<<endl;

	  char YEAR[200], sM_PARA[200], pM_PARA[200], eff_PARA[200];

       
        for(int i = 0+1; i<p_sM.getSize();i++){

        const char *sM_para = p_sM[i]->GetName();

        //strcpy(YEAR, "");

        strcpy(YEAR, Form("_%d", year));

        strcpy(sM_PARA, sM_para);

        strcat (sM_PARA, YEAR);

        p_sM[i]->SetName(sM_PARA);

        w_sigMDCB->var(sM_PARA)->setConstant(kTRUE) ;
       }

        for(int i = 0+1; i<p_pM.getSize();i++){

        const char *pM_para = p_pM[i]->GetName();

        strcpy(YEAR, Form("_%d", year));

        //strcpy(YEAR, "");

        strcpy(pM_PARA, pM_para);

        strcat (pM_PARA, YEAR);

        p_pM[i]->SetName(pM_PARA);

        w_MDCB_Kstar->var(pM_PARA)->setConstant(kTRUE) ;
       }


        for(int i = 0+4; i<p_eff.getSize();i++){

        const char *eff_para = p_eff[i]->GetName();

        strcpy(YEAR, Form("_%d", year));

        //strcpy(YEAR, "");

        strcpy(eff_PARA, eff_para);

        strcat (eff_PARA, YEAR);

        p_eff[i]->SetName(eff_PARA);

        w_sig2D->var(eff_PARA)->setConstant(kTRUE) ;
       }


       RooRealVar nSig("nSig","signal events",20,10,500);
       RooRealVar nBkgComb("nBkgComb","number of exponential bkg events",20,10,500);

       RooConstVar PeakFrac("PeakFrac","Kstar Yield Constraint", peak_frac);
       RooProduct nBkgPeak("nBkgPeak","Peaking Bkg events",RooArgList(nSig, PeakFrac));

       RooAbsPdf *sigMDCB = w_sigMDCB->pdf("f_sigMDCB");

       RooAbsPdf *sig2D = w_sig2D->pdf("f_sig2D_ts");

       RooAbsPdf *bkgCombM = w_bkgCombM->pdf("bkgCombM");

       RooAbsPdf *MDCB_Kstar = w_MDCB_Kstar->pdf("f_bkgM_KStar");

       RooProdPdf sig3D("sig3D","sig3D", RooArgList(*sigMDCB, *sig2D));    //2016 3D Signal Model
       
       RooAddPdf model("model","S+Exp+peakBkg", RooArgList(sig3D, *bkgCombM, *MDCB_Kstar), RooArgList(nSig,nBkgComb,nBkgPeak));    ///Total Model

       model.Print("t");

   // for debugging, check model tree
       model.printCompactTree();
       model.graphVizTree("model.dot");

   // turn off some messages
   RooMsgService::instance().setStreamStatus(0, kFALSE);
   RooMsgService::instance().setStreamStatus(1, kFALSE);
   RooMsgService::instance().setStreamStatus(2, kFALSE);

   // needed to remove some sporious messages from RooAddPdf when generating events
   //model.fixCoefNormalization(variables);

   RooFitResult* fitres = model.fitTo(*Data, Extended(kTRUE), Save(kTRUE), Minimizer("Minuit") , Minos(kTRUE)) ;

   fitres->Print("v"); 

   TCanvas *c = new TCanvas("c","c", 1000,800);
   TPad *pad1   = new TPad("pad1","pad1", 0.1, 0.25, 0.995, 0.97);
   TPad *pad2   = new TPad("pad2","pad2", 0.1, 0.02, 0.995, 0.24);
   pad1->Draw();
   pad2->Draw();

   pad1->cd();

  RooPlot *xframe = Bmass->frame(Title("B invariant mass"),Bins(20)); //50

  xframe->SetTitle("");
  xframe->SetMinimum(0.);
  xframe->GetXaxis()->SetTitleOffset(1.2);
  xframe->GetXaxis()->SetTitle("#bf{m}(#bf{K^{+}K^{-}}#mu^{+}#mu^{-}) #bf{[GeV]}");
  xframe->GetYaxis()->SetTitle("");
  xframe->GetYaxis()->SetRangeUser(0.0, 1.5*xframe->GetMaximum());

  Data->plotOn(xframe,RooFit::Name("data"));

  model.plotOn(xframe,RooFit::Name("pdf"),LineStyle(1), LineColor(4),LineWidth(5)); ///pdf for both sig+bkg

  RooHist* hpull = xframe->pullHist() ;

  double chi2 = xframe->chiSquare();
  cout << "chi2/dof = " << chi2 << endl;

  model.plotOn(xframe, Components(*sigMDCB), FillStyle(3013), FillColor(kBlue), LineColor(3), DrawOption("F") ); //pdf for only signal

  model.plotOn(xframe, Components(*bkgCombM), LineColor(2) ); //pdf for Comb. Bkg

  model.plotOn(xframe, Components(*MDCB_Kstar), /*FillStyle(3013),*/ FillColor(kOrange), LineColor(kOrange), DrawOption("F") ); //pdf for only signal

  xframe->Draw();

  TPaveText* paveText = new TPaveText(0.50,0.50,0.90,0.80,"NDC");
  paveText->SetBorderSize(0.0);
  paveText->SetFillColor(kWhite);
  paveText->SetFillStyle(0);
  paveText->SetTextSize(0.04);//0.015
  paveText->AddText(Form("N_{Sig}  = %.0f #pm %.0f", nSig.getVal(), nSig.getError()));
  paveText->AddText(Form("N_{Comb Bkg} = %.0f #pm %.0f", nBkgComb.getVal(), nBkgComb.getError()));
  paveText->AddText(Form("N_{Peaking Bkg} = %.0f #pm %.000f", nBkgPeak.getVal(), nSig.getError()*PeakFrac.getVal()));
 //paveText->AddText(Form("#lambda  = %.4f #pm %.4f", lambda.getVal(), lambda.getError()));
  paveText->AddText(Form("#chi^{2}/ndf = %.2f",chi2));
  paveText->Draw();

 pad2->cd();

 RooPlot* pull = Bmass->frame(Title("Pull Distribution")) ;
 pull->addPlotable(hpull,"BX");
 hpull->SetFillColor(kBlue);
 pull->SetMinimum(-3);
 pull->SetMaximum(+3);
 pull->GetYaxis()->SetNdivisions(5);
 pull->SetTitle("#bf{Pull Distribution}");
 pull->GetYaxis()->SetTitle( "Pull" );
 pull->GetXaxis()->SetTitle("");
 pull->GetXaxis()->SetLabelOffset(1);
 pull->GetXaxis()->SetLimits( Bmass->getMin(), Bmass->getMax() );
 pull->Draw("P");

c->Update();

c->SaveAs(Form("./1D_data_Fit_%d_%s_test.png",year,bin));

}