#include "TFile.h"
#include "TGraphAsymmErrors.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TAxis.h"
#include <iostream>
#include <sstream>
#include <iomanip>
#include "AtlasUtils.h"
#include "AtlasStyle.C"
#include "AtlasLabels.C"

#include "ComputeTotalCrossSection.C"

#define saveAsRoot      true


void MPIsystFromData(std::string hname="j_bjet_Pt_binned")
{

  if (hname=="XXX"){
    std::cout<<"No histogram selected - possible choices are: \n";
    std::cout<<"\"j_bjet_Pt_binned\""<<std::endl;
    std::cout<<"\"j_bjet_Y_binned\""<<std::endl;
    std::cout<<"\"j_bjet_Yboost_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZY_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZYf_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZdR620_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZdPhiF20_binned\""<<std::endl;
    std::cout<<"\"j_bjet_DY20_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZY_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPtv3_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_dR_binned\""<<std::endl; 
    std::cout<<"\"j_2bjet_Mbb_binned\""<<std::endl;

    return;
  }

  //     TFile* file0 = new TFile("../aMCatNLO/Zbb_all.v3.Rivet.root");
  //     TFile* file0 = new TFile("../aMCatNLO/lecceProd/histo_oldSetNewCuts_MCFMscales_5fZbMPI_1_6samples.v3.Rivet.root"); /// MPI to be scaled with 0.5
  //  TFile* file0 = new TFile("../aMCatNLO/lecceProd/histo_oldSetNewCuts_MCFMscales_5fZj_1_4samples.v3.Rivet.root"); /// MPI to be scaled with 0.5 x (Z_noj - Zj)/Zinclusive = 668.68/868.17
  //  TFile* file0 = new TFile("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales.v3.Rivet.root"); /// MPi to be scaled by 0.5 x (Z_noj - Zj)/Zinclusive = 668.68/868.17
  //     TFile* file0 = new TFile("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_1to5.v3.Rivet.root"); /// MPi to be scaled by 0.5 x (Z_noj - Zj)/Zinclusive = 668.68/868.17
     TFile* file0 = new TFile("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_1to10.v3.Rivet.root"); /// MPi to be scaled by 0.5 x (Z_noj - Zj)/Zinclusive = 668.68/868.17
  //   TFile* file0 = new TFile("../aMCatNLO/LouvainProd/histo_Louvain_4fZbb_MCFMscales.v3.Rivet.root");     /// MPI to be scaled by 0.5 
     TFile* file1 = new TFile("../aMCatNLO/MPI_v3_fromGavin_Aug22.root");
     TFile* fMCFM;
     if (hname.substr(0,3)=="j_b"){
       fMCFM = new TFile("../MCFMall/ZbMCFMgraphMSTW2008_aug13.root");
     }
     else if (hname.substr(0,3)=="j_2"){
       fMCFM  = new TFile("../MCFMall/ZbbMCFMgraphMSTW2008_aug13.root");
     }
     //
//   KEY: TH1D	zb_Z_pT_MPI;1	Z pT
//   KEY: TH1D	zb_Z_y_MPI;1	Z rapidity
//   KEY: TH1D	zb_dR_MPI;1	#Delta R
//   KEY: TH1D	zb_dphi_MPI;1	#Delta #phi
//   KEY: TH1D	zb_dy_MPI;1	#Delta y
//   KEY: TH1D	zb_yboost_MPI;1	y_{boost}
//   KEY: TH1D	zb_jet_pT_MPI;1	jet pT
//   KEY: TH1D	zb_jet_y_MPI;1	jet rapidity
//   KEY: TH1D	zbb_M_MPI;1	Mass(b,b)
//   KEY: TH1D	zbb_Z_pT_MPI;1	Z pT
//   KEY: TH1D	zbb_Z_y_MPI;1	Z rapidity
//   KEY: TH1D	zbb_dR_MPI;1	#Delta R(b,b)
     //

     Double_t scaleF = 0.5;
     //     Double_t scaleF = 0.5 *  668.68/868.17;

     TGraphAsymmErrors* gr0 = (TGraphAsymmErrors*)file0->Get(hname.c_str());
     std::string hnameMPI="XXX";
     std::string hnameMCFM="XXX";
     if (hname=="j_bjet_Pt_binned")          {
       hnameMPI="zb_jet_pT_MPI";
       hnameMCFM= "Pt_gae_";

     }
     else if (hname=="j_bjet_Y_binned")      {
	 hnameMPI="zb_jet_y_MPI";
	 scaleF = 2.*scaleF;
	 hnameMCFM= "Y_gae_";
     }
     else if (hname=="j_bjet_Yboost_binned") {
       hnameMPI="zb_yboost_MPI";
       hnameMCFM="Yboost_gae_";
     }
     //
     else if (hname=="j_bjet_DY20_binned")     {
       hnameMPI="zb_dy_MPI";
       hnameMCFM = "ZdY_gae_";
     }
     else if (hname=="j_bjet_ZdPhiF20_binned") {
	 hnameMPI="zb_dphi_MPI";
	 scaleF = 2.*scaleF;
	 hnameMCFM = "ZdPhi_gae_";
     }
     else if (hname=="j_bjet_ZdR620_binned")   {
       hnameMPI="zb_dR_MPI";
       hnameMCFM = "ZdR_gae_";
     }
     //
     else if (hname=="j_bjet_ZPt_binned")      {
       hnameMPI="zb_Z_pT_MPI";
       hnameMCFM = "ZPt_gae_";
     }
     else if (hname=="j_bjet_ZY_binned")       {
	 scaleF = 2.*scaleF;
	 hnameMPI="zb_Z_y_MPI";
	 hnameMCFM = "ZY_gae_";
     }
     //
     else if (hname=="j_2bjet_ZPtv3_binned")   {
       hnameMPI="zbb_Z_pT_MPI";
       hnameMCFM = "ZPT_gae_";
     }
     //else if (hname=="j_2bjet_ZPt_binned")     hnameMPI="zbb_Z_pT_MPI";
     else if (hname=="j_2bjet_ZY_binned")      {
	 scaleF = 2.*scaleF;
	 hnameMPI="zbb_Z_y_MPI";
	 hnameMCFM = "ZY_gae_";
     }
     else if (hname=="j_2bjet_dR_binned")      {
       hnameMPI="zbb_dR_MPI";
       hnameMCFM = "bbdR_gae_";
     }
     else if (hname=="j_2bjet_Mbb_binned")     {
       hnameMPI="zbb_M_MPI";
       hnameMCFM = "bbM_gae_";
     }

     std::cout<<"**************************** Trying to add "<<hname<<" and "<<hnameMPI<<std::endl;
     std::cout<<" MPI from Alpgen samples scaled by "<<scaleF<<std::endl;
     
     TGraphAsymmErrors* gr1;
     TH1F* htmp = (TH1F*)file1->Get(hnameMPI.c_str());
     if (htmp != NULL) gr1 = new TGraphAsymmErrors(htmp);
     else gr1=NULL;
     if (gr0==NULL || gr1==NULL) {
	 std::cout<<" some histogram / graph not found for hname=<"<<hname<<"> hnameMPI=<"<<hnameMPI<<">"<<std::endl;
	 file1->ls();
	 return;
     }
     hnameMCFM=hnameMCFM+"stat";
     std::cout<<" looking for gr with name <"<<hnameMCFM<<">"<<std::endl;
     TGraphAsymmErrors* grMCFM   = (TGraphAsymmErrors*)fMCFM->Get(hnameMCFM.c_str());
     if (grMCFM==NULL) return;
     
     Double_t yn[100];
     Double_t eupyn[100];
     Double_t edwyn[100];
     Double_t ys[100];
     Double_t eys[100];
     
     int Npoints=gr0->GetN();
     for (int k=0; k<Npoints; ++k)
     {
	 Double_t dx = 2.*(gr0->GetX()[k]);
	 yn[k]  = 0.*gr0->GetY()[k];
	 yn[k] += 1000.*scaleF*gr1->GetY()[k];
	 Double_t ey0=gr0->GetEYhigh()[k];
	 Double_t ey1=scaleF*gr1->GetEYhigh()[k];
	 eupyn[k]=1000.*ey1;
	 ey1=scaleF*gr1->GetEYlow()[k];
	 edwyn[k]=1000.*ey1;

	 ys[k]  = 0.5 * yn[k]/(grMCFM->GetY()[k]);
	 eys[k] = 0.;

     }


     TFile* foutRoot;
     if (saveAsRoot) 
       {
	 std::string fn = "MergedHistos.root";
	 foutRoot = new TFile(fn.c_str(),"UPDATE");
       }



     TGraphAsymmErrors* grnMCFM = new TGraphAsymmErrors(Npoints, 
							gr0->GetX(), grMCFM->GetY(), 
						        gr0->GetEXlow(), gr0->GetEXhigh(),
							grMCFM->GetEYlow(), grMCFM->GetEYhigh());
     TGraphAsymmErrors* mergedGr = new TGraphAsymmErrors(Npoints, 
							 gr0->GetX(), yn, 
							 gr0->GetEXlow(), gr0->GetEXhigh(),
							 edwyn, eupyn);


     TGraphAsymmErrors* systMPIDataComp = new TGraphAsymmErrors(Npoints, 
							 gr0->GetX(), ys, 
							 gr0->GetEXlow(), gr0->GetEXhigh(),
							 eys, eys);

     std::cout<<"Merged "<<std::endl;
     Double_t xsmerged = ComputeTotalCrossSectionWithOpt(mergedGr, 0);
     
     std::cout<<"W/O MPI "<<std::endl;
     Double_t xsorig = ComputeTotalCrossSectionWithOpt(gr0, 0);
     
     std::cout<<"MPI "<<std::endl;
     Double_t xsmpi = ComputeTotalCrossSectionWithOpt(gr1, 0);
     
     std::cout<<"MPI syst"<<std::endl;
     Double_t xssyst = ComputeTotalCrossSectionWithOpt(systMPIDataComp, 0);
     
     std::cout<<"MCFM xsec "<<std::endl;
     Double_t xsMFM = ComputeTotalCrossSectionWithOpt(grnMCFM,0);

     TCanvas* c1 = new TCanvas("MCFM","MCFM MPI");
     grnMCFM->Draw("AC");
     TCanvas* c2 = new TCanvas("merging mpi","Mering MPI");
     mergedGr->Draw("APE");
     //     gr0->Draw("PECSAME");
     //     gr1->Draw("PESAME");
     TCanvas* c3 = new TCanvas("mpi syst"," MPI syst");
     systMPIDataComp->Draw("APE");

     if (saveAsRoot)
       {
	   mergedGr->Write(hname.c_str());
	   systMPIDataComp->Write(("MPIsyst_"+hname).c_str());
       }
     
     if (saveAsRoot) 
       {
	 foutRoot->Close();
       }

}