#define BCR_PDF_cxx

#include <BCR_PDF.h>
TString chNum[10]={"0","1","2","3","4","5","6","7","8","9"};
bool debug=kFALSE;


BCR_PDF::BCR_PDF(const BCR_PDF &other) : TNamed(other)
{
  for (int lg=0; lg<n_grid; lg++) {
    for (int lx=0; lx<n_bin_x; lx++) 
      xl_bin[lg][lx] = other.xl_bin[lg][lx];
    for (int lq2=0; lq2<n_bin_q2; lq2++) 
      q2l_bin[lg][lq2] = other.q2l_bin[lg][lq2];
    
    for (int lx=0; lx<n_bin_x; lx++) 
      for (int lq2=0; lq2<n_bin_q2; lq2++) 
	f[lg][lx][lq2] = other.f[lg][lx][lq2];    
  }
}

BCR_PDF & BCR_PDF::operator=(const BCR_PDF & other)
{
  if (this == &other) {
    cout << " assignment to itself ?" << endl;
    return *this;
  }

  TNamed::operator=(other); 
  for (int lg=0; lg<n_grid; lg++) {
    for (int lx=0; lx<n_bin_x; lx++) 
      xl_bin[lg][lx] = other.xl_bin[lg][lx];
    for (int lq2=0; lq2<n_bin_q2; lq2++) 
      q2l_bin[lg][lq2] = other.q2l_bin[lg][lq2];
    
    for (int lx=0; lx<n_bin_x; lx++) 
      for (int lq2=0; lq2<n_bin_q2; lq2++) 
	f[lg][lx][lq2] = other.f[lg][lx][lq2];    
  }
  
  return *this;
}

void BCR_PDF::Reset()
{
  for (int lg=0; lg<n_grid; lg++) {
    for (int lx=0;lx<n_bin_x;lx++)
      xl_bin[lg][lx] = 0.;
    
    for (int lq2=0;lq2<n_bin_q2;lq2++)
      q2l_bin[lg][lq2] = 0;
    
    for (int lq2=0;lq2<n_bin_q2;lq2++)
      for (int lx=0;lx<n_bin_x;lx++) 
	f[lg][lx][lq2] = 0.;      
  }
}
void BCR_PDF::ReadGrids(TString chPDF,TString chPath,Int_t iOff)
{
  TString chFile;
  for (int lg=0; lg<n_grid; lg++) {

    if (lg==0) {
      Int_t icycle = iOff*25+lg+1;
      TString chCycle = GetCycleString(icycle);
      chFile = chPath +  "/" +chPDF + chCycle + ".data";
    } else if (lg==1) {
      Int_t icycle = iOff*25+lg;
      TString chCycle = GetCycleString(icycle);
      chFile = chPath + "/" + "b_" +chPDF + chCycle + ".data";
    } else {
      Int_t icycle = iOff*25+lg;
      TString chCycle = GetCycleString(icycle);
      chFile = chPath + "/" + chPDF + chCycle + ".data";
    } 
	
    if (debug) cout << chFile << endl;

    ifstream fin;
    fin.open(chFile);
    for (int lx=0;lx<n_bin_x;lx++) {
      fin >> xl_bin[lg][lx];
      //    cout << lx << ": " << xl_bin[lx] << endl;
    }
    for (int lq2=0;lq2<n_bin_q2;lq2++) {
      fin >> q2l_bin[lg][lq2];
      //    cout << lq2 << ": " << q2l_bin[lq2] << endl;
      q2l_bin[lg][lq2] = TMath::Log(q2l_bin[lg][lq2]);
    }
    
    for (int lq2=0;lq2<n_bin_q2;lq2++) {
      for (int lx=0;lx<n_bin_x-1;lx++) {
	fin >> f[lg][lx][lq2];      
	//      cout << lx << ", " << lq2 << ": " << f[lx][lq2] << endl;
      }
    }
    fin.close();
    
  }
}


void BCR_PDF::ReadZeusKN(TString chPDF,TString chPath)
{
  TString chFile = chPath 
    +  "/kn_zeusxs4bcr." +chPDF + ".out";

  Float_t xold=-1., q2old=-1.;
  Int_t ic=-1,ix=-1,iq2=-1;
  Float_t q2=0.,xx=0.,ffff=0.,fmax=0.,fmin=0.;
  FILE * fp = fopen(chFile,"r");
  while (fscanf(fp,"%f %f %f %f %f",
		&q2,&xx,&ffff,&fmax,&fmin) != EOF) {
    ic++;
    if (q2!=q2old) {
      q2old = q2;
      ix=-1;
      iq2++;
    }
    if (xx!=xold) {
      xold = xx;
      ix++;
    }
    cout << setw(6) << ic
	 << setw(4) << iq2
	 << setw(4) << ix
	 << ": " 
	 << q2 << "   "
      	 << xx << "   "
	 << ffff << "   "
	 << fmax << "   "
	 << fmin << "   "
	 <<endl;


    if ( (ix  >= 0 && ix  < n_bin_x) 
	 &&
	 (iq2 >= 0 && iq2  < n_bin_q2)) {
      
      xl_bin[0][ix] = TMath::Log(xx);
      xl_bin[1][ix] = TMath::Log(xx);
      xl_bin[2][ix] = TMath::Log(xx);
      
      q2l_bin[0][iq2] = TMath::Log(q2);
      q2l_bin[1][iq2] = TMath::Log(q2);
      q2l_bin[2][iq2] = TMath::Log(q2);
      
      f[0][ix][iq2] = ffff;
      f[1][ix][iq2] = fmax;
      f[2][ix][iq2] = fmin;
      cout << setw(6) << ic
	   << setw(4) << iq2
	   << setw(4) << ix
	   << ": " 
	   << TMath::Exp(q2l_bin[0][iq2]) << "   "
	   << TMath::Exp(xl_bin[0][ix])   << "   "
	   << f[0][ix][iq2] << "   "
	   << f[1][ix][iq2] << "   "
	   << f[2][ix][iq2] << "   "
	   <<endl;


    } else {
      cout << "BCR_PDF::ReadZeusKN index out of bound:  "
	   << ix << "  " << iq2 << endl;
    }    
  }
  fclose(fp);
  fp = NULL;
  
}

void BCR_PDF::ReadHector(TString chPDF,TString chPath)
{
  TString chFile = chPath + "/" +chPDF + ".dat";

  Float_t xold=-1., q2old=-1.;
  Int_t ic=-1,ix=-1,iq2=-1;
  Float_t q2=0.,xx=0.,fun=0.,xsec=0.;
  FILE * fp = fopen(chFile,"r");
  while (fscanf(fp,"%f %f %f %f",
		&q2,&xx,&xsec,&fun) != EOF) {
    ic++;
    if (q2 != q2old) {
      q2old = q2;
      ix=-1;
      iq2++;
    }
    if (xx-xold) {
      xold = xx;
      ix++;
    }


    if ( (ix  >= 0 && ix  < n_bin_x) 
	 &&
	 (iq2 >= 0 && iq2  < n_bin_q2)) {
      
      xl_bin[0][ix] = TMath::Log(xx);
      
      q2l_bin[0][iq2] = TMath::Log(q2);
      
      f[0][ix][iq2] = fun;



    } else {
      cout << "BCR_PDF::ReadHector index out of bound:  "
	   << ix << "  " << iq2 << "  "
	   << q2 << "  " << q2old << "   " 
	   << xx << "  " << xold << "   " 
	   << endl;
    }    
  }
  fclose(fp);
  fp = NULL;
  
}



Double_t BCR_PDF::Eval(const Double_t xin, const Double_t q2in, const Int_t ig) const
{
  Double_t x = TMath::Log(xin);
  Double_t y = TMath::Log(q2in);
  
  if (ig<0||ig>n_grid) return 0.;

  Int_t i = n_bin_x-1;
  for (int lx=1; lx <n_bin_x-1; lx++) {
    if ( x < xl_bin[ig][lx] ) {i=lx;break;}
    if (xl_bin[ig][lx] >= 0.) {i=lx;break;}
  }
  i = i-1; //  FORTRAN to C++ conversion

  Int_t j = n_bin_q2-1;
  for (int lq2=1; lq2 <n_bin_q2-1; lq2++) {
    if ( y < q2l_bin[ig][lq2] ) {j=lq2;break;}
  }
  j = j-1; //  FORTRAN to C++ conversion




  Double_t dx = xl_bin[ig][i+1]-xl_bin[ig][i];
  Double_t xd = (x-xl_bin[ig][i])/dx;
  Double_t dy = q2l_bin[ig][j+1]-q2l_bin[ig][j];
  Double_t yd = (y-q2l_bin[ig][j])/dy;

  Double_t rc = (f[ig][i][j]
		 + xd*(f[ig][i+1][j]-f[ig][i][j])
		 + yd*(f[ig][i][j+1]-f[ig][i][j])
		 + xd*yd*(f[ig][i+1][j+1]+f[ig][i][j]-
			  f[ig][i+1][j]-f[ig][i][j+1])
		 );

  if (debug)
    cout << xin << "  " 
	 << q2in << "  "
	 << i << "  " 
	 << j << "  " 
	 << rc << "  " 
	 << endl;
  return rc;
}

Double_t BCR_PDF::Value(const Double_t xin, const Double_t q2in) const
{
  return Eval(xin,q2in,0);
}
Double_t BCR_PDF::ErrExp(const Double_t xin, const Double_t q2in) const 
{
  return Eval(xin,q2in,1);
}

Double_t BCR_PDF::ErrModP(const Double_t xin, const Double_t q2in) const
{
  Double_t  err2sum=0;
  for (int ll=2;ll<10;ll++) {
    Double_t val  = Value(xin,q2in);
    Double_t err1 = TMath::Max(Eval(xin,q2in,ll*2)-val,0.);
    Double_t err2 = TMath::Max(Eval(xin,q2in,ll*2+1)-val,0.);
    Double_t err  = TMath::Max(err1,err2);
    err2sum = err2sum+err*err;
  }
  return TMath::Sqrt(err2sum);
}
Double_t BCR_PDF::ErrModM(const Double_t xin, const Double_t q2in) const 
{
  Double_t  err2sum=0;
  for (int ll=2;ll<10;ll++) {
    Double_t val  = Value(xin,q2in);
    Double_t err1 = TMath::Max(val-Eval(xin,q2in,ll*2),0.);
    Double_t err2 = TMath::Max(val-Eval(xin,q2in,ll*2+1),0.);
    Double_t err  = TMath::Max(err1,err2);
    err2sum = err2sum+err*err;
  }
  return TMath::Sqrt(err2sum);
}

Double_t BCR_PDF::ErrTotP(const Double_t xin, const Double_t q2in) const
{
  Double_t exp = ErrExp(xin,q2in);
  Double_t mod = ErrModP(xin,q2in);
  return TMath::Sqrt(exp*exp+mod*mod);
}
Double_t BCR_PDF::ErrTotM(const Double_t xin, const Double_t q2in) const
{
  Double_t exp = ErrExp(xin,q2in);
  Double_t mod = ErrModM(xin,q2in);
  return TMath::Sqrt(exp*exp+mod*mod);
}


TString BCR_PDF::GetCycleString(UInt_t iCycle)
{
  iCycle=iCycle%1000;

  TString chCycle;  
  if      (iCycle<10) 
    chCycle = ("_00"+chNum[iCycle]);
  else if (iCycle<100)
    chCycle = ("_0"+chNum[iCycle/10]+chNum[iCycle%10]);
  else if (iCycle<1000)
    chCycle = ("_"+chNum[iCycle/100]+chNum[(iCycle/10)%10]+chNum[iCycle%10]);

  if (debug) cout << setw(5) << iCycle << ":   " << chCycle << endl;
  return chCycle;
}

void BCR_PDF::Print(Double_t x,Double_t q2) const
{
  cout << "BCR_PDF::Print  " << GetName() << " " << GetTitle() << endl;
  cout << Value(x,q2) 
       << " +- " << ErrExp(x,q2) << " (stat) "
       << " + "  << ErrModP(x,q2) 
       << " - "  << ErrModM(x,q2) 
       << endl;
  cout << Value(x,q2) 
       << " + "  << ErrTotP(x,q2) 
       << " - "  << ErrTotM(x,q2) 
       << endl;
}
void BCR_PDF::SelfTest() 
{
  cout << "Executing SelfTest on BCR_PDF: "  << endl;
  cout << "Name:  " << GetName() << endl;
  cout << "Title: " << GetTitle() << endl;;
  for (int lx=0; lx < n_bin_x; lx++) {
    cout << lx << ":  " << xl_bin[0][lx] << endl;
  }

}

TGraphAsymmErrors * BCR_PDF::ErrBand(TString chMode,
				     TString asFunOf,
				     Double_t rMin,
				     Double_t rMax,
				     TString  forFixed,
				     Double_t rFix,
				     Double_t s)
{
  Int_t nP=100;
  Double_t vxax[nP],vval[nP],vexl[nP],vexh[nP],veyl[nP],veyh[nP];
  Double_t vx[nP], vq2[nP], vy[nP];

  for (int lp=0;lp<nP;lp++) {
    if        (asFunOf=="x" && forFixed=="Q2") {      
      vq2[lp] = rFix;      
      Double_t rlmin = TMath::Log10(rMin);
      Double_t rlmax = TMath::Log10(rMax);
      vx[lp] = TMath::Power(10.,rlmin+lp*(rlmax-rlmin)/nP);
      vy[lp] = vq2[lp]/s/vx[lp]; 
      vxax[lp] = vx[lp];
      
    } else if (asFunOf=="Q2" && forFixed=="x") {      
      vx[lp] = rFix;      
      Double_t rlmin = TMath::Log10(rMin);
      Double_t rlmax = TMath::Log10(rMax);
      vq2[lp] = TMath::Power(10.,rlmin+lp/nP*(rlmax-rlmin));
      vy[lp] = vq2[lp]/s/vx[lp]; 
      vxax[lp] = vq2[lp];
    }
  
    vval[lp] = Value(vx[lp],vq2[lp]);
    vexl[lp]=0.0;
    vexh[lp]=0.0;
    if        (chMode=="Exp") {
      veyh[lp] = ErrExp(vx[lp],vq2[lp]);
      veyl[lp] = veyh[lp];
    } else if (chMode=="Mod") {
      veyh[lp] = ErrModP(vx[lp],vq2[lp]);
      veyl[lp] = ErrModM(vx[lp],vq2[lp]);
    } else if (chMode=="Tot") {
      veyh[lp] = ErrTotP(vx[lp],vq2[lp]);
      veyl[lp] = ErrTotM(vx[lp],vq2[lp]);
    } 
    
    if (debug) cout << lp << ": "  
		    << " x = " << vx[lp]
		    << " Q2 = " << vq2[lp]
		    << "  : "
		    << "vxax " << vxax[lp] << " " 
		    << "vval " << vval[lp] << " " 
		    << "vexl " << vexl[lp] << " " 
		    << "vexh " << vexh[lp] << " " 
		    << "veyl " << veyl[lp] << " " 
		    << "veyh " << veyh[lp] << " " 
		    << endl;

  } // end of loop over points in band

  TGraphAsymmErrors * theBand =
    new TGraphAsymmErrors(nP,vxax,vval,vexl,vexh,veyl,veyh);

  return theBand;
}

void BCR_PDF::FullDump(){
  for (int lq2=0;  lq2<n_bin_q2; lq2++) {
    for (int lx=0; lx<n_bin_x;   lx++) {
      cout << setw(4) << lq2
	   << setw(4) << lx << "  " 
	   << TMath::Exp(q2l_bin[0][lq2]) << "  " 
	   << TMath::Exp(xl_bin[0][lx]) << "  " 
	   << f[0][lx][lq2] << "  " 
	   << f[1][lx][lq2] << "  " 
	   << f[2][lx][lq2]
	   << endl;
    }
  }

}