#ifndef DFitFunctors_Extended_h
#define DFitFunctors_Extended_h

#include "DFitFunctor.h"

using namespace std;

class DFitFunctor_1OverX : public DFitFunctor
{
	public:
		DFitFunctor_1OverX(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("1OverX", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}
		}

		DFitFunctor_1OverX(TF1* locFunc) : DFitFunctor("1OverX", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			return locParamArray[0]/(locParamArray[1] + pow(locX[0], locParamArray[2]));
		}

   ClassDef(DFitFunctor_1OverX, 1)
};

class DFitFunctor_ModifiedTanH : public DFitFunctor
{
	public:
		DFitFunctor_ModifiedTanH(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("ModifiedTanH", locInitParams, locFuncColor), dExponentialScaleFactor(1.0), dMultiplierScaleFactor(1.0)
		{
			if(locInitParams.size() != 4)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Modified tanh #it{#alpha}";  dParamNames[1] = "Modified tanh #it{#beta}";
			dParamNames[2] = "Modified tanh #it{A}";  dParamNames[3] = "Modified tanh #it{B}";

			dParamLimits[0] = pair<double, double>(0.0, 100.0*dInitParams[0]); // bad terms for exponential powers: won't converge properly
			// bad terms for multipliers: function will either be negative or will contain singularity
			dParamLimits[2] = pair<double, double>(0.0, 100.0*dInitParams[2]);
			dParamLimits[3] = pair<double, double>(0.0, 100.0*dInitParams[3]);
		}

		DFitFunctor_ModifiedTanH(TF1* locFunc) : DFitFunctor("ModifiedTanH", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			if(locParamArray[0] > locParamArray[1])
				return 0.0; // bad terms for exponential powers: won't converge properly

			double locNumerator = locParamArray[2]*TMath::Exp(locParamArray[0]*locX[0]/dExponentialScaleFactor)/dMultiplierScaleFactor; //: F_t = a/(b+1)
			double locDenominator = locParamArray[3]*TMath::Exp(locParamArray[1]*locX[0]/dExponentialScaleFactor)/dMultiplierScaleFactor + 1.0;
			if(!(fabs(locDenominator) > 0.0))
				return 0.0;
			return locNumerator/locDenominator;
		}

		static void Find_TanHInitGuessValues(double locM, double locT, double locF_t, double locF_m, double& locBeta, double& locAlpha, double& locB, double& locA)
		{
			//alpha sqrt term is = sqrt(1.0 - 6.0*locG + locG*locG);
				//where locG = TMath::Exp(-1.0*locBeta*(locM - locT)); //locG is locF on paper
			//roots are locG = 3 +/- 2*sqrt(2)
				//beta > 0, locM > locT, so G = e^(-Ax) with A & x > 0: G always < 1
				//so crosses 0 at locG = 3 - 2*sqrt(2)
					//3 - 2*sqrt(2) = e^(-1.0*locBeta*(locM - locT))
					//ln(3 - 2*sqrt(2)) = -1.0*locBeta*(locM - locT)
					//crosses 0 at beta = -1.0*ln(3 - 2*sqrt(2))/(locM - locT)
			double locRangeMin = -1.0*log(3.0 - 2.0*sqrt(2.0))/(locM - locT);
			double locRangeMax = 10.0*locRangeMin;
			unsigned int locNumXPoints = 10000.0;

			TF1* locFunc = new TF1("BetaFinder", DFitFunctor_ModifiedTanH::Calc_CrazyTanH, locRangeMin, locRangeMax, 3);
			locFunc->SetParameters(locM, locT, locF_t);
			locFunc->SetNpx(locNumXPoints);
			locBeta = locFunc->GetX(locF_m, locRangeMin, locRangeMax);//, 1.E-10, 100, true)
			delete locFunc;

			double locG = TMath::Exp(-1.0*locBeta*(locM - locT)); //locG is locF on paper
			double loc1MinusG = 1.0 - locG;

			double locSqRootTerm = sqrt(1.0 - 6.0*locG + locG*locG);
			double locAlpha1 = locBeta*(loc1MinusG + locSqRootTerm)/(2.0*loc1MinusG);
			double locAlpha2 = locBeta*(loc1MinusG - locSqRootTerm)/(2.0*loc1MinusG);

			locAlpha = locAlpha1;
			if(!((locAlpha > -1.0) || (locAlpha < 1.0)))
				locAlpha = locAlpha2; //alpha1 is NaN
			else if((locAlpha1 > locBeta) || (locAlpha1 <= 0.0))
				locAlpha = locAlpha2;

			locB = locAlpha*TMath::Exp(-1.0*locBeta*(locM - locT))/(locBeta - locAlpha); //force it to have a maximum at x = locM
			locA = (locB + 1.0)*locF_t;
		}

		static double Calc_CrazyTanH(double *locX, double *locParamArray)
		{
			double locBeta = locX[0];
			if(locBeta <= 0.0)
				return 0.0; //bad term for exponential power: won't converge properly
			double locM = locParamArray[0];
			double locT = locParamArray[1];
			double locF_t = locParamArray[2];

			double locG = TMath::Exp(-1.0*locBeta*(locM - locT)); //locG is locF on paper
			double loc1MinusG = 1.0 - locG;

			double locSqRootTerm = sqrt(1.0 - 6.0*locG + locG*locG);
			double locAlpha1 = locBeta*(loc1MinusG + locSqRootTerm)/(2.0*loc1MinusG);
			double locAlpha2 = locBeta*(loc1MinusG - locSqRootTerm)/(2.0*loc1MinusG);

			double locAlpha = locAlpha1;
			if(!((locAlpha > -1.0) || (locAlpha < 1.0)))
				locAlpha = locAlpha2; //alpha1 is NaN
			else if((locAlpha1 > locBeta) || (locAlpha1 <= 0.0))
				locAlpha = locAlpha2;

			double locB = locAlpha*TMath::Exp(-1.0*locBeta*(locM - locT))/(locBeta - locAlpha); //force it to have a maximum at x = locM
			double locA = (locB + 1.0)*locF_t;
			double locF_m = locA*TMath::Exp(locAlpha*(locM - locT))/(locB*TMath::Exp(locBeta*(locM - locT)) + 1.0);

			//cout << "b, a1, a2, a, B, A, F_m = " << locBeta << ", " << locAlpha1 << ", " << locAlpha2 << ", " << locAlpha << ", " << locB << ", " << locA << ", " << locF_m << endl;
			if(!((locAlpha > -1.0) || (locAlpha < 1.0)))
				return 0.0; //NaN

			if((locAlpha <= 0.0) || (locAlpha > locBeta))
				return 0.0; // bad terms for exponential power: won't converge properly
			if((locA <= 0.0) || (locB <= 0.0))
				return 0.0; // bad terms for multipliers: function will either be negative or will contain singularity

			return locF_m;
		}

		//to make sure that the error matrix converges properly
		double dExponentialScaleFactor;
		double dMultiplierScaleFactor;

   ClassDef(DFitFunctor_ModifiedTanH, 1)
};

class DFitFunctor_ModifiedTanHOneExponentTerm : public DFitFunctor
{
	public:
		DFitFunctor_ModifiedTanHOneExponentTerm(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("ModifiedTanHOneExponentTerm", locInitParams, locFuncColor), dExponentialScaleFactor(1.0), dMultiplierScaleFactor(1.0)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Modified tanh #it{#alpha}";
			dParamNames[1] = "Modified tanh #it{A}";  dParamNames[2] = "Modified tanh #it{B}";

			dParamLimits[0] = pair<double, double>(0.0, 100.0*dInitParams[0]); // bad terms for exponential power: won't converge properly
			// bad terms for multipliers: function will either be negative or will contain singularity
			dParamLimits[1] = pair<double, double>(0.0, 100.0*dInitParams[1]);
			dParamLimits[2] = pair<double, double>(0.0, 100.0*dInitParams[2]);
		}

		DFitFunctor_ModifiedTanHOneExponentTerm(TF1* locFunc) : DFitFunctor("ModifiedTanHOneExponentTerm", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			double locNumerator = locParamArray[1]/dMultiplierScaleFactor;
			double locDenominator = locParamArray[2]/dMultiplierScaleFactor + TMath::Exp(-1.0*locParamArray[0]*locX[0]/dExponentialScaleFactor);

			if(!(fabs(locDenominator) > 0.0))
				return 0.0;
			return locNumerator/locDenominator;
		}

		//to make sure that the error matrix converges properly
		double dExponentialScaleFactor;
		double dMultiplierScaleFactor;

   ClassDef(DFitFunctor_ModifiedTanHOneExponentTerm, 1)
};

#endif