// test fit of fumili using a gauss function 
// with and without providing derivatives

#include "Math/IParamFunction.h"
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "TF1.h"
#include "TH1.h"
#include "TVirtualFitter.h"
#include "HFitInterface.h"
#include "Math/MinimizerOptions.h"
#include "TMinuit.h"
	
class GaussModelFunction : public ROOT::Math::IParamMultiGradFunction { 

public: 

	GaussModelFunction(double amp = 1, double mean = 0, double sigma = 1)  { 
		fPar[0] = amp;
		fPar[1] = mean;
		fPar[2] = sigma; 
	}
	GaussModelFunction(const double *p) { 
		fPar[0] = p[0];
		fPar[1] = p[1];
		fPar[2] = p[2]; 
	}

	
	unsigned int NPar() const { return 3; }
	unsigned int NDim() const { return 1; }

	void SetParameters(const double * p) {
		std::copy(p,p+3,fPar);
	}

	void SetParameters(double amp, double mean, double sigma)  { 
		fPar[0] = amp;
		fPar[1] = mean;
		fPar[2] = sigma; 
	}
	
	const double * Parameters() const { return fPar; }

		
	std::string ParameterName(unsigned int i) const { 
		if (i == 0) return "Amp";
		if (i == 1) return "Mean";
		if (i == 2) return "Sigma";
		return "Invalid";
	}
	
	ROOT::Math::IBaseFunctionMultiDim * Clone() const { 
		return new GaussModelFunction(fPar[0],fPar[1],fPar[2] );
	} 

	void ParameterGradient(const double * x , const double * p, double * grad ) const { 
cerr << "ParameterGradient, params: " << p[0] << " " << p[1] << " " << p[2] << endl;
		if (p == 0) { 
			ParameterGradient(x,fPar,grad); 
			return; 
		}
		double amp =	p[0]; 
		double mean =  p[1];
		double sigma = p[2];
		double y = (x[0]-mean)/sigma;
		double y1 = (x[0]-(mean+0.02))/sigma;
		double y2 = (x[0]-(mean-0.02))/sigma;
		grad[0] =  exp(-0.5*y*y);
		grad[1] = amp*grad[0]*y/sigma;
//		grad[1] = (amp * std::exp( -0.5 * y1 * y1)-amp * std::exp( -0.5 * y2 * y2))/0.04;
		grad[2] = grad[1] * y; 
	}


private:

	double DoEvalPar(const double * x, const double *p) const { 
		if (p == 0) { 
			return DoEvalPar(x,fPar); 
		}
		double amp =	p[0]; 
		double mean =  p[1];
		double sigma = p[2];
		double z =  (x[0]- mean)/ sigma; 
		return amp * std::exp( -0.5 * z * z);
	}

	double DoParameterDerivative(const double *x, const double * p, unsigned int ipar) const { 
cerr << "DoParameterDerivative" << endl;
		double g[3];
		ParameterGradient(x,p,g);
		return g[ipar];
	}

	double fPar[3];

};
	
	

void exampleGradFit(std::string fitterName="Minuit",std::string algoname="Migrad")  { 

	TH1D * h1 = new TH1D("h1","h1",200,-5,5); 
	h1->FillRandom("gaus"); 

	GaussModelFunction  gaus;

	ROOT::Math::MinimizerOptions::SetDefaultMinimizer(fitterName.c_str(),algoname.c_str()); 
	ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.01); 

	ROOT::Fit::BinData d; 
	ROOT::Fit::FillData(d,h1,0); 
	ROOT::Fit::Fitter fitter; 

	gaus.SetParameters(100,0,1);



	std::cout <<"\n\n**********************************************************\n";
	std::cout <<"Chi2 Fit using analytical gradient of Gaussian function \n\n";

	// set first the funciton (needed if parameter limits must be set)
	fitter.SetFunction(gaus);  
	// set limit on sigma 
	fitter.Config().ParSettings(2).SetLimits(0,100);
	fitter.Config().ParSettings(0).SetLimits(0, 200);
	fitter.Config().ParSettings(1).SetLimits(-1, 1);
	//fitter.Config().MinimizerOptions().SetPrintLevel(3);
	fitter.Fit(d);
	fitter.Result().Print(std::cout);

	// to  draw need to create a TF1 from pointer to fit function
	// ponter needs to be alive as TF1 so needs to clone
	if (fitter.Result().FittedFunction() ) {
		ROOT::Math::IParamMultiFunction * ff = dynamic_cast<ROOT::Math::IParamMultiFunction *>
			(fitter.Result().FittedFunction()->Clone() );
		TF1 * f = new TF1("fitFunc",ff,-5,5,0);
		h1->Draw();
		f->Draw("same");
	}

}