#ifndef DFitFunctors_Miscellaneous_h
#define DFitFunctors_Miscellaneous_h

#include <vector>
#include <string>
#include <sstream>

#include "TMath.h"
#include "TObject.h"
#include "TRandom3.h"
#include "TH1D.h"
#include "TNtuple.h"
#include "TMatrixDSym.h"
#include "TF1.h"

#include "DFitFunctor.h"

using namespace std;

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

			for(size_t loc_i = 0; loc_i < dParamNames.size(); ++loc_i)
			{
				ostringstream locParamName;
				locParamName << "p_{" << loc_i << "}";
				dParamNames[loc_i] = locParamName.str();
			}
		}

		DFitFunctor_Polynomial(TF1* locFunc) : DFitFunctor("Polynomial", locFunc),
		dPolynomialOrder(locFunc->GetNpar() - 1){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: [0] + [1]*x + [2]*x*x + ...
			double locValue = 0.0;
			for(unsigned int loc_i = 0; loc_i <= dPolynomialOrder; ++loc_i)
			{
				double locTerm = locParamArray[loc_i];
				for(unsigned int loc_j = 0; loc_j < loc_i; ++loc_j)
					locTerm *= locX[0];
				locValue += locTerm;
			}
			return locValue;
		}

	private:
		 unsigned int dPolynomialOrder; //0 for flat, 1 for linear, etc.

	ClassDef(DFitFunctor_Polynomial, 1)
};

class DFitFunctor_LinearAboveZero : public DFitFunctor
{
	public:
		//input vector is intercept, slope
		//fit params are f(locFitRangeMin), f(locFitRangeMax)
		DFitFunctor_LinearAboveZero(const vector<double>& locInitParams, double locFitRangeMin, double locFitRangeMax, int locFuncColor = kBlack) : 
		DFitFunctor("LinearAboveZero", locInitParams, locFuncColor), dFitRangeMin(locFitRangeMin), dFitRangeMax(locFitRangeMax)
		{
			if(locInitParams.size() != 2)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			//convert to new coordinate system
			//y = mx + b
			dInitParams[0] = locInitParams[1]*locFitRangeMin + locInitParams[0];
			if(dInitParams[0] < 0.0)
				dInitParams[0] = 0.0;

			dInitParams[1] = locInitParams[1]*locFitRangeMax + locInitParams[0];
			if(dInitParams[1] < 0.0)
				dInitParams[1] = 0.0;

			for(size_t loc_i = 0; loc_i < dParamNames.size(); ++loc_i)
			{
				ostringstream locParamName;
				locParamName << "y_{" << loc_i << "}";
				dParamNames[loc_i] = locParamName.str();
			}

			//Force above zero!
			double locMaxLimit = (dInitParams[0] > 0.0) ? 100.0*dInitParams[0] : 100.0;
			dParamLimits[0] = pair<double, double>(0.0, locMaxLimit);
			locMaxLimit = (dInitParams[1] > 0.0) ? 100.0*dInitParams[1] : 100.0;
			dParamLimits[1] = pair<double, double>(0.0, locMaxLimit);
		}

		DFitFunctor_LinearAboveZero(TF1* locFunc) : DFitFunctor("LinearAboveZero", locFunc)
		{
			dFunc->GetRange(dFitRangeMin, dFitRangeMax);
		}

		virtual double operator()(double* locX, double* locParamArray)
		{
			double locIntercept = 0.0, locSlope = 0.0;
			Calc_PolyParms(locParamArray[0], locParamArray[1], locIntercept, locSlope);
			return (locSlope*locX[0] + locIntercept);
		}

		void Calc_PolyParms(double& locIntercept, double& locSlope)
		{
			if(dFunc == NULL)
				return;
			Calc_PolyParms(dFunc->GetParameter(0), dFunc->GetParameter(1), locIntercept, locSlope);
		}

		void Calc_PolyParms(double locParam1, double locParam2, double& locIntercept, double& locSlope)
		{
			locSlope = (locParam2 - locParam1)/(dFitRangeMax - dFitRangeMin);
			locIntercept = locParam2 - locSlope*dFitRangeMax; //y = mx + b, b = y - mx
		}

	private:

		double dFitRangeMin;
		double dFitRangeMax;

	ClassDef(DFitFunctor_LinearAboveZero, 1)
};

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

			dParamNames[0] = "Voigtian Integral";  dParamNames[1] = "Voigtian #mu";
			dParamNames[2] = "Gaussian #sigma";  dParamNames[3] = "Breit-Wigner #Gamma (FWHM)";

			dParamLimits[0] = pair<double, double>(0.0, 100.0*dInitParams[0]);
			dParamLimits[2] = pair<double, double>(0.0, 100.0*dInitParams[2]);
			dParamLimits[3] = pair<double, double>(0.0, 100.0*dInitParams[3]);
		}

		DFitFunctor_Voigtian(TF1* locFunc) : DFitFunctor("Voigtian", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Integral, Peak Center, Gaussian Sigma, Gamma (Full Width at Half Maximum)
			return locParamArray[0]*TMath::Voigt(locX[0] - locParamArray[1], locParamArray[2], locParamArray[3]);
		}

		void Calc_CutRange(double locNumGammas, double& locCutRangeMin, double& locCutRangeMax) const
		{
			//THIS IS INACCURATE: assumes Breit-Wigner dominates
			double locMean = dFunc->GetParameter(1);
			double locGamma = dFunc->GetParameter(3);

//    double locFactor = 0.5*locGamma*tan(-0.5*locSignalPercentage*TMath::Pi());
//    locCutRangeMin = locMean + locFactor;
//    locCutRangeMax = locMean - locFactor;

			locCutRangeMin = locMean - locNumGammas*locGamma;
			locCutRangeMax = locMean + locNumGammas*locGamma;
		}
	ClassDef(DFitFunctor_Voigtian, 1)
};

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

			dParamNames[0] = "Landau Height";  dParamNames[1] = "Landau #mu";  dParamNames[2] = "Landau #sigma";

			dParamLimits[0] = pair<double, double>(0.0, 100.0*dInitParams[0]);
			dParamLimits[2] = pair<double, double>(0.0, 100.0*dInitParams[2]);
		}

		DFitFunctor_Landau(TF1* locFunc) : DFitFunctor("Landau", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Height, Mean, Sigma
			return locParamArray[0]*TMath::Landau(locX[0], locParamArray[1], locParamArray[2]);
		}

	ClassDef(DFitFunctor_Landau, 1)
};

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

			dParamNames[0] = "Exponential Scale Factor";  dParamNames[1] = "Exponential Exponent Factor";
		}

		DFitFunctor_Exponential(TF1* locFunc) : DFitFunctor("Exponential", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			return locParamArray[0]*TMath::Exp(locParamArray[1]*locX[0]);
		}

	ClassDef(DFitFunctor_Exponential, 1)
};

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

			dParamNames[0] = "Scale Factor";  dParamNames[1] = "Mean Lifetime";
		}

		DFitFunctor_Lifetime(TF1* locFunc) : DFitFunctor("Lifetime", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			return locParamArray[0]*TMath::Exp(-1.0*locX[0]/locParamArray[1]);
		}

	ClassDef(DFitFunctor_Lifetime, 1)
};

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

			dParamNames[0] = "Square Root Scale Factor";  dParamNames[1] = "Square Root Stretch Factor";  dParamNames[2] = "Square Root Offset";
		}

		DFitFunctor_SquareRoot(TF1* locFunc) : DFitFunctor("SquareRoot", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			if((locParamArray[2] > locX[0]) || (locParamArray[1] < 0.0))
				return 0.0;
			return locParamArray[0]*sqrt(locParamArray[1]*(locX[0] - locParamArray[2]));
		}

	ClassDef(DFitFunctor_SquareRoot, 1)
};

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

		DFitFunctor_SphericalHarmonic(TF1* locFunc) : DFitFunctor("SphericalHarmonic", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//http://mathworld.wolfram.com/SphericalHarmonic.html
			//cross terms are zero: functions are orthogonal, and at a given L, integrating over phi (0, 2pi) zeros them. 
			double locCosSq = locX[0]*locX[0];
			double locSinSq = 1.0 - locCosSq;
			//Y_L^(Lz) i.e. Y_L^M

			if(dL == 0) // |[0]^(1/2)*Y_0^0|^2
				return locParamArray[0]*locParamArray[0]/2.0;

			else if(dL == 1) // |([0]^(1/2)*(Y_1^-1 + Y_1^1) + [1]^(1/2)*Y_1^0)|^2
				return 0.75*2.0*locParamArray[0]*locSinSq + 1.5*locParamArray[1]*locCosSq;

			else if(dL == 2) // |([0]^(1/2)*(Y_2^-2 + Y_2^2) + [1]^(1/2)*(Y_2^-1 + Y_2^1) + [2]^(1/2)*Y_2^0)|^2
			{
				double locValue = 15.0*2.0*locParamArray[0]*locSinSq*locSinSq/16.0 + 15.0*locParamArray[1]*locCosSq*locSinSq/4.0;
				return locValue + 5.0*locParamArray[2]*(3.0*locCosSq - 1.0)*(3.0*locCosSq - 1.0)/8.0;
			}

			else if(dL == 3) // |([0]^(1/2)*(Y_3^-3 + Y_3^3) + [1]^(1/2)*(Y_3^-2 + Y_3^2) + [2]^(1/2)*(Y_3^-1 + Y_3^1) + [3]^(1/2)*Y_3^0)|^2
			{
				double locValue = 35.0*2.0*locParamArray[0]*locSinSq*locSinSq*locSinSq/32.0 + 105.0*locParamArray[1]*locCosSq*locSinSq*locSinSq/16.0;
				locValue += 21.0*2.0*locParamArray[2]*locSinSq*(5.0*locCosSq - 1.0)*(5.0*locCosSq - 1.0)/32.0;
				locValue += 7.0*locParamArray[3]*(25.0*locCosSq*locCosSq*locCosSq - 30.0*locCosSq*locCosSq + 9.0*locCosSq)/8.0;
				return locValue;
			}

			else
				return 1.0;
		}

		unsigned int dL;

	ClassDef(DFitFunctor_SphericalHarmonic, 1)
};

class DFitFunctor_Cosine : public DFitFunctor
{
	public:
		DFitFunctor_Cosine(const vector<double>& locInitParams, bool locDegreesFlag = true, int locFuncColor = kBlack) : 
		DFitFunctor("Cosine", locInitParams, locFuncColor), dDegreesFlag(locDegreesFlag)
		{
			if(locInitParams.size() != 2)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Amplitude";
			if(locDegreesFlag)
				dParamNames[1] = "Phase (Degrees)";
			else
				dParamNames[1] = "Phase (Radians)";
		}

		DFitFunctor_Cosine(TF1* locFunc, bool locDegreesFlag = true) : DFitFunctor("Cosine", locFunc),
		dDegreesFlag(locDegreesFlag){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			double locXAngle = dDegreesFlag ? locX[0]*TMath::Pi()/180.0 : locX[0];
			return locParamArray[0]*cos(2.0*(locXAngle - locParamArray[1]));
		}

	private:
		bool dDegreesFlag;

	ClassDef(DFitFunctor_Cosine, 1)
};

#endif