#ifndef ETA_INTERPOLATION_BASE_H
#define ETA_INTERPOLATION_BASE_H

#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2D.h>
#include <TH2F.h>
#endif
#include <cmath> 
#include "slsInterpolation.h"
#include "tiffIO.h"

class etaInterpolationBase : public slsInterpolation {
  
 public:
 
 etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nsy=25, int nb=-1, int nby=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns,nsy), hhx(NULL), hhy(NULL), heta(NULL), nbetaX(nb), nbetaY(nby), etamin(emin), etamax(emax) {
    // cout << "eb " << nb << " " << emin << " " << emax << endl;  
    // cout << nb << " " << etamin << " " << etamax << endl;
    if (nbetaX<=0) {
      //cout << "aaa:" <<endl;
      nbetaX=nSubPixelsX*10;
    } 
    if (nbetaY<=0) {
      //cout << "aaa:" <<endl;
      nbetaY=nSubPixelsY*10;
    }
    if (etamin>=etamax) {
      etamin=-1;
      etamax=2;
    }
    etastepX=(etamax-etamin)/nbetaX;
    etastepY=(etamax-etamin)/nbetaY;
    heta=new int[nbetaX*nbetaY];
    hhx=new float[nbetaX*nbetaY];
    hhy=new float[nbetaX*nbetaY];
    rangeMin=etamin;
    rangeMax=etamax;
    flat= new double[nSubPixelsX*nSubPixelsY];
    hintcorr=new int [nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY];
    
  };
  
 etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
   nbetaX=orig->nbetaX;
   nbetaY=orig->nbetaY;
   etamin=orig->etamin;
   etamax=orig->etamax;
   rangeMin=orig->rangeMin;
   rangeMax=orig->rangeMax;


    etastepX=(etamax-etamin)/nbetaX;
    etastepY=(etamax-etamin)/nbetaY;
    heta=new int[nbetaX*nbetaY];
    memcpy(heta,orig->heta,nbetaX*nbetaY*sizeof(int));
    hhx=new float[nbetaX*nbetaY];
    memcpy(hhx,orig->hhx,nbetaX*nbetaY*sizeof(float));
    hhy=new float[nbetaX*nbetaY];
    memcpy(hhy,orig->hhy,nbetaX*nbetaY*sizeof(float));
    hintcorr=new int [nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY];

 };

  



  virtual void resetFlatField() {
    for (int ibx=0; ibx<nbetaX*nbetaY; ibx++) {
      heta[ibx]=0;
      hhx[ibx]=0;
      hhy[ibx]=0;
    }


  };
    
  int *setEta(int *h, int nb=-1, int nby=-1, double emin=1, double emax=0)
  {  
    if (h) {
      if (heta) delete [] heta;
      heta=h;
      nbetaX=nb;
      nbetaY=nby;
      if (nbetaX<=0) nbetaX=nSubPixelsX*10;
      if (nbetaY<=0) nbetaY=nSubPixelsY*10;
      
      etamin=emin;
      etamax=emax;
      if (etamin>=etamax) {
	etamin=-1;
	etamax=2;
      }
      rangeMin=etamin;
      rangeMax=etamax;
      etastepX=(etamax-etamin)/nbetaX;
      etastepY=(etamax-etamin)/nbetaY;
    }
    return heta;
  };
  
   int *setFlatField(int *h, int nb=-1, int nby=-1, double emin=1, double emax=0)
  {  
    return setEta(h, nb, nby, emin, emax);
  };



  int *getFlatField(){return setEta(NULL);};
  
  int *getFlatField(int &nb, int &nby, double &emin, double &emax){
    nb=nbetaX; 
    nby=nbetaY; 
    emin=etamin; 
    emax=etamax; 
    return getFlatField();
  }; 
  
  
  void *writeFlatField(const char * imgname) {
    float *gm=NULL;
    gm=new float[nbetaX*nbetaY];
    for (int ix=0; ix<nbetaX; ix++) {
      for (int iy=0; iy<nbetaY; iy++) {
	gm[iy*nbetaX+ix]=heta[iy*nbetaX+ix];
      }
    } 
    WriteToTiff(gm, imgname, nbetaX, nbetaY);   
    delete [] gm;
    return NULL; 
  };
  
  int readFlatField(const char * imgname, double emin=1, double emax=0) {
    if (emax>=1) etamax=emax;
    if (emin<=0) etamin=emin;   
   
    if (etamin>=etamax) {
      etamin=-1;
      etamax=2;
    }
    
    etastepX=(etamax-etamin)/nbetaX;
    etastepY=(etamax-etamin)/nbetaY;
    uint32 nnx;
    uint32 nny;
    float *gm=ReadFromTiff(imgname, nnx, nny);
    /* if (nnx!=nny) { */
    /*   cout << "different number of bins in x " << nnx << "  and y " << nny<< " !"<< endl; */
    /*   cout << "Aborting read"<< endl; */
    /*   return 0; */
    /* } */
    nbetaX=nnx;
    nbetaY=nny;
    if (gm) {
      if (heta) {
	delete [] heta;
	delete [] hhx;
	delete [] hhy;
      }
      
      heta=new int[nbetaX*nbetaY];
      hhx=new float[nbetaX*nbetaY];
      hhy=new float[nbetaX*nbetaY];
      
      for (int ix=0; ix<nbetaX; ix++) {
	for (int iy=0; iy<nbetaY; iy++) {
	  heta[iy*nbetaX+ix]=gm[iy*nbetaX+ix];
	}
      }
      delete [] gm;
      return 1;
    }
    return 0;
  };
  
    
  
    

float *gethhx()
  {
    // hhx->Scale((double)nSubPixels);
    return hhx;
  };
  
  float *gethhy()
  {
    // hhy->Scale((double)nSubPixels);
    return hhy;
    };
  virtual int addToFlatField(double etax, double etay){
    int ex,ey; 
    ex=(etax-etamin)/etastepX;
    ey=(etay-etamin)/etastepY;
    if (ey<nbetaY && ex<nbetaX && ex>=0 && ey>=0)
      heta[ey*nbetaX+ex]++; 
    return 0;    
  };
  
  
  // virtual void prepareInterpolation(int &ok)=0;


  void debugSaveAll(int ind=0) {
    int  ibx, iby;
   char tit[10000];
   
   float tot_eta=0;

  float *etah=new float[nbetaX*nbetaY];
  // int etabins=nbeta;
  int ibb=0;

  for (int ii=0; ii<nbetaX*nbetaY; ii++) {
   
      etah[ii]=heta[ii];
      tot_eta+=heta[ii];
  }
  sprintf(tit,"/scratch/eta_%d.tiff",ind);
  WriteToTiff(etah, tit, nbetaX, nbetaY);


  for (int ii=0; ii<nbetaX*nbetaY; ii++) {
    ibb=(hhx[ii]*nSubPixelsX);
    etah[ii]=ibb;
  }
  sprintf(tit,"/scratch/eta_hhx_%d.tiff",ind);
  WriteToTiff(etah, tit, nbetaX, nbetaY);
	  
  for (int ii=0; ii<nbetaX*nbetaY; ii++) {
    ibb=hhy[ii]*nSubPixelsY;
    etah[ii]=ibb;
  }
  sprintf(tit,"/scratch/eta_hhy_%d.tiff",ind);
  WriteToTiff(etah, tit, nbetaX, nbetaY);
	  
  
  float *ftest=new float[nSubPixelsX*nSubPixelsY];

  for (int ib=0; ib<nSubPixelsX*nSubPixelsY; ib++) ftest[ib]=0;
  

  //int ibx=0, iby=0;
  
  for (int ii=0; ii<nbetaX*nbetaY; ii++) {
    
    ibx=nSubPixelsX*hhx[ii];
    iby=nSubPixelsY*hhy[ii];
    if (ibx<0) ibx=0;
    if (iby<0) iby=0;
    if (ibx>=nSubPixelsX) ibx=nSubPixelsX-1;
    if (iby>=nSubPixelsY) iby=nSubPixelsY-1;
    

    if (ibx>=0 && ibx<nSubPixelsX && iby>=0 && iby<nSubPixelsY) {
      //
      // if (ibx>0 && iby>0) cout << ibx << " " << iby << " " << ii << endl;
      ftest[ibx+iby*nSubPixelsX]+=heta[ii];
    } else
      cout << "Bad interpolation "<< ii << " " << ibx << " " << iby<< endl; 
    
  }

  sprintf(tit,"/scratch/ftest_%d.tiff",ind);
  WriteToTiff(ftest, tit, nSubPixelsX, nSubPixelsY);

  //int ibx=0, iby=0;
  tot_eta/=nSubPixelsX*nSubPixelsY;
  int nbad=0;
  for (int ii=0; ii<nbetaX*nbetaY; ii++) {
    ibx=nSubPixelsX*hhx[ii];
    iby=nSubPixelsY*hhy[ii];
    if (ftest[ibx+iby*nSubPixelsX]<tot_eta*0.5) {
      etah[ii]=1;
      nbad++;
    } else if(ftest[ibx+iby*nSubPixelsX]>tot_eta*2.){
      etah[ii]=2;
      nbad++;
    } else
      etah[ii]=0;
  }
  sprintf(tit,"/scratch/eta_bad_%d.tiff",ind);
  WriteToTiff(etah, tit, nbetaX, nbetaY);
  // cout << "Index: " << ind << "\t Bad bins: "<< nbad << endl;
  //int ibx=0, iby=0;

  delete [] ftest;
  delete [] etah;

  }


 protected:
  

  double calcDiff(double avg, float *hx, float *hy) {
    //double p_tot=0;
    double diff=0, d;
    //double bsize=1./nSubPixels;
    int nbad=0;
    double p_tot_x[nSubPixelsX], p_tot_y[nSubPixelsY], p_tot[nSubPixelsX*nSubPixelsY];
    double maxdiff=0, mindiff=avg*nSubPixelsX*nSubPixelsY;
   
    int ipx, ipy;
    for (ipy=0; ipy<nSubPixelsY; ipy++) {
      for (ipx=0; ipx<nSubPixelsX; ipx++) {
	p_tot[ipx+ipy*nSubPixelsX]=0;
      }
      p_tot_y[ipy]=0;
      p_tot_x[ipy]=0;
    }
   
    for (int ibx=0; ibx<nbetaX; ibx++) {
      for (int iby=0; iby<nbetaY; iby++) {
	ipx=hx[ibx+iby*nbetaX]*nSubPixelsX;
	if (ipx<0) ipx=0;
	if (ipx>=nSubPixelsX) ipx=nSubPixelsX-1;
	
	ipy=hy[ibx+iby*nbetaX]*nSubPixelsY;
	if (ipy<0) ipy=0;
	if (ipy>=nSubPixelsY) ipy=nSubPixelsY-1;
	
	p_tot[ipx+ipy*nSubPixelsX]+=heta[ibx+iby*nbetaX];
	p_tot_y[ipy]+=heta[ibx+iby*nbetaX];
	p_tot_x[ipx]+=heta[ibx+iby*nbetaX];
      }
    }
    
    
    //  cout << endl << endl;
    for (ipy=0; ipy<nSubPixelsY; ipy++) { 
     cout.width(5);
     //flat_y[ipy]=p_tot_y[ipy];//avg/nSubPixels;
      for (ipx=0; ipx<nSubPixelsX; ipx++) {

	//	flat_x[ipx]=p_tot_x[ipx];///avg/nSubPixels;
	flat[ipx+nSubPixelsX*ipy]=p_tot[ipx+nSubPixelsX*ipy];///avg;
	d=p_tot[ipx+nSubPixelsX*ipy]-avg;
	if (d<0) d*=-1.;
	if (d>5*sqrt(avg) )
	  nbad++;
	diff+=d*d;
	if (d<mindiff) mindiff=d;
	if (d>maxdiff) maxdiff=d;
	//	cout << setprecision(4) << p_tot[ipx+nSubPixels*ipy] << "   ";
      }

      /* cout << "** "  << setprecision(4) <<  flat_y[ipy]; */
      //cout << "\n";
    }
    /* cout << "**" << endl; cout.width(5); */
    /* for (ipx=0; ipx<nSubPixels; ipx++) { */
    /*   cout  << setprecision(4) <<  flat_x[ipx] << " "; */
    /* } */
    //cout << "**" << endl; cout.width(5);
    //cout << "Min diff: " << mindiff/sqrt(avg) << " Max diff: " << maxdiff/sqrt(avg) << " Nbad: " << nbad << endl;
    


    //   cout << "Bad pixels: " << 100.*(float)nbad/((float)(nSubPixels*nSubPixels)) << " %" << endl;
    return sqrt(diff);
  }
  
  float *hhx;
  float *hhy;
  int *heta;
  int nbetaX, nbetaY;
  double etamin, etamax, etastepX, etastepY;
  double rangeMin, rangeMax;



  double *flat;
  int *hintcorr;
  

};

#endif