#ifndef DFitFunctors_BreitWigners_h
#define DFitFunctors_BreitWigners_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"

using namespace std;

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

class DFitFunctor_NonRelBreitWigner : public DFitFunctor
{
	//input param [0] should be height at centroid; will set different value in constructor though
	public:
		DFitFunctor_NonRelBreitWigner(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("NonRelBreitWigner", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

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

			dInitParams[0] *= dInitParams[2]*TMath::Pi()/2.0; //A = h*gamma*pi/2

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

		DFitFunctor_NonRelBreitWigner(TF1* locFunc) : DFitFunctor("NonRelBreitWigner", locFunc){}

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

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

			double locMean = dFunc->GetParameter(1);
			double locGamma = dFunc->GetParameter(2);
/*
//    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;
		}

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

		using DFitFunctor::Calc_Integral;

	protected:
		DFitFunctor_NonRelBreitWigner(string locFitType, const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor(locFitType, locInitParams, locFuncColor){}

		DFitFunctor_NonRelBreitWigner(string locFitType, TF1* locFunc) : DFitFunctor(locFitType, locFunc){}

	ClassDef(DFitFunctor_NonRelBreitWigner, 1)
};

class DFitFunctor_RelBreitWigner : public DFitFunctor
{
	//input param [0] should be height at centroid; will set different value in constructor though
	public:
		DFitFunctor_RelBreitWigner(const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor("RelBreitWigner", locInitParams, locFuncColor)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Breit-Wigner Height";  dParamNames[1] = "Breit-Wigner #mu";  dParamNames[2] = "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]);
		}

		DFitFunctor_RelBreitWigner(TF1* locFunc) : DFitFunctor("RelBreitWigner", locFunc){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Height, Peak Center, Gamma (Full Width at Half Maximum)
			double locDeltaMassSq = locX[0]*locX[0] - locParamArray[1]*locParamArray[1];
			double locMSqGammaSq = locParamArray[1]*locParamArray[1]*locParamArray[2]*locParamArray[2];
			return locParamArray[0]*locMSqGammaSq/(locDeltaMassSq*locDeltaMassSq + locMSqGammaSq);
		}

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

			double locMean = dFunc->GetParameter(1);
			double locGamma = dFunc->GetParameter(2);
			locCutRangeMin = locMean - locNumGammas*locGamma;
			locCutRangeMax = locMean + locNumGammas*locGamma;
		}

	protected:
		DFitFunctor_RelBreitWigner(string locFitType, const vector<double>& locInitParams, int locFuncColor = kBlack) : 
		DFitFunctor(locFitType, locInitParams, locFuncColor){}

		DFitFunctor_RelBreitWigner(string locFitType, TF1* locFunc) : DFitFunctor(locFitType, locFunc){}

	ClassDef(DFitFunctor_RelBreitWigner, 1)
};

class DFitFunctor_RelBreitWignerGammaDepends : public DFitFunctor
{
	//input param [0] should be height at centroid; will set different value in constructor though
	public:
		//dProductMasses is masses of products (if A -> B, C: Masses are B, C)
		DFitFunctor_RelBreitWignerGammaDepends(const vector<double>& locInitParams, unsigned int locL, const vector<double>& locProductMasses, int locFuncColor = kBlack) : 
		DFitFunctor("RelBreitWignerGammaDepends", locInitParams, locFuncColor), dL(locL), dProductMasses(locProductMasses), dDecayProductRadius(1.0)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Breit-Wigner Height";  dParamNames[1] = "Breit-Wigner #mu";  dParamNames[2] = "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]);
		}

		DFitFunctor_RelBreitWignerGammaDepends(TF1* locFunc, unsigned int locL, const vector<double>& locProductMasses) : 
		DFitFunctor("RelBreitWignerGammaDepends", locFunc), dL(locL), dProductMasses(locProductMasses), dDecayProductRadius(1.0){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Height, Peak Center, Gamma (Full Width at Half Maximum)
			//breit-wigner
			double locGamma = Calc_Gamma(locParamArray[2], locX[0], locParamArray[1], dProductMasses[0], dProductMasses[1]);

			double locValue = locParamArray[0]*locParamArray[1]*locParamArray[1]*locParamArray[2]*locParamArray[2];

			double locDeltaMassSq = locX[0]*locX[0] - locParamArray[1]*locParamArray[1];
			locValue /= locDeltaMassSq*locDeltaMassSq + locParamArray[1]*locParamArray[1]*locGamma*locGamma;
			return locValue;
		}

		double Calc_Gamma(double locGamma0, double locM, double locM0, double locProductMass0, double locProductMass1) const
		{
			double locBarrierFactorSquared = Calc_BarrierFactorSquared(locM, locM0, locProductMass0, locProductMass1);
			double locGamma = locGamma0*locBarrierFactorSquared;

			double locPhaseSpaceSqrtTerm = Calc_PhaseSpaceSqrtTerm(locProductMass0, locProductMass1, locM);
			locPhaseSpaceSqrtTerm /= Calc_PhaseSpaceSqrtTerm(locProductMass0, locProductMass1, locM0);

			locGamma *= pow(locPhaseSpaceSqrtTerm, double(dL) + 0.5);
			locGamma *= pow(locM0/locM, double(2*dL + 2));
			return locGamma;
		}

		double Calc_BarrierFactorSquared(double locM, double locM0, double locProduct1Mass, double locProduct2Mass) const
		{
			//B'_L
			if(dL == 0)
				return 1.0;

			double locD = dDecayProductRadius/0.197327; //the impact parameter (meson radius) //default value is 1 fm, but is converted to units of GeV^(-1)
			double locQ = Calc_PhaseSpaceTerm(locProduct1Mass, locProduct2Mass, locM);
			double locQ0 = Calc_PhaseSpaceTerm(locProduct1Mass, locProduct2Mass, locM0);
			double locZ = locD*locD*locQ*locQ;
			double locZ0 = locD*locD*locQ0*locQ0;

			if(dL == 1)
				return (1.0 + locZ0)/(1.0 + locZ); //Bl' for L = 1
//				return 2.0*locZ/(1.0 + locZ); //Bl for L = 1

			double locZ0Sq = locZ0*locZ0;
			double locZSq = locZ*locZ;

			if(dL == 2)
				return (locZ0Sq + 3.0*locZ0 + 9.0)/(locZSq + 3.0*locZ + 9.0);
//				return 13.0*locZ*locZ/((locZ - 3.0)*(locZ - 3.0) + 9.0*locZ); //Bl for L = 2

			if(dL == 3)
				return (locZ0Sq*locZ0 + 6.0*locZ0Sq + 45.0*locZ0 + 225.0)/(locZSq*locZ + 6.0*locZSq + 45.0*locZ + 225.0);
//				return 277.0*locZ*locZ*locZ/(locZ*(locZ - 15.0)*(locZ - 15.0) + 9.0*(2.0*locZ - 5.0)*(2.0*locZ - 5.0)); //Bl for L = 3

			return 1.0;
		}

		double Calc_PhaseSpaceTerm(double locProduct1Mass, double locProduct2Mass, double locInvariantMass) const
		{
			double locPhaseSpaceSqrtTerm = Calc_PhaseSpaceSqrtTerm(locProduct1Mass, locProduct2Mass, locInvariantMass);
			return sqrt(locPhaseSpaceSqrtTerm)/(2.0*locInvariantMass);
		}

		double Calc_PhaseSpaceSqrtTerm(double locProduct1Mass, double locProduct2Mass, double locInvariantMass) const
		{
			if(locInvariantMass <= (locProduct1Mass + locProduct2Mass))
				return 0.0;

			double locInvariantMassSq = locInvariantMass*locInvariantMass;
			double locMassSum = locProduct1Mass + locProduct2Mass;
			double locDeltaMass = locProduct1Mass - locProduct2Mass;
			double locSqrtTerm = (locInvariantMassSq - locMassSum*locMassSum)*(locInvariantMassSq - locDeltaMass*locDeltaMass);

			if(locSqrtTerm < 0.0)
				return 0.0;
			return locSqrtTerm;
		}

		unsigned int dL; //orbital angular momentum (s-wave = 0, p-wave = 1, etc.)
		vector<double> dProductMasses;
		double dDecayProductRadius; //default is 1 fm //for barrier factor calculation

	protected:
		DFitFunctor_RelBreitWignerGammaDepends(string locFitType, const vector<double>& locInitParams, unsigned int locL, const vector<double>& locProductMasses, int locFuncColor = kBlack) : 
		DFitFunctor(locFitType, locInitParams, locFuncColor), dL(locL), dProductMasses(locProductMasses), dDecayProductRadius(1.0){}

		DFitFunctor_RelBreitWignerGammaDepends(string locFitType, TF1* locFunc, unsigned int locL, const vector<double>& locProductMasses) : 
		DFitFunctor(locFitType, locFunc), dL(locL), dProductMasses(locProductMasses), dDecayProductRadius(1.0){}

	ClassDef(DFitFunctor_RelBreitWignerGammaDepends, 1)
};

class DFitFunctor_PhaseSpaceRelBreitWignerGammaDepends : public DFitFunctor_RelBreitWignerGammaDepends
{
	//input param [0] should be height at centroid; will set different value in constructor though
	public:
		//only supports locNumBodies = 2 or 3:
			//2 is g + A -> B -> C + D
			//3 is g + A -> B + C, C -> D + F
			//in both cases, target is A
		//input param [0] should be height at centroid; will set different value in constructor though
		//dProductMasses is masses of products (2 body: C, D) (3 body: B, C, D, F)
		DFitFunctor_PhaseSpaceRelBreitWignerGammaDepends(const vector<double>& locInitParams, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass, int locFuncColor = kBlack) : 
		DFitFunctor_RelBreitWignerGammaDepends("PhaseSpaceRelBreitWignerGammaDepends", locInitParams, locL, locProductMasses, locFuncColor), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Breit-Wigner Scale Factor";  dParamNames[1] = "Breit-Wigner #mu";  dParamNames[2] = "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]);

			//set init param [0]
			if(dNumBodies == 2)
			{
				double locWSq = dInitParams[1]*dInitParams[1];
				dInitParams[0] = locInitParams[0]*locWSq*(locWSq - dTargetMass*dTargetMass);
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], dInitParams[1]));
			}
			else if(dNumBodies == 3)
			{
				double locWSq = dW*dW;
				double locPhaseSpaceFactor = sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[2], dProductMasses[3], dProductMasses[1]));
				locPhaseSpaceFactor *= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], dW));
				locPhaseSpaceFactor /= 2.0*dProductMasses[1]*locWSq*(locWSq - dTargetMass*dTargetMass);
				dInitParams[0] = locInitParams[0]/locPhaseSpaceFactor;
//cout << "input 0 value, init 0 value, init phase space factor (f) = " << locInitParams[0] << ", " << dInitParams[0] << ", " << locPhaseSpaceFactor << endl;
			}
		}

		DFitFunctor_PhaseSpaceRelBreitWignerGammaDepends(TF1* locFunc, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass) : 
		DFitFunctor_RelBreitWignerGammaDepends("PhaseSpaceRelBreitWignerGammaDepends", locFunc, locL, locProductMasses), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Scale-factor, Peak Center, Gamma (Full Width at Half Maximum)
			//breit-wigner
//cout << "num bodies, W, masses 0/1/2/3 = " << dNumBodies << ", " << dW << ", " << dProductMasses[0] << ", " << dProductMasses[1] << ", " << dProductMasses[2] << ", " << dProductMasses[3] << endl;
			//gamma
			double locGamma = 0.0;
			if(dNumBodies == 2)
				locGamma = Calc_Gamma(locParamArray[2], locX[0], locParamArray[1], dProductMasses[0], dProductMasses[1]);
			else if(dNumBodies == 3)
				locGamma = Calc_Gamma(locParamArray[2], locX[0], locParamArray[1], dProductMasses[2], dProductMasses[3]);
			else
				return 0.0;
//cout << "gamma = " << locGamma << endl;
			//phase space:
			double locValue = 0.0;
			if(dNumBodies == 2)
			{
				double locXSq = locX[0]*locX[0];
				locValue = sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], locX[0]));
				locValue /= locXSq*(locXSq - dTargetMass*dTargetMass);
			}
			else if(dNumBodies == 3)
			{
				double locWSq = dW*dW;
				locValue = sqrt(Calc_PhaseSpaceSqrtTerm(locX[0], dProductMasses[0], dW));
				locValue *= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[2], dProductMasses[3], locX[0]));
				locValue /= 2.0*locX[0]*locWSq*(locWSq - dTargetMass*dTargetMass);
			}
			else
				return 0.0;
//cout << "phase space = " << locValue << endl;
//cout << "params = " << locParamArray[0] << ", " << locParamArray[1] << ", " << locParamArray[2] << endl;
			//b-w
			locValue *= locParamArray[0]*locParamArray[1]*locParamArray[1]*locParamArray[2]*locParamArray[2];
			double locDeltaMassSq = locX[0]*locX[0] - locParamArray[1]*locParamArray[1];
			locValue /= locDeltaMassSq*locDeltaMassSq + locParamArray[1]*locParamArray[1]*locGamma*locGamma;
//cout << "value = " << locValue << endl;
			return locValue;
		}

	protected:
		unsigned int dNumBodies;
		double dW;
		double dTargetMass;

		DFitFunctor_PhaseSpaceRelBreitWignerGammaDepends(string locFitType, const vector<double>& locInitParams, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass, int locFuncColor = kBlack) : 
		DFitFunctor_RelBreitWignerGammaDepends(locFitType, locInitParams, locL, locProductMasses, locFuncColor), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass){}

		DFitFunctor_PhaseSpaceRelBreitWignerGammaDepends(string locFitType, TF1* locFunc, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass) : 
		DFitFunctor_RelBreitWignerGammaDepends(locFitType, locFunc, locL, locProductMasses), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass){}

	ClassDef(DFitFunctor_PhaseSpaceRelBreitWignerGammaDepends, 1)
};

class DFitFunctor_PhaseSpaceRelBreitWignerGammaDependsMassSq : public DFitFunctor_RelBreitWignerGammaDepends
{
	//input param [0] should be height at centroid; will set different value in constructor though
	public:
		//only supports locNumBodies = 2 or 3:
			//2 is g + A -> B -> C + D
			//3 is g + A -> B + C, C -> D + F
			//in both cases, target is A
		//input param [0] should be height at centroid; will set different value in constructor though
		//dProductMasses is masses of products (2 body: C, D) (3 body: B, C, D, F)
		DFitFunctor_PhaseSpaceRelBreitWignerGammaDependsMassSq(const vector<double>& locInitParams, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass, int locFuncColor = kBlack) : 
		DFitFunctor_RelBreitWignerGammaDepends("PhaseSpaceRelBreitWignerGammaDependsMassSq", locInitParams, locL, locProductMasses, locFuncColor), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Breit-Wigner Scale Factor";  dParamNames[1] = "Breit-Wigner #mu";  dParamNames[2] = "Breit-Wigner #Gamma (FWHM)";

			//set init param [0]
			if(dNumBodies == 2)
			{
				double locWSq = dInitParams[1]*dInitParams[1];
				dInitParams[0] = locInitParams[0]*locWSq*(locWSq - dTargetMass*dTargetMass);
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], dInitParams[1]));
			}
			else if(dNumBodies == 3)
			{
				double locWSq = dW*dW;
				dInitParams[0] = locInitParams[0]*2.0*dProductMasses[1]*locWSq*(locWSq - dTargetMass*dTargetMass);
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], dW));
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[2], dProductMasses[3], dProductMasses[1]));
			}

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

		DFitFunctor_PhaseSpaceRelBreitWignerGammaDependsMassSq(TF1* locFunc, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass) : 
		DFitFunctor_RelBreitWignerGammaDepends("PhaseSpaceRelBreitWignerGammaDependsMassSq", locFunc, locL, locProductMasses), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass){}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//locParamArray Indices are: Scale-factor, Peak Center, Gamma (Full Width at Half Maximum)
			//breit-wigner

			double locSqrtX = sqrt(locX[0]);

			//gamma
			double locGamma = 0.0;
			if(dNumBodies == 2)
				locGamma = Calc_Gamma(locParamArray[2], locSqrtX, locParamArray[1], dProductMasses[0], dProductMasses[1]);
			else if(dNumBodies == 3)
				locGamma = Calc_Gamma(locParamArray[2], locSqrtX, locParamArray[1], dProductMasses[2], dProductMasses[3]);
			else
				return 0.0;

			//phase space:
			double locValue = 0.0;
			if(dNumBodies == 2)
			{
				locValue = sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], locSqrtX));
				locValue /= locX[0]*(locX[0] - dTargetMass*dTargetMass);
			}
			else if(dNumBodies == 3)
			{
				double locWSq = dW*dW;
				locValue = sqrt(Calc_PhaseSpaceSqrtTerm(locSqrtX, dProductMasses[0], dW));
				locValue *= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[2], dProductMasses[3], locSqrtX));
				locValue /= 2.0*locSqrtX*locWSq*(locWSq - dTargetMass*dTargetMass);
			}
			else
				return 0.0;

			//b-w
			locValue *= locParamArray[0]*locParamArray[1]*locParamArray[1]*locParamArray[2]*locParamArray[2];
			double locDeltaMassSq = locX[0] - locParamArray[1]*locParamArray[1];
			locValue /= locDeltaMassSq*locDeltaMassSq + locParamArray[1]*locParamArray[1]*locGamma*locGamma;

			return locValue;
		}

	protected:
		unsigned int dNumBodies;
		double dW;
		double dTargetMass;

		DFitFunctor_PhaseSpaceRelBreitWignerGammaDependsMassSq(string locFitType, const vector<double>& locInitParams, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass, int locFuncColor = kBlack) : 
		DFitFunctor_RelBreitWignerGammaDepends(locFitType, locInitParams, locL, locProductMasses, locFuncColor), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass){}

		DFitFunctor_PhaseSpaceRelBreitWignerGammaDependsMassSq(string locFitType, TF1* locFunc, unsigned int locL, const vector<double>& locProductMasses, unsigned int locNumBodies, double locW, double locTargetMass) : 
		DFitFunctor_RelBreitWignerGammaDepends(locFitType, locFunc, locL, locProductMasses), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass){}

	ClassDef(DFitFunctor_PhaseSpaceRelBreitWignerGammaDependsMassSq, 1)
};

class DFitFunctor_PhaseSpaceNonRelBreitWigner : public DFitFunctor_NonRelBreitWigner
{
	public:
		//only supports locNumBodies = 2 or 3:
			//2 is g + A -> B -> C + D
			//3 is g + A -> B + C, C -> D + F
			//in both cases, target is A
		//input param [0] should be height at centroid; will set different value in constructor though
		//dProductMasses is masses of products (2 body: C, D) (3 body: B, C, D, F)
		DFitFunctor_PhaseSpaceNonRelBreitWigner(const vector<double>& locInitParams, double locW, double locTargetMass, unsigned int locNumBodies, const vector<double>& locProductMasses, int locFuncColor = kBlack) : 
		DFitFunctor_NonRelBreitWigner("PhaseSpaceNonRelBreitWigner", locInitParams, locFuncColor), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass), dProductMasses(locProductMasses)
		{
			if(locInitParams.size() != 3)
			{
				cout << "ERROR: INCORRECT NUMBER OF PARAMETERS IN DFitFunctor CONSTRUCTOR. EXITING" << endl;
				abort();
			}

			dParamNames[0] = "Breit-Wigner Scale Factor";  dParamNames[1] = "Breit-Wigner #mu";  dParamNames[2] = "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]);

			//set init param [0]
			double locWSq = dW*dW;
			if(dNumBodies == 2)
			{
				dInitParams[0] = locInitParams[0]*locWSq*(locWSq - dTargetMass*dTargetMass);
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], dW));
			}
			else if(dNumBodies == 3)
			{
				dInitParams[0] = locInitParams[0]*2.0*dProductMasses[1]*locWSq*(locWSq - dTargetMass*dTargetMass);
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], dW));
				dInitParams[0] /= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[2], dProductMasses[3], dProductMasses[1]));
			}
		}

		DFitFunctor_PhaseSpaceNonRelBreitWigner(TF1* locFunc, double locW, double locTargetMass, unsigned int locNumBodies, const vector<double>& locProductMasses) : 
		DFitFunctor_NonRelBreitWigner("PhaseSpaceNonRelBreitWigner", locFunc), 
		dNumBodies(locNumBodies), dW(locW), dTargetMass(locTargetMass), dProductMasses(locProductMasses){}

		double Calc_PhaseSpaceSqrtTerm(double locProduct1Mass, double locProduct2Mass, double locInvariantMass) const
		{
			if(locInvariantMass <= (locProduct1Mass + locProduct2Mass))
				return 0.0;

			double locInvariantMassSq = locInvariantMass*locInvariantMass;
			double locMassSum = locProduct1Mass + locProduct2Mass;
			double locDeltaMass = locProduct1Mass - locProduct2Mass;
			double locSqrtTerm = (locInvariantMassSq - locMassSum*locMassSum)*(locInvariantMassSq - locDeltaMass*locDeltaMass);

			if(locSqrtTerm < 0.0)
				return 0.0;
			return locSqrtTerm;
		}

		virtual double operator()(double* locX, double* locParamArray)
		{
			//resonance portion
			double locGammaSq = locParamArray[2]*locParamArray[2];
			double locDeltaMass = locX[0] - locParamArray[1];
			double locValue = locParamArray[0]*locGammaSq/(4.0*locDeltaMass*locDeltaMass + locGammaSq);

			//phase space portion
			double locXSq = locX[0]*locX[0];
			if(dNumBodies == 2)
			{
				locValue *= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[0], dProductMasses[1], locX[0]));
				locValue /= locXSq*(locXSq - dTargetMass*dTargetMass);
			}
			else if(dNumBodies == 3)
			{
				double locWSq = dW*dW;
				locValue *= sqrt(Calc_PhaseSpaceSqrtTerm(locX[0], dProductMasses[0], dW));
				locValue *= sqrt(Calc_PhaseSpaceSqrtTerm(dProductMasses[2], dProductMasses[3], locX[0]));
				locValue /= 2.0*locX[0]*locWSq*(locWSq - dTargetMass*dTargetMass);
			}
			else
				return 0.0;
			return locValue;
		}

		virtual void Calc_Integral(double& locIntegral, double& locIntegralError) const
		{
			DFitFunctor::Calc_Integral(locIntegral, locIntegralError);
		}

		using DFitFunctor::Calc_Integral;

		unsigned int dNumBodies;
		double dW;
		double dTargetMass;
		vector<double> dProductMasses;

	ClassDef(DFitFunctor_PhaseSpaceNonRelBreitWigner, 1)
};

#endif