//Gaussian Fits
void Create_dGaussianEdgeGraphs(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locLowEdge, TGraphErrors*& locHighEdge, double locNumSigmas, unsigned int locNumBinsPerGroup, double locMinimumCounts, double locFitRangeFindHeightCutoff_Low = 0.1, double locFitRangeFindHeighMidSpot = -2.0, double locFitRangeFindHeightCutoff_High = 0.1);
void Create_dGaussianEdgeGraphs(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locLowEdge, TGraphErrors*& locHighEdge, TGraphErrors*& locMeans, TGraphErrors*& locSigmas, double locNumSigmas, unsigned int locNumBinsPerGroup, double locMinimumCounts, double locFitRangeFindHeightCutoff_Low = 0.1, double locFitRangeFindHeighMidSpot = -2.0, double locFitRangeFindHeightCutoff_High = 0.1);
void Fit_dGaussians(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locMeans, TGraphErrors*& locSigmas, unsigned int locNumBinsPerGroup, double locMinimumCounts, double locFitRangeFindHeightCutoff_Low = 0.1, double locFitRangeFindHeighMidSpot = -2.0, double locFitRangeFindHeightCutoff_High = 0.1);

//Means without Gaussian fits
void Calc_dMeanValue(TH1* loc1DHist, double locCutOffFraction, double& locMean, double& locUncertainty);
void Find_dBinEdges(TH1* loc1DHist, double locCutOffFraction, unsigned int& locMinBinToUse, unsigned int& locMaxBinToUse);
TGraphErrors* Calc_dYMeanValue(TH2* locInput2DHist, double locCutOffFraction, int locNumBinsPerGroup = 1, double locMinimumCounts = 0.0);

//Filter Hists
TH1* Filter_dHist(TH1* locHist);
TH2* Filter_dHist(TH2* locHist);

//Print Hist Contents to text file
void PrintToFile_dHist(TH1* locHist, string locOutputFileName, bool locIncludeErrorsFlag = true);

//MISCELLANEOUS
TObjArray* Get_Hists(TObjArray* locFileArray, string locDirPath, string locHistName);

void Create_dGaussianEdgeGraphs(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locLowEdge, TGraphErrors*& locHighEdge, double locNumSigmas, unsigned int locNumBinsPerGroup, double locMinimumCounts, double locFitRangeFindHeightCutoff_Low, double locFitRangeFindHeighMidSpot, double locFitRangeFindHeightCutoff_High)
{
	TGraphErrors* locMeans;
	TGraphErrors* locSigmas;
	Create_dGaussianEdgeGraphs(locInput2DHist, locHistogramSlices, locLowEdge, locHighEdge, locMeans, locSigmas, locNumSigmas, locNumBinsPerGroup, locMinimumCounts, locFitRangeFindHeightCutoff_Low, locFitRangeFindHeighMidSpot, locFitRangeFindHeightCutoff_High);
}

void Create_dGaussianEdgeGraphs(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locLowEdge, TGraphErrors*& locHighEdge, TGraphErrors*& locMeans, TGraphErrors*& locSigmas, double locNumSigmas, unsigned int locNumBinsPerGroup, double locMinimumCounts, double locFitRangeFindHeightCutoff_Low, double locFitRangeFindHeighMidSpot, double locFitRangeFindHeightCutoff_High)
{
	Fit_dGaussians(locInput2DHist, locHistogramSlices, locMeans, locSigmas, locNumBinsPerGroup, locMinimumCounts, locFitRangeFindHeightCutoff_Low, locFitRangeFindHeighMidSpot, locFitRangeFindHeightCutoff_High);

	unsigned int locNumPoints = locMeans->GetN();
	double* locLowValues = new double[locNumPoints];
	double* locHighValues = new double[locNumPoints];
	double* locUncertainties = new double[locNumPoints];
	for(unsigned int loc_i = 0; loc_i < locNumPoints; ++loc_i)
	{
		locLowValues[loc_i] = (locMeans->GetY())[loc_i] - locNumSigmas*((locSigmas->GetY())[loc_i]);
		locHighValues[loc_i] = (locMeans->GetY())[loc_i] + locNumSigmas*((locSigmas->GetY())[loc_i]);
		locUncertainties[loc_i] = sqrt((locMeans->GetEY())[loc_i]*(locMeans->GetEY())[loc_i] + locNumSigmas*locNumSigmas*(locSigmas->GetEY())[loc_i]*(locSigmas->GetEY())[loc_i]);
	}

	ostringstream locNumSigmaStream;
	locNumSigmaStream << locNumSigmas;

	locLowEdge = new TGraphErrors(locNumPoints, locMeans->GetX(), locLowValues, locMeans->GetEX(), locUncertainties);
	locLowEdge->SetName((string(locInput2DHist->GetName()) + string("_PeakLowEdge")).c_str());
	locLowEdge->SetTitle(locInput2DHist->GetTitle());
	locLowEdge->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locLowEdge->GetYaxis()->SetTitle((string("#mu - ") + locNumSigmaStream.str() + string("#sigma ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());

	locHighEdge = new TGraphErrors(locNumPoints, locMeans->GetX(), locHighValues, locMeans->GetEX(), locUncertainties);
	locHighEdge->SetName((string(locInput2DHist->GetName()) + string("_PeakHighEdge")).c_str());
	locHighEdge->SetTitle(locInput2DHist->GetTitle());
	locHighEdge->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locHighEdge->GetYaxis()->SetTitle((string("#mu + ") + locNumSigmaStream.str() + string("#sigma ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());
}

void Fit_dGaussians(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locMeans, TGraphErrors*& locSigmas, unsigned int locNumBinsPerGroup, double locMinimumCounts, double locFitRangeFindHeightCutoff_Low, double locFitRangeFindHeighMidSpot, double locFitRangeFindHeightCutoff_High)
{
	//assumes no background
	unsigned int locGraphArrayIndex = 0;

	unsigned int locMaxArraySize = locInput2DHist->GetNbinsX()/locNumBinsPerGroup + 1;
	double* locXArray = new double[locMaxArraySize];
	double* locXArrayUncertainty = new double[locMaxArraySize];
	double* locMeanArray = new double[locMaxArraySize];
	double* locMeanUncertainties = new double[locMaxArraySize];
	double* locSigmaArray = new double[locMaxArraySize];
	double* locSigmaUncertainties = new double[locMaxArraySize];

	DPlotFitter locPlotFitter;

	locHistogramSlices = new TObjArray();
	ostringstream locHistTitle;
	for(int loc_i = 1; loc_i <= locInput2DHist->GetNbinsX(); loc_i += locNumBinsPerGroup)
	{
		ostringstream locBinStream;
		locBinStream << loc_i;
		string locHistName = string(locInput2DHist->GetName()) + string("_XBin") + locBinStream.str();
		TH1* loc1DHist = locInput2DHist->ProjectionY(locHistName.c_str(), loc_i, loc_i + locNumBinsPerGroup - 1);

		if(loc1DHist->GetEntries() < locMinimumCounts)
			continue;
		locHistTitle.str("");
		locHistTitle << locInput2DHist->GetTitle() << ", " << locInput2DHist->GetXaxis()->GetBinLowEdge(loc_i) << " <= ";
		locHistTitle << locInput2DHist->GetXaxis()->GetTitle() << " < " << locInput2DHist->GetXaxis()->GetBinUpEdge(loc_i + locNumBinsPerGroup - 1);
		loc1DHist->SetTitle(locHistTitle.str().c_str());
		double locBinCenter = (locNumBinsPerGroup%2 == 1) ? locInput2DHist->GetXaxis()->GetBinCenter(loc_i + locNumBinsPerGroup/2) : locInput2DHist->GetXaxis()->GetBinLowEdge(loc_i + locNumBinsPerGroup/2);

		unsigned int locMinBinToUse, locMaxBinToUse;
		double locFitRangeFindHeightCutoff = (locBinCenter < locFitRangeFindHeighMidSpot) ? locFitRangeFindHeightCutoff_Low : locFitRangeFindHeightCutoff_High;
		Find_dBinEdges(loc1DHist, locFitRangeFindHeightCutoff, locMinBinToUse, locMaxBinToUse);

		double locFitRangeMin = loc1DHist->GetXaxis()->GetBinCenter(locMinBinToUse);
		double locFitRangeMax = loc1DHist->GetXaxis()->GetBinCenter(locMaxBinToUse);

		//setup fit params
		double locGaussMean = 0.5*(locFitRangeMin + locFitRangeMax);
		double locGaussWidth = (locGaussMean - locFitRangeMin)/2.5; //fit range to edge ~2.5 sigma
		int locPeakBin = loc1DHist->GetBinContent(loc1DHist->GetXaxis()->FindBin(locGaussMean));
		double locGaussHeight = 0.9*loc1DHist->GetBinContent(locPeakBin);

		//setup fit functors
		vector<double> locInitParams(3, 0.0);
		locInitParams[0] = locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = locGaussWidth;
		DFitFunctor_Gaussian* locFitFunctor_Gaussian = new DFitFunctor_Gaussian(locInitParams, kBlack);

		//fit gauss
		DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_Gaussian);
		locPlotFitter.Fit_Hist(loc1DHist, locFitControl);
		locHistogramSlices->AddAtAndExpand(loc1DHist, locGraphArrayIndex);

		//get fit params
		locXArray[locGraphArrayIndex] = locBinCenter;
		locXArrayUncertainty[locGraphArrayIndex] = 0.0;
		locMeanArray[locGraphArrayIndex] = locFitFunctor_Gaussian->dFunc->GetParameter(1);
		locMeanUncertainties[locGraphArrayIndex] = locFitFunctor_Gaussian->dFunc->GetParError(1);
		locSigmaArray[locGraphArrayIndex] = locFitFunctor_Gaussian->dFunc->GetParameter(2);
		locSigmaUncertainties[locGraphArrayIndex] = locFitFunctor_Gaussian->dFunc->GetParError(2);
		++locGraphArrayIndex;

/*
		double locMeanShift = 0.0;
		if((loc1DHist->GetMean() - locGaussMean) > 0.005)
			locMeanShift = 0.03;
		else if((locGaussMean - loc1DHist->GetMean()) > 0.005)
			locMeanShift = -0.03;

		//setup fit functors
		vector<DFitFunctor*> locSignalFitFunctors;
		vector<double> locInitParams(3, 0.0);
		locInitParams[0] = locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = locGaussWidth;
		DFitFunctor_Gaussian* locFitFunctor_Gaussian = new DFitFunctor_Gaussian(locInitParams, kBlue, true);
		locSignalFitFunctors.push_back(locFitFunctor_Gaussian);
		locInitParams[0] = 0.05*locGaussHeight;  locInitParams[1] = locGaussMean + locMeanShift;  locInitParams[2] = 2.0*locGaussWidth;
		DFitFunctor_Gaussian* locFitFunctor_Gaussian2 = new DFitFunctor_Gaussian(locInitParams, kViolet, true);
		locSignalFitFunctors.push_back(locFitFunctor_Gaussian2);

		//linear background
		vector<DFitFunctor*> locBackgroundFitFunctors;
		locInitParams.resize(2);
		locInitParams[0] = 0.01*locGaussHeight;  locInitParams[1] = 0.0;//  locInitParams[2] = 0.0;
		DFitFunctor_Polynomial* locFitFunctor_Polynomial = new DFitFunctor_Polynomial(locInitParams, kRed, true);
		locBackgroundFitFunctors.push_back(locFitFunctor_Polynomial);
*/
/*
		//fit gauss
		locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locSignalFitFunctors, locBackgroundFitFunctors, true);
		locPlotFitter.Fit_Hist(loc1DHist, locFitControl);
		locHistogramSlices->AddAtAndExpand(loc1DHist, locGraphArrayIndex);

		//find which was narrower
		DFitFunctor_Gaussian* locFitFunctor_NarrowerGaussian = locFitFunctor_Gaussian;
		if(locFitFunctor_Gaussian2->dFunc->GetParameter(2) < locFitFunctor_Gaussian->dFunc->GetParameter(2))
			locFitFunctor_NarrowerGaussian = locFitFunctor_Gaussian2;

		//get fit params
		locXArray[locGraphArrayIndex] = locBinCenter;
		locXArrayUncertainty[locGraphArrayIndex] = 0.0;
		locMeanArray[locGraphArrayIndex] = locFitFunctor_NarrowerGaussian->dFunc->GetParameter(1);
		locMeanUncertainties[locGraphArrayIndex] = locFitFunctor_NarrowerGaussian->dFunc->GetParError(1);
		locSigmaArray[locGraphArrayIndex] = locFitFunctor_NarrowerGaussian->dFunc->GetParameter(2);
		locSigmaUncertainties[locGraphArrayIndex] = locFitFunctor_NarrowerGaussian->dFunc->GetParError(2);
		++locGraphArrayIndex;
*/
	}

	locMeans = new TGraphErrors(locGraphArrayIndex, locXArray, locMeanArray, locXArrayUncertainty, locMeanUncertainties);
	locMeans->SetName((string(locInput2DHist->GetName()) + string("_Means")).c_str());
	locMeans->SetTitle(locInput2DHist->GetTitle());
	locMeans->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locMeans->GetYaxis()->SetTitle((string("Mean ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());

	locSigmas = new TGraphErrors(locGraphArrayIndex, locXArray, locSigmaArray, locXArrayUncertainty, locSigmaUncertainties);
	locSigmas->SetName((string(locInput2DHist->GetName()) + string("_Sigmas")).c_str());
	locSigmas->SetTitle(locInput2DHist->GetTitle());
	locSigmas->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locSigmas->GetYaxis()->SetTitle((string("#sigma ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());
}


void Calc_dMeanValue(TH1* loc1DHist, double locCutOffFraction, double& locMean, double& locUncertainty)
{
	//find the weighted mean & uncertainty
	unsigned int locMinBinToUse, locMaxBinToUse;
	Find_dBinEdges(loc1DHist, locCutOffFraction, locMinBinToUse, locMaxBinToUse);
	loc1DHist->GetXaxis()->SetRange(locMinBinToUse, locMaxBinToUse);
	locMean = loc1DHist->GetMean();
	locUncertainty = loc1DHist->GetMeanError();
}

void Find_dBinEdges(TH1* loc1DHist, double locCutOffFraction, unsigned int& locMinBinToUse, unsigned int& locMaxBinToUse)
{
	//locCutOffFraction: ignore bins with less counts than the max * this
	int locMaxBin = loc1DHist->GetMaximumBin();
	double locMaxValue = loc1DHist->GetBinContent(locMaxBin);

	//find min/max bins to use
	int locNumBins = loc1DHist->GetNbinsX();
	locMinBinToUse = locMaxBin;
	//min
	for(int loc_i = locMaxBin - 1; loc_i > 0; --loc_i)
	{
		if(loc1DHist->GetBinContent(loc_i) < (locMaxValue*locCutOffFraction))
			break;
		locMinBinToUse = loc_i;
	}
	//max
	locMaxBinToUse = locMaxBin;
	for(int loc_i = locMaxBin + 1; loc_i <= locNumBins; ++loc_i)
	{
		if(loc1DHist->GetBinContent(loc_i) < (locMaxValue*locCutOffFraction))
			break;
		locMaxBinToUse = loc_i;
	}
}

TGraphErrors* Calc_dYMeanValue(TH2* locInput2DHist, double locCutOffFraction, int locNumBinsPerGroup, double locMinimumCounts)
{
	//locCutOffFraction: ignore bins with less counts than the max * this

	unsigned int locNumBins = locInput2DHist->GetNbinsX();
	double* locXArray = new double[locNumBins];
	double* locYArray = new double[locNumBins];
	double* locXArrayUncertainty = new double[locNumBins];
	double* locYArrayUncertainty = new double[locNumBins];

	unsigned int locGraphArrayIndex = 0;
	for(int loc_i = 1; loc_i <= locNumBins; loc_i += locNumBinsPerGroup)
	{
		ostringstream locBinStream;
		locBinStream << loc_i;
		string locHistName = string(locInput2DHist->GetName()) + string("_XBin") + locBinStream.str();
		TH1D* loc1DHist = locInput2DHist->ProjectionY(locHistName.c_str(), loc_i, loc_i + locNumBinsPerGroup - 1);

		if(loc1DHist->GetEntries() < locMinimumCounts)
			continue;
		double locMean, locUncertainty;
		Calc_dMeanValue(loc1DHist, locCutOffFraction, locMean, locUncertainty);
		if(!(fabs(locMean) > 0.0))
			continue;

		locXArrayUncertainty[locGraphArrayIndex] = 0.0;
		double locBinCenter = (locNumBinsPerGroup%2 == 1) ? locInput2DHist->GetXaxis()->GetBinCenter(loc_i + locNumBinsPerGroup/2) : locInput2DHist->GetXaxis()->GetBinLowEdge(loc_i + locNumBinsPerGroup/2);
		locXArray[locGraphArrayIndex] = locBinCenter;
		locYArray[locGraphArrayIndex] = locMean;
		locYArrayUncertainty[locGraphArrayIndex] = locUncertainty;
		++locGraphArrayIndex;
	}

	TGraphErrors* locGraphErrors = new TGraphErrors(locGraphArrayIndex, locXArray, locYArray, locXArrayUncertainty, locYArrayUncertainty);
	locGraphErrors->SetName((string(locInput2DHist->GetName()) + string("_Means")).c_str());
	locGraphErrors->SetTitle(locInput2DHist->GetTitle());
	locGraphErrors->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locGraphErrors->GetYaxis()->SetTitle((string("Mean ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());

	return locGraphErrors;
}

TH1* Filter_dHist(TH1* locHist)
{
	double locMaxDiffFraction = 1.3;
	vector<Int_t> locBinsToClear;
	for(Int_t locBin = 2; locBin < locHist->GetNbinsX(); ++locBin)
	{
		double locLowBinContent = locHist->GetBinContent(locBin - 1);
		double locHighBinContent = locHist->GetBinContent(locBin + 1);
		if((!(locLowBinContent > 0.0)) || (!(locHighBinContent > 0.0)))
		{
			//zero all around it, clear it
			locBinsToClear.push_back(locBin);
			continue;
		}
		double locAdjacentAvg = 0.5*(locLowBinContent + locHighBinContent);
		double locFractionDifference = locHist->GetBinContent(locBin)/locAdjacentAvg;
		if(locFractionDifference > locMaxDiffFraction)
			locBinsToClear.push_back(locBin);
	}
	for(size_t loc_i = 0; loc_i < locBinsToClear.size(); ++loc_i)
	{
		locHist->SetBinContent(locBinsToClear[loc_i], 0.0);
		locHist->SetBinError(locBinsToClear[loc_i], 0.0);
	}
	return locHist;
}

TH2* Filter_dHist(TH2* locHist)
{
	double locMaxDiffFraction = 1.3;
	vector<pair<Int_t, Int_t> > locBinsToClear;
	for(Int_t locXBin = 2; locXBin < locHist->GetNbinsX(); ++locXBin)
	{
		for(Int_t locYBin = 2; locYBin < locHist->GetNbinsY(); ++locYBin)
		{
			double locLowBinContent = locHist->GetBinContent(locXBin, locYBin - 1);
			double locHighBinContent = locHist->GetBinContent(locXBin, locYBin + 1);

			double locLeftBinContent = locHist->GetBinContent(locXBin - 1, locYBin);
			double locRightBinContent = locHist->GetBinContent(locXBin + 1, locYBin);

			if((!(locLowBinContent > 0.0)) && (!(locHighBinContent > 0.0)) && (!(locLowBinContent > 0.0)) && (!(locHighBinContent > 0.0)))
			{
				//zero all around it, clear it
				locBinsToClear.push_back(pair<Int_t, Int_t>(locXBin, locYBin));
				continue;
			}
			if((!(locLowBinContent > 0.0)) || (!(locHighBinContent> 0.0)))
				continue;
			if((!(locLeftBinContent > 0.0)) || (!(locRightBinContent> 0.0)))
				continue;

			double locAdjacentAvg = 0.25*(locLowBinContent + locHighBinContent + locLeftBinContent + locRightBinContent);
			double locFractionDifference = locHist->GetBinContent(locXBin, locYBin)/locAdjacentAvg;
			if(locFractionDifference > locMaxDiffFraction)
				locBinsToClear.push_back(pair<Int_t, Int_t>(locXBin, locYBin));
		}
	}
	for(size_t loc_i = 0; loc_i < locBinsToClear.size(); ++loc_i)
	{
		locHist->SetBinContent(locBinsToClear[loc_i].first, locBinsToClear[loc_i].second, 0.0);
		locHist->SetBinError(locBinsToClear[loc_i].first, locBinsToClear[loc_i].second, 0.0);
	}
	return locHist;
}

void PrintToFile_dHist(TH1* locHist, string locOutputFileName, bool locIncludeErrorsFlag)
{
	ofstream locOutputFileStream;
	locOutputFileStream.open(locOutputFileName.c_str());

	//header
	locOutputFileStream << "# " << locHist->GetNbinsX() << " Histogram Bins, X-min = " << locHist->GetXaxis()->GetXmin() << ", X-max = " << locHist->GetXaxis()->GetXmax() << endl;
	locOutputFileStream << "#Format: Bin-index (starting from 1, up THROUGH N-bins), Bin-content";
	if(locIncludeErrorsFlag)
		locOutputFileStream << ", Bin-content-error";
	locOutputFileStream << endl;

	for(int locBin = 1; locBin <= locHist->GetNbinsX(); ++locBin)
		locOutputFileStream << locBin << " " << locHist->GetBinContent(locBin) << " " << locHist->GetBinError(locBin) << endl;

	locOutputFileStream.close();
}

TObjArray* Get_Hists(TObjArray* locFileArray, string locDirPath, string locHistName)
{
	TObjArray* locHistArray = new TObjArray(locFileArray->GetEntriesFast());
	for(int loc_i = 0; loc_i < locFileArray->GetEntriesFast(); ++loc_i)
	{
		TFile* locFile = (TFile*)locFileArray->At(loc_i);
		locFile->cd();
		gDirectory->cd(locDirPath.c_str());
		TObject* locHist = gDirectory->Get(locHistName.c_str());
		locHistArray->AddAt(locHist, loc_i);
	}
	return locHistArray;
}