#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 addMPI4amcnlo(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");
     //
//   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";
     if (hname=="j_bjet_Pt_binned")          hnameMPI="zb_jet_pT_MPI";
     else if (hname=="j_bjet_Y_binned")      {
	 hnameMPI="zb_jet_y_MPI";
	 scaleF = 2.*scaleF;
     }
     else if (hname=="j_bjet_Yboost_binned") hnameMPI="zb_yboost_MPI";
     //
     else if (hname=="j_bjet_DY20_binned")     hnameMPI="zb_dy_MPI";
     else if (hname=="j_bjet_ZdPhiF20_binned") {
	 hnameMPI="zb_dphi_MPI";
	 scaleF = 2.*scaleF;	
     }
     else if (hname=="j_bjet_ZdR620_binned")   hnameMPI="zb_dR_MPI";
     //
     else if (hname=="j_bjet_ZPt_binned")      hnameMPI="zb_Z_pT_MPI";
     else if (hname=="j_bjet_ZY_binned")       {
	 scaleF = 2.*scaleF;
	 hnameMPI="zb_Z_y_MPI";
     }
     //
     else if (hname=="j_2bjet_ZPtv3_binned")   hnameMPI="zbb_Z_pT_MPI";
     //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";
     }
     else if (hname=="j_2bjet_dR_binned")      hnameMPI="zbb_dR_MPI";
     else if (hname=="j_2bjet_Mbb_binned")     hnameMPI="zbb_M_MPI";

     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;
     }
     
     Double_t yn[100];
     Double_t eupyn[100];
     Double_t edwyn[100];
     
     int Npoints=gr0->GetN();
     for (int k=0; k<Npoints; ++k)
     {
	 Double_t dx = 2.*(gr0->GetX()[k]);
	 yn[k]  = gr0->GetY()[k];
	 yn[k] += scaleF*gr1->GetY()[k];
	 Double_t ey0=gr0->GetEYhigh()[k];
	 Double_t ey1=scaleF*gr1->GetEYhigh()[k];
	 eupyn[k]=sqrt(ey0*ey0+ey1*ey1);
	 ey0=gr0->GetEYlow()[k];
	 ey1=scaleF*gr1->GetEYlow()[k];
	 edwyn[k]=sqrt(ey0*ey0+ey1*ey1);
     }


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



     TGraphAsymmErrors* mergedGr = new TGraphAsymmErrors(Npoints, 
							 gr0->GetX(), yn, 
							 gr0->GetEXlow(), gr0->GetEXhigh(),
							 edwyn, eupyn);
     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);
     
     TCanvas* c1 = new TCanvas("merging mpi","Mering MPI");
     mergedGr->Draw("APEC");
     gr0->Draw("PECSAME");
     gr1->Draw("PESAME");
     

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

}