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

Double_t show_integrate(std::string gName, std::string gavinn, TFile* iFile, TFile* jFile,  TFile* kFile,  
			TFile* iShTotFile,  TFile* iShZbFile,  TFile* iShnoZbFile)
{
  
    TH1D* gg = (TH1D*)kFile->Get(gavinn.c_str()); 
    TGraphAsymmErrors* gMPInoZbJimmyOr = new TGraphAsymmErrors(gg);
  Double_t yTot[1000];
  Double_t EyTot[1000];

  double scaleF = 1.;
  if (gName=="j_bjet_ZdPhiF20_binned" || gName=="j_bjet_Y_binned" || gName=="j_bjet_ZY_binned" || gName=="j_2bjet_ZY_binned")
  {
      scaleF = 2.;
  }
  for (int k=0;k<gMPInoZbJimmyOr->GetN(); ++k)
  {
      yTot[k]=scaleF*gMPInoZbJimmyOr->GetY()[k];
      EyTot[k]=scaleF*gMPInoZbJimmyOr->GetEYlow()[k];
  }
  TGraphAsymmErrors* gMPInoZbJimmy = new TGraphAsymmErrors(gMPInoZbJimmyOr->GetN(), gMPInoZbJimmyOr->GetX(), yTot, 
					 gMPInoZbJimmyOr->GetEXlow(), gMPInoZbJimmyOr->GetEXhigh(), 
					 EyTot, EyTot);
  

  TGraphAsymmErrors* gMCNLO_MPI   = (TGraphAsymmErrors*)iFile->Get(gName.c_str());  
  TGraphAsymmErrors* gMCNLO_noMPI = (TGraphAsymmErrors*)jFile->Get(gName.c_str());  
  TGraphAsymmErrors* gMPIZbHerwigPP; 
  TGraphAsymmErrors* gMPITotHerJimmy; 
  Double_t yZbHerwigPP[1000];
  Double_t EyZbHerwigPP[1000];
  Double_t yDiffTot[1000];
  Double_t EyDiffTot[1000];
  Double_t yDiffZb[1000];
  Double_t EyDiffZb[1000];
  Double_t yDiffnoZb[1000];
  Double_t EyDiffnoZb[1000];


  TGraphAsymmErrors* gMPI_Sh_noZb = (TGraphAsymmErrors*)iShnoZbFile->Get(gName.c_str());  
  TGraphAsymmErrors* gMPI_Sh_tot = (TGraphAsymmErrors*)iShTotFile->Get(gName.c_str());  
  TGraphAsymmErrors* gMPI_Sh_Zb  = (TGraphAsymmErrors*)iShZbFile->Get(gName.c_str());  

  if (gMPI_Sh_tot==NULL) 
  {
      std::cout<<" not found "<<std::endl;
      
      return;
  }
  
	  
  for (int k=0;k<gMCNLO_MPI->GetN(); ++k)
  {
      yZbHerwigPP[k]  = gMCNLO_MPI->GetY()[k]-gMCNLO_noMPI->GetY()[k];
      EyZbHerwigPP[k] = sqrt(   (gMCNLO_MPI->GetEYlow()[k])*(gMCNLO_MPI->GetEYlow()[k])
         			+(gMCNLO_noMPI->GetEYlow()[k])*(gMCNLO_noMPI->GetEYlow()[k])  );
      yTot[k] = yZbHerwigPP[k]+gMPInoZbJimmy->GetY()[k];
      EyTot[k] = sqrt( EyZbHerwigPP[k]*EyZbHerwigPP[k] + (gMPInoZbJimmy->GetEYlow()[k])*(gMPInoZbJimmy->GetEYlow()[k]) );


      yDiffTot[k] = yTot[k]-gMPI_Sh_tot->GetY()[k];
      EyDiffTot[k] = sqrt(EyTot[k]*EyTot[k] - (gMPI_Sh_tot->GetEYlow()[k])*(gMPI_Sh_tot->GetEYlow()[k]));

      yDiffZb[k] = yZbHerwigPP[k]-gMPI_Sh_Zb->GetY()[k];
      EyDiffZb[k] = sqrt(EyZbHerwigPP[k]*EyZbHerwigPP[k] - (gMPI_Sh_Zb->GetEYlow()[k])*(gMPI_Sh_Zb->GetEYlow()[k]));
      
      yDiffnoZb[k] = gMPInoZbJimmy->GetY()[k]-gMPI_Sh_Zb->GetY()[k];
      EyDiffnoZb[k] = sqrt((gMPInoZbJimmy->GetEYlow()[k])*(gMPInoZbJimmy->GetEYlow()[k]) + (gMPI_Sh_noZb->GetEYlow()[k])*(gMPI_Sh_noZb->GetEYlow()[k]));
      

  }
  gMPIZbHerwigPP = new TGraphAsymmErrors(gMCNLO_MPI->GetN(), gMCNLO_MPI->GetX(), yZbHerwigPP, 
					 gMCNLO_MPI->GetEXlow(), gMCNLO_MPI->GetEXhigh(), 
					 EyZbHerwigPP, EyZbHerwigPP);
  
  gMPITotHerJimmy = new TGraphAsymmErrors(gMCNLO_MPI->GetN(), gMCNLO_MPI->GetX(), yTot, 
					 gMCNLO_MPI->GetEXlow(), gMCNLO_MPI->GetEXhigh(), 
					 EyTot, EyTot);
  TGraphAsymmErrors* gDiffTot = new TGraphAsymmErrors(gMCNLO_MPI->GetN(), gMCNLO_MPI->GetX(), yDiffTot, 
						      gMCNLO_MPI->GetEXlow(), gMCNLO_MPI->GetEXhigh(), 
						      EyDiffTot, EyDiffTot);
  TGraphAsymmErrors* gDiffZb  = new TGraphAsymmErrors(gMCNLO_MPI->GetN(), gMCNLO_MPI->GetX(), yDiffZb, 
						      gMCNLO_MPI->GetEXlow(), gMCNLO_MPI->GetEXhigh(), 
						      EyDiffZb, EyDiffZb);
  TGraphAsymmErrors* gDiffnoZb = new TGraphAsymmErrors(gMCNLO_MPI->GetN(), gMCNLO_MPI->GetX(), yDiffnoZb, 
						      gMCNLO_MPI->GetEXlow(), gMCNLO_MPI->GetEXhigh(), 
						      EyDiffnoZb, EyDiffnoZb);

  TCanvas* c = new TCanvas(gName.c_str(),gName.c_str());
  c->Divide(2,2);
  
  c->cd(1);
  gMCNLO_MPI->Draw("AEPC");
  gMCNLO_noMPI->SetLineColor(2);
  gMCNLO_noMPI->Draw("EPCSAME");
  
  gMPI_Sh_tot->SetLineColor(3);  
  gMPI_Sh_Zb->SetLineColor(3);  
  gMPI_Sh_noZb->SetLineColor(3);  
  gMPI_Sh_tot->SetMarkerColor(3);  
  gMPI_Sh_Zb->SetMarkerColor(3);  
  gMPI_Sh_noZb->SetMarkerColor(3);  

  c->cd(2);
  gMPITotHerJimmy->Draw("AEPC");
  gMPI_Sh_tot->Draw("EPCSAME");
  myText(0.6,0.8,1,"Total MPI");
  
  c->cd(3);	
  gMPInoZbJimmy->Draw("AEPC");
  gMPI_Sh_noZb->Draw("EPCSAME");
  myText(0.6,0.8,1,"Z+X MPI");

  c->cd(4);	
  gMPIZbHerwigPP->Draw("AEPC");
  gMPI_Sh_Zb->Draw("EPCSAME");
  myText(0.6,0.8,1,"Zb+X MPI");

  c->SaveAs(("MPI_"+gName+".eps").c_str());
  //  c->SaveAs(("MPI_"+gName+".eps").c_str());

  std::cout<<"****************************** "<<gName<<" *********TOTAL MPI*****************"<<std::endl;
  std::cout<<"*****ALPGEN************************* "<<std::endl;
  Double_t xsAl  = ComputeTotalCrossSection(gMPITotHerJimmy);
  std::cout<<"*****SHERPA************************* "<<std::endl;
  Double_t xsSh  = ComputeTotalCrossSection(gMPI_Sh_tot);
  std::cout<<"*****ALPGEN-SHERPA****************** "<<std::endl;	
  Double_t xsDiffTot  = ComputeTotalCrossSection(gDiffTot);
 
  std::cout<<"****************************** "<<gName<<" *********nonZb MPI****************"<<std::endl;
  std::cout<<"*****ALPGEN************************* "<<std::endl;
  Double_t xsAlnoZb   = ComputeTotalCrossSection(gMPInoZbJimmy);
  std::cout<<"*****SHERPA************************* "<<std::endl;
  Double_t xsShnoZb   = ComputeTotalCrossSection(gMPI_Sh_noZb);
  std::cout<<"*****ALPGEN-SHERPA****************** "<<std::endl;	
  Double_t xsDiffnoZb = ComputeTotalCrossSection(gDiffnoZb);

  std::cout<<"****************************** "<<gName<<" ********* Zb MPI *****************"<<std::endl;
  std::cout<<"*****ALPGEN************************* "<<std::endl;
  Double_t xsAlZb     = ComputeTotalCrossSection(gMPIZbHerwigPP);
  std::cout<<"*****SHERPA************************* "<<std::endl;
  Double_t xsShZb   = ComputeTotalCrossSection(gMPI_Sh_Zb);
  std::cout<<"*****ALPGEN-SHERPA****************** "<<std::endl;	
  Double_t xsDiffZb   = ComputeTotalCrossSection(gDiffZb);

  
}

void CompareMPI(std::string aa="pippo")
{
    TFile* kFile = new TFile("../aMCatNLO/MPI_v3_fromGavin_Aug22.root");
    TFile* iFile = new TFile("../aMCatNLO/lecceProd/histo_oldSetNewCuts_MCFMscales_5fZbMPI_1_6samples.v3.Rivet.root");
    TFile* jFile = new TFile("../aMCatNLO/lecceProd/histo_oldSetNewCuts_MCFMscales_5fZbnoMPI_1_6samples.v3.Rivet.root");

    TFile* iShnoZbFile  = new TFile("../SherpaPythia/hMPI_Sherpa_ZlandZcME.root");
    TFile* iShTotFile = new TFile("../SherpaPythia/hMPI_Sherpa_ZlandZcandZbME.root");
    TFile* iShZbFile  = new TFile("../SherpaPythia/hMPI_Sherpa_ZbME.root");
    
//    TFile* jFile = new TFile("/afs/le.infn.it/project/atlasewk/Generators/aMC_at_NLO/lecceProd/oldsetup13Aug/histo_oldSet5fZj.v3.Rivet.root");
    
    kFile->ls();
    
    gROOT->LoadMacro("AtlasUtils.C");
  std::string gname; 
  std::string gavinn;
  Double_t xs = 0.; 
  Double_t xs1 = 0.;
  TH1D* gg;

  std::cout<<"First  histo file = "<<iFile->GetName()<<std::endl;
  std::cout<<"Second histo file = "<<jFile->GetName()<<std::endl;
  
  gname = "j_bjet_Pt_binned";
  gavinn = "zb_jet_pT_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);
  //xs = show_integrate(gname.c_str(), jFile);
  //xs = show_integrate(gname.c_str(), kFile);
  //xs = ComputeTotalCrossSection(ggg);
  //return;
  

  gname = "j_bjet_Y_binned";
  gavinn = "zb_jet_y_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);
  //return;
  

  gname = "j_bjet_Yboost_binned";
  gavinn = "zb_yboost_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);
  //return;
  
  
  gname = "j_bjet_ZPt_binned";
  gavinn = "zb_Z_pT_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);

  gname = "j_bjet_ZY_binned";
  gavinn = "zb_Z_y_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);
  
  gname = "j_bjet_ZdPhiF20_binned";
  gavinn = "zb_dphi_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);

  gname = "j_bjet_ZdR620_binned";
  gavinn = "zb_dR_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);

  gname = "j_bjet_DY20_binned";
  gavinn = "zb_dy_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);

  gname = "j_2bjet_ZPtv3_binned";
  gavinn = "zbb_Z_pT_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);
  
  gname = "j_2bjet_ZY_binned";
  gavinn = "zbb_Z_y_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);

  gname = "j_2bjet_dR_binned";
  gavinn = "zbb_dR_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);

  gname = "j_2bjet_Mbb_binned";
  gavinn = "zbb_M_MPI";
  xs = show_integrate(gname, gavinn, iFile, jFile, kFile, iShTotFile, iShZbFile, iShnoZbFile);


  return;  
}