#include <TMath.h>
#include <TGraph2D.h>
#include <TRandom.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TF2.h>
#include <TH1.h>
#include "RooSimultaneous.h"



void kFactorsFit()
{
  //  TCanvas *c1 = new TCanvas("c1","comparison",100,10,1100,800);
 //   TPaveText *pt = new TPaveText(.05,.1,.95,.8);
//    c1->SetWindowSize(700, 800);
    char nl[26];
    char nl3[26];
     char nl2[26];
    gStyle->SetOptStat(0);
gStyle->SetOptTitle(0);
  

    TFile *f02_0 = new TFile("/Users/Anton/rootf/DY0_v19.root");
    TFile *f02_1 = new TFile("/Users/Anton/rootf/DY1_v19.root");
    TFile *f02_2 = new TFile("/Users/Anton/rootf/DY2_v19.root");
    
    TFile *f02_ttbar = new TFile("/Users/Anton/rootf/TTbar_v19.root");
    TFile *f02_STt = new TFile("/Users/Anton/rootf/STt_v19.root");
    TFile *f02_STta = new TFile("/Users/Anton/rootf/STta_v19.root");
    TFile *f02_STtW = new TFile("/Users/Anton/rootf/STtW_v19.root");
    TFile *f02_STtWa = new TFile("/Users/Anton/rootf/STtWa_v19.root");
    TFile *f02_STs = new TFile("/Users/Anton/rootf/STs_v19.root");
    
    TFile *f02_WW = new TFile("/Users/Anton/rootf/WW_v19.root");
    TFile *f02_ZZ = new TFile("/Users/Anton/rootf/ZZ_v19.root");
    TFile *f02_WZ = new TFile("/Users/Anton/rootf/WZ_v19.root");
    
    TFile *f02_data = new TFile("/Users/Anton/rootf/DATA_v19.root");
    
    double totalLumi = 5440.447  + 2385.032 + 4254.329 + 4043.236 + 3105.456 + 7544.016 + 8526.350 + 216.594;
    double scale0 = (45338390 - 4558628)/4754 ; // nlo 49953989
    double scale1 = (36422242 - 13666818)/888.9 ; // nlo 49940088
    double scale2 = (27326526 - 15005687)/348.8 ; // nlo
    
    double scaleTTbar = (29820504 - 14487223)/831.76 ; // nlo
    double scaleSTt = (5960044)/136*0.35 ;
    double scaleSTta = (3999090)/81*0.35 ;
    double scaleSTtW = (6944090)/35.6 ;
    double scaleSTtWa = (6925238)/35.6 ;
    
    double scaleSTs = (811495 - 188505)/96.74 ;
    
    double scaleWW = (1979000)/118.7 ;
    double scaleZZ = (971600)/16.5 ;
    double scaleWZ = (1000000)/47.13  ;
    
    
    TH1F *heta_ttbar_notag = (TH1F*)f02_ttbar->Get("hYstar");
    TH1F *heta_STt_notag = (TH1F*)f02_STt->Get("hYstar");
    TH1F *heta_STs_notag = (TH1F*)f02_STs->Get("hYstar");
    TH1F *heta_STta_notag = (TH1F*)f02_STta->Get("hYstar");
    TH1F *heta_STtW_notag = (TH1F*)f02_STtW->Get("hYstar");
    TH1F *heta_STtWa_notag = (TH1F*)f02_STtWa->Get("hYstar");
    TH1F *heta_WW_notag = (TH1F*)f02_WW->Get("hYstar");
    TH1F *heta_ZZ_notag = (TH1F*)f02_ZZ->Get("hYstar");
    TH1F *heta_WZ_notag = (TH1F*)f02_WZ->Get("hYstar");
    TH1F *heta_ttbar_ctag = (TH1F*)f02_ttbar->Get("hYstar_ctag");
    TH1F *heta_STt_ctag = (TH1F*)f02_STt->Get("hYstar_ctag");
    TH1F *heta_STs_ctag = (TH1F*)f02_STs->Get("hYstar_ctag");
    TH1F *heta_STta_ctag = (TH1F*)f02_STta->Get("hYstar_ctag");
    TH1F *heta_STtW_ctag = (TH1F*)f02_STtW->Get("hYstar_ctag");
    TH1F *heta_STtWa_ctag = (TH1F*)f02_STtWa->Get("hYstar_ctag");
    TH1F *heta_WW_ctag = (TH1F*)f02_WW->Get("hYstar_ctag");
    TH1F *heta_ZZ_ctag = (TH1F*)f02_ZZ->Get("hYstar_ctag");
    TH1F *heta_WZ_ctag = (TH1F*)f02_WZ->Get("hYstar_ctag");
    TH1F *heta_ttbar_btag = (TH1F*)f02_ttbar->Get("hYstar_btag");
    TH1F *heta_STt_btag = (TH1F*)f02_STt->Get("hYstar_btag");
    TH1F *heta_STs_btag = (TH1F*)f02_STs->Get("hYstar_btag");
    TH1F *heta_STta_btag = (TH1F*)f02_STta->Get("hYstar_btag");
    TH1F *heta_STtW_btag = (TH1F*)f02_STtW->Get("hYstar_btag");
    TH1F *heta_STtWa_btag = (TH1F*)f02_STtWa->Get("hYstar_btag");
    TH1F *heta_WW_btag = (TH1F*)f02_WW->Get("hYstar_btag");
    TH1F *heta_ZZ_btag = (TH1F*)f02_ZZ->Get("hYstar_btag");
    TH1F *heta_WZ_btag = (TH1F*)f02_WZ->Get("hYstar_btag");
    
    TH1F *heta_DY0_bottom_notag = (TH1F*)f02_0->Get("hYstar_bottom_notag");
    TH1F *heta_DY1_bottom_notag = (TH1F*)f02_1->Get("hYstar_bottom_notag");
    TH1F *heta_DY2_bottom_notag = (TH1F*)f02_2->Get("hYstar_bottom_notag");
    TH1F *heta_DY0_charm_notag = (TH1F*)f02_0->Get("hYstar_charm_notag");
    TH1F *heta_DY1_charm_notag = (TH1F*)f02_1->Get("hYstar_charm_notag");
    TH1F *heta_DY2_charm_notag = (TH1F*)f02_2->Get("hYstar_charm_notag");
    TH1F *heta_DY0_light_notag = (TH1F*)f02_0->Get("hYstar_light_notag");
    TH1F *heta_DY1_light_notag = (TH1F*)f02_1->Get("hYstar_light_notag");
    TH1F *heta_DY2_light_notag = (TH1F*)f02_2->Get("hYstar_light_notag");
    TH1F *heta_DY0_bottom_ctag = (TH1F*)f02_0->Get("hYstar_bottom_ctag");
    TH1F *heta_DY1_bottom_ctag = (TH1F*)f02_1->Get("hYstar_bottom_ctag");
    TH1F *heta_DY2_bottom_ctag = (TH1F*)f02_2->Get("hYstar_bottom_ctag");
    TH1F *heta_DY0_charm_ctag = (TH1F*)f02_0->Get("hYstar_charm_ctag");
    TH1F *heta_DY1_charm_ctag = (TH1F*)f02_1->Get("hYstar_charm_ctag");
    TH1F *heta_DY2_charm_ctag = (TH1F*)f02_2->Get("hYstar_charm_ctag");
    TH1F *heta_DY0_light_ctag = (TH1F*)f02_0->Get("hYstar_light_ctag");
    TH1F *heta_DY1_light_ctag = (TH1F*)f02_1->Get("hYstar_light_ctag");
    TH1F *heta_DY2_light_ctag = (TH1F*)f02_2->Get("hYstar_light_ctag");
    TH1F *heta_DY0_bottom_btag = (TH1F*)f02_0->Get("hYstar_bottom_btag");
    TH1F *heta_DY1_bottom_btag = (TH1F*)f02_1->Get("hYstar_bottom_btag");
    TH1F *heta_DY2_bottom_btag = (TH1F*)f02_2->Get("hYstar_bottom_btag");
    TH1F *heta_DY0_charm_btag = (TH1F*)f02_0->Get("hYstar_charm_btag");
    TH1F *heta_DY1_charm_btag = (TH1F*)f02_1->Get("hYstar_charm_btag");
    TH1F *heta_DY2_charm_btag = (TH1F*)f02_2->Get("hYstar_charm_btag");
    TH1F *heta_DY0_light_btag = (TH1F*)f02_0->Get("hYstar_light_btag");
    TH1F *heta_DY1_light_btag = (TH1F*)f02_1->Get("hYstar_light_btag");
    TH1F *heta_DY2_light_btag = (TH1F*)f02_2->Get("hYstar_light_btag");
    
    heta_DY0_bottom_notag->Scale((totalLumi   )/scale0);
    heta_DY1_bottom_notag->Scale((totalLumi   )/scale1);
    heta_DY2_bottom_notag->Scale((totalLumi   )/scale2);
    heta_DY0_light_notag->Scale((totalLumi   )/scale0);
    heta_DY1_light_notag->Scale((totalLumi   )/scale1);
    heta_DY2_light_notag->Scale((totalLumi   )/scale2);
    heta_DY0_charm_notag->Scale((totalLumi   )/scale0);
    heta_DY1_charm_notag->Scale((totalLumi   )/scale1);
    heta_DY2_charm_notag->Scale((totalLumi   )/scale2);
    heta_DY0_bottom_ctag->Scale((totalLumi   )/scale0);
    heta_DY1_bottom_ctag->Scale((totalLumi   )/scale1);
    heta_DY2_bottom_ctag->Scale((totalLumi   )/scale2);
    heta_DY0_light_ctag->Scale((totalLumi   )/scale0);
    heta_DY1_light_ctag->Scale((totalLumi   )/scale1);
    heta_DY2_light_ctag->Scale((totalLumi   )/scale2);
    heta_DY0_charm_ctag->Scale((totalLumi   )/scale0);
    heta_DY1_charm_ctag->Scale((totalLumi   )/scale1);
    heta_DY2_charm_ctag->Scale((totalLumi   )/scale2);
    heta_DY0_bottom_btag->Scale((totalLumi   )/scale0);
    heta_DY1_bottom_btag->Scale((totalLumi   )/scale1);
    heta_DY2_bottom_btag->Scale((totalLumi   )/scale2);
    heta_DY0_light_btag->Scale((totalLumi   )/scale0);
    heta_DY1_light_btag->Scale((totalLumi   )/scale1);
    heta_DY2_light_btag->Scale((totalLumi   )/scale2);
    heta_DY0_charm_btag->Scale((totalLumi   )/scale0);
    heta_DY1_charm_btag->Scale((totalLumi   )/scale1);
    heta_DY2_charm_btag->Scale((totalLumi   )/scale2);
    
    
    heta_DY0_bottom_notag->Add(heta_DY1_bottom_notag);
    heta_DY0_bottom_notag->Add(heta_DY2_bottom_notag);
    heta_DY0_light_notag->Add(heta_DY1_light_notag);
    heta_DY0_light_notag->Add(heta_DY2_light_notag);
    heta_DY0_charm_notag->Add(heta_DY1_charm_notag);
    heta_DY0_charm_notag->Add(heta_DY2_charm_notag);
    
    heta_DY0_bottom_ctag->Add(heta_DY1_bottom_ctag);
    heta_DY0_bottom_ctag->Add(heta_DY2_bottom_ctag);
    heta_DY0_light_ctag->Add(heta_DY1_light_ctag);
    heta_DY0_light_ctag->Add(heta_DY2_light_ctag);
    heta_DY0_charm_ctag->Add(heta_DY1_charm_ctag);
    heta_DY0_charm_ctag->Add(heta_DY2_charm_ctag);
    
    heta_DY0_bottom_btag->Add(heta_DY1_bottom_btag);
    heta_DY0_bottom_btag->Add(heta_DY2_bottom_btag);
    heta_DY0_light_btag->Add(heta_DY1_light_btag);
    heta_DY0_light_btag->Add(heta_DY2_light_btag);
    heta_DY0_charm_btag->Add(heta_DY1_charm_btag);
    heta_DY0_charm_btag->Add(heta_DY2_charm_btag);
    
    
    TH1F *heta_data_notag = (TH1F*)f02_data->Get("hYstar");
    TH1F *heta_data_ctag = (TH1F*)f02_data->Get("hYstar_ctag");
    TH1F *heta_data_btag = (TH1F*)f02_data->Get("hYstar_btag");
    
    cout<<"scale ttbar, stt, sts, stta, sttw, sttwa, ww, zz, wz:  "<<scaleTTbar<<" , "<<scaleSTt<<" , "<<scaleSTs<<" , "<<scaleSTta<<" , "<<scaleSTtW<<" , "<<scaleSTtWa<<" , "<<scaleWW<<" , "<<scaleZZ<<" , "<<scaleWZ<<endl;
    
    heta_ttbar_notag->Scale((totalLumi  )/scaleTTbar);
    heta_STt_notag->Scale((totalLumi  )/scaleSTt);
    heta_STs_notag->Scale((totalLumi  )/scaleSTs);
    heta_STta_notag->Scale((totalLumi  )/scaleSTta);
    heta_STtW_notag->Scale((totalLumi  )/scaleSTtW);
    heta_STtWa_notag->Scale((totalLumi  )/scaleSTtWa);
    heta_WW_notag->Scale((totalLumi  )/scaleWW);
    heta_ZZ_notag->Scale((totalLumi  )/scaleZZ);
    heta_WZ_notag->Scale((totalLumi  )/scaleWZ);
    heta_ttbar_ctag->Scale((totalLumi  )/scaleTTbar);
    heta_STt_ctag->Scale((totalLumi  )/scaleSTt);
    heta_STs_ctag->Scale((totalLumi  )/scaleSTs);
    heta_STta_ctag->Scale((totalLumi  )/scaleSTta);
    heta_STtW_ctag->Scale((totalLumi  )/scaleSTtW);
    heta_STtWa_ctag->Scale((totalLumi  )/scaleSTtWa);
    heta_WW_ctag->Scale((totalLumi  )/scaleWW);
    heta_ZZ_ctag->Scale((totalLumi  )/scaleZZ);
    heta_WZ_ctag->Scale((totalLumi  )/scaleWZ);
    heta_ttbar_btag->Scale((totalLumi  )/scaleTTbar);
    heta_STt_btag->Scale((totalLumi  )/scaleSTt);
    heta_STs_btag->Scale((totalLumi  )/scaleSTs);
    heta_STta_btag->Scale((totalLumi  )/scaleSTta);
    heta_STtW_btag->Scale((totalLumi  )/scaleSTtW);
    heta_STtWa_btag->Scale((totalLumi  )/scaleSTtWa);
    heta_WW_btag->Scale((totalLumi  )/scaleWW);
    heta_ZZ_btag->Scale((totalLumi  )/scaleZZ);
    heta_WZ_btag->Scale((totalLumi  )/scaleWZ);
    heta_WW_notag->Add(heta_ZZ_notag);
    heta_WW_notag->Add(heta_WZ_notag);
    heta_ttbar_notag->Add(heta_STt_notag);
    heta_ttbar_notag->Add(heta_STs_notag);
    heta_ttbar_notag->Add(heta_STta_notag);
    heta_ttbar_notag->Add(heta_STtW_notag);
    heta_ttbar_notag->Add(heta_STtWa_notag);
    heta_ttbar_notag->Add(heta_WW_notag);
    heta_data_notag->Add(heta_ttbar_notag, -1);
    heta_WW_ctag->Add(heta_ZZ_ctag);
    heta_WW_ctag->Add(heta_WZ_ctag);
    heta_ttbar_ctag->Add(heta_STt_ctag);
    heta_ttbar_ctag->Add(heta_STs_ctag);
    heta_ttbar_ctag->Add(heta_STta_ctag);
    heta_ttbar_ctag->Add(heta_STtW_ctag);
    heta_ttbar_ctag->Add(heta_STtWa_ctag);
    heta_ttbar_ctag->Add(heta_WW_ctag);
    heta_data_ctag->Add(heta_ttbar_ctag, -1);
    heta_WW_btag->Add(heta_ZZ_btag);
    heta_WW_btag->Add(heta_WZ_btag);
    heta_ttbar_btag->Add(heta_STt_btag);
    heta_ttbar_btag->Add(heta_STs_btag);
    heta_ttbar_btag->Add(heta_STta_btag);
    heta_ttbar_btag->Add(heta_STtW_btag);
    heta_ttbar_btag->Add(heta_STtWa_btag);
    heta_ttbar_btag->Add(heta_WW_btag);
    heta_data_btag->Add(heta_ttbar_btag, -1);
    
    gSystem->Load("libRooFit");
    using namespace RooFit;
    RooWorkspace w("w",kTRUE) ;
    
    RooRealVar x("x","x",0.0,5.0);
    
  //  heta_data_notag->Scale(0.1);
  //  heta_DY0_bottom_notag->Scale(0.1);
  //  heta_DY0_charm_notag->Scale(0.1);
  //  heta_DY0_light_notag->Scale(0.1);
    
    RooDataHist data("data","data",x,heta_data_btag);
    RooDataHist mc1("mc1","scaled mc1",x,heta_DY0_bottom_btag);
    RooDataHist mc2("mc2","scaled mc2",x,heta_DY0_charm_btag);
    RooDataHist mc3("mc3","scaled mc3",x,heta_DY0_light_btag);
    
    RooDataHist datac("datac","datac",x,heta_data_ctag);
    RooDataHist mc1c("mc1c","scaled mc1",x,heta_DY0_bottom_ctag);
    RooDataHist mc2c("mc2c","scaled mc2",x,heta_DY0_charm_ctag);
    RooDataHist mc3c("mc3c","scaled mc3",x,heta_DY0_light_ctag);
    
    RooDataHist datano("datano","datano",x,heta_data_notag);
    RooDataHist mc1no("mc1no","scaled mc1",x,heta_DY0_bottom_notag);
    RooDataHist mc2no("mc2no","scaled mc2",x,heta_DY0_charm_notag);
    RooDataHist mc3no("mc3no","scaled mc3",x,heta_DY0_light_notag);
    
    
  
    
    RooHistPdf modelmc1("modelmc1","modelmc1",x, mc1);
    RooHistPdf modelmc2("modelmc2","modelmc2",x, mc2);
    RooHistPdf modelmc3("modelmc3","modelmc3",x, mc3);
    
    RooHistPdf modelmc1c("modelmc1c","modelmc1c",x, mc1c);
    RooHistPdf modelmc2c("modelmc2c","modelmc2c",x, mc2c);
    RooHistPdf modelmc3c("modelmc3c","modelmc3c",x, mc3c);
    
    RooHistPdf modelmc1no("modelmc1no","modelmc1no",x, mc1no);
    RooHistPdf modelmc2no("modelmc2no","modelmc2no",x, mc2no);
    RooHistPdf modelmc3no("modelmc3no","modelmc3no",x, mc3no);
    
    RooRealVar mean1_yield("mean1_yield","fraction of background 1 ",10000,0.0,9000000) ;
     RooRealVar mean2_yield("mean2_yield","fraction of background 1 ",10000,0.0,9000000) ;
     RooRealVar mean3_yield("mean3_yield","fraction of background 1 ",10000,0.0,9000000) ;
    
    RooRealVar mean1("mean1","fraction of background 1 ",1.0,0.7,1.3) ;
    RooFormulaVar mc1_yield_n("mean_1","74197.0*mean1",
                       RooArgList(mean1)) ;
    RooRealVar mean2("mean2","fraction of background 2 ",1.0,0.7,1.3) ;
    RooFormulaVar mc2_yield_n("mean_2","7006.0*mean2",
                              RooArgList(mean2)) ;
    RooRealVar mean3("mean3","fraction of background 3 ",1.0,0.99,1.01) ;
    RooFormulaVar mc3_yield_n("mean_3","2040.0*mean3",
                              RooArgList(mean3)) ;
    
    RooFormulaVar mc1_yield_nc("mean_1c","20824.0*mean1",
                              RooArgList(mean1)) ;
    RooFormulaVar mc2_yield_nc("mean_2c","49472.0*mean2",
                              RooArgList(mean2)) ;
    RooFormulaVar mc3_yield_nc("mean_3c","25327.0*mean3",
                              RooArgList(mean3)) ;
    
    RooFormulaVar mc1_yield_nno("mean_1no","154667.0*mean1",
                               RooArgList(mean1)) ;
    RooFormulaVar mc2_yield_nno("mean_2no","258050.0*mean2",
                               RooArgList(mean2)) ;
    RooFormulaVar mc3_yield_nno("mean_3no","2086300.0*mean3",
                               RooArgList(mean3)) ;
    
    RooAddPdf model1("model1","model1",RooArgList(modelmc1,modelmc2,modelmc3),RooArgList(mc1_yield_n, mc2_yield_n, mc3_yield_n));
    RooAddPdf model2("model2","model2",RooArgList(modelmc1c,modelmc2c,modelmc3c),RooArgList(mc1_yield_nc, mc2_yield_nc, mc3_yield_nc));
    RooAddPdf model3("model3","model3",RooArgList(modelmc1no,modelmc2no,modelmc3no),RooArgList(mc1_yield_nno, mc2_yield_nno, mc3_yield_nno));
    
  TCanvas* c22 = new TCanvas("X","Y",1);
   // heta_data_notag->Draw();
    RooPlot* dframe = x.frame(Title("Data"));
    datano.plotOn(dframe, RooFit::LineColor(kRed));
  
   // model.fitTo(data,Extended(),Save());
 //   model2.fitTo(datac,Extended(),Save());
 //   model2.paramOn(dframe);
 //   cout<<"data, m1, m2, m3:  "<<heta_data_notag->Integral()<<" , "<<heta_DY0_bottom_notag->Integral()<<" , "<<heta_DY0_charm_notag->Integral()<<" , "<<heta_DY0_light_notag->Integral()<<endl;
 //   cout<<"y1, 2, 3:  "<<mc1_yield.getValV()/heta_DY0_bottom_btag->Integral()<<" , "<<mc2_yield.getValV()/heta_DY0_charm_btag->Integral()<<" , "<<mc3_yield.getValV()/heta_DY0_light_btag->Integral()<<" , "<<endl;
    dframe->Draw();
    
   // heta_data_ctag->Draw();
  
   
    
    RooCategory sample("sample","sample") ;
    sample.defineType("sample1") ;
    sample.defineType("sample2") ;
    sample.defineType("sample3") ;
    RooDataHist combData("combData","combined data",x,Index(sample),Import("sample1",data),Import("sample2",datac));
    
    RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
    simPdf.addPdf(model1,"sample1") ;
    simPdf.addPdf(model2,"sample2") ;
  // simPdf.addPdf(model3,"sample3") ;
    simPdf.fitTo(combData) ;
  //  model3.fitTo(datano);
    simPdf.getParameters(combData)->Print("v") ;
    cout<<"sdfas:  "<<model1.coefList()<<endl;
    cout<<"model2:  "<<mean1.getValV()<<" , "<<mean2.getValV()<<" , "<<mean3.getValV()<<endl;
    cout<<"btag model2 res:  "<<mean1.getValV()*heta_DY0_bottom_btag->Integral() + mean2.getValV()*heta_DY0_charm_btag->Integral() + mean3.getValV()*heta_DY0_light_btag->Integral()<<"  data :  "<<heta_data_btag->Integral()<<endl;
    cout<<"ctag model2 res:  "<<mean1.getValV()*heta_DY0_bottom_ctag->Integral() + mean2.getValV()*heta_DY0_charm_ctag->Integral() + mean3.getValV()*heta_DY0_light_ctag->Integral()<<"  data :  "<<heta_data_ctag->Integral()<<endl;
    cout<<"notag model2 res:  "<<mean1.getValV()*heta_DY0_bottom_notag->Integral() + mean2.getValV()*heta_DY0_charm_notag->Integral() + mean3.getValV()*heta_DY0_light_notag->Integral()<<"  data :  "<<heta_data_notag->Integral()<<endl;
    
    cout<<"notag uncorrected b, c , light   "<<heta_DY0_bottom_notag->Integral()<<" , "<<heta_DY0_charm_notag->Integral()<<" , "<<heta_DY0_light_notag->Integral()<<endl;
    // Associate model with the physics state and model_ctl with the control state
  //  simPdf.addPdf(model,"physics") ;
  //  simPdf.addPdf(model_ctl,"control") ;
    
        }