#ifndef DFitFunctors_Gaussians_h
#define DFitFunctors_Gaussians_h

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

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

#include "DFitFunctor.h"
#include "DFitUtilities.h"
#include "../PlotData/DStatUtilities.h"

using namespace std;

//NOTE: FIT TYPES CANNOT HAVE UNDERSCORES!!!
//THEY ALSO CANNOT BE "Total" "Signal" OR "Background"

class DFitFunctor_Gaussian : public DFitFunctor
{
	public:
		//input param [0] should be height, will be converted to area
		DFitFunctor_Gaussian(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("Gaussian", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Gaussian Integral";  dParamNames[1] = "Gaussian #mu";  dParamNames[2] = "Gaussian #sigma";

			dInitParams[0] *= dInitParams[2]*sqrt(2.0*TMath::Pi()); //convert from height to area

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

		DFitFunctor_Gaussian(TF1* locFunc) : DFitFunctor("Gaussian", locFunc){}

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

		void Calc_CutRange(double locNumSigmas, double& locCutRangeMin, double& locCutRangeMax) const
		{
			if(dFunc == NULL)
			{
				locCutRangeMin = -9.9E9;
				locCutRangeMax = -9.9E9;
				return;
			}
			double locMean = dFunc->GetParameter(1);
			double locWidth = dFunc->GetParameter(2);
			locCutRangeMin = locMean - locNumSigmas*locWidth;
			locCutRangeMax = locMean + locNumSigmas*locWidth;
		}

		double Get_FitGaussianHeight(void) const
		{
			if(dFunc == NULL)
				return 0.0;
			return dFunc->GetParameter(0)/(dFunc->GetParameter(2)*sqrt(2.0*TMath::Pi()));
		}

		virtual void Calc_Integral(double& locIntegral, double& locIntegralError) const
		{
			locIntegral = dFunc->GetParameter(0);
			locIntegralError = dCovarianceMatrix(0, 0);
		}

		using DFitFunctor::Calc_Integral;

	ClassDef(DFitFunctor_Gaussian, 1)
};

class DFitFunctor_DoubleGaussian : public DFitFunctor
{
	public:
		//input params [0] & [3] should be heights, will be converted to areas
		DFitFunctor_DoubleGaussian(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("DoubleGaussian", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 5)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Gaussian 1 Integral";  dParamNames[1] = "Gaussian #mu";  dParamNames[2] = "Gaussian 1 #sigma";
			dParamNames[3] = "Gaussian 2 Integral";  dParamNames[4] = "Gaussian 2 #sigma";

			dInitParams[0] *= dInitParams[2]*sqrt(2.0*TMath::Pi()); //convert from height to area
			dInitParams[3] *= dInitParams[4]*sqrt(2.0*TMath::Pi()); //convert from height to area

			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]);
			dParamLimits[4] = pair<double, double>(0.0, 100.0*dInitParams[4]);
		}

		DFitFunctor_DoubleGaussian(TF1* locFunc) : DFitFunctor("DoubleGaussian", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Integral1, Mean1, Sigma1, Integral2, Sigma2
			double locValue = locParamArray[0]*TMath::Gaus(locX[0], locParamArray[1], locParamArray[2], kTRUE);
			locValue += locParamArray[3]*TMath::Gaus(locX[0], locParamArray[1], locParamArray[4], kTRUE);
			return locValue;
		}

		pair<double, double> Calc_StandardDeviation(void) const
		{
			//variance = 2nd moment of probability distribution function: Gaussians are summed so variances are summed
				//with multiplicative amplitude factors
			if(dFunc == NULL)
				return pair<double, double>(0.0, 0.0);
			if((dFunc->GetParameter(0) < 0.0) || (dFunc->GetParameter(3) < 0.0))
				return pair<double, double>(0.0, 0.0);

			if(!(dFunc->GetParameter(0) > 0.0) && !(dFunc->GetParameter(3) > 0.0))
				return pair<double, double>(0.0, 0.0);
			if(!(dFunc->GetParameter(0) > 0.0))
				return pair<double, double>(dFunc->GetParameter(4), dFunc->GetParError(4));
			if(!(dFunc->GetParameter(3) > 0.0))
				return pair<double, double>(dFunc->GetParameter(2), dFunc->GetParError(2));

			//locIntegral1*locSigma1*locSigma1
			vector<pair<double, double> > locTerms;
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(0), dFunc->GetParError(0)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(2), dFunc->GetParError(2)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(2), dFunc->GetParError(2)));
			pair<double, double> locNumTerm1 = DStatUtilities::Calc_Product(locTerms, 1.0);

			//locIntegral2*locSigma2*locSigma2
			locTerms.clear();
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(3), dFunc->GetParError(3)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(4), dFunc->GetParError(4)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(4), dFunc->GetParError(4)));
			pair<double, double> locNumTerm2 = DStatUtilities::Calc_Product(locTerms, 1.0);

			//numerator
			pair<double, double> locNumerator = DStatUtilities::Calc_Sum(locNumTerm1, locNumTerm2);

			//denominator
			locTerms.clear();
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(0), dFunc->GetParError(0)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(3), dFunc->GetParError(3)));
			pair<double, double> locDenominator = DStatUtilities::Calc_Sum(locTerms);

			pair<double, double> locStandardDeviation = DStatUtilities::Calc_Ratio(locNumerator, locDenominator);
			locStandardDeviation.first = sqrt(locStandardDeviation.first);
			locStandardDeviation.second /= 2.0*locStandardDeviation.first;

//a = b^1/2
//b = a^2
//da/db = 1/2*b^(-1/2) = 1/(2a)
//sig_a = sig_b*da/db
			return locStandardDeviation;
		}

		void Calc_CutRange(double locNumSigmas, double& locCutRangeMin, double& locCutRangeMax) const
		{
			if(dFunc == NULL)
			{
				locCutRangeMin = -9.9E9;
				locCutRangeMax = -9.9E9;
				return;
			}

			double locMean = dFunc->GetParameter(1);
			pair<double, double> locStandardDeviation = Calc_StandardDeviation();
			locCutRangeMin = locMean - locNumSigmas*locStandardDeviation.first;
			locCutRangeMax = locMean + locNumSigmas*locStandardDeviation.first;
		}

		void Get_FitGaussianHeights(double& locGaussianHeight1, double& locGaussianHeight2) const
		{
			if(dFunc == NULL)
			{
				locGaussianHeight1 = 0.0;
				locGaussianHeight2 = 0.0;
				return;
			}
			locGaussianHeight1 = dFunc->GetParameter(0)/(dFunc->GetParameter(2)*sqrt(2.0*TMath::Pi()));
			locGaussianHeight2 = dFunc->GetParameter(3)/(dFunc->GetParameter(4)*sqrt(2.0*TMath::Pi()));
		}

		virtual void Calc_Integral(double& locIntegral, double& locIntegralError) const
		{
			locIntegral = dFunc->GetParameter(0) + dFunc->GetParameter(3);
			locIntegralError = dCovarianceMatrix(0, 0) + dCovarianceMatrix(3, 3) + 2.0*dCovarianceMatrix(0, 3);
		}

		using DFitFunctor::Calc_Integral;

	ClassDef(DFitFunctor_DoubleGaussian, 1)
};

class DFitFunctor_DoubleGaussian_Fixed : public DFitFunctor
{
	//[0] is not fixed, the rest are
	public:
		//input param [0] should be height, will be converted to area
		DFitFunctor_DoubleGaussian_Fixed(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("DoubleGaussian_Fixed", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 5)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Total Integral";  dParamNames[1] = "Gaussian #mu";  dParamNames[2] = "Gaussian 1 #sigma";
			dParamNames[3] = "Gaussian 1 / 2 Integral Ratio";  dParamNames[4] = "Gaussian 2 #sigma";

			//input is height of first gaussian
			//Integral 1 = height_1*width*sqrt(2*pi)
			//Integral 2 = Integral_1 / dInitParams[3]
			//Total = Integral 1 + Integral 2
			//Total = Integral_1*(1 + 1.0 / dInitParams[3])
			//Total = height_1*width_1*sqrt(2*pi)*(1 + 1.0 / dInitParams[3])
			dInitParams[0] *= dInitParams[2]*sqrt(2.0*TMath::Pi())*(1.0 + 1.0/dInitParams[3]); //convert from height to area
			dParamLimits[0] = pair<double, double>(0.0, 100.0*dInitParams[0]);

			Fix_Parameter(1);
			Fix_Parameter(2);
			Fix_Parameter(3);
			Fix_Parameter(4);
		}

		DFitFunctor_DoubleGaussian_Fixed(TF1* locFunc) : DFitFunctor("DoubleGaussian", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Integral1, Mean1, Sigma1, Integral2, Sigma2
			//Int1 + Int2 = Total
			//Int1 / Int2 = Ratio //Int1 = Ratio * Int2
			//Int2 = Int1/Ratio

			//Ratio * Int2 + Int2 = Total //Int2*(1 + Ratio) = Total //Int2 = Total / (1 + Ratio)
			double locIntegral2 = locParamArray[0] / (1.0 + locParamArray[3]);
			double locIntegral1 = locParamArray[0] - locIntegral2;

			double locValue = locIntegral1*TMath::Gaus(locX[0], locParamArray[1], locParamArray[2], kTRUE);
			locValue += locIntegral2*TMath::Gaus(locX[0], locParamArray[1], locParamArray[4], kTRUE);
			return locValue;
		}

		pair<double, double> Calc_StandardDeviation(void) const
		{
			//variance = 2nd moment of probability distribution function: Gaussians are summed so variances are summed
				//with multiplicative amplitude factors

			if(dFunc == NULL)
				return pair<double, double>(-1.0, 0.0);
			if((dFunc->GetParameter(0) < 0.0) || (dFunc->GetParameter(3) < 0.0))
				return pair<double, double>(-1.0, 0.0); //should be impossible ...

			vector<pair<double, double> > locTerms;
			locTerms.clear();
			locTerms.push_back(pair<double, double>(1.0, 0.0));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(3), dFunc->GetParError(3)));
			pair<double, double> locIntegralTotal(dFunc->GetParameter(0), dFunc->GetParError(0));

			pair<double, double> locIntegral2 = DStatUtilities::Calc_Ratio(locIntegralTotal, locTerms); // [0]/(1 + [3])
			pair<double, double> locIntegral1 = DStatUtilities::Calc_Sum(locIntegralTotal, locIntegral2, 1.0, -1.0); // [0] - Int2

			//locIntegral1*locSigma1*locSigma1
			locTerms.clear();
			locTerms.push_back(locIntegral1);
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(2), dFunc->GetParError(2)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(2), dFunc->GetParError(2)));
			pair<double, double> locNumTerm1 = DStatUtilities::Calc_Product(locTerms, 1.0);

			//locIntegral2*locSigma2*locSigma2
			locTerms.clear();
			locTerms.push_back(locIntegral2);
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(4), dFunc->GetParError(4)));
			locTerms.push_back(pair<double, double>(dFunc->GetParameter(4), dFunc->GetParError(4)));
			pair<double, double> locNumTerm2 = DStatUtilities::Calc_Product(locTerms, 1.0);

			//numerator
			pair<double, double> locNumerator = DStatUtilities::Calc_Sum(locNumTerm1, locNumTerm2);

			//denominator
			locTerms.clear();
			locTerms.push_back(locIntegral1);
			locTerms.push_back(locIntegral2);
			pair<double, double> locDenominator = DStatUtilities::Calc_Sum(locTerms);

			pair<double, double> locStandardDeviation = DStatUtilities::Calc_Ratio(locNumerator, locDenominator);
			locStandardDeviation.first = sqrt(locStandardDeviation.first);
			locStandardDeviation.second /= 2.0*locStandardDeviation.first;

//a = b^1/2
//b = a^2
//da/db = 1/2*b^(-1/2) = 1/(2a)
//sig_a = sig_b*da/db
			return locStandardDeviation;
		}

		void Calc_CutRange(double locNumSigmas, double& locCutRangeMin, double& locCutRangeMax) const
		{
			if(dFunc == NULL)
			{
				locCutRangeMin = -9.9E9;
				locCutRangeMax = -9.9E9;
				return;
			}

			double locMean = dFunc->GetParameter(1);
			pair<double, double> locStandardDeviation = Calc_StandardDeviation();
			locCutRangeMin = locMean - locNumSigmas*locStandardDeviation.first;
			locCutRangeMax = locMean + locNumSigmas*locStandardDeviation.first;
		}

		void Get_FitGaussianHeights(double& locGaussianHeight1, double& locGaussianHeight2) const
		{
			if(dFunc == NULL)
			{
				locGaussianHeight1 = 0.0;
				locGaussianHeight2 = 0.0;
				return;
			}
			locGaussianHeight1 = dFunc->GetParameter(0)/(dFunc->GetParameter(2)*sqrt(2.0*TMath::Pi()));
			locGaussianHeight2 = dFunc->GetParameter(3)/(dFunc->GetParameter(4)*sqrt(2.0*TMath::Pi()));
		}

		virtual void Calc_Integral(double& locIntegral, double& locIntegralError) const
		{
			locIntegral = dFunc->GetParameter(0) + dFunc->GetParameter(3);
			locIntegralError = dCovarianceMatrix(0, 0) + dCovarianceMatrix(3, 3) + 2.0*dCovarianceMatrix(0, 3);
		}

		using DFitFunctor::Calc_Integral;

	ClassDef(DFitFunctor_DoubleGaussian_Fixed, 1)
};

class DFitFunctor_TripleGaussian : public DFitFunctor
{
	public:
		//input params [0] & [3] should be heights, will be converted to areas
		DFitFunctor_TripleGaussian(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("TripleGaussian", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 7)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Gaussian 1 Integral";  dParamNames[1] = "Gaussian #mu";  dParamNames[2] = "Gaussian 1 #sigma";
			dParamNames[3] = "Gaussian 2 Integral";  dParamNames[4] = "Gaussian 2 #sigma";
			dParamNames[5] = "Gaussian 3 Integral";  dParamNames[6] = "Gaussian 3 #sigma";

			dInitParams[0] *= dInitParams[2]*sqrt(2.0*TMath::Pi()); //convert from height to area
			dInitParams[3] *= dInitParams[4]*sqrt(2.0*TMath::Pi()); //convert from height to area
			dInitParams[5] *= dInitParams[6]*sqrt(2.0*TMath::Pi()); //convert from height to area

			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]);
			dParamLimits[4] = pair<double, double>(0.0, 100.0*dInitParams[4]);
			dParamLimits[5] = pair<double, double>(0.0, 100.0*dInitParams[5]);
			dParamLimits[6] = pair<double, double>(0.0, 100.0*dInitParams[6]);
		}

		DFitFunctor_TripleGaussian(TF1* locFunc) : DFitFunctor("TripleGaussian", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Integral1, Mean1, Sigma1, Integral2, Sigma2
			double locValue = locParamArray[0]*TMath::Gaus(locX[0], locParamArray[1], locParamArray[2], kTRUE);
			locValue += locParamArray[3]*TMath::Gaus(locX[0], locParamArray[1], locParamArray[4], kTRUE);
			locValue += locParamArray[5]*TMath::Gaus(locX[0], locParamArray[1], locParamArray[6], kTRUE);
			return locValue;
		}

		double Calc_StandardDeviation(void) const
		{
			//variance = 2nd moment of probability distribution function: Gaussians are summed so variances are summed
				//with multiplicative amplitude factors

			if(dFunc == NULL)
				return -9.9E9;

			double locIntegral1 = dFunc->GetParameter(0);
			double locIntegral2 = dFunc->GetParameter(3);
			double locIntegral3 = dFunc->GetParameter(5);
			double locSigma1 = dFunc->GetParameter(2);
			double locSigma2 = dFunc->GetParameter(4);
			double locSigma3 = dFunc->GetParameter(6);

			double locTerm = locIntegral1*locSigma1*locSigma1 + locIntegral2*locSigma2*locSigma2 + locIntegral3*locSigma3*locSigma3;
			return sqrt(locTerm/(locIntegral1 + locIntegral2 + locIntegral3));
		}

		void Calc_CutRange(double locNumSigmas, double& locCutRangeMin, double& locCutRangeMax) const
		{
			if(dFunc == NULL)
			{
				locCutRangeMin = -9.9E9;
				locCutRangeMax = -9.9E9;
				return;
			}

			double locMean = dFunc->GetParameter(1);
			double locStandardDeviation = Calc_StandardDeviation();
			locCutRangeMin = locMean - locNumSigmas*locStandardDeviation;
			locCutRangeMax = locMean + locNumSigmas*locStandardDeviation;
		}

		void Get_FitGaussianHeights(double& locGaussianHeight1, double& locGaussianHeight2, double& locGaussianHeight3) const
		{
			if(dFunc == NULL)
			{
				locGaussianHeight1 = 0.0;
				locGaussianHeight2 = 0.0;
				locGaussianHeight3 = 0.0;
				return;
			}
			locGaussianHeight1 = dFunc->GetParameter(0)/(dFunc->GetParameter(2)*sqrt(2.0*TMath::Pi()));
			locGaussianHeight2 = dFunc->GetParameter(3)/(dFunc->GetParameter(4)*sqrt(2.0*TMath::Pi()));
			locGaussianHeight3 = dFunc->GetParameter(5)/(dFunc->GetParameter(6)*sqrt(2.0*TMath::Pi()));
		}

		virtual void Calc_Integral(double& locIntegral, double& locIntegralError) const
		{
			locIntegral = dFunc->GetParameter(0) + dFunc->GetParameter(3) + dFunc->GetParameter(5);
			locIntegralError = dCovarianceMatrix(0, 0) + dCovarianceMatrix(3, 3) + dCovarianceMatrix(5, 5);
			locIntegralError += 2.0*dCovarianceMatrix(0, 3) + 2.0*dCovarianceMatrix(0, 5) + 2.0*dCovarianceMatrix(3, 5);
		}

		using DFitFunctor::Calc_Integral;

	ClassDef(DFitFunctor_TripleGaussian, 1)
};

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

			dParamNames[0] = "Gaussian Height";  dParamNames[1] = "Gaussian #mu";
			dParamNames[2] = "Gaussian Low-Side #sigma";  dParamNames[3] = "Gaussian High-Side #sigma";

			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_AsymmetricGaussian(TF1* locFunc) : DFitFunctor("AsymmetricGaussian", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Height, Mean, Low-Side Sigma, High-Side Sigma
			double locValue = locParamArray[0]*(locX[0] < locParamArray[1])*TMath::Gaus(locX[0], locParamArray[1], locParamArray[2]);
			locValue += locParamArray[0]*(locX[0] >= locParamArray[1])*TMath::Gaus(locX[0], locParamArray[1], locParamArray[3]);
			return locValue;
		}

		void Calc_CutRange(double locSignalFraction, double& locCutRangeMin, double& locCutRangeMax) const
		{
			if(dFunc == NULL)
			{
				locCutRangeMin = -9.9E9;
				locCutRangeMax = -9.9E9;
				return;
			}

			double locMean = dFunc->GetParameter(1);
			double locSigma1 = dFunc->GetParameter(2);
			double locSigma2 = dFunc->GetParameter(3);
			double locMaxSigma = (locSigma2 > locSigma1) ? locSigma2 : locSigma1;

			DFitUtilities::Calc_CutRange_SignalFraction(dFunc, locSignalFraction, locCutRangeMin, locCutRangeMax, locMean, locMaxSigma);
		}

	ClassDef(DFitFunctor_AsymmetricGaussian, 1)
};

#endif