/*#if !defined (__CINT__) || defined (__MAKECINT__)
#include "Rtypes.h"
#endif*/


int NDATA=1000;
int NMEAS=1000;
double SIGMAVEL=1.;
double WPEN=0.2;
double Y[1000];
double T;
double UBIAS;
double TH;
double YMEAS[1000],VELMEAS[1000];

double EFIT[1000];
double EFITINTERP[1000];
double VFIT[1000];
double vscale=0.;


#include <Math/Interpolator.h>
//#include <TMath.h>
#include <TMinuit.h>
#include <vector>

/*extern "C" {
	Bool_t interp1spline(UInt_t,UInt_t, Double_t*, Double_t*, Double_t*, Double_t*);
	Double_t Mobility(Int_t,Double_t,Double_t);
	void fcn(int& , double*, double& , double* ,int);
	void calculateEfield(UInt_t, Double_t*, Double_t*, UInt_t, Double_t, Double_t,    
	Double_t*, Double_t*, Double_t*, Double_t, Double_t, Double_t, UInt_t,
	UInt_t, Double_t*,Double_t , Double_t,
	Double_t ,Double_t);	
}*/	
	

Bool_t interp1spline(UInt_t N1,UInt_t N2,Double_t *x1,Double_t *y1,Double_t *x2,Double_t *y2i){
	/* -Introduction:
	 * - Materials:
	 * 	- N1,N2: arrays length.
	 * 	- x1,y1: x,y of the curve to be interpolated on the points x2. x1,x2 ODERED!!!
	 * 	- x2: point where to interpolate y1.
	 * 	- y2i: output! output values of y1 on x2.
	*/
	std::vector<double> X1;
	std::vector<double> Y1;
	for (UInt_t i=0;i<N1;i++){
		X1.push_back(x1[i]);
		Y1.push_back(y1[i]);
	}
	ROOT::Math::Interpolator inter(X1,Y1);
	// Serarch min and max to check extremes
	Double_t x1min=x1[0];
	Double_t x1max=x1[N1-1];
	for (UInt_t i=0;i<N1;i++){
		if (x1[i]<x1min)
			x1min=x1[i];
		if (x1[i]>x1max)
			x1max=x1[i];
	}
	for (UInt_t i=0;i<N2;i++){
		if (x2[i]<x1min){
	//		printf("Warning: interp1spline: Extrapolation lower extreme.\n");
			y2i[i]=inter.Eval(x1min)+inter.Deriv(x1min)*(x2[i]-x1min)+inter.Deriv2(x1min)*pow(x2[i]-x1min,2);
			continue;
		}
		if (x2[i]>x1max){
	//		printf("Warning: interp1spline: Extrapolation higher extreme.\n");
			y2i[i]=inter.Eval(x1max)+inter.Deriv(x1max)*(x2[i]-x1min)+inter.Deriv2(x1max)*pow(x2[i]-x1min,2);
			continue;
		}
		y2i[i]=inter.Eval(x2[i]);

	}

		return true;
}



// Mobility
Double_t Mobility(Int_t Charge,Double_t T,Double_t E){
  Double_t lfm=0,hfm=0;
  Double_t vsatn,vsatp,vsat;
  Double_t betap,betan;
  Double_t alpha;
	      Double_t bb,cc,E0;
  E=E*1e-2;
    if(Charge>0)
    {
     E0=2970*pow(T/300,5.63);
     bb=9.57e-8*pow(T/300,-0.155);
     cc=-3.24e-13;
     lfm=457*pow(T/300,-2.80);
     if(E>E0) hfm=1./(1/lfm+bb*(E-E0)+cc*pow(E-E0,2)); else hfm=lfm;
    }
  else
    {
      E0=2970*pow(T/300,5.63);
      lfm=1430*pow(T/300,-1.99);
      vsatn=1.05e7*pow(T/300,-3.02);
      if(E>E0) hfm=1./(1/lfm+1/vsatn*(E-E0)); else hfm=lfm;
    }
return hfm/1e4;
}

////////// Minuit
void fcn(int &npar, double *gin, double &f, double *par,int iflag)
{
	int k,i;
	double uk=0.;
	double chi2 = 0.0;
	double mue,muh;
	double Ei[NMEAS];
	double dy=0.;
	for(k=0;k<NMEAS; k++)
	{
		interp1spline(NDATA,NMEAS,Y,par,YMEAS,Ei);
	}
	vscale=0.;
	for(k=0;k<NMEAS-1; k++)
	{
		muh=Mobility(+1.,T,Ei[k]);
		mue=Mobility(-1.,T,Ei[k]);
		dy=YMEAS[k+1]-YMEAS[k];
		vscale+=1./UBIAS*VELMEAS[k]/(mue+muh)*dy;
	}
		muh=Mobility(+1.,T,Ei[NMEAS-1]);
		mue=Mobility(-1.,T,Ei[NMEAS-1]);
	vscale+=1./UBIAS*VELMEAS[NMEAS-1]/(mue+muh)*dy; // Last step

	for(k=0;k<NMEAS; k++)
	{
		muh=Mobility(+1.,T,Ei[k]);
		mue=Mobility(-1.,T,Ei[k]);
		uk=(muh+mue)*Ei[k];
		chi2+=1./pow(SIGMAVEL,2)*pow(1.-VELMEAS[k]/(vscale*uk),2);
	//	printf("velk = %e \tvscal*uk=%e\t uk=%e \tvscale=%e \n",velmeas[k],vscale*uk,uk,vscale);

	}
	for(i=1;i<NDATA-1; i++)
	{
		uk=(muh+mue)*par[i];
		chi2+=WPEN*pow((0.5*(par[i-1]+par[i+1])-par[i])/(par[i]),2);
	}
	//printf("chi2= %e\n",chi2);
	for(i=0;i<NDATA; i++)
			EFIT[i]=par[i];
	f = chi2;
}

void calculateEfield(UInt_t nmeas, Double_t *ymeas,Double_t *velmeas, UInt_t ndata, Double_t Offsety, Double_t Dy,    
	Double_t *Efit, Double_t *Efitinterp,Double_t *Vfit,Double_t U,Double_t th=300e-6,Double_t T0=273.15,UInt_t Nfilt=1,
	UInt_t NITERFIT=100000,Double_t *Efit0=NULL,Double_t Efitmin=0., Double_t Efitmax=10e7,
	Double_t sigmavel=1.,Double_t wpen=0.2){




	// Set global values
	NDATA=ndata;
	NMEAS=nmeas;
	T=T0;
	UBIAS=U;
	TH=th;
	WPEN=wpen;
	SIGMAVEL=sigmavel;
	for (UInt_t i=0;i<nmeas;i++){
		YMEAS[i]=ymeas[i];
		VELMEAS[i]=velmeas[i];
	}
	for (UInt_t i=0;i<ndata;i++){
		Efit[i]=EFIT[i];
	}
	for (UInt_t i=0;i<ndata;i++){
		Y[i]=Offsety+i*Dy;
	}

	TMinuit minuit(NDATA);
	minuit.SetFCN(fcn);
	Double_t arglist[10];
	Int_t ierflg = 0;
	Double_t mue,muh;
	   Double_t amin,edm,errdef;
	   Int_t nvpar,nparx,icstat;
	double step[NDATA];
	double start[NDATA];

	for (UInt_t i=0;i<NMEAS-Nfilt;i++){
		Double_t temp=0.;
		for (UInt_t j=0;j<Nfilt;j++)
			temp+=velmeas[i+j]/Nfilt;
		velmeas[i]=temp;
	}

	// Minuit
	if (Efit0==NULL){
		for (UInt_t i=0;i<NDATA;i++){
			start[i]=(2*U/th)-i/NDATA*(2*U/th);
		}
	}
	else{
		for (UInt_t i=0;i<NDATA;i++){
			start[i]=Efit0[i];
		}
	}
	for (UInt_t i=0;i<NDATA;i++){
		step[i]=U/th/10.;
	}

	arglist[0] = 1;
	//TMinuit *gMinuit = new TMinuit(3);
	printf("\n Teste dich! \n");
	gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
	
	for (UInt_t i=0;i<NDATA;i++){
		Char_t str[10];
		sprintf(str,"Ek%d",i);
		minuit.mnparm(i,str, start[i], step[i], Efitmin,Efitmax, ierflg);
	}
		printf("\n Teste dich! \n");
		gMinuit->mnexcm("CALL FCN", arglist ,1,ierflg);
		printf("\n Teste dich! \n");
	// Now ready for minimization step
	   arglist[0] = NITERFIT;
	   arglist[1] = 1.;
	   gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
	 printf("\n Teste dich! \n");
	// Print results
	// Test interp
	
	interp1spline(NDATA,NMEAS,Y,EFIT,YMEAS,EFITINTERP);
	// Velocity
	for (UInt_t i=0;i<NMEAS;i++){
		muh=Mobility(+1.,T,EFITINTERP[i]);
		mue=Mobility(-1.,T,EFITINTERP[i]);
		VFIT[i]=(muh+mue)*EFITINTERP[i];
	//	printf("mue=%e,  muh=%e,   E=%e,   v=%e \n",mue,muh,Efitinterp[i],vfit[i]);
	}
	for (UInt_t i=0;i<NMEAS;i++){
		VFIT[i]=VFIT[i]*vscale;//vfitmax;
	}


	// Return values of interest
	for (UInt_t i=0;i<nmeas;i++){
		Efitinterp[i]=EFITINTERP[i];
		Vfit[i]=VFIT[i];
	}
	for (UInt_t i=0;i<ndata;i++){
		Efit[i]=EFIT[i];
	}

}