#include <iostream>
#include <sstream>
#include <string>
#include <fstream>
#include <map>
#include <set>
#include <string>
#include <vector>

#include <TFile.h>
#include <TH2.h>
#include <TH1.h>
#include <TGraph2D.h>
#include <TGraphErrors.h>

#include "PlotFitter/DPlotFitter.h"
#include "PlotFitter/DFitControl.h"
#include "PlotFitter/DFitUtilities.h"
#include "PlotFitter/DFitFunctors_Miscellaneous.h"
#include "PlotFitter/DFitFunctors_Gaussians.h"

#include "PlotDrawer/DPlotDrawer.h"
#include "PlotDrawer/DCanvasInstructions.h"

#include "particleType.h"

//User-callable functions
void Fit_Resolutions(string locInputFileName, bool locFirstPassFlag);
void Smooth_Resolutions(string locInputFileName);
void Plot_Resoultions(string locInputFileName);
void Draw_Resolutions(string locInputFileName);

//Fit resolutions
void Fit_Resolution(TH2* loc2DHist, double locYSignalEdge, int locBaseYRebinFactor, int locMinYRebin);
TH1* Get_HistProjection(TH2* loc2DHist, int& locStartXBin, double locYSignalEdge, int locBaseYRebinFactor, int locMinYRebin, double& locAverageX);
DFitControl* Fit_Histogram(TH1* locHist, double locYSignalEdge, double locAverageX);
void Fit_ToLinears(TGraphErrors* locGraph);
void Cleanup_Graph(TGraphErrors* locGraph, int locNumFitPoints);

//fit control
void Find_InitGaussParams(TH1* locHist, double locWidestYSignalEdge, double& locGaussHeight, double& locGaussMean, double& locGaussWidth, double& locApproximateIntegral);
DFitControl* Attempt_DoubleGaussFit(TH1* locHist, double locGaussHeight, double locGaussMean, double locInitGaussWidth, double& locFitRangeMin, double& locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);
DFitControl* Attempt_GaussFit(TH1* locHist, double locGaussHeight, double locGaussMean, double locInitGaussWidth, double& locFitRangeMin, double& locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);
bool Eval_Fit(TH1* locHist, DFitControl* locFitControl);
bool Eval_Fit_Gaussian(TH1* locHist, DFitControl* locFitControl);

//perform fits
DFitControl* Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);
DFitControl* Fit_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);

//create fit functors
DFitFunctor* Build_PolyBackground_Flat(TH1* locHist, double locFitRangeMin, double locFitRangeMax);
DFitFunctor* Build_PolyBackground_Linear(TH1* locHist, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax);
DFitFunctor_DoubleGaussian* Build_DoubleGaussian(double locGaussHeight, double locGaussMean, double locGaussWidth);
DFitFunctor_Gaussian* Build_Gaussian(double locGaussHeight, double locGaussMean, double locGaussWidth);

//utilities
double* Convert_VectorToArray(const vector<double>& locVector);
double Calc_InitGaussHeight(TH1* locHist, double locGaussMean, double locGaussWidth, double locNumSigmas);

bool gDebugFlag = false;

/*********************************************************** FIT RESOLUTIONS ***********************************************************/

/*
 From ~/Scratch/gluex/
.x ../Load.C
.L ../scripts_trackeff/Fit_Resolutions.C+
//DON'T FORGET THE BOOL ARGUMENT!!
Fit_Resolutions("trackeff_Proton_2pi.root");
Fit_Resolutions("trackeff_Pi-_2pi.root");
Fit_Resolutions("trackeff_Pi+_2pi.root");

Fit_Resolutions("trackeff_Pi-_3pi.root");
Fit_Resolutions("trackeff_Pi+_3pi.root");

Fit_Resolutions("trackeff_Proton_4pi.root");
Fit_Resolutions("trackeff_Pi-_4pi.root");
Fit_Resolutions("trackeff_Pi+_4pi.root");
*/
void Fit_Resolutions(string locInputFileName, bool locFirstPassFlag)
{
	//1st pass: delta-phi & delta-p
	//2nd pass: delta-theta //background is too high without above cuts
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	string locOutputFileName = locInputFileName + ".fit.root";
	TFile* locOutputFile = new TFile(locOutputFileName.c_str(), "RECREATE");
	gDirectory->mkdir("Resolution", "Resolution");
	gDirectory->cd("Resolution");

	vector<string> locVsTypes = {"VsP", "VsTheta", "VsPhi"};

	//setup fit/binning controls
	map<string, int> locYRebinMap;
	locYRebinMap["DeltaPOverP"] = 5;
	locYRebinMap["DeltaTheta"] = 5;
	locYRebinMap["DeltaPhi"] = 8;
	map<string, int> locMinYRebinMap;
	locMinYRebinMap["DeltaPOverP"] = 2;
	locMinYRebinMap["DeltaTheta"] = 4;
	locMinYRebinMap["DeltaPhi"] = 4;
	map<string, double> locYSignalEdgeMap; //should be the WIDEST that it gets to be
	locYSignalEdgeMap["DeltaPOverP"] = 0.2;
	locYSignalEdgeMap["DeltaTheta"] = 20.0;
	locYSignalEdgeMap["DeltaPhi"] = 20.0;

	//loop over energy bins
	int locEBin = 0;
	while(true)
	{
		cout << "ebin = " << locEBin << endl;
		//make energy bin dir name
		ostringstream locEBinDir;
		locEBinDir << "EnergyBin_" << locEBin;
		string locDirName = string("Resolution/") + locEBinDir.str();

		//see if dir exists
		TObject* locObject = locInputFile->Get(locDirName.c_str());
		if(locObject == NULL)
			break; //at the end

		//make output dir
		gDirectory->mkdir(locEBinDir.str().c_str(), locEBinDir.str().c_str());
		gDirectory->cd(locEBinDir.str().c_str());

		//loop over delta-types
		for(auto locYSignalEdgePair : locYSignalEdgeMap)
		{
			if(locFirstPassFlag == (locYSignalEdgePair.first == "DeltaTheta"))
				continue; //if 1st pass, not delta-theta. if 2nd pass, only delta-theta

			double locYSignalEdge = locYSignalEdgePair.second;
			int locBaseYRebinFactor = locYRebinMap[locYSignalEdgePair.first];
			int locMinYRebin = locMinYRebinMap[locYSignalEdgePair.first];

			//loop over vs-types
			for(auto locVsString : locVsTypes)
			{
				//get histogram & fit
				string locHistPath = locDirName + string("/") + locYSignalEdgePair.first + locVsString;
				//cout << "locHistPath = " << locHistPath << endl;
				TH2* loc2DHist = (TH2*)locInputFile->Get(locHistPath.c_str());
				Fit_Resolution(loc2DHist, locYSignalEdge, locBaseYRebinFactor, locMinYRebin);
			}
		}

		//get ready for next energy bin
		gDirectory->cd("..");
		++locEBin;
	}

	//save results
	locOutputFile->Write();
	locOutputFile->Close();
}

void Fit_Resolution(TH2* loc2DHist, double locYSignalEdge, int locBaseYRebinFactor, int locMinYRebin)
{
	int locStartXBin = 1;
	vector<double> locXs, locMeans, locMeanUncertainties, locSigmas, locSigmaUncertainties;

	//loop over x bins, making projections & fitting
	while(true)
	{
		//get next projection
		double locAverageX = 0.0;
		TH1* locProjHist = Get_HistProjection(loc2DHist, locStartXBin, locYSignalEdge, locBaseYRebinFactor, locMinYRebin, locAverageX);
		if(locProjHist == nullptr)
			break; //none left to fit

		//fit it
		DFitControl* locFitControl = Fit_Histogram(locProjHist, locYSignalEdge, locAverageX);
		if(locFitControl == nullptr)
			continue;

		//save fit results
		locXs.push_back(locAverageX);

		//get fit functions
		DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = dynamic_cast<DFitFunctor_DoubleGaussian*>(locFitControl->dSignalFitGroup->dFitFunctors[0]);
		if(locFitFunctor_DoubleGaussian != nullptr)
		{
			locMeans.push_back(locFitFunctor_DoubleGaussian->dFunc->GetParameter(1));
			locMeanUncertainties.push_back(locFitFunctor_DoubleGaussian->dFunc->GetParError(1));
			pair<double, double> locSigmaPair = locFitFunctor_DoubleGaussian->Calc_StandardDeviation();
			locSigmas.push_back(locSigmaPair.first);
			locSigmaUncertainties.push_back(locSigmaPair.second);
		}
		else
		{
			DFitFunctor_Gaussian* locFitFunctor_Gaussian = dynamic_cast<DFitFunctor_Gaussian*>(locFitControl->dSignalFitGroup->dFitFunctors[0]);
			locMeans.push_back(locFitFunctor_Gaussian->dFunc->GetParameter(1));
			locMeanUncertainties.push_back(locFitFunctor_Gaussian->dFunc->GetParError(1));
			locSigmas.push_back(locFitFunctor_Gaussian->dFunc->GetParameter(2));
			locSigmaUncertainties.push_back(locFitFunctor_Gaussian->dFunc->GetParError(2));
		}
	}

	//convert vectors to arrays
	unsigned int locArraySize = locMeans.size();
	double* locXArray = Convert_VectorToArray(locXs);
	double* locMeanArray = Convert_VectorToArray(locMeans);
	double* locMeanUncertaintyArray = Convert_VectorToArray(locMeanUncertainties);
	double* locSigmaArray = Convert_VectorToArray(locSigmas);
	double* locSigmaUncertaintyArray = Convert_VectorToArray(locSigmaUncertainties);

	//create graph of means
	TGraphErrors* locMeanGraph = new TGraphErrors(locMeans.size(), locXArray, locMeanArray, NULL, locMeanUncertaintyArray);
	locMeanGraph->SetName((string(loc2DHist->GetName()) + string("_Means")).c_str());
	locMeanGraph->SetTitle(loc2DHist->GetTitle());
	locMeanGraph->GetXaxis()->SetTitle(loc2DHist->GetXaxis()->GetTitle());
	locMeanGraph->GetYaxis()->SetTitle((string("Mean ") + string(loc2DHist->GetYaxis()->GetTitle())).c_str());
	locMeanGraph->Write();

	//create graph of sigmas
	TGraphErrors* locSigmaGraph = new TGraphErrors(locMeans.size(), locXArray, locSigmaArray, NULL, locSigmaUncertaintyArray);
	locSigmaGraph->SetName((string(loc2DHist->GetName()) + string("_Sigmas")).c_str());
	locSigmaGraph->SetTitle(loc2DHist->GetTitle());
	locSigmaGraph->GetXaxis()->SetTitle(loc2DHist->GetXaxis()->GetTitle());
	locSigmaGraph->GetYaxis()->SetTitle((string("#sigma ") + string(loc2DHist->GetYaxis()->GetTitle())).c_str());
	locSigmaGraph->Write();
}

TH1* Get_HistProjection(TH2* loc2DHist, int& locStartXBin, double locYSignalEdge, int locBaseYRebinFactor, int locMinYRebin, double& locAverageX)
{
	//y-range over which to evaluate integral (the signal range)
	int locYRangeLowBin = loc2DHist->GetYaxis()->FindBin(-1.0*locYSignalEdge);
	int locYRangeHighBin = loc2DHist->GetYaxis()->FindBin(locYSignalEdge);

	//initialize loop parameters
	int locInitStartXBin = locStartXBin;
	int locEndXRangeBin = locStartXBin - 1; //will increment at beginning of loop
	double locTotalSignalIntegral = 0.0;

	//target signal integral
	string loc2DHistName = loc2DHist->GetName();
	double locTargetSignalIntegral = 4000.0;
	if(loc2DHistName.find("DeltaPVs") != string::npos)
		locTargetSignalIntegral = 4000.0;
	else if(loc2DHistName.find("VsTheta") != string::npos)
		locTargetSignalIntegral = 4000.0;
	else if(loc2DHistName.find("VsPhi") != string::npos)
		locTargetSignalIntegral = 6000.0;
	else //vs p
		locTargetSignalIntegral = 4000.0;

	//find x-range of bins to use
	//cout << "get hist projection, start bin = " << locStartXBin << endl;
	vector<double> locSignalIntegrals;
	while(locTotalSignalIntegral < locTargetSignalIntegral)
	{
		//increment bin
		++locEndXRangeBin;
		if(locEndXRangeBin > loc2DHist->GetNbinsX())
			return nullptr; //not enough statistics to fit

		//get bin integral
		double locBinIntegral = loc2DHist->Integral(locEndXRangeBin, locEndXRangeBin, locYRangeLowBin, locYRangeHighBin);

		//compute background integral
		int locBackgroundSearchNumBins = 4;
		double locAverageBackground = loc2DHist->Integral(locEndXRangeBin, locEndXRangeBin, locYRangeLowBin - 1 - locBackgroundSearchNumBins, locYRangeLowBin - 1)/double(locBackgroundSearchNumBins);
		locAverageBackground += loc2DHist->Integral(locEndXRangeBin, locEndXRangeBin, locYRangeHighBin + 1, locYRangeHighBin + 1 + locBackgroundSearchNumBins)/double(locBackgroundSearchNumBins);
		locAverageBackground /= 2.0;
		double locBackgroundIntegral = locAverageBackground*(locYRangeHighBin - locYRangeLowBin + 1);

		//increase total signal integral
		double locSignalIntegral = locBinIntegral - locBackgroundIntegral;
		if((locInitStartXBin == 1) && locSignalIntegrals.empty() && (locSignalIntegral < 0.05*locTargetSignalIntegral))
		{
			//at the beginning of the search for this histogram, signal fraction is so small, just ignore the bin
			++locStartXBin;
			continue;
		}

		//save the result
		locSignalIntegrals.push_back(locSignalIntegral);
		locTotalSignalIntegral += locSignalIntegral;
	}

	//compute average x:
	locAverageX = 0.0;
	for(size_t loc_i = 0; loc_i < locSignalIntegrals.size(); ++loc_i)
	{
		double locX = loc2DHist->GetXaxis()->GetBinCenter(locStartXBin + loc_i);
		locAverageX += locX*locSignalIntegrals[loc_i];
	}
	locAverageX /= locTotalSignalIntegral;

	//Make histogram projection
	ostringstream locHistName;
	locHistName << loc2DHist->GetName() << "_" << locStartXBin << "_" << locEndXRangeBin;
	TH1* locProjHist = (TH1*)loc2DHist->ProjectionY(locHistName.str().c_str(), locStartXBin, locEndXRangeBin);

	//set title
	ostringstream locHistTitle;
	locHistTitle << loc2DHist->GetTitle() << ", " << loc2DHist->GetXaxis()->GetBinLowEdge(locStartXBin) << " #leq ";
	locHistTitle << loc2DHist->GetXaxis()->GetTitle() << " < " << loc2DHist->GetXaxis()->GetBinUpEdge(locEndXRangeBin);
	locProjHist->SetTitle(locHistTitle.str().c_str());

	//Rebin histogram as needed
	int locRebinTerm = locBaseYRebinFactor;
	locRebinTerm /= (locTotalSignalIntegral/locTargetSignalIntegral);
	if(locRebinTerm < locMinYRebin)
		locRebinTerm = locMinYRebin;
	while((locProjHist->GetNbinsX() % locRebinTerm) != 0)
		++locRebinTerm; //make sure it rebins evenly
	if(locRebinTerm > 1)
		locProjHist->Rebin(locRebinTerm);

	//zero bins with negative content
	for(int locBin = 1; locBin <= locProjHist->GetNbinsX(); ++locBin)
	{
		if(locProjHist->GetBinContent(locBin) < 0.0)
			locProjHist->SetBinContent(locBin, 0.0);
	}

	locStartXBin = locEndXRangeBin + 1; //for next time
	return locProjHist;
}

void Find_InitGaussParams(TH1* locHist, double locWidestYSignalEdge, double& locGaussHeight, double& locGaussMean, double& locGaussWidth, double& locApproximateIntegral)
{
	string locHistName = locHist->GetName();

	//get gauss mean
	int locWidestYHighBin = locHist->FindBin(locWidestYSignalEdge);
	int locMeanSearchLowBin = locHist->FindBin(-0.2*locWidestYSignalEdge);
	int locMeanSearchHighBin = locHist->FindBin(0.2*locWidestYSignalEdge);
	locHist->GetXaxis()->SetRange(locMeanSearchLowBin, locMeanSearchHighBin);
	int locMaxBin = locHist->GetMaximumBin();
	locGaussMean = locHist->GetBinCenter(locMaxBin);
	locHist->GetXaxis()->SetRange(1, locHist->GetNbinsX()); //restore

	//get width: the input edge is assumed to be the WIDEST it can be: see if it is narrower
	//step to the side until hist height becomes relatively flat, or increases: fit linear slopes to rolling range of 4 bins
	int locFitRangeNumBins = 7;
	//check the high side
	int locFitRangeMinBin = locHist->FindBin(locGaussMean) - 1; //-1: will increase below
	double locInitialSlope = 1.0;
	while((locFitRangeMinBin + locFitRangeNumBins) < locWidestYHighBin)
	{
		++locFitRangeMinBin;
		double locFitRangeMin = locHist->GetBinLowEdge(locFitRangeMinBin);
		double locFitRangeMax = locHist->GetBinLowEdge(locFitRangeMinBin + locFitRangeNumBins);
		double locIntercept = 0.0, locSlope = 0.0;
		TF1* locLinearFunc = DFitUtilities::Guess_Params_Linear_Fit(locHist, locFitRangeMin, locFitRangeMax, locIntercept, locSlope);
		if(locLinearFunc == nullptr)
			continue;
		delete locLinearFunc;
		if(locInitialSlope > 0.0) //not set yet: set it now
			locInitialSlope = locSlope;
		if(fabs(locSlope) < fabs(0.05*locInitialSlope))
			break; //slope is significantly reduced: we are now at the edge
		if(locHistName.find("DeltaTheta") != string::npos) //also does vs phi
		{
			if(fabs(locSlope) < fabs(0.2*locInitialSlope))
				break; //slope is significantly reduced: we are now at the edge
		}
	}
	locGaussWidth = (locHist->GetBinLowEdge(locFitRangeMinBin + locFitRangeNumBins) - locGaussMean)/3.0;

//	gDebugFlag = (locHistName == "DeltaThetaVsTheta_5_5");
	if(gDebugFlag)
		cout << "init width = " << locGaussWidth << endl;

	//get init height guess
	locGaussHeight = Calc_InitGaussHeight(locHist, locGaussMean, locGaussWidth, 5.0);
	if(locHistName.find("DeltaTheta") != string::npos)
		return; //integral is ignored for these

	//initial fit (gauss + linear): get second guess for width
	double locInitFitRangeMin = locGaussMean - 5.0*locGaussWidth;
	double locInitFitRangeMax = locGaussMean + 5.0*locGaussWidth;
	double locCutRangeMin, locCutRangeMax;
	DFitControl* locFitControl = Fit_Gaussian(locHist, locGaussHeight, locGaussMean, locGaussWidth, locInitFitRangeMin, locInitFitRangeMax, locCutRangeMin, locCutRangeMax);
	if(locFitControl != nullptr)
	{
		DFitFunctor_Gaussian* locFitFunctor_Gaussian = (DFitFunctor_Gaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];
		locGaussWidth = locFitFunctor_Gaussian->dFunc->GetParameter(2);
		locApproximateIntegral = locFitFunctor_Gaussian->dFunc->GetParameter(0);

		//reset the fit
		locFitControl->Remove_FitFunctions(locHist);
		delete locFitControl;
	}
	else
		locApproximateIntegral = locGaussHeight*locGaussWidth*sqrt(2.0*TMath::Pi());
}

DFitControl* Fit_Histogram(TH1* locHist, double locWidestYSignalEdge, double locAverageX)
{
	string locHistName = locHist->GetName();
	double locGaussHeight = 0.0, locGaussMean = 0.0, locGaussWidth = 0.0, locApproximateIntegral = 0.0;
	Find_InitGaussParams(locHist, locWidestYSignalEdge, locGaussHeight, locGaussMean, locGaussWidth, locApproximateIntegral);


	/*
	//delta theta vs theta
	if(locHistName.find("DeltaThetaVsTheta") != string::npos)
	{
		//only 1 fit: entire fit range, except, at low-theta, truncate fit range (asymmetric due to boundary)
		double locFitRangeMin = locHist->GetXaxis()->GetBinLowEdge(0);
		double locFitRangeMax = locHist->GetXaxis()->GetBinUpEdge(locHist->GetNbinsX());

		//adjust fit range min
		int locFirstBinWithData = 1;
		for(; locFirstBinWithData <= locHist->GetNbinsX(); ++locFirstBinWithData)
		{
			if(locHist->GetBinContent(locFirstBinWithData) > 0.0)
				break;
		}
		double locLowEdge = locHist->GetXaxis()->GetBinLowEdge(locFirstBinWithData) + 5.0; //+ 5 degrees
		if(locLowEdge > locFitRangeMin)
			locFitRangeMin = locLowEdge;
		if(locFitRangeMin > 0.0)
			locFitRangeMin = 0.0;

		//fit
		double locCutRangeMin, locCutRangeMax;
		return Attempt_DoubleGaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	}
*/
	//delta phi
	if(locHistName.find("DeltaPhi") != string::npos)
	{
		//only 1 fit: entire fit range
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locHist->GetXaxis()->GetBinLowEdge(0);
		double locFitRangeMax = locHist->GetXaxis()->GetBinUpEdge(locHist->GetNbinsX());
		return Attempt_DoubleGaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	}

	if((locHistName.find("DeltaThetaVsTheta") != string::npos) && (locAverageX < 5.0))
	{
		//only 1 fit: very short range
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locGaussMean;
		double locFitRangeMax = locGaussMean + 0.5*locAverageX;
		locGaussWidth = (locFitRangeMax - locFitRangeMin)/3.0;
		return Attempt_GaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	}

	if((locHistName.find("DeltaThetaVsTheta") != string::npos) && (locAverageX < 20.0))
	{
		//only 1 fit: very short range
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locGaussMean - 3.0*locGaussWidth;
		double locFitRangeMax = locGaussMean + 3.0*locGaussWidth;
		return Attempt_GaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	}

	//delta theta vs p & phi
	if((locHistName.find("DeltaThetaVsP") != string::npos) || (locHistName.find("DeltaThetaVsPhi") != string::npos))
	{
		//only 1 fit: very short range
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locGaussMean - 3.0*locGaussWidth;
		double locFitRangeMax = locGaussMean + 3.0*locGaussWidth;
		return Attempt_GaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	}

	/*	else if((locHistName.find("DeltaThetaVsP") != string::npos) || (locHistName.find("DeltaThetaVsPhi") != string::npos))
	{
		//distributions are horribly distorted: only 1 fit: +/- 4 sigma
		//find the first x-bin with data
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locGaussMean - 4.0*locGaussWidth;
		double locFitRangeMax = locGaussMean + 4.0*locGaussWidth;
		return Attempt_DoubleGaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	}
*/
	//delta-theta vs theta
		//low angles: would benefit from
	//delta-theta vs p & phi:

	//delta-p and delta-theta vs theta

	//sometimes when there's a lot of background, the fit will encapsulate the background in the signal peak
		//to avoid this, first get a reasonable guess for what the signal integral is: do a +/- 5 sigma fit
		//then, for all fits, the signal integral will be compared to this value. If it is much larger, the result is rejected, regardless of the chisq
			//it is allowed to be significantly smaller, in case this fit gives a bad reference
		//if the 5 sigma fit fails, then 1.2*integral_from_3sigma_fit will be used
	locApproximateIntegral *= 1.2; //3sigma
	if(locHistName.find("DeltaTheta") == string::npos)
	{
		//attempt fit
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locGaussMean - 5.0*locGaussWidth;
		double locFitRangeMax = locGaussMean + 5.0*locGaussWidth;
		DFitControl* locFitControl = Attempt_DoubleGaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
		if(locFitControl != nullptr)
		{
			DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = (DFitFunctor_DoubleGaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];
			double locIntegralError = 0.0;
			locFitFunctor_DoubleGaussian->Calc_Integral(locApproximateIntegral, locIntegralError);

			locFitControl->Remove_FitFunctions(locHist);
			delete locFitControl;
		}
	}

	//ok, will have a different routine based on whether there is a background peak that dominates the signal peak or not
	int locMaxBin = locHist->GetMaximumBin();
	int locMeanBin = locHist->GetXaxis()->FindBin(locGaussMean);
	bool locSignalDominatesFlag = (abs(locMaxBin - locMeanBin) < 0.05*locHist->GetNbinsX());

//gDebugFlag = (locHistName == "DeltaThetaVsTheta_5_5");

	//if signal is dominant peak: start with a wide fit range, slowly decrease, stop once the chisq is reasonable
	//otherwise, start with narrow range, slowly increase, stop once the chisq is reasonable
	double locFitRangeNumSigma = locSignalDominatesFlag ? 9.5 : 4.0; //will start at 9, 4.5
	double locMaxChiSqPerDF = 5.0;
	DFitControl* locFitControl = nullptr;
	while(true)
	{
		if(locSignalDominatesFlag)
		{
			//reduce the range
			locFitRangeNumSigma -= 0.5;
			if(locFitRangeNumSigma < 2.9)
				break;
		}
		else
		{
			//increase fit range
			locFitRangeNumSigma += 0.5;
			if(locFitRangeNumSigma > 9.1)
				break;
		}

		//get fit range
		double locFitRangeMin = locGaussMean - locFitRangeNumSigma*locGaussWidth;
		double locFitRangeMax = locGaussMean + locFitRangeNumSigma*locGaussWidth;
		int locFirstBinWithData = 1;
		if(locHistName.find("DeltaTheta") != string::npos)
		{
			for(; locFirstBinWithData <= locHist->GetNbinsX(); ++locFirstBinWithData)
			{
				if(locHist->GetBinContent(locFirstBinWithData) > 10.0)
					break;
			}
			double locLowEdge = locHist->GetXaxis()->GetBinLowEdge(locFirstBinWithData);
			if(locHistName.find("DeltaThetaVsTheta") != string::npos)
				locLowEdge += 5.0;
			if(locLowEdge > locFitRangeMin)
				locFitRangeMin = locLowEdge;
			if(locFitRangeMin > locGaussMean)
				locFitRangeMin = locGaussMean;
		}

if(gDebugFlag)
	cout << "#sigma, fit range: " << locFitRangeNumSigma << ", " << locFitRangeMin << ", " << locFitRangeMax << endl;

		//attempt fit
		double locCutRangeMin, locCutRangeMax;
		locFitControl = Attempt_DoubleGaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);

if(gDebugFlag)
	cout << "fit range #sigma, init gauss width = " << locFitRangeNumSigma << ", " << locGaussWidth << endl;
		if(locFitControl == nullptr)
			continue;

		//check integral
		DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = (DFitFunctor_DoubleGaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];
		double locIntegral = 0.0, locIntegralError = 0.0;
		locFitFunctor_DoubleGaussian->Calc_Integral(locIntegral, locIntegralError);
		if(((locIntegral - locApproximateIntegral)/locApproximateIntegral > 0.33) && (locHistName.find("DeltaTheta") == string::npos))
		{
			locFitControl->Remove_FitFunctions(locHist);
			delete locFitControl;
			locFitControl = nullptr;
			continue; //huge increase in integral: likely a bad fit!
		}

		Double_t locChiSq = locFitControl->dTotalFunc->GetChisquare();
		Int_t locNDF = locFitControl->dTotalFunc->GetNDF();
		double locChiSqPerDF = locChiSq/(double(locNDF));
		if(locChiSqPerDF <= locMaxChiSqPerDF)
			break;

		//chisq not good enough, reduce fit range and try again
		locFitControl->Remove_FitFunctions(locHist);
		delete locFitControl;
		locFitControl = nullptr;
	}

	return locFitControl;
}

/*********************************************************** SMOOTH RESOLUTIONS ***********************************************************/

/*
 From ~/Scratch/gluex/
.x ../Load.C
.L ../scripts_trackeff/Fit_Resolutions.C+
Smooth_Resolutions("trackeff_Proton.root.fit.root");
Smooth_Resolutions("trackeff_Pi-_2pi.root.fit.root");
Smooth_Resolutions("trackeff_Pi+_2pi.root.fit.root");
Smooth_Resolutions("trackeff_Pi-_4pi.root.fit.root");
Smooth_Resolutions("trackeff_Pi+_4pi.root.fit.root");
*/
void Smooth_Resolutions(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	string locOutputFileName = locInputFileName + ".smoothed.root";
	TFile* locOutputFile = new TFile(locOutputFileName.c_str(), "RECREATE");
	gDirectory->mkdir("Resolution", "Resolution");
	gDirectory->cd("Resolution");

	vector<string> locHistTypes = {"DeltaPOverPVsP", "DeltaPOverPVsTheta", "DeltaPOverPVsPhi", "DeltaThetaVsP", "DeltaThetaVsTheta",
			"DeltaThetaVsPhi", "DeltaPhiVsP", "DeltaPhiVsTheta", "DeltaPhiVsPhi"};

	//loop over energy bins
	int locEBin = 0;
	while(true)
	{
		cout << "ebin = " << locEBin << endl;
		ostringstream locEBinDir;
		locEBinDir << "EnergyBin_" << locEBin;
		string locDirName = string("Resolution/") + locEBinDir.str();

		//see if dir exists
		TObject* locObject = locInputFile->Get(locDirName.c_str());
		if(locObject == NULL)
			break; //at the end

		//make output dir
		gDirectory->mkdir(locEBinDir.str().c_str(), locEBinDir.str().c_str());
		gDirectory->cd(locEBinDir.str().c_str());

		//loop over graph-types
		for(auto locHistName : locHistTypes)
		{
//			if(locHistName != "DeltaThetaVsP")
//				continue;

			//get histogram & fit
			string locGraphPath = locDirName + string("/") + locHistName + string("_Means");
			cout << "locGraphPath = " << locGraphPath << endl;
			TGraphErrors* locMeanGraph = (TGraphErrors*)locInputFile->Get(locGraphPath.c_str());
			if(locMeanGraph != nullptr)
			{
				locMeanGraph = (TGraphErrors*)locMeanGraph->Clone(locMeanGraph->GetName());
				Fit_ToLinears(locMeanGraph);
				locMeanGraph->Write();
			}

			locGraphPath = locDirName + string("/") + locHistName + string("_Sigmas");
			cout << "locGraphPath = " << locGraphPath << endl;
			TGraphErrors* locSigmaGraph = (TGraphErrors*)locInputFile->Get(locGraphPath.c_str());
			if(locSigmaGraph != nullptr)
			{
				locSigmaGraph = (TGraphErrors*)locSigmaGraph->Clone(locSigmaGraph->GetName());
				Fit_ToLinears(locSigmaGraph);
				locSigmaGraph->Write();
			}
		}

		gDirectory->cd("..");
		++locEBin;
	}

	locOutputFile->Write();
	locOutputFile->Close();
}

void Fit_ToLinears(TGraphErrors* locGraph)
{
	//strategy: do following steps until the chisq/df is good.
		//first fit a flat function.
		//then fit a linear function (skip this step if vs. phi)
		//then try two linear functions ... then 3, 4, etc. until fit is good
		//the linear functions are constrained to be piecewise continuous
	//for vs phi functions, the function is constrained to be the same at -180 as 180
		//so, a linear fit makes no sense (would be constrained to be flat

	Cleanup_Graph(locGraph, 10);
	Cleanup_Graph(locGraph, 5);

	size_t locGraphN = locGraph->GetN();
	if(locGraphN == 0)
		return;
	double* locGraphX = locGraph->GetX();
	double* locGraphY = locGraph->GetY();
	double* locGraphEY = locGraph->GetEY();
	double locTolerance = 0.001;
	double locMaxChiSqPerDF = (locGraphN > 20) ? 5.0 : 3.0;
	double locFitRangeMin = locGraphX[0] - locTolerance;
	double locFitRangeMax = locGraphX[locGraphN - 1] + locTolerance;

	TF1* locFitFunc = nullptr;
	string locGraphName = locGraph->GetName();
	string locFuncName = locGraphName + string("_Func");

	//first try flat fit
	locFitFunc = new TF1(locFuncName.c_str(), "[0]", locFitRangeMin, locFitRangeMax);
	locFitFunc->SetParameter(0, 0.0);
	locFitFunc->SetLineColor(kRed);
	locGraph->Fit(locFitFunc, "QRN+");
	if(locFitFunc->GetNDF() == 0)
		return; //yikes. we're done
	double locChiSqPerDF = locFitFunc->GetChisquare()/double(locFitFunc->GetNDF());
	if(locChiSqPerDF < locMaxChiSqPerDF)
	{
		locGraph->GetListOfFunctions()->Add(locFitFunc);
		return;
	}

	//nope, try piecewise-linear fits
	ostringstream locFitFuncStream;
	bool locVsPhiFlag = (locGraphName.find("VsPhi") != string::npos);
	bool locFirstTimeFlag = true;
	size_t locNumFitPieces = locVsPhiFlag ? 1 : 0; //incremented at beginning of loop //for vs. phi minimum is 2
//	if(locVsPhiFlag)
//		return;
	vector<double> locFitParams; //initial values
	while(locChiSqPerDF > locMaxChiSqPerDF)
	{
		++locNumFitPieces; //# linear functions
		if(locNumFitPieces > 50)
		{
			//ok, this is getting ridiculous. just do a simple fit and get out of here
			//simple fit: flat for vs phi, linear otherwise
			delete locFitFunc;
			if(locVsPhiFlag)
			{
				locFitFunc = new TF1(locFuncName.c_str(), "[0]", locFitRangeMin, locFitRangeMax);
				locFitFunc->SetParameter(0, 0.0);
			}
			else
			{
				locFitFunc = new TF1(locFuncName.c_str(), "[0] + [1]*x", locFitRangeMin, locFitRangeMax);
				locFitFunc->SetParameters(0.0, 0.0);
			}
			locFitFunc->SetLineColor(kRed);
			locGraph->Fit(locFitFunc, "QRN+");
			locGraph->GetListOfFunctions()->Add(locFitFunc);
			return;
		}

		if(locVsPhiFlag)
		{
			locFitFuncStream.str("");
			locFitFuncStream << "([0] + [1]*x)";
			size_t locPieceCount = 1;
			while(locPieceCount < (locNumFitPieces - 1))
			{
				size_t locInterceptIndex = 2*locPieceCount;
				ostringstream locIntersectionStream;
//				locIntersectionStream << "(([2] - [0])/([1] - [3]))"; // [0] + [1]*x = [2] + [3]*x, x = ([2] - [0])/([1] - [3])
				locIntersectionStream << "(([" << locInterceptIndex << "] - [" << locInterceptIndex - 2 << "])/([" << locInterceptIndex - 1 << "] - [" << locInterceptIndex + 1 << "]))"; // [0] + [1]*x = [2] + [3]*x, x = ([2] - [0])/([1] - [3])
				locFitFuncStream << "*(x < " << locIntersectionStream.str() << ") + ([" << locInterceptIndex << "] + [" << locInterceptIndex + 1 << "]*x)*(x >= " << locIntersectionStream.str() << ")";
				++locPieceCount;
			}

			// last function has to pass through previous intercept (x is fit param) and new intercept:
			// slope = (y_end - y_previous)/(x_end - x_previous)
			// slope = (y_begin - y_previous)/(x_begin + 360.0 - x_previous)

			// plug in y_begin
			// y_begin = [0] + [1]*x_begin
			// slope = ([0] + [1]*x_begin - y_previous)/(x_begin + 360.0 - x_previous)

			// plug in x_begin
			// x_begin = -180.0;
			// slope = ([0] - [1]*180.0 - y_previous)/(180.0 - x_previous)

			size_t locIntersectionIndex = 2*locPieceCount;

			// y_previous = [I - 1]*x_previous + [I - 2]
			// x_previous = [I] // locInterceptIndex
			ostringstream locYPrevious;
			locYPrevious << "([" << locIntersectionIndex - 1 << "]*[" << locIntersectionIndex << "] + [" << locIntersectionIndex - 2 << "])";

			// slope = ([0] - [1]*180.0 - y_previous)/(180.0 - x_previous)
			ostringstream locSlopeStream;
			locSlopeStream << "([0] - [1]*180.0 - " << locYPrevious.str() << ")/(180.0 - [" << locIntersectionIndex << "])";

			// equation:
			// eq = intercept + slope*x
			// intercept = y_previous - slope*x_previous
			// eq = y_previous - slope*x_previous + slope*x
			// eq = y_previous + (x - x_previous)*slope

			ostringstream locEquation;
			locEquation << "(" << locYPrevious.str() << " + (x - [" << locIntersectionIndex << "])*" << locSlopeStream.str() << ")";
			locFitFuncStream << "*(x < [" << locIntersectionIndex << "]) + " << locEquation.str() << "*(x >= [" << locIntersectionIndex << "])";

			//fit func = ([0] + [1]*x)*(x < [2]) + ([0] + [1]*180.0 + (x - 180.0)*([0] - [1]*180.0 - ([1]*[2] + [0]))/(180.0 - [2]))*(x >= [2])
		}
		else if(locFirstTimeFlag)
			locFitFuncStream << "([0] + [1]*x)";
		else
		{
			//apply max limit to previous function (e.g. x < [2]), add new function, apply minimum to it (e.g. x >= 2)
			size_t locInterceptIndex = 2*locNumFitPieces - 2;
			ostringstream locIntersectionStream;
			locIntersectionStream << "(([" << locInterceptIndex << "] - [" << locInterceptIndex - 2 << "])/([" << locInterceptIndex - 1 << "] - [" << locInterceptIndex + 1 << "]))"; // [0] + [1]*x = [2] + [3]*x, x = ([2] - [0])/([1] - [3])
			locFitFuncStream << "*(x < " << locIntersectionStream.str() << ") + ([" << locInterceptIndex << "] + [" << locInterceptIndex + 1 << "]*x)*(x >= " << locIntersectionStream.str() << ")";
		}

		//setup fit params
		if(locFitParams.empty())
		{
			if(locVsPhiFlag)
			{
				double locSlope = (locGraphY[locGraphN/2] - locGraphY[0])/(locGraphX[locGraphN/2] - locGraphX[0]);
				double locIntercept = locGraphY[0] - locSlope*locGraphX[0]; //b = y - mx
				locFitParams.resize(3);
				locFitParams[0] = locIntercept;  locFitParams[1] = locSlope;  locFitParams[2] = locGraphX[locGraphN/2];
			}
			else
			{
				double locSlope = (locGraphY[locGraphN - 1] - locGraphY[0])/(locGraphX[locGraphN - 1] - locGraphX[0]);
				double locIntercept = locGraphY[0] - locSlope*locGraphX[0]; //b = y - mx
				locFitParams.resize(2);
				locFitParams[0] = locIntercept;  locFitParams[1] = locSlope;
			}
		}
		else
		{
			//reset fit params to those of init prior to last fit, so can eval starting func to prepare for next iteration
			for(size_t loc_i = 0; loc_i < locFitParams.size(); ++loc_i)
			{
				locFitFunc->SetParameter(loc_i, locFitParams[loc_i]);
				if(gDebugFlag)
					cout << "i, value = " << loc_i << ", " << locFitParams[loc_i] << endl;
			}

			//new node: location of point with largest discrepancy from previous function //excepting the edges
			double locLargestNumSigma = -1.0;
			int locNewNodeIndex = 1;
			for(size_t loc_j = 1; loc_j < (locGraphN - 1); ++loc_j)
			{
				if(!(locGraphEY[loc_j] > 0.0))
					continue;
				double locNumSigma = fabs(locGraphY[loc_j] - locFitFunc->Eval(locGraphX[loc_j]))/locGraphEY[loc_j];
				if(locNumSigma < locLargestNumSigma)
					continue;
				locLargestNumSigma = locNumSigma;
				locNewNodeIndex = loc_j;
				if(gDebugFlag)
					cout << "x, y, func y, ey, #sigma, largest #sigma = " << locGraphX[loc_j] << ", " << locGraphY[loc_j] << ", " << locFitFunc->Eval(locGraphX[loc_j]) << ", " << locGraphEY[loc_j] << ", " << locNumSigma << ", " << locLargestNumSigma << endl;
			}
			double locNewNodeX = locGraphX[locNewNodeIndex];
			double locNewNodeY = locGraphY[locNewNodeIndex];
			if(gDebugFlag)
				cout << "new node x, y, #sigma = " << locNewNodeX << ", " << locNewNodeY << ", " << locLargestNumSigma << endl;

			//ok, now need to find which line segment to break up, and to insert the new node instead
			double locPreviousIntersectionX = locGraphX[0];
			double locIntersectionX = locGraphX[locGraphN - 1];
			size_t locParamReplaceIndex = 0;

			//loop through linears, find the intersection AFTER the new node
			bool locIntersectionFoundFlag = false;
			for(size_t loc_i = 2; loc_i < locFitParams.size(); loc_i += 2)
			{
				//get intersection between this segment (loc_i) and the previous segment (loc_i - 2)
				if(locVsPhiFlag && (loc_i == locFitParams.size() - 1))
					locIntersectionX = locFitParams[loc_i];
				else
					locIntersectionX = (locFitParams[loc_i] - locFitParams[loc_i - 2])/(locFitParams[loc_i - 1] - locFitParams[loc_i + 1]);

				//is it before or after the new node?
				if(locIntersectionX < locNewNodeX) //before
					locPreviousIntersectionX = locIntersectionX;
				else //after
				{
					locParamReplaceIndex = loc_i - 2;
					locIntersectionFoundFlag = true;
					break;
				}
			}
			//[0] + [1]*x, [2] + [3]*x, [4]
			if(!locIntersectionFoundFlag)
			{
				locIntersectionX = locGraphX[locGraphN - 1];
				locParamReplaceIndex = locVsPhiFlag ? locFitParams.size() - 1 : locFitParams.size() - 2;
			}
			if(gDebugFlag)
				cout << "intersection, previous = " << locIntersectionX << ", " << locPreviousIntersectionX << endl;

			//The above intersection is the node AFTER the new node, and is between this segment (loc_i - 2) and the previous . Break up the previous line segment
			double locIntersectionY = locFitFunc->Eval(locIntersectionX);
			double locPreviousIntersectionY = locFitFunc->Eval(locPreviousIntersectionX);
			if(gDebugFlag)
				cout << "intersection y, previous y = " << locIntersectionY << ", " << locPreviousIntersectionY << endl;

			//first segment
			double locSlope_First = (locNewNodeY - locPreviousIntersectionY)/(locNewNodeX - locPreviousIntersectionX);
			double locIntercept_First = locNewNodeY - locSlope_First*locNewNodeX; //b = y - mx

			//set first segment
			if(gDebugFlag)
			{
				cout << "fit param size, replace index = " << locFitParams.size() << ", " << locParamReplaceIndex << endl;
				cout << "first fit params: " << locIntercept_First << ", " << locSlope_First << endl;
			}
			locFitParams[locParamReplaceIndex] = locIntercept_First;
			//[0] + [1]*x, [2] + [3]*x, [4]

			if(locVsPhiFlag && (locParamReplaceIndex == locFitParams.size() - 1))
			{
				locFitParams.push_back(locSlope_First);
				locFitParams.push_back(locNewNodeX);
			}
			else
			{
				//finish first segment
				locFitParams[locParamReplaceIndex + 1] = locSlope_First;

				//second segment
				double locSlope_Second = (locIntersectionY - locNewNodeY)/(locIntersectionX - locNewNodeX);
				double locIntercept_Second = locNewNodeY - locSlope_Second*locNewNodeX; //b = y - mx

				//save/insert new parameters
				if(gDebugFlag)
					cout << "second fit params: " << locIntercept_Second << ", " << locSlope_Second << endl;
				if(locFitParams.size() == (locParamReplaceIndex + 2))
				{
					locFitParams.push_back(locIntercept_Second);
					locFitParams.push_back(locSlope_Second);
				}
				else
				{
					locFitParams.insert(locFitParams.begin() + locParamReplaceIndex + 2, locIntercept_Second);
					locFitParams.insert(locFitParams.begin() + locParamReplaceIndex + 3, locSlope_Second);
				}
			}
		}

		//delete function for next fit
		delete locFitFunc;

		// create new function
		locFitFunc = new TF1(locFuncName.c_str(), locFitFuncStream.str().c_str(), locFitRangeMin, locFitRangeMax);
		locFitFunc->SetNpx(1000);
		locFitFunc->SetLineColor(kRed);
		if(gDebugFlag)
			cout << "fit func = " << locFitFuncStream.str() << endl;

		//set fit params
		for(size_t loc_i = 0; loc_i < locFitParams.size(); ++loc_i)
		{
			locFitFunc->SetParameter(loc_i, locFitParams[loc_i]);
			if(gDebugFlag)
				cout << "i, value = " << loc_i << ", " << locFitParams[loc_i] << endl;
		}

//		if(locNumFitPieces == 3)
//			break;
//		if(locFitParams.size() > locGraphN)
//			break;
		locGraph->Fit(locFitFunc, "QRN+");

		//check chisq/df
		if(gDebugFlag)
			cout << "chisq, ndf = " << locFitFunc->GetChisquare() << ", " << locFitFunc->GetNDF() << endl;
		if(locFitFunc->GetNDF() <= 0)
			break; //yikes. we're done
		locChiSqPerDF = locFitFunc->GetChisquare()/double(locFitFunc->GetNDF());
		if(gDebugFlag)
			cout << "chisqperdf, max = " << locChiSqPerDF << ", " << locMaxChiSqPerDF << endl;
		locFirstTimeFlag = false;
	}

	if(locFitFunc != nullptr)
		locGraph->GetListOfFunctions()->Add(locFitFunc);
}

void Cleanup_Graph(TGraphErrors* locGraph, int locNumFitPoints)
{
	//fit to linear
	//then, loop through points
	//if a point is way-far-away from the fit, but the adjacent points are close, then disregard the point and refit
	//repeat until done

	int locGraphN = locGraph->GetN();
	if(locGraphN == 0)
		return;
	double* locGraphX = locGraph->GetX();
	double* locGraphY = locGraph->GetY();
	double* locGraphEY = locGraph->GetEY();

	string locGraphName = locGraph->GetName();
	string locFuncName = locGraphName + string("_Func");

	//loop over range, fitting in pieces
	int locMinFitIndex = 0;
	set<size_t> locPointsToRemove;
	while(true)
	{
		int locMaxFitIndex = locMinFitIndex + locNumFitPoints - 1;
		if(locMaxFitIndex >= locGraphN)
			break;

		double locTolerance = 0.001;
		double locFitRangeMin = locGraphX[locMinFitIndex] - locTolerance;
		double locFitRangeMax = locGraphX[locMaxFitIndex] + locTolerance;

		//linear fit
		TF1* locFitFunc = new TF1(locFuncName.c_str(), "[0] + [1]*x", locFitRangeMin, locFitRangeMax);
		locFitFunc->SetParameters(0.0, 0.0);
		locGraph->Fit(locFitFunc, "QRN+");

		//find points to remove, checking only the central bins of the fit range
		int locPlusMinusToCheck = (locNumFitPoints > 8) ? 3 : 1;
		for(int loc_i = locMinFitIndex + locPlusMinusToCheck; loc_i <= locMaxFitIndex - locPlusMinusToCheck; ++loc_i)
		{
			double locFitY = locFitFunc->Eval(locGraphX[loc_i]);
			double locDeltaY = locGraphY[loc_i] - locFitY;
			double locNumSigma = fabs(locDeltaY)/locGraphEY[loc_i];
	//		if(locPrintFlag)
	//			cout << "y, #sigma = " << locGraphY[loc_i] << ", " << locNumSigma << endl;
			if(locNumSigma < 3.0)
				continue; //point not very far from fit: is ok

			//check nearby points //using sigma of the original point, to compare like distances
			int locNumCloseToFitPoints = 0;
			int locPlusMinusNearby = (locNumFitPoints > 8) ? 2 : 1;
			double locNearbyThresholdSigma = (locNumFitPoints > 8) ? 2.0 : 3.0;
			for(int loc_j = loc_i - locPlusMinusNearby; loc_j <= loc_i + locPlusMinusNearby; ++loc_j)
			{
				double locNearbyFitY = locFitFunc->Eval(locGraphX[loc_j]);
				double locNearbyDeltaY = locGraphY[loc_j] - locNearbyFitY;
				double locNearbyNumSigma = fabs(locNearbyDeltaY)/locGraphEY[loc_i];
				if(locNearbyNumSigma < locNearbyThresholdSigma)
					++locNumCloseToFitPoints;
			}

			//if many other points close, reject this point
			int locNumCloseTightThreshold = (locPlusMinusNearby == 2) ? 3 : 2;
			if((locNumSigma > 3.0) && (locNumCloseToFitPoints >= locNumCloseTightThreshold))
				locPointsToRemove.insert(loc_i);
			else if((locNumSigma > 4.0) && (locNumCloseToFitPoints >= (locNumCloseTightThreshold - 1)))
				locPointsToRemove.insert(loc_i);
		}

		delete locFitFunc;
		++locMinFitIndex;
	}

//	cout << "remove points: " << locGraph->GetName() << ": ";
	for(auto locIterator = locPointsToRemove.rbegin(); locIterator != locPointsToRemove.rend(); ++locIterator)
	{
//		cout << locGraphX[*locIterator] << ", ";
		locGraph->RemovePoint(*locIterator);
	}
//	cout << endl;
}

/*********************************************************** PLOT RESOLUTIONS ***********************************************************/

/*
 From ~/Scratch/gluex/
//MERGE pass1 & pass2 results first!!
.x ../Load.C
.L ../scripts_trackeff/Fit_Resolutions.C+
Plot_Resoultions("trackeff_Proton.root.fit.root.smoothed.root");
Plot_Resoultions("trackeff_Pi-_2pi.root.fit.root.smoothed.root");
Plot_Resoultions("trackeff_Pi-_4pi.root.fit.root.smoothed.root");
Plot_Resoultions("trackeff_Pi+_2pi.root.fit.root.smoothed.root");
Plot_Resoultions("trackeff_Pi+_4pi.root.fit.root.smoothed.root");
*/
void Plot_Resoultions(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	vector<string> locHistTypes = {"DeltaPOverPVsP", "DeltaPOverPVsTheta", "DeltaPOverPVsPhi", "DeltaThetaVsP", "DeltaThetaVsTheta",
			"DeltaThetaVsPhi", "DeltaPhiVsP", "DeltaPhiVsTheta", "DeltaPhiVsPhi"};

	map<string, int> locTotalNumPointsMap_Means, locTotalNumPointsMap_Sigmas;
	for(auto locHistName : locHistTypes)
	{
		locTotalNumPointsMap_Means[locHistName] = 0;
		locTotalNumPointsMap_Sigmas[locHistName] = 0;
	}

	//loop over energy bins
	int locEBin = 0;
	int locNumEBins = 0;
	map<string, vector<TGraphErrors*> > locMeanGraphMap, locSigmaGraphMap;
	while(true)
	{
		cout << "ebin = " << locEBin << endl;
		ostringstream locEBinDir;
		locEBinDir << "EnergyBin_" << locEBin;
		string locDirName = string("Resolution/") + locEBinDir.str();

		//see if dir exists
		TObject* locObject = locInputFile->Get(locDirName.c_str());
		if(locObject == NULL)
			break; //at the end
		++locNumEBins;

		//loop over graph-types
		for(auto locHistName : locHistTypes)
		{
			//get histogram & fit
			string locGraphPath = locDirName + string("/") + locHistName + string("_Means");
			TGraphErrors* locGraph = (TGraphErrors*)locInputFile->Get(locGraphPath.c_str());
			locMeanGraphMap[locHistName].push_back(locGraph);
			if(locGraph != nullptr)
				locTotalNumPointsMap_Means[locHistName] += locGraph->GetN();

			locGraphPath = locDirName + string("/") + locHistName + string("_Sigmas");
			locGraph = (TGraphErrors*)locInputFile->Get(locGraphPath.c_str());
			locSigmaGraphMap[locHistName].push_back(locGraph);
			if(locGraph != nullptr)
				locTotalNumPointsMap_Sigmas[locHistName] += locGraph->GetN();
		}

		++locEBin;
	}

	string locOutputFileName = locInputFileName + ".plots.root";
	TFile* locOutputFile = new TFile(locOutputFileName.c_str(), "RECREATE");
	gDirectory->mkdir("Resolution", "Resolution");
	gDirectory->cd("Resolution");

	//Create TGraph2D's
	double locMinBeamEnergy = 3.0;
	double locMaxBeamEnergy = 12.0;
	double locDeltaEnergyPerBin = (locMaxBeamEnergy - locMinBeamEnergy)/double(locNumEBins);
	for(auto locHistName : locHistTypes)
	{
		int locNumPoints_Means = locTotalNumPointsMap_Means[locHistName];
		int locNumPoints_Sigmas = locTotalNumPointsMap_Sigmas[locHistName];
		cout << "hist name = " << locHistName << endl;
		if((locNumPoints_Means == 0) || (locNumPoints_Sigmas == 0))
			continue;

		//means
		int locNShift = 0;
		double* loc2DXArray_Means = new double[locNumPoints_Means]; //p, theta, phi
		double* loc2DYArray_Means = new double[locNumPoints_Means]; //beam energy
		double* loc2DZArray_Means = new double[locNumPoints_Means]; //the mean or sigma
		auto& locMeanGraphs = locMeanGraphMap[locHistName];

		//loop
		double locBeamEnergy = locMinBeamEnergy + 0.5*locDeltaEnergyPerBin;
		for(auto locGraph : locMeanGraphs)
		{
			string locGraphName = locGraph->GetName();
			string locFuncName = locGraphName + string("_Func");
			TF1* locFunc = (TF1*)locGraph->GetListOfFunctions()->FindObject(locFuncName.c_str());

			int locGraphN = locGraph->GetN();
			double* locGraphX = locGraph->GetX();
			double* locGraphY = locGraph->GetY();
			for(int loc_i = 0; loc_i < locGraphN; ++loc_i)
			{
				loc2DXArray_Means[locNShift + loc_i] = locGraphX[loc_i];
				loc2DYArray_Means[locNShift + loc_i] = locBeamEnergy;
				if(locFunc == nullptr)
					loc2DZArray_Means[locNShift + loc_i] = locGraphY[loc_i];
				else
					loc2DZArray_Means[locNShift + loc_i] = locFunc->Eval(locGraphX[loc_i]);
			}
			locNShift += locGraphN;
			locBeamEnergy += locDeltaEnergyPerBin;
		}

		//sigmas
		locNShift = 0;
		double* loc2DXArray_Sigmas = new double[locNumPoints_Sigmas]; //p, theta, phi
		double* loc2DYArray_Sigmas = new double[locNumPoints_Sigmas]; //beam energy
		double* loc2DZArray_Sigmas = new double[locNumPoints_Sigmas]; //the mean or sigma
		auto& locSigmaGraphs = locSigmaGraphMap[locHistName];

		//loop
		locBeamEnergy = locMinBeamEnergy + 0.5*locDeltaEnergyPerBin;
		for(auto locGraph : locSigmaGraphs)
		{
			string locGraphName = locGraph->GetName();
			string locFuncName = locGraphName + string("_Func");
			TF1* locFunc = (TF1*)locGraph->GetListOfFunctions()->FindObject(locFuncName.c_str());

			int locGraphN = locGraph->GetN();
			double* locGraphX = locGraph->GetX();
			double* locGraphY = locGraph->GetY();
			for(int loc_i = 0; loc_i < locGraphN; ++loc_i)
			{
				loc2DXArray_Sigmas[locNShift + loc_i] = locGraphX[loc_i];
				loc2DYArray_Sigmas[locNShift + loc_i] = locBeamEnergy;
				if(locFunc == nullptr)
					loc2DZArray_Sigmas[locNShift + loc_i] = locGraphY[loc_i];
				else
					loc2DZArray_Sigmas[locNShift + loc_i] = locFunc->Eval(locGraphX[loc_i]);
			}
			locNShift += locGraphN;
			locBeamEnergy += locDeltaEnergyPerBin;
		}

		TGraph2D* locGraph2D_Means = new TGraph2D(locNumPoints_Means, loc2DXArray_Means, loc2DYArray_Means, loc2DZArray_Means);
		string locGraphName = locHistName + "_Means";
		string locGraphTitle = string(";") + string(locMeanGraphs.back()->GetXaxis()->GetTitle()) + string(";Beam Energy (GeV);") + string(locMeanGraphs.back()->GetYaxis()->GetTitle());
		locGraph2D_Means->SetName(locGraphName.c_str());
		locGraph2D_Means->SetTitle(locGraphTitle.c_str());
		locGraph2D_Means->Write();

		TGraph2D* locGraph2D_Sigmas = new TGraph2D(locNumPoints_Sigmas, loc2DXArray_Sigmas, loc2DYArray_Sigmas, loc2DZArray_Sigmas);
		locGraphName = locHistName + "_Sigmas";
		locGraphTitle = string(";") + string(locSigmaGraphs.back()->GetXaxis()->GetTitle()) + string(";Beam Energy (GeV);") + string(locSigmaGraphs.back()->GetYaxis()->GetTitle());
		locGraph2D_Sigmas->SetName(locGraphName.c_str());
		locGraph2D_Sigmas->SetTitle(locGraphTitle.c_str());
		locGraph2D_Sigmas->Write();
	}

	//save results
	locOutputFile->Write();
	locOutputFile->Close();
}

/*
.x ../Load.C
.L ../scripts_trackeff/Fit_Resolutions.C+
Draw_Resolutions("trackeff_Proton.root.fit.root.smoothed.root.plots.root");
Draw_Resolutions("trackeff_Pi-_2pi.root.fit.root.smoothed.root.plots.root");
Draw_Resolutions("trackeff_Pi-_4pi.root.fit.root.smoothed.root.plots.root");
Draw_Resolutions("trackeff_Pi+_2pi.root.fit.root.smoothed.root.plots.root");
Draw_Resolutions("trackeff_Pi+_4pi.root.fit.root.smoothed.root.plots.root");
*/
void Draw_Resolutions(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	vector<string> locHistTypes = {"DeltaPOverPVsP", "DeltaPOverPVsTheta", "DeltaPOverPVsPhi", "DeltaThetaVsP", "DeltaThetaVsTheta",
			"DeltaThetaVsPhi", "DeltaPhiVsP", "DeltaPhiVsTheta", "DeltaPhiVsPhi"};

	//loop over graph-types
	string locDrawArgument = "TRI1Z";
	for(auto locHistName : locHistTypes)
	{
		//get histogram & fit
		string locGraphPath = string("Resolution/") + locHistName + string("_Means");
		TGraph2D* locGraph2D = (TGraph2D*)locInputFile->Get(locGraphPath.c_str());
		if(locGraph2D != nullptr)
		{
			new TCanvas();
			locGraph2D->Draw(locDrawArgument.c_str());
		}

		locGraphPath = string("Resolution/") + locHistName + string("_Sigmas");
		locGraph2D = (TGraph2D*)locInputFile->Get(locGraphPath.c_str());
		if(locGraph2D != nullptr)
		{
			new TCanvas();
			locGraph2D->Draw(locDrawArgument.c_str());
		}
	}
}

/*********************************************************** FIT UTILITIES ***********************************************************/

DFitControl* Attempt_DoubleGaussFit(TH1* locHist, double locGaussHeight, double locGaussMean, double locInitGaussWidth, double& locFitRangeMin, double& locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax)
{
	if(locFitRangeMin < locHist->GetXaxis()->GetBinLowEdge(1))
		locFitRangeMin = locHist->GetXaxis()->GetBinLowEdge(1);
	if(locFitRangeMax > locHist->GetXaxis()->GetBinUpEdge(locHist->GetNbinsX()))
		locFitRangeMax = locHist->GetXaxis()->GetBinUpEdge(locHist->GetNbinsX());

	//fit
	if(fabs(locFitRangeMax - locFitRangeMin) <= locHist->GetBinWidth(1))
		return nullptr;

//	if(string(locHist->GetName()) == "DeltaPhiVsP_24_24")
//		cout << "height, mean, width, #sigmas, fit min, fit max = " << locGaussHeight << ", " << locGaussMean << ", " << locInitGaussWidth << ", " << locNumSigmas << ", " << locFitRangeMin << ", " << locFitRangeMax << endl;
	DFitControl* locFitControl = Fit_DoubleGauss(locHist, locGaussHeight, locGaussMean, locInitGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);

	//eval background
	if(!Eval_Fit(locHist, locFitControl))
	{
		locFitControl->Remove_FitFunctions(locHist);
		delete locFitControl;
		return nullptr;
	}

	return locFitControl;
}

DFitControl* Attempt_GaussFit(TH1* locHist, double locGaussHeight, double locGaussMean, double locInitGaussWidth, double& locFitRangeMin, double& locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax)
{
	if(locFitRangeMin < locHist->GetXaxis()->GetBinLowEdge(1))
		locFitRangeMin = locHist->GetXaxis()->GetBinLowEdge(1);
	if(locFitRangeMax > locHist->GetXaxis()->GetBinUpEdge(locHist->GetNbinsX()))
		locFitRangeMax = locHist->GetXaxis()->GetBinUpEdge(locHist->GetNbinsX());

	//fit
	if(fabs(locFitRangeMax - locFitRangeMin) <= locHist->GetBinWidth(1))
		return nullptr;

//	if(string(locHist->GetName()) == "DeltaPhiVsP_24_24")
//		cout << "height, mean, width, #sigmas, fit min, fit max = " << locGaussHeight << ", " << locGaussMean << ", " << locInitGaussWidth << ", " << locNumSigmas << ", " << locFitRangeMin << ", " << locFitRangeMax << endl;
	DFitControl* locFitControl = Fit_Gaussian(locHist, locGaussHeight, locGaussMean, locInitGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);

	//eval background
	if(!Eval_Fit_Gaussian(locHist, locFitControl))
	{
		locFitControl->Remove_FitFunctions(locHist);
		delete locFitControl;
		return nullptr;
	}

	return locFitControl;
}

bool Eval_Fit(TH1* locHist, DFitControl* locFitControl)
{
	DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = (DFitFunctor_DoubleGaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];
	DFitFunctor_Polynomial* locFitFunctor_Polynomial = (DFitFunctor_Polynomial*)locFitControl->dBackgroundFitGroup->dFitFunctors[0];

	double locGaussMean = locFitFunctor_DoubleGaussian->dFunc->GetParameter(1);
	double locGaussSigma = locFitFunctor_DoubleGaussian->Calc_StandardDeviation().first;

	double locFitRangeMin = locFitControl->dFitRangeMin;
	double locFitRangeMax = locFitControl->dFitRangeMax;

	//check to see if the fit mean is not close to the center of the histogram
	double locFitRangeCenter = 0.5*(locFitRangeMin + locFitRangeMax);
	double locFracDistToCenter = 0.0;
	if(locGaussMean > locFitRangeCenter)
		locFracDistToCenter = (locGaussMean - locFitRangeCenter)/(locFitRangeMax - locFitRangeCenter); //5 / 20: 1
	else
		locFracDistToCenter = (locFitRangeCenter - locGaussMean)/(locFitRangeCenter - locFitRangeMin); //20

	string locHistName = locHist->GetName();
	if((locFracDistToCenter > 0.75) && (locHistName.find("DeltaThetaVsTheta") == string::npos))
		return false;

	double locFirstSigma = locFitFunctor_DoubleGaussian->dFunc->GetParameter(2);
	double locSecondSigma = locFitFunctor_DoubleGaussian->dFunc->GetParameter(4);

	double locWiderSigma = (locSecondSigma > locFirstSigma) ? locSecondSigma : locFirstSigma;
	double locNarrowerSigma = (locSecondSigma > locFirstSigma) ? locFirstSigma : locSecondSigma;

	double locFirstArea = locFitFunctor_DoubleGaussian->dFunc->GetParameter(0);
	double locSecondArea = locFitFunctor_DoubleGaussian->dFunc->GetParameter(3);
	double locWiderArea = (locSecondSigma > locFirstSigma) ? locSecondArea : locFirstArea;
	double locNarrowerArea = (locSecondSigma > locFirstSigma) ? locFirstArea : locSecondArea;

//	if(string(locHist->GetName()) == "DeltaThetaVsTheta_28_28")
//		cout << "wide/narrow sigmas/areas: " << locWiderSigma << ", " << locNarrowerSigma << ", " << locWiderArea << ", " << locNarrowerArea << endl;
//	if((locWiderSigma/locNarrowerSigma > 3.0) && (locWiderArea/locNarrowerArea < 1.0))
//		return false; //if wide sigma is very wide, but area is small

	return true;
}

bool Eval_Fit_Gaussian(TH1* locHist, DFitControl* locFitControl)
{
	DFitFunctor_Gaussian* locFitFunctor_Gaussian = (DFitFunctor_Gaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];
	DFitFunctor_Polynomial* locFitFunctor_Polynomial = (DFitFunctor_Polynomial*)locFitControl->dBackgroundFitGroup->dFitFunctors[0];

	double locGaussMean = locFitFunctor_Gaussian->dFunc->GetParameter(1);
	double locGaussSigma = locFitFunctor_Gaussian->dFunc->GetParameter(2);

	double locFitRangeMin = locFitControl->dFitRangeMin;
	double locFitRangeMax = locFitControl->dFitRangeMax;

	//check to see if the fit mean is not close to the center of the histogram
	double locFitRangeCenter = 0.5*(locFitRangeMin + locFitRangeMax);
	double locFracDistToCenter = 0.0;
	if(locGaussMean > locFitRangeCenter)
		locFracDistToCenter = (locGaussMean - locFitRangeCenter)/(locFitRangeMax - locFitRangeCenter); //5 / 20: 1
	else
		locFracDistToCenter = (locFitRangeCenter - locGaussMean)/(locFitRangeCenter - locFitRangeMin); //20

	string locHistName = locHist->GetName();
	if((locFracDistToCenter > 0.75) && (locHistName.find("DeltaThetaVsTheta") == string::npos))
		return false;

	return true;
}

DFitControl* Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax)
{
	//Build Functors
	DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = Build_DoubleGaussian(locGaussHeight, locGaussMean, locGaussWidth);
	DFitFunctor* locFitFunctor_Background = nullptr; //Build_PolyBackground(locHist, locFitRangeMin, locFitRangeMax);
	locFitFunctor_Background = Build_PolyBackground_Linear(locHist, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax);

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_DoubleGaussian, locFitFunctor_Background);
	locFitControl->dLogLikelihoodFlag = false;

	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);
//	cout << "FIT STATUS = " << locFitStatus << endl;

	//cut range
	locFitFunctor_DoubleGaussian->Calc_CutRange(3.0, locCutRangeMin, locCutRangeMax);

	return locFitControl;
}

DFitControl* Fit_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax)
{
	//Build Functors
	DFitFunctor_Gaussian* locFitFunctor_Gaussian = Build_Gaussian(locGaussHeight, locGaussMean, locGaussWidth);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground_Linear(locHist, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax);

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_Gaussian, locFitFunctor_Background);
	locFitControl->dLogLikelihoodFlag = false;

	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);
//	cout << "FIT STATUS = " << locFitStatus << endl;

	//cut range
	locFitFunctor_Gaussian->Calc_CutRange(3.0, locCutRangeMin, locCutRangeMax);

	return locFitControl;
}

double Calc_InitGaussHeight(TH1* locHist, double locGaussMean, double locGaussWidth, double locNumSigmas)
{
	int locLowSideBin = locHist->GetXaxis()->FindBin(locGaussMean - locNumSigmas*locGaussWidth);
	if(locLowSideBin < 1)
		locLowSideBin = 1;
	int locHighSideBin = locHist->GetXaxis()->FindBin(locGaussMean + locNumSigmas*locGaussWidth);
	if(locHighSideBin > locHist->GetNbinsX())
		locHighSideBin = locHist->GetNbinsX();
	double locAvgBackground = 0.5*(locHist->GetBinContent(locLowSideBin) + locHist->GetBinContent(locHighSideBin));

	int locPeakBin = locHist->GetXaxis()->FindBin(locGaussMean);
	return locHist->GetBinContent(locPeakBin) - locAvgBackground;
}

DFitFunctor* Build_PolyBackground_Flat(TH1* locHist, double locFitRangeMin, double locFitRangeMax)
{
	double locIntercept = 0.0;
	DFitUtilities::Guess_Params_Flat(locHist, locFitRangeMin, locFitRangeMax, locIntercept);

	vector<double> locInitParams(1, locIntercept);
	return (new DFitFunctor_Polynomial(locInitParams, kMagenta));
}

DFitFunctor* Build_PolyBackground_Linear(TH1* locHist, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax)
{
	double locSlope = 0.0, locIntercept = 0.0;

	if((locGaussMean - 3.0*locGaussWidth) < locFitRangeMin)
	{
		locIntercept = locHist->GetBinContent(locHist->GetXaxis()->FindBin(locFitRangeMax));
//		return (new DFitFunctor_Polynomial(vector<double>(1, locIntercept), kMagenta)); //flat!
	}
	else if((locGaussMean + 3.0*locGaussWidth) > locFitRangeMax)
	{
		locIntercept = locHist->GetBinContent(locHist->GetXaxis()->FindBin(locFitRangeMin));
//		return (new DFitFunctor_Polynomial(vector<double>(1, locIntercept), kMagenta)); //flat!
	}
	else
		DFitUtilities::Guess_Params_Linear(locHist, locFitRangeMin, locFitRangeMax, locIntercept, locSlope, true);

	vector<double> locInitParams;
	locInitParams.resize(2);
	locInitParams[0] = locIntercept;  locInitParams[1] = locSlope;
	return (new DFitFunctor_LinearAboveZero(locInitParams, locFitRangeMin, locFitRangeMax, kMagenta));
}

DFitFunctor_DoubleGaussian* Build_DoubleGaussian(double locGaussHeight, double locGaussMean, double locGaussWidth)
{
	vector<double> locInitParams(5, 0.0);
	//peak gaussian
	locInitParams[0] = 0.85*locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = 0.9*locGaussWidth;
	//non-gaussian-tail gaussian
	locInitParams[3] = 0.15*locGaussHeight;  locInitParams[4] = 1.5*locGaussWidth;
	DFitFunctor_DoubleGaussian* locDoubleGaussian = new DFitFunctor_DoubleGaussian(locInitParams, kBlue);
	return locDoubleGaussian;
}

DFitFunctor_Gaussian* Build_Gaussian(double locGaussHeight, double locGaussMean, double locGaussWidth)
{
	vector<double> locInitParams(3, 0.0);
	locInitParams[0] = locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = locGaussWidth;

	DFitFunctor_Gaussian* locGaussian = new DFitFunctor_Gaussian(locInitParams, kBlue);
	return locGaussian;
}

double* Convert_VectorToArray(const vector<double>& locVector)
{
	double* locArray = new double[locVector.size()];
	for(size_t loc_i = 0; loc_i < locVector.size(); ++loc_i)
		locArray[loc_i] = locVector[loc_i];
	return locArray;
}