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

#include "tpc/spacecharge/CubSpline.h"
#include "tpc/spacecharge/BSpline.h"
#include <iostream>
#include <vector>
#include "TF1.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "TRandom.h"

// helper class ------------------------------------------------
class Function {
public:
  
  void setCSP(CubSpline* csp) {
    _csp = csp;
    _knots = _csp->getKnots();
    _length = _knots->size();
    _coeffs = new std::vector<double>(_length,0);
    std::cout<<"ncoeffs="<<_coeffs->size()<<std::endl;
  }
   
  double eval(double *x, double *p) {
    //user is responsible for proper length of parameter array!
    for(int i=0; i<_length; i++) 
      (*_coeffs)[i]=p[i];  //array to vector conversion

    _csp->setCoeffs(_coeffs);
    double result = _csp->eval(x[0]);
    return result;
  }

private:
  CubSpline* _csp;
  std::vector<double>* _knots;
  std::vector<double>* _coeffs;
  int _length;

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



void testCubSpline() {

int Nknots = 15;
int Nlamda = Nknots+8;
std::vector<double>* knots = new std::vector<double>;
std::vector<double>* coeff = new std::vector<double>;
double lL = 0;
double uL = 15;
double step = (uL - lL)/Nknots;

double k;

for(int i=0; i<Nlamda; ++i){
  k = lL + (i-4)*step;
  knots->push_back(k); 
  coeff->push_back(0);  //alternating 1's and 0's
}

std::cout<<"\n\nFirst 10 knots:"<<std::endl;
for(int k=0; k<10; k++)
  std::cout<<knots->at(k)<<"  ";
std::cout<<std::endl;

// B-Spline testing --------------------------------------------

BSpline* Spline = new BSpline(knots, 10);
Spline->print();
for(int j=0; j<5; j++) {
  std::cout<<"\nValue at x = "<<(knots->at(10)) - j*step<<": "
	   <<Spline->eval(knots->at(10) - j*step);
}

std::cout<<"\nInitialized BSpline\n"<<std::endl;

// -------------------------------------------------------------


TCanvas* c1 = new TCanvas();



// Cub-Spline testing ------------------------------------------
CubSpline* csp = new CubSpline(knots, coeff);

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

//TF1* sin = new TF1("blub", "sin(x)",0,10);
TH1D* hist = new TH1D("blub","blub",30 ,0,15);
for ( Int_t i = 0; i<1000; i++ ) hist->Fill(gRandom->Gaus(6,1));
for ( Int_t i = 0; i<1000; i++ ) hist->Fill(gRandom->Gaus(9,1));
TF1* f1 =new TF1("f1",func,&Function::eval,3,15,23,"Function", "eval");
Double_t param[23] = {0,0,0,0,0,0,0,10,10,10,10,10,10,0,0,0,0,0,0,0,0,0,0};

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

// f1->Draw();


}