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

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

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

#include "PlotData/DStatUtilities.h"

#include "particleType.h"

//MAIN
void Fit_MassPeaks(string locInputFileName, Particle_t locPID);
void Smooth_MassCuts(string locInputFileName);

//fit control
DFitControl* Fit_MassPeak(TH1* locHist, Particle_t locPID);
bool Eval_Fit(TH1* locHist, DFitControl* locFitControl, double locNumSigmas);
DFitControl* Attempt_DoubleGaussFit(TH1* locHist, double locGaussHeight, double locGaussMean, double locInitGaussWidth, double& locFitRangeMin, double& locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);

//UTILITY
void Rebin_Histogram(TH1* locHist);
void Find_InitGaussParams(TH1* locHist, double& locGaussHeight, double& locGaussMean, double& locGaussWidth, double& locApproximateIntegral);

//smooth
TF1* Fit_ToLinears(TGraphErrors* locGraph);
void Cleanup_Graph(TGraphErrors* locGraph);

//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);

double gMinPeakYield = 2000.0;
bool gDebugFlag = false;

/*********************************************************** SUBTRACT RF SIDEBANDS AND FIT ***********************************************************/

/*
.x ../Load.C
.L ../scripts_trackeff/Fit_MassCuts.C+
Fit_MassPeaks("trackeff_Proton.root", Proton);
Fit_MassPeaks("trackeff_Pi-_2pi.root", PiMinus);
Fit_MassPeaks("trackeff_Pi+_2pi.root", PiPlus);
Fit_MassPeaks("trackeff_Pi-_4pi.root", PiMinus);
Fit_MassPeaks("trackeff_Pi+_4pi.root", PiPlus);
*/
void Fit_MassPeaks(string locInputFileName, Particle_t locPID)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

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

	// Found & Missing
	vector<TH2*> loc2DMassHists;
	loc2DMassHists.push_back((TH2*)locInputFile->Get("MissingMass/MissingMassVsBeamEnergy_Found"));
	loc2DMassHists.push_back((TH2*)locInputFile->Get("MissingMass/MissingMassVsBeamEnergy_Missing"));

	// Reconned
//	loc2DMassHists.push_back((TH2*)locInputFile->Get("MissingMass/ReconMissingMassSquaredVsBeamEnergy"));


	//also, recombine found/missing into one histogram
	//can also combine kinematics for track found/missing
	//then, can use this to sideband subtract, giving kinematics in the peak
		//dunno what fraction are found & missing though
	TH2* locCombinedMassHist = (TH2*)loc2DMassHists[0]->Clone("MissingMassVsBeamEnergy");
	locCombinedMassHist->Add(loc2DMassHists[1]);
	loc2DMassHists.push_back(locCombinedMassHist);

	//Now, split up into beam energy bins

	//Find # beam energy bins
	int locEBin = 0; //will be #bins at the end
	while(true)
	{
		//make energy bin dir name
		ostringstream locEBinDir;
		locEBinDir << "EnergyBin_" << locEBin;
		string locDirName = string("Resolution/") + locEBinDir.str();
		if(locInputFile->Get(locDirName.c_str()) == NULL)
			break; //at the end
		++locEBin;
	}

	//TH2*: found, missing
	map<TH2*, vector<double> > locXs, locMeans, locMeanUncertainties, locSigmas, locSigmaUncertainties;

	//project & fit
	int locNumBeamBinsToMerge = loc2DMassHists[0]->GetNbinsX()/locEBin;
	for(auto& loc2DHist : loc2DMassHists) //loop over found/missing
	{
		int locStartXBin = 1;
		for(int loc_i = 0; loc_i < locEBin; ++loc_i) //loop over beam E bins
		{
			int locEndXRangeBin = locStartXBin + locNumBeamBinsToMerge - 1;
			double locMinBeamEnergy = loc2DHist->GetXaxis()->GetBinLowEdge(locStartXBin);
			double locMaxBeamEnergy = loc2DHist->GetXaxis()->GetBinUpEdge(locEndXRangeBin);

			//project
			ostringstream locHistName;
			locHistName << loc2DHist->GetName() << "_" << loc_i;
			TH1* locProjHist = (TH1*)loc2DHist->ProjectionY(locHistName.str().c_str(), locStartXBin, locEndXRangeBin);
			cout << "hist name = " << locHistName.str() << endl;

			//set title
			ostringstream locHistTitle;
			locHistTitle << locProjHist->GetTitle() << ", " << locMinBeamEnergy << " < E_{#gamma} (GeV) < " << locMaxBeamEnergy;
			locProjHist->SetTitle(locHistTitle.str().c_str());

			//Rebin
			Rebin_Histogram(locProjHist);

			//Fit to double-Gaussian + linear background bound above zero
			DFitControl* locFitControl = Fit_MassPeak(locProjHist, locPID);
			if(locFitControl == nullptr)
			{
				locStartXBin += locNumBeamBinsToMerge;
				continue;
			}

			//get fit functions
			DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = (DFitFunctor_DoubleGaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];
			DFitFunctor* locFitFunctor_Background = locFitControl->dBackgroundFitGroup->dFitFunctors[0];

			//save fit results
			double locAvgBeamEnergy = 0.5*(locMinBeamEnergy + locMaxBeamEnergy);
			locXs[loc2DHist].push_back(locAvgBeamEnergy);
			locMeans[loc2DHist].push_back(locFitFunctor_DoubleGaussian->dFunc->GetParameter(1));
			locMeanUncertainties[loc2DHist].push_back(locFitFunctor_DoubleGaussian->dFunc->GetParError(1));

			pair<double, double> locSigmaPair = locFitFunctor_DoubleGaussian->Calc_StandardDeviation();
			locSigmas[loc2DHist].push_back(locSigmaPair.first);
			locSigmaUncertainties[loc2DHist].push_back(locSigmaPair.second);

			//prepare for next iteration
			locStartXBin += locNumBeamBinsToMerge;
		}
	}

	//make graphs
	for(auto& locMapPairs : locXs)
	{
		TH2* loc2DHist = locMapPairs.first;

		string locHistName = loc2DHist->GetName();
		string locHistTitle = loc2DHist->GetTitle();

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

		//create graph of means
		TGraphErrors* locMeanGraph = new TGraphErrors(locArraySize, locXArray, locMeanArray, NULL, locMeanUncertaintyArray);
		locMeanGraph->SetName((locHistName + string("_Means")).c_str());
		locMeanGraph->SetTitle(locHistTitle.c_str());
		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(locArraySize, locXArray, locSigmaArray, NULL, locSigmaUncertaintyArray);
		locSigmaGraph->SetName((locHistName + string("_Sigmas")).c_str());
		locSigmaGraph->SetTitle(locHistTitle.c_str());
		locSigmaGraph->GetXaxis()->SetTitle(loc2DHist->GetXaxis()->GetTitle());
		locSigmaGraph->GetYaxis()->SetTitle((string("#sigma ") + string(loc2DHist->GetYaxis()->GetTitle())).c_str());
		locSigmaGraph->Write();
	}

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

void Rebin_Histogram(TH1* locHist)
{
	//compute average jaggedness of histogram
	double locAverageDifferenceFraction = 0.0;
	for(int loc_i = 1; loc_i < locHist->GetNbinsX(); ++loc_i)
	{
		double locDifference = locHist->GetBinContent(loc_i + 1) - locHist->GetBinContent(loc_i);
		double locDifferenceFraction = locDifference/locHist->GetBinContent(loc_i);
		locAverageDifferenceFraction += fabs(locDifferenceFraction);
	}
	locAverageDifferenceFraction /= double(locHist->GetNbinsX() - 1);

	int locRebinTerm = int(locAverageDifferenceFraction/0.02 + 1.5);
	while((locHist->GetNbinsX() % locRebinTerm) != 0)
		++locRebinTerm; //make sure it rebins evenly
	if(locRebinTerm > 1)
		locHist->Rebin(locRebinTerm);
//	/cout << "avg diff frac, rebin term = " << locAverageDifferenceFraction << ", " << locRebinTerm << endl;

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

DFitControl* Fit_MassPeak(TH1* locHist, Particle_t locPID)
{
	//assume it's clean: no messiness
	double locGaussHeight = 0.0, locGaussMean = ParticleMass(locPID), locGaussWidth = 0.0, locApproximateIntegral = 0.0;
	Find_InitGaussParams(locHist, locGaussHeight, locGaussMean, locGaussWidth, locApproximateIntegral);

	//found
	string locHistName = locHist->GetName();
	if(locHistName.find("Found") != 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);
	}

	//start with a wide fit range, slowly decrease, stop once the chisq is reasonable
	double locFitRangeNumSigma = 9.5; //will start at 9
	double locMaxChiSqPerDF = 5.0;
	DFitControl* locFitControl = nullptr;
	while(true)
	{
		//reduce the range
		locFitRangeNumSigma -= 0.5;
		if(locFitRangeNumSigma < 2.9)
			break;

		//attempt fit
		double locCutRangeMin, locCutRangeMax;
		double locFitRangeMin = locGaussMean - locFitRangeNumSigma*locGaussWidth;
		double locFitRangeMax = locGaussMean + locFitRangeNumSigma*locGaussWidth;
		locFitControl = Attempt_DoubleGaussFit(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
		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)
		{
			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;
}

void Find_InitGaussParams(TH1* locHist, double& locGaussHeight, double& locGaussMean, double& locGaussWidth, double& locApproximateIntegral)
{
	//get slope at the end of the histogram
	double locLinearFitMin = locHist->GetBinLowEdge(1);
	double locLinearFitMax = locHist->GetBinLowEdge(locHist->GetNbinsX()/20);
	double locLinFitIntercept = 0.0, locLinFitSlope = 0.0;
	TF1* locLinearFunc = DFitUtilities::Guess_Params_Linear_Fit(locHist, locLinearFitMin, locLinearFitMax, locLinFitIntercept, locLinFitSlope);

	//get width:
	//step to the side until hist height becomes relatively flat, or increases: fit linear slopes to rolling range of 4 bins
	int locFitRangeNumBins = 7;
	//first the low side
	int locFitRangeMinBin = locHist->FindBin(locGaussMean) - locFitRangeNumBins + 1; //+1: will decrease below
	while(locFitRangeMinBin > 1)
	{
		--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(fabs(locSlope - locLinFitSlope)/locLinFitSlope < 0.05)
			break; //slope is significantly reduced: we are now at the edge
	}
	locGaussWidth = (locGaussMean - locHist->GetBinLowEdge(locFitRangeMinBin))/3.0;

	//get init height guess
	locGaussHeight = Calc_InitGaussHeight(locHist, locGaussMean, locGaussWidth, 5.0);

	//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());
}

/*********************************************************** SMOOTH MASS CUTS ***********************************************************/

/*
 From ~/Scratch/gluex/
.x ../Load.C
.L ../scripts_trackeff/Fit_MassCuts.C+
Smooth_MassCuts("trackeff_Proton.root.mass_fit.root");
Smooth_MassCuts("trackeff_Pi-_2pi.root.mass_fit.root");
Smooth_MassCuts("trackeff_Pi+_2pi.root.mass_fit.root");
Smooth_MassCuts("trackeff_Pi-_4pi.root.mass_fit.root");
Smooth_MassCuts("trackeff_Pi+_4pi.root.mass_fit.root");
*/
void Smooth_MassCuts(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("MissingMass", "MissingMass");
	gDirectory->cd("MissingMass");

	vector<string> locHistTypes = {"MissingMassVsBeamEnergy_Found", "MissingMassVsBeamEnergy_Missing"};

	//bool: found flag
	map<bool, vector<double> > locAvgBeamEnergies, locFactors, locFactorUncertainties, locSidebandRangeFractions;

	//loop over graph-types
	for(size_t loc_i = 0; loc_i < locHistTypes.size(); ++loc_i)
	{
		string locHistName = locHistTypes[loc_i];
		bool locFoundFlag = (loc_i == 0);

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

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

		//for mass sideband subtraction:
		//need integrals: #total events in sidebands, background under signal,
		if((locMeanFitFunc == nullptr) || (locSigmaFitFunc == nullptr))
			continue;

		//compute the sideband subtraction factor (is -1* the factor)
		int locEBin = 0;
		while(true)
		{
			ostringstream locMassHistName;
			locMassHistName << "MissingMass/" << locHistName << "_" << locEBin;
			TH1* locMassHist = (TH1*)locInputFile->Get(locMassHistName.str().c_str());
			if(locMassHist == NULL)
				break; //at the end

			string locFuncName = string(locMassHist->GetName()) + string("_Total");
			TF1* locTotalFitFunc = locMassHist->GetFunction(locFuncName.c_str());
			if(locTotalFitFunc == NULL) //fit not performed
				continue;

			//get the fit
			locFuncName = string(locMassHist->GetName()) + string("_Signal_0");
			TF1* locSignalFitFunc = locMassHist->GetFunction(locFuncName.c_str());
			DFitFunctor_DoubleGaussian* locSignalFitFunctor = new DFitFunctor_DoubleGaussian(locSignalFitFunc);

			locFuncName = string(locMassHist->GetName()) + string("_Background_0");
			TF1* locBackgroundFitFunc = locMassHist->GetFunction(locFuncName.c_str());
			DFitFunctor_LinearAboveZero* locBackgroundFitFunctor = new DFitFunctor_LinearAboveZero(locBackgroundFitFunc);

			//fit control
			DFitControl* locFitControl = new DFitControl(locTotalFitFunc, nullptr, locSignalFitFunctor, locBackgroundFitFunctor);

			//compute 3sigma cut values
			double locAvgBeamEnergy = 3.25 + 0.5*locEBin; //HARDCODED: BAD!!!
			locAvgBeamEnergies[locFoundFlag].push_back(locAvgBeamEnergy);
			double locMean = locMeanFitFunc->Eval(locAvgBeamEnergy);
			double loc3Sigma = 3.0*locSigmaFitFunc->Eval(locAvgBeamEnergy);

			double locCutRangeMin = locMean - loc3Sigma;
			double locCutRangeMax = locMean + loc3Sigma;

			// calc yield trusting background in a given range
			DFitControl::DYieldData locYieldData;
			locFitControl->Calc_Yield(locMassHist, false, locCutRangeMin, locCutRangeMax, locYieldData);
			double locBackgroundUnderPeak = locBackgroundFitFunctor->dFunc->Integral(locCutRangeMin, locCutRangeMax);

			double locBinWidth = locMassHist->GetBinLowEdge(2) - locMassHist->GetBinLowEdge(1);
			double locTotalIntegral = locMassHist->Integral(1, locMassHist->GetNbinsX(), "width");
			double locSidebandIntegral = locTotalIntegral - locYieldData.dYield*locBinWidth - locBackgroundUnderPeak;
			double locSidebandScaleFraction = locBackgroundUnderPeak/locSidebandIntegral;
			locFactors[locFoundFlag].push_back(locSidebandScaleFraction);
			locFactorUncertainties[locFoundFlag].push_back(0.0);

			double locHistRange = locMassHist->GetXaxis()->GetBinUpEdge(locMassHist->GetNbinsX()) - locMassHist->GetBinLowEdge(1);
			double locSidebandRangeFraction = (locHistRange - (locCutRangeMax - locCutRangeMin))/locHistRange;
			locSidebandRangeFractions[locFoundFlag].push_back(locSidebandRangeFraction);

			++locEBin;
		}
	}

	for(auto& locMapPair : locAvgBeamEnergies)
	{
		bool locFoundFlag = locMapPair.first;

		//convert vectors to arrays
		unsigned int locArraySize = locAvgBeamEnergies[locFoundFlag].size();
		double* locXArray = Convert_VectorToArray(locAvgBeamEnergies[locFoundFlag]);
		double* locFactorArray = Convert_VectorToArray(locFactors[locFoundFlag]);
		double* locFactorUncertaintyArray = Convert_VectorToArray(locFactorUncertainties[locFoundFlag]);
		double* locSidebandRangeFractionArray = Convert_VectorToArray(locSidebandRangeFractions[locFoundFlag]);

		//create graph of factors
		ostringstream locGraphName;
		locGraphName << "SidebandSubtractionFactor_" << (locFoundFlag ? "Found" : "Missing");
		TGraph* locFactorGraph = new TGraph(locArraySize, locXArray, locFactorArray);
		locFactorGraph->SetName(locGraphName.str().c_str());
		locFactorGraph->GetXaxis()->SetTitle("Beam Energy (GeV)");
		locFactorGraph->GetYaxis()->SetTitle("Sideband Subtraction Factor");
		locFactorGraph->Write();

		//create graph of ranges
		locGraphName.str("");
		locGraphName << "SidebandRangeFraction_" << (locFoundFlag ? "Found" : "Missing");
		TGraph* locRangeGraph = new TGraph(locArraySize, locXArray, locSidebandRangeFractionArray);
		locRangeGraph->SetName(locGraphName.str().c_str());
		locRangeGraph->GetXaxis()->SetTitle("Beam Energy (GeV)");
		locRangeGraph->GetYaxis()->SetTitle("Sideband Range Fraction");
		locRangeGraph->Write();
	}

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

TF1* Fit_ToLinears(TGraphErrors* locGraph)
{
	//strategy: do following steps until the chisq/df is good.
		//first fit a flat function.
		//then fit a linear function
		//then try two linear functions ... then 3, 4, etc. until fit is good
		//the linear functions are constrained to be piecewise continuous

	Cleanup_Graph(locGraph);

	size_t locGraphN = locGraph->GetN();
	if(locGraphN == 0)
		return nullptr;
	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 locFitFunc; //yikes. we're done
	double locChiSqPerDF = locFitFunc->GetChisquare()/double(locFitFunc->GetNDF());
	if(locChiSqPerDF < locMaxChiSqPerDF)
	{
		locGraph->GetListOfFunctions()->Add(locFitFunc);
		return locFitFunc;
	}

	//nope, try piecewise-linear fits
	ostringstream locFitFuncStream;
	bool locFirstTimeFlag = true;
	size_t locNumFitPieces = 0; //incremented at beginning of loop
	vector<double> locFitParams; //initial values
	while(locChiSqPerDF > locMaxChiSqPerDF)
	{
		++locNumFitPieces; //# linear functions
//if(locNumFitPieces == 5)
//	break;


		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())
		{
			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)
				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 = 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]

			//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);

	return locFitFunc;
}

void Cleanup_Graph(TGraphErrors* locGraph)
{
	//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 locNumFitPoints = 10;
	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
		for(int loc_i = locMinFitIndex + locNumFitPoints/3; loc_i <= locMaxFitIndex - locNumFitPoints/3; ++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;
			for(int loc_j = loc_i - 2; loc_j <= loc_i + 2; ++loc_j)
			{
				double locNearbyFitY = locFitFunc->Eval(locGraphX[loc_j]);
				double locNearbyDeltaY = locGraphY[loc_j] - locNearbyFitY;
				double locNearbyNumSigma = fabs(locNearbyDeltaY)/locGraphEY[loc_i];
				if(locNearbyNumSigma < 2.0)
					++locNumCloseToFitPoints;
			}

			//if many other points close, reject this point
			if((locNumSigma > 4.0) && (locNumCloseToFitPoints > 2))
				locPointsToRemove.insert(loc_i);
			else if((locNumSigma > 3.0) && (locNumCloseToFitPoints > 3))
				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;
}

/*********************************************************** 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, 3.0))
	{
		locFitControl->Remove_FitFunctions(locHist);
		delete locFitControl;
		return nullptr;
	}

	return locFitControl;
}

bool Eval_Fit(TH1* locHist, DFitControl* locFitControl, double locNumSigmas)
{
	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);
	if((locFirstArea < 0.0) || (locSecondArea < 0.0))
		return false;
	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
	if(locWiderSigma/locNarrowerSigma > 20.0)
		return false; //ridiculous

	DFitControl::DYieldData locYieldData;
	locFitControl->Calc_Yield(locHist, true, false, locYieldData); //trust signal, yield computed from fit function
	cout << "yield = " << locYieldData.dYield << endl;
	if(locYieldData.dYield < gMinPeakYield)
		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;
}