#include "TFile.h"
#include "TGraphAsymmErrors.h"

grAdd4MCNLO(std::string hname="XXX")
{
  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_ZdR620_binned\""<<std::endl;
    std::cout<<"\"j_bjet_DY20_binned\""<<std::endl;
    std::cout<<"\"j_bjet_ZdPhiF20_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZPtv3_binned\""<<std::endl;
    std::cout<<"\"j_2bjet_ZY_binned\""<<std::endl; 
    std::cout<<"\"j_2bjet_dR_binned\""<<std::endl; 
    std::cout<<"\"j_2bjet_Mbb_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_Pt_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_Y_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_Yboost_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_ZY_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_ZdR620_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_DY20_binned\""<<std::endl;
    std::cout<<"\"j2m_bjet_ZdPhiF20_binned\""<<std::endl;
    std::cout<<"\"j2m_2bjet_ZPt_binned\""<<std::endl;
    std::cout<<"\"j2m_2bjet_ZPtv3_binned\""<<std::endl;
    std::cout<<"\"j2m_2bjet_ZY_binned\""<<std::endl; 
    std::cout<<"\"j2m_2bjet_dR_binned\""<<std::endl; 
    std::cout<<"\"j2m_2bjet_Mbb_binned\""<<std::endl;

    return -1;
  }
  

// TFile f0("Rivet_23471000.root");
// TFile f1("Rivet_23471234.root");
// TFile f2("Rivet_23471774.root");
// TFile f3("Rivet_23471884.root");
// TFile f4("Rivet_23471999.root");
// TFile f5("Rivet_23474444.root");
// TFile f6("Rivet_23478944.root");
// TFile f7("Rivet_89441234.root");
// TFile f8("Rivet_23472347.root");
// TFile f9("Rivet_89551234.root");
//////////////////// 
// TFile f0("/einstein2/chiodini/atlas/rivet/g2_aug_13_23471000/Rivet.root");
// TFile f1("/einstein2/chiodini/atlas/rivet/g2_aug_13_23471234/Rivet.root");
// TFile f2("/einstein2/chiodini/atlas/rivet/g2_aug_13_23471774/Rivet.root");
// TFile f3("/einstein2/chiodini/atlas/rivet/g2_aug_13_23471884/Rivet.root");
// TFile f4("/einstein2/chiodini/atlas/rivet/g2_aug_13_23471999/Rivet.root");
// TFile f5("/einstein2/chiodini/atlas/rivet/g2_aug_13_23472347/Rivet.root");
// TFile f6("/einstein2/chiodini/atlas/rivet/g2_aug_13_23474444/Rivet.root");
// TFile f7("/einstein2/chiodini/atlas/rivet/g2_aug_13_23478944/Rivet.root");
// TFile f8("/einstein2/chiodini/atlas/rivet/g2_aug_13_89441234/Rivet.root");
// TFile f9("/einstein2/chiodini/atlas/rivet/g2_aug_13_89551234/Rivet.root");

/*
TFile f0("/einstein2/chiodini/atlas/rivet/g3_aug_13_23471000/Rivet.root");
TFile f1("/einstein2/chiodini/atlas/rivet/g3_aug_13_23471234/Rivet.root");
TFile f2("/einstein2/chiodini/atlas/rivet/g3_aug_13_23471774/Rivet.root");
TFile f3("/einstein2/chiodini/atlas/rivet/g3_aug_13_23471884/Rivet.root");
TFile f4("/einstein2/chiodini/atlas/rivet/g3_aug_13_23471999/Rivet.root");
TFile f5("/einstein2/chiodini/atlas/rivet/g3_aug_13_23472347/Rivet.root");
TFile f6("/einstein2/chiodini/atlas/rivet/g3_aug_13_23474444/Rivet.root");
TFile f7("/einstein2/chiodini/atlas/rivet/g3_aug_13_23478944/Rivet.root");
TFile f8("/einstein2/chiodini/atlas/rivet/g3_aug_13_89441234/Rivet.root");
TFile f9("/einstein2/chiodini/atlas/rivet/g3_aug_13_89551234/Rivet.root");
*/
TFile f0("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales.v3.Rivet.root");
TFile f1("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_02.v3.Rivet.root");
TFile f2("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_03.v3.Rivet.root");
TFile f3("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_04.v3.Rivet.root");
TFile f4("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_05.v3.Rivet.root");
TFile f5("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_06.v3.Rivet.root");
TFile f6("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_07.v3.Rivet.root");
TFile f7("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_08.v3.Rivet.root");
TFile f8("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_09.v3.Rivet.root");
TFile f9("../aMCatNLO/LouvainProd/histo_Louvain_5fZj_MCFMscales_10.v3.Rivet.root");

/*
/einstein2/chiodini/atlas/rivet/

e si chiamano

   -rw-r--r-- 1 chiodini gruppo1 55969 Sep 11 09:18 g3_aug_13_23471000/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56089 Sep 11 09:18 g3_aug_13_23471234/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56073 Sep 11 09:18 g3_aug_13_23471774/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56030 Sep 11 09:18 g3_aug_13_23471884/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56114 Sep 11 09:18 g3_aug_13_23471999/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56090 Sep 11 09:18 g3_aug_13_23472347/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56065 Sep 11 09:18 g3_aug_13_23474444/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 55957 Sep 11 09:19 g3_aug_13_23478944/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 56101 Sep 11 09:19 g3_aug_13_89441234/Rivet.root
   -rw-r--r-- 1 chiodini gruppo1 55946 Sep 11 09:19 g3_aug_13_89551234/Rivet.root
*/

TGraphAsymmErrors *g0 = (TGraphAsymmErrors*)f0.Get(hname.c_str());
TGraphAsymmErrors *g1 = (TGraphAsymmErrors*)f1.Get(hname.c_str());
TGraphAsymmErrors *g2 = (TGraphAsymmErrors*)f2.Get(hname.c_str());
TGraphAsymmErrors *g3 = (TGraphAsymmErrors*)f3.Get(hname.c_str());
TGraphAsymmErrors *g4 = (TGraphAsymmErrors*)f4.Get(hname.c_str());
TGraphAsymmErrors *g5 = (TGraphAsymmErrors*)f5.Get(hname.c_str());
TGraphAsymmErrors *g6 = (TGraphAsymmErrors*)f6.Get(hname.c_str());
TGraphAsymmErrors *g7 = (TGraphAsymmErrors*)f7.Get(hname.c_str());
TGraphAsymmErrors *g8 = (TGraphAsymmErrors*)f8.Get(hname.c_str());
TGraphAsymmErrors *g9 = (TGraphAsymmErrors*)f9.Get(hname.c_str());
// TGraphAsymmErrors *g5 = (TGraphAsymmErrors*)f5.Get(hname.c_str());
// TGraphAsymmErrors *g6 = (TGraphAsymmErrors*)f6.Get(hname.c_str());
// TGraphAsymmErrors *g7 = (TGraphAsymmErrors*)f7.Get(hname.c_str());
// TGraphAsymmErrors *g8 = (TGraphAsymmErrors*)f8.Get(hname.c_str());
// TGraphAsymmErrors *g9 = (TGraphAsymmErrors*)f9.Get(hname.c_str());
Double_t  Ygr[100];
Double_t EYgr[100];
int Npoints = g0->GetN();

std::cout<<" ////////////////////////////// averaging distribution <"<<hname<<">"<<std::endl;

for (int k=0; k<Npoints; ++k)
    {
	Ygr[k]   = 0.;
	EYgr[k]  = 0.;	
	std::cout<<" ********************************** bin k "<<k<<std::endl;
	
	//	for (int j=0, j<5, ++j)
	//	{
	    // media pesata 
	    double sum_weight = 0.;
	    double weight = 0.;

	    weight = 1./(g0->GetEYlow()[k]*g0->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g0->GetY()[k]);
	    sum_weight = sum_weight + weight;
	    std::cout<<"sameple j=0  ----  y, ey = "<<g0->GetY()[k]<<" - "<<g0->GetEYlow()[k]<<" + "<<g0->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g1->GetEYlow()[k]*g1->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g1->GetY()[k]);
	    sum_weight = sum_weight + weight;
	    std::cout<<"sameple j=1  ----  y, ey = "<<g1->GetY()[k]<<" - "<<g1->GetEYlow()[k]<<" + "<<g1->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g2->GetEYlow()[k]*g2->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g2->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=2  ----  y, ey = "<<g2->GetY()[k]<<" - "<<g2->GetEYlow()[k]<<" + "<<g2->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g3->GetEYlow()[k]*g3->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g3->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=3  ----  y, ey = "<<g3->GetY()[k]<<" - "<<g3->GetEYlow()[k]<<" + "<<g3->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g4->GetEYlow()[k]*g4->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g4->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=4  ----  y, ey = "<<g4->GetY()[k]<<" - "<<g4->GetEYlow()[k]<<" + "<<g4->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g5->GetEYlow()[k]*g5->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g5->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=5  ----  y, ey = "<<g5->GetY()[k]<<" - "<<g5->GetEYlow()[k]<<" + "<<g5->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g6->GetEYlow()[k]*g6->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g6->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=6  ----  y, ey = "<<g6->GetY()[k]<<" - "<<g6->GetEYlow()[k]<<" + "<<g6->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g7->GetEYlow()[k]*g7->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g7->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=7  ----  y, ey = "<<g7->GetY()[k]<<" - "<<g7->GetEYlow()[k]<<" + "<<g7->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g8->GetEYlow()[k]*g8->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g8->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=8  ----  y, ey = "<<g8->GetY()[k]<<" - "<<g8->GetEYlow()[k]<<" + "<<g8->GetEYhigh()[k]<<std::endl;

	    weight = 1./(g9->GetEYlow()[k]*g9->GetEYlow()[k]);	    
	    Ygr[k]  = Ygr[k] + weight * (g9->GetY()[k]);
	    sum_weight = sum_weight + weight;	    
	    std::cout<<"sameple j=9  ----  y, ey = "<<g9->GetY()[k]<<" - "<<g9->GetEYlow()[k]<<" + "<<g9->GetEYhigh()[k]<<std::endl;

	    Ygr[k]  = Ygr[k]  / sum_weight;
	    EYgr[k] = sqrt(1./sum_weight);
	    std::cout<<"AVERAGE  ----  y, ey = "<<Ygr[k]<<" +/- "<<EYgr[k]<<std::endl;

    }
TGraphAsymmErrors *gNew = new TGraphAsymmErrors(Npoints, g0->GetX(), Ygr, g0->GetEXhigh(), g0->GetEXlow(), EYgr, EYgr);

gNew->SetLineColor(2);

gNew->Draw("APE");
g0->Draw("PESAME");
g1->Draw("PESAME");
g2->Draw("PESAME");
g3->Draw("PESAME");
g4->Draw("PESAME");
gNew->Draw("PESAME");



 TFile *f = new TFile("hTotal.root","UPDATE");
 gNew->Write(hname.c_str());
 f->Close();
 
}