#include "Math/IFunction.h"
#include "Math/IParamFunction.h"
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TH1.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "Math/FitMethodFunction.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGraph2DErrors.h"
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <vector>

using namespace std;

vector<double> dMassFit;
vector<double> dVolumeFit;
vector<double> dChargeFit;
vector<double> dErrorFit;

void read_fit_file(ostringstream& ossFit, int nEntriesFit){

	string Value; 
	char GetLine[1000];

	ifstream file;
	file.open(ossFit.str().c_str());

	if(file.good() == 0){
		cout << "File does not exist, terminating ..." << endl;
		exit(0);
	}

	int nCounter = 0;

	for(int i = 0; i < 3; i++){
		file.getline(GetLine,1000);
	}

	while((nCounter < nEntriesFit)){
		file >> Value;
		dMassFit.push_back(atof(Value.c_str()));	
		file >> Value;
		dErrorFit.push_back(atof(Value.c_str()));	
		file >> Value;
		dVolumeFit.push_back(atof(Value.c_str()));	
		file >> Value;
		dChargeFit.push_back(atoi(Value.c_str()));					

		for(int i = 0; i < 5; i++){
		file >> Value;
		}


		nCounter++;
	}

	file.close();

}

void DisplayFitData(){

	for(int i = 0; i < 8; i++){

		cout << dMassFit.at(i) << " " << dVolumeFit.at(i) << " " << dChargeFit.at(i) << " " << dErrorFit.at(i) << endl;
	}
}


class MyParametricFunction: public ROOT::Math::IParametricFunctionMultiDim
{
	private:
		const double* pars;

	public:
		double DoEvalPar(const double* x, const double* p) const
		{
			return (p[0] + p[1]*p[2]/x[0] + p[1]*p[2]*p[2]*x[1]*x[1]/(x[0]*x[0]));		
		}

		unsigned int NDim() const
		{
			return 2;
		}

		ROOT::Math::IParametricFunctionMultiDim* Clone() const
		{
			return new MyParametricFunction();
		}

		const double* Parameters() const 
		{
			return pars;
		}

		void SetParameters(const double* p)
		{
			pars = p;
		}

		unsigned int NPar() const
		{
			return 2;
		}

};


void NonLinLSFit() { 

  
  	int nEntriesFit = 8;

	ostringstream ossFit;
	ossFit << "fit_mass_table_beta_1_qmass_0.01_tc.dat";

	read_fit_file(ossFit, nEntriesFit);
  	
  	TGraph2DErrors *gr = new TGraph2DErrors(nEntriesFit);

	DisplayFitData();

	//Setting the points for the graph
	for(int i = 0; i < nEntriesFit; i++){

		gr->SetPoint(i,dVolumeFit[i],dChargeFit[i],dMassFit[i]);
		gr->SetPointError(i,0,0,dErrorFit[i]);
	}
	
	MyParametricFunction* fobj = new MyParametricFunction();

	int npar = 3;
	//TF1 * MyFunc = new TF1("MyFunc",fobj,&MyParamFunc::Extrapolation,0,1,npar,"MyParamFunction","Extrapolation");
	TF1 * MyFunc = new TF1("MyFunc",fobj,0,1,npar,"MyParametricFunction");

	ROOT::Math::WrappedMultiTF1 EP(*MyFunc,2);	

	ROOT::Fit::DataOptions opt;
	ROOT::Fit::DataRange rangeB;
	
	rangeB.SetRange(0,0.007);

	ROOT::Fit::BinData dataEP(opt,rangeB);
	ROOT::Fit::FillData(dataEP, gr);
	
	ROOT::Fit::Chi2Function EPChi2(dataEP, EP);

	ROOT::Fit::Fitter EPfitter;

	const int EPparams = 3;
	double FitterParams[EPparams] = {0.6,-0.3,0.3}; 

	EPfitter.Config().SetParamsSettings(3,FitterParams);

	EPfitter.Config().ParSettings(0).SetLimits(0.5,1.0);
	EPfitter.Config().ParSettings(1).SetLimits(0.,-0.5);
	EPfitter.Config().ParSettings(2).SetLimits(0.5,-0.5);
	
	EPfitter.Config().ParSettings(0).SetStepSize(0.0001);
	EPfitter.Config().ParSettings(1).SetStepSize(0.0001);
	EPfitter.Config().ParSettings(2).SetStepSize(0.0001);
	
	EPfitter.Config().MinimizerOptions().SetPrintLevel(1);
	EPfitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");

	EPfitter.FitFCN(EPChi2,0,dataEP.Size(),true);
	ROOT::Fit::FitResult result = EPfitter.Result();

	result.Print(std::cout);

}