#ifndef DFitUtilities_h
#define DFitUtilities_h

#include <math.h>
#include <iostream>

#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TDirectory.h"

using namespace std;

class DFitUtilities : public TObject
{
	public:
		inline static bool Check_ValidFitRange(TH1* locHist, double locFitRangeMin, double locFitRangeMax)
		{
			if((locFitRangeMax <= locHist->GetXaxis()->GetXmin()) || (locFitRangeMin >= locHist->GetXaxis()->GetXmax()))
			{
				cout << "ERROR: FIT RANGE OUTSIDE OF DATA RANGE FOR OBJECT: " << locHist->GetName() << ". FIT NOT PERFORMED." << endl;
				return false;
			}
/*
			int locFirstFitBin = locHist->FindBin(locFitRangeMin);
			int locLastFitBin = locHist->FindBin(locFitRangeMax);
			if(!(locHist->Integral(locFirstFitBin, locLastFitBin) > 0.0))
			{
				cout << "ERROR: NO DATA IN FIT RANGE FOR HIST: " << locHist->GetName() << ". FIT NOT PERFORMED." << endl;
				return false;
			}
*/
			return true;
		}

		inline static TF1* Guess_Params_Linear_Fit(TH1* locHist, double locFitRangeMin, double locFitRangeMax, double& locIntercept, double& locSlope)
		{
			if(!Check_ValidFitRange(locHist, locFitRangeMin, locFitRangeMax))
				return NULL;
			TF1* locFunc = new TF1("temp", "[0] + [1]*x", locFitRangeMin, locFitRangeMax);
			locHist->Fit(locFunc, "QRN+");
			locIntercept = locFunc->GetParameter(0);
			locSlope = locFunc->GetParameter(1);
			return locFunc;
		}

		inline static void Guess_Params_Linear(TH1* locHist, double locFitRangeMin, double locFitRangeMax, double& locIntercept, double& locSlope, bool locForceGreaterThanZeroFlag = false)
		{
			if(!Check_ValidFitRange(locHist, locFitRangeMin, locFitRangeMax))
				return;

			//use avg bin content within +/- 2 bins
			int locXMinBin = locHist->FindBin(locFitRangeMin);
			double locYAtXMin = 0.0;
			int locMinBinToSearch = ((locXMinBin - 2) <= 0) ? 1 : locXMinBin - 2;
			for(int loc_i = locMinBinToSearch; loc_i <= locXMinBin + 2; ++loc_i)
				locYAtXMin += locHist->GetBinContent(loc_i);
			locYAtXMin /= 5.0;

			//use avg bin content within +/- 2 bins
			int locXMaxBin = locHist->FindBin(locFitRangeMax);
			double locYAtXMax = 0.0;
			int locMaxBinToSearch = ((locXMaxBin + 2) > locHist->GetNbinsX()) ? locHist->GetNbinsX() : locXMaxBin + 2;
			for(int loc_i = locXMaxBin - 2; loc_i <= locMaxBinToSearch; ++loc_i)
				locYAtXMax += locHist->GetBinContent(loc_i);
			locYAtXMax /= 5.0;

			if(locForceGreaterThanZeroFlag && (locYAtXMin <= 0.0))
				locYAtXMin = 0.01;
			if(locForceGreaterThanZeroFlag && (locYAtXMax <= 0.0))
				locYAtXMax = 0.01;

			locSlope = (locYAtXMax - locYAtXMin)/(locFitRangeMax - locFitRangeMin);
			locIntercept = locYAtXMax - locSlope*locFitRangeMax;
		}

		inline static TF1* Guess_Params_Flat_Fit(TH1* locHist, double locFitRangeMin, double locFitRangeMax, double& locIntercept)
		{
			if(!Check_ValidFitRange(locHist, locFitRangeMin, locFitRangeMax))
				return NULL;
			TF1* locFunc = new TF1("temp", "[0]", locFitRangeMin, locFitRangeMax);
			locHist->Fit(locFunc, "QRN+");
			locIntercept = locFunc->GetParameter(0);
			return locFunc;
		}

		inline static void Guess_Params_Flat(TH1* locHist, double locFitRangeMin, double locFitRangeMax, double& locIntercept)
		{
			if(!Check_ValidFitRange(locHist, locFitRangeMin, locFitRangeMax))
				return;
			double locYAtXMin = locHist->GetBinContent(locHist->FindBin(locFitRangeMin));
			double locYAtXMax = locHist->GetBinContent(locHist->FindBin(locFitRangeMax));
			locIntercept = 0.5*(locYAtXMax + locYAtXMin);
		}

		inline static TF1* Guess_Params_Exponential(TH1* locHist, double locFitRangeMin, double locFitRangeMax)
		{
			if(!Check_ValidFitRange(locHist, locFitRangeMin, locFitRangeMax))
				return NULL;
			TF1* locFunc = new TF1("temp", "[0]*exp([1]*x)", locFitRangeMin, locFitRangeMax);
			locHist->Fit(locFunc, "QRN+");
			return locFunc;
		}

		inline static void Calc_ParabolaParams(double locCurve_VerticalChange, double locCurve_HorizontalRangeOfVerticalChange, double locCriticalPointX, double locCriticalPointY, double& locFlatParam, double& locLinearParam, double& locQuadParam)
		{
			//locCurve_VerticalChange is CHANGE, STARTING AT THE CRITICAL POINT, AND GOING IN EITHER DIRECTION
			//locCurve_HorizontalRangeOfVerticalChange CAN BE ANY SIGN
			//y = a*x^2 + b*x + c
			//start at origin, a = 1, b = c = 0
			//set a: vertical v, horizontal h
			//y = a*x*x, v = a*h*h, a = v/(h*h)
			//move up (c)
			//move center to m, x -> (x - m)
			//y = a*x^2 + b*x + c, y = a*(x - m)^2 + b*(x - m) + c, y = a*x^2 - 2*a*m*x + a*m*m + b*x - b*m + c
			//y = (a)*x^2 + (b - 2*a*m)*x + (c + a*m*m - b*m)
			//y = a*x^2 + f*x + g, where f = b - 2*a*m, and g = c + a*m*m - b*m
			locQuadParam = locCurve_VerticalChange/(locCurve_HorizontalRangeOfVerticalChange*locCurve_HorizontalRangeOfVerticalChange);
			locLinearParam = -2.0*locQuadParam*locCriticalPointX;
			locFlatParam = locCriticalPointY + locQuadParam*locCriticalPointX*locCriticalPointX;
		}

		static double Calc_EntriesInRange(const TH1* locHist, double locFitRangeMin, double locFitRangeMax, double& locAverageBinWidth);

		static void Calc_CutRange_SignalFraction(TF1* locFunc, double locSignalFraction, double& locCutRangeMin, double& locCutRangeMax, double locMean, double locSigma);

		static void Fit_Slices_Gaussian(TH2 *loc2DHist, TH1*& locMeanHist, TH1*& locSigmaHist, TH1*& locChiSqHist);

	ClassDef(DFitUtilities, 1)
};

#endif