#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>
#include <iomanip>

using namespace std;

int det_entries(ostringstream& oss){
	
	char string[256];
	char GetLine[1000];

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

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

	int nCounter = -1;
	int nNoLines;

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

	while(file.good() && !file.eof()){
		for(int i = 0; i < 9; i++){
		file >> string;
		cout << string << " ";
		}
		cout << endl; 
		nCounter++;
	}

	nNoLines = nCounter;
	
	cout << "file has " << nNoLines << " lines"  << endl;

	file.close();

	return nNoLines;
}

void read_fit_file(ostringstream& ossFit, int nEntriesFit, vector<double>& dMassFit,vector<double>& dVolumeFit,vector<double>& dChargeFit,vector<double>& dErrorFit){

	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(vector<double>& dMassFit,vector<double>& dVolumeFit,vector<double>& dChargeFit,vector<double>& dErrorFit, int nEntries){

	for(int i = 0; i < nEntries; 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[0]*x[0])*x[1]*x[1]);		
		}

		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 3;
		}

};


// definition of shared parameter
// background function 
int iparObsA[3] = {0,1,2};

// signal + background function 
int iparObsB[3] = {3,4,2};

class GlobalChi2 : public ROOT::Math::FitMethodFunction { 

public:
   GlobalChi2( int dim, int npoints, ROOT::Math::FitMethodFunction & f1,  ROOT::Math::FitMethodFunction & f2) :
      ROOT::Math::FitMethodFunction(dim,npoints),
      fChi2_1(&f1), fChi2_2(&f2) {}


   ROOT::Math::IMultiGenFunction * Clone() const { 
      // copy using default copy-ctor
      // i.e. function pointer will be copied (and not the functions)
      return new GlobalChi2(*this);   
   }

   double  DataElement(const double *par, unsigned int ipoint, double *g = 0) const { 
      // implement evaluation of single chi2 element
      double p1[3];
      for (int i = 0; i < 3; ++i) p1[i] = par[iparObsA[i] ];

      double p2[3]; 
      for (int i = 0; i < 3; ++i) p2[i] = par[iparObsB[i] ];

      double g1[3]; 
      double g2[3];
      double value = 0; 
      if (g != 0) { 
         for (int i = 0; i < 6; ++i) g[i] = 0; 
      }

      if (ipoint < fChi2_1->NPoints() ) {
         value = fChi2_1->DataElement(p1, ipoint, g1);
         // update gradient values
         if (g != 0) { 
            for (int i= 0; i < 3; ++i) g[iparObsA[i]] = g1[i];         
         }
      }
      else { 
         unsigned int jpoint = ipoint- fChi2_1->NPoints();
         assert (jpoint < fChi2_2->NPoints() );
         value =  fChi2_2->DataElement(p2, jpoint, g2);
         // update gradient values
         if ( g != 0) { 
            for (int i= 0; i < 3; ++i) g[iparObsB[i]] = g2[i];
         }
      }
        
      return value; 
      
   }


   // needed if want to use Fumili or Fumili2
   virtual Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; }

private:
   // parameter vector is first background (in common 1 and 2) and then is signal (only in 2)
   virtual double DoEval (const double *par) const {
      double p1[3];
      for (int i = 0; i < 3; ++i) p1[i] = par[iparObsA[i] ];

      double p2[3]; 
      for (int i = 0; i < 3; ++i) p2[i] = par[iparObsB[i] ];

      return (*fChi2_1)(p1) + (*fChi2_2)(p2);
   } 

   const  ROOT::Math::FitMethodFunction * fChi2_1;
   const  ROOT::Math::FitMethodFunction * fChi2_2;
};


void NonLinLSFit() { 

	vector<double> dMassFitA;
	vector<double> dVolumeFitA;
	vector<double> dChargeFitA;
	vector<double> dErrorFitA;

	vector<double> dMassFitB;
	vector<double> dVolumeFitB;
	vector<double> dChargeFitB;
	vector<double> dErrorFitB;

	int beta = 5;
	double qmass = 0.03;
	
	string ObsA = "Pion";
	string ObsB = "WL_exp_1";

	ostringstream ossFitA;
	ossFitA << "fit_mass_table_beta_" << beta << "_qmass_" << fixed << setprecision(2) << qmass << "_tc_" << ObsA.c_str() << ".dat";

	ostringstream ossFitB;
	ossFitB << "fit_mass_table_beta_" << beta << "_qmass_" << fixed << setprecision(2) << qmass << "_tc_" << ObsB.c_str() << ".dat";

	int nEntriesFitA = det_entries(ossFitA);
	int nEntriesFitB = det_entries(ossFitB);

	read_fit_file(ossFitA, nEntriesFitA, dMassFitA, dVolumeFitA, dChargeFitA, dErrorFitA);
	read_fit_file(ossFitB, nEntriesFitB, dMassFitB, dVolumeFitB, dChargeFitB, dErrorFitB);
  	
		cout << "Generating Graphs" << endl;
  	//TGraph2DErrors *grA = new TGraph2DErrors(nEntriesFitA);
  	//TGraph2DErrors *grB = new TGraph2DErrors(nEntriesFitB);
	TGraph2DErrors** gr = new TGraph2DErrors*[2]();

	gr[0] = new TGraph2DErrors();
	gr[1] = new TGraph2DErrors();
		cout << "success..." << endl;

	cout << "First data set:" << endl;
	DisplayFitData(dMassFitA, dVolumeFitA, dChargeFitA, dErrorFitA, nEntriesFitA);

	cout << "Second data set:" << endl;
	DisplayFitData(dMassFitB, dVolumeFitB, dChargeFitB, dErrorFitB, nEntriesFitB);

	//Setting the points for the graph grA
	cout << "Setting data points for grA:" << endl;
	for(int i = 0; i < nEntriesFitA; i++){

		gr[0]->SetPoint(i,dVolumeFitA[i],dChargeFitA[i],dMassFitA[i]);
		gr[0]->SetPointError(i,0,0,dErrorFitB[i]);

		cout << "grA->SetPoint(" << i << "," << dVolumeFitA.at(i)  << "," << dChargeFitA.at(i)<< "," << dMassFitA.at(i) << "),  grA->SetPointError(" << i << ",0,0," << dErrorFitA.at(i) << ")" << endl;
	}
	//Setting the points for the graph grB
	cout << endl << "Setting data points for grB:" << endl;
	for(int i = 0; i < nEntriesFitB; i++){

		//grB->SetPoint(i,dVolumeFitB[i],dChargeFitB[i],dMassFitB[i]);
		//grB->SetPointError(i,0,0,dErrorFitB[i]);
		gr[1]->SetPoint(i,dVolumeFitB[i],dChargeFitB[i],dMassFitB[i]);
		gr[1]->SetPointError(i,0,0,dErrorFitB[i]);

		cout << "grB->SetPoint(" << i << "," << dVolumeFitB.at(i)  << "," << dChargeFitB.at(i)<< "," << dMassFitB.at(i) << "),  grA->SetPointError(" << i << ",0,0," << dErrorFitB.at(i) << ")" << endl;
	}
	
	MyParametricFunction* fobjA = new MyParametricFunction();
	MyParametricFunction* fobjB = new MyParametricFunction();

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

	ROOT::Math::WrappedMultiTF1 EPA(*MyFuncA,2);	
	ROOT::Math::WrappedMultiTF1 EPB(*MyFuncB,2);	

	ROOT::Fit::DataOptions optA;
	ROOT::Fit::DataOptions optB;

	ROOT::Fit::DataRange rangeA;
	ROOT::Fit::DataRange rangeB;

	ROOT::Fit::BinData dataEPA(optA,rangeA);
	ROOT::Fit::BinData dataEPB(optB,rangeB);
	ROOT::Fit::FillData(dataEPA, gr[0]);
	ROOT::Fit::FillData(dataEPB, gr[1]);
	
	ROOT::Fit::Chi2Function EPAChi2(dataEPA, EPA);
	ROOT::Fit::Chi2Function EPBChi2(dataEPB, EPB);

	ROOT::Fit::Fitter EPfitter;

	const int EPparams = 5;
	double FitterParams[EPparams] = {0.35,-0.03,0.03,0.1,-0.009}; 

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

	EPfitter.Config().ParSettings(0).SetLimits(0.,2.0);
	EPfitter.Config().ParSettings(1).SetLimits(-0.5,0.5);
	EPfitter.Config().ParSettings(2).SetLimits(0.,0.5);
	EPfitter.Config().ParSettings(3).SetLimits(0,2.0);
	EPfitter.Config().ParSettings(4).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().ParSettings(3).SetStepSize(0.0001);
	EPfitter.Config().ParSettings(4).SetStepSize(0.0001);

	GlobalChi2 globalChi2(EPparams,dataEPA.Size()+dataEPB.Size(),EPAChi2, EPBChi2);

	
	EPfitter.Config().MinimizerOptions().SetPrintLevel(1);
	EPfitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
	//EPfitter.Config().SetMinimizer("Minuit2","Migrad");

	EPfitter.FitFCN(globalChi2,0,dataEPA.Size()+dataEPB.Size(),true);

	ROOT::Fit::FitResult result = EPfitter.Result();

	cout << "Size of fitted data: " << dataEPA.Size()+dataEPB.Size() << endl;

	result.Print(std::cout);

}