#include "TGraphAsymmErrors.h"
#include <iostream>
#include <iomanip>

Double_t getErr(TGraphAsymmErrors* gr, std::string upORdw="UP", int optErr=0){
     
     Double_t* grY    = gr->GetY();
     Double_t* grEYl  = gr->GetEYlow();
     Double_t* grEYh  = gr->GetEYhigh();
     Double_t* grX    = gr->GetX();
     Double_t* grEXl  = gr->GetEXlow();
     Double_t* grEXh  = gr->GetEXhigh();


     int Npoints = gr->GetN();
     Double_t dX;
     Double_t errh;
     Double_t errl;
     Double_t xsec = 0.;
     Double_t exsech = 0.;
     Double_t exsecl = 0.;

     double sgrEYl = 0.;
     for (int k=0; k<Npoints; ++k)
       {

	 dX = grEXl[k]+grEXh[k];
	 //std::cout<<" do we really come here ?? dx, y, eyl, eyh "<<dX<<" "<<grY[k]<<" "<<grEYl[k]<<" "<<grEYh[k]<<std::endl;
	 
	 
	 xsec   = xsec   + dX * grY  [k];
	 if (optErr==1) 
	 {
	     //std::cout<<" adding linearly for correlated errors"<<std::endl;
	     
	     exsech = exsech + dX * grEYh[k];
	     exsecl = exsecl + dX * grEYl[k];
	     //std::cout<<" do we come here "<<xsec<<" "<< exsecl<<" "<<exsech<<std::endl;
	 }
	 else if (optErr==0)
	 {
	     //std::cout<<" adding in quadrature "<<std::endl;
	     
	     exsech = exsech + dX * grEYh[k] * dX * grEYh[k];
	     exsecl = exsecl + dX * grEYl[k] * dX * grEYl[k];
	     sgrEYl = sgrEYl + grEYl[k];
	 }
	 
	 //std::cout<<" computing err: index k = "<<k<<" dx, grEYl/h = "<<dX<<" "<<grEYl  [k]<<"/"<<grEYh  [k]<<std::endl;
	 //std::cout<<" computing err: index k = "<<k<<" quadr sum EYl = "<<" "<<sgrEYl<<" my error low "<<exsecl<<std::endl;

       }
      if (optErr==0) 
      {
	  exsech = sqrt(exsech);
	  exsecl = sqrt(exsecl);	
	  sgrEYl = sqrt(sgrEYl);
      }
      
      //std::cout<<"END of computing err: quadr sum EYl = "<<sgrEYl<<" my error low "<<exsecl<<std::endl;
      //std::cout<<"===Err=== "<<xsec <<" - "<<exsecl<<" + "<<exsech<<std::endl;
      if      (upORdw=="UP") return exsech;
      return exsecl;
    
}
Double_t getSystErr(TGraphAsymmErrors* gr, std::string upORdw="UP", int optErr=1)
{
    return getErr(gr, upORdw, optErr);
}
Double_t getStatErr(TGraphAsymmErrors* gr, std::string upORdw="UP", int optErr=0)
{
    return getErr(gr, upORdw, optErr);
}


Double_t ComputeTotalCrossSection(TGraphAsymmErrors* gr){
     
    return ComputeTotalCrossSectionWithOpt(gr, 0);
}

Double_t ComputeTotalCrossSectionWithOpt(TGraphAsymmErrors* gr, int optErr=0){
     

     Double_t* grY    = gr->GetY();
     Double_t* grEYl  = gr->GetEYlow();
     Double_t* grEYh  = gr->GetEYhigh();
     Double_t* grX    = gr->GetX();
     Double_t* grEXl  = gr->GetEXlow();
     Double_t* grEXh  = gr->GetEXhigh();


     int Npoints = gr->GetN();
     Double_t dX;
     Double_t errh;
     Double_t errl;
     Double_t xsec = 0.;
     Double_t exsech = 0.;
     Double_t exsecl = 0.;

      for (int k=0; k<Npoints; ++k)
       {

	 dX = grEXl[k]+grEXh[k];
	 
	 xsec   = xsec   + dX * grY  [k];
	 if (optErr==1) 
	 {
	     //std::cout<<" adding linearly for correlated errors"<<std::endl;
	     
	     exsech = exsech + dX * grEYh[k];
	     exsecl = exsecl + dX * grEYl[k];
	 }
	 else if (optErr==0)
	 {
	     //std::cout<<" adding in quadrature "<<std::endl;
	     
	     exsech = exsech + dX * grEYh[k] * dX * grEYh[k];
	     exsecl = exsecl + dX * grEYl[k] * dX * grEYl[k];
	 }
	 
	 //std::cout<<" k = "<<k<<" dx, grY = "<<dX<<" "<<grY  [k]<<std::endl;

       }
      if (optErr==0) 
      {
	  exsech = sqrt(exsech);
	  exsecl = sqrt(exsecl);	  
      }
      
      std::cout<<"=== Xsec (fb) = "<<xsec <<" - "<<exsecl<<" + "<<exsech<<std::endl;
      return xsec;
}