//testing script for the spline stuff - needs to be compiled
//Felix

#include "tpc/spacecharge/BiCubSpline.h"
#include "tpc/spacecharge/BSpline.h"
#include <iostream>
#include <vector>
#include "TF2.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TRandom.h"

// helper class ------------------------------------------------
class Function {
public:
  
  void setCSP(BiCubSpline* csp) {
    _csp = csp;
    _kx = _csp->getKx();
    _ky = _csp->getKy();
    _xl = _kx->size();
    _yl = _ky->size();
    _coeffs = new std::vector<std::vector<double>*>;
    std::cout<<"ncoeffs="<<_coeffs->size()<<std::endl;
  }
   
  double eval(double *x, double *p) {
    //user is responsible for proper length of parameter array!
    _csp->setCoeffsByArray(p);
    
    double result = _csp->eval(x[0], x[1]);
    return result;
  }

private:
  BiCubSpline* _csp;
  std::vector<double>* _kx;
  std::vector<double>* _ky;
  std::vector<std::vector<double>*>* _coeffs;
  int _xl, _yl;

};
// -------------------------------------------------------------



void testBiCubSpline() {

int NknotsX = 4;
int NknotsY = 4; 
int NlamdaX = NknotsX+8;
int NlamdaY = NknotsY+8;
std::vector<double>* _kx = new std::vector<double>;
std::vector<double>* _ky = new std::vector<double>;
std::vector<std::vector<double>*>* coeff=new std::vector<std::vector<double>*>;
double lLx = 0;
double uLx = 15;
double lLy = 0;
double uLy = 15;
double stepx = (uLx - lLx)/NknotsX;
double stepy = (uLy - lLy)/NknotsY;

double kx,ky;

for(int i=0; i<NlamdaX; ++i){
  kx = lLx + (i-4)*stepx;
  _kx->push_back(kx); 
}
for(int j=0; j<NlamdaY; ++j){
  ky = lLy + (j-4)*stepy;
  _ky->push_back(ky); 
}

for(int i=0; i<NlamdaY; i++) {
  coeff->push_back(new std::vector<double>(NlamdaX, 0));
}

TCanvas* c1 = new TCanvas();



// Cub-Spline testing ------------------------------------------
BiCubSpline* csp = new BiCubSpline(_kx, _ky, coeff);

Function* func = new Function();   
func->setCSP(csp);


//TF1* sin = new TF1("blub", "sin(x)",0,10);
TH2D* hist = new TH2D("blub","blub",40 ,0,15, 30,0,15);
for ( Int_t i = 0; i<1000; i++ ) hist->Fill(gRandom->Gaus(5,1), 
					    gRandom->Gaus(3,1));
//for ( Int_t i = 0; i<1000; i++ ) hist->Fill(gRandom->Gaus(9,1),
//					    gRandom->Gaus(10,1));
TF2* f2 =new TF2("f2",func,&Function::eval,0,15,0,15,NlamdaX*NlamdaY,
		 "Function", "eval");
Double_t param[14*14] = {0};
for(int i=5; i<7;i++)
   param[i]=10;
for(int i=17; i<19;i++)
   param[i]=10;


f1->SetParameters(param);
f1->DrawClone();
hist->Fit("f1","R");
hist->Draw("SURF");

// f1->Draw();


}