#include "DDataConversion.h"

ClassImp(DDataConversion)

int gDebugLevel = 0;

TH1D* DDataConversion::Convert_GraphToHist(const TGraphErrors* locGraphErrors)
{
	//assumes no bins are missing!!
	unsigned int locNumGraphBins = locGraphErrors->GetN();
	double* locXArray = locGraphErrors->GetX();
	double* locYArray = locGraphErrors->GetY();
	double* locYUncertaintiesArray = locGraphErrors->GetEY();

	double* locBinEdges = new double[locNumGraphBins + 1];
	for(unsigned int loc_i = 0; loc_i < locNumGraphBins; ++loc_i)
	{
		if(loc_i == 0)
		{
			double locBinWidth = locXArray[loc_i + 1] - locXArray[loc_i];
			double locNextBinEdge = 0.5*(locXArray[loc_i] + locXArray[loc_i + 1]);
			locBinEdges[loc_i] = locNextBinEdge - locBinWidth;
			locBinEdges[loc_i + 1] = locNextBinEdge;
		}
		else if(loc_i == (locNumGraphBins - 1))
		{
			double locBinWidth = locXArray[loc_i] - locXArray[loc_i - 1];
			double locPreviousBinEdge = 0.5*(locXArray[loc_i] + locXArray[loc_i - 1]);
			locBinEdges[loc_i + 1] = locPreviousBinEdge + locBinWidth;
		}
		else
		{
			double locNextBinEdge = 0.5*(locXArray[loc_i] + locXArray[loc_i + 1]);
			locBinEdges[loc_i + 1] = locNextBinEdge;
		}
	}

	string locHistName = string(locGraphErrors->GetName()) + string("_Hist");
	string locHistTitle = string(locGraphErrors->GetTitle()) + string(locGraphErrors->GetXaxis()->GetTitle()) + string(locGraphErrors->GetYaxis()->GetTitle());
	TH1D* locHist = new TH1D(locHistName.c_str(), locHistTitle.c_str(), locNumGraphBins, locBinEdges);
	for(unsigned int loc_i = 0; loc_i < locNumGraphBins; ++loc_i)
	{
		locHist->SetBinContent(loc_i, locYArray[loc_i]);
		locHist->SetBinError(loc_i, locYUncertaintiesArray[loc_i]);
	}
	return locHist;
}

TGraphErrors* DDataConversion::Convert_HistToGraph(TH1D *locHist, bool locExcludeZeroesFlag)
{
	unsigned int locNumNonZeroBins = 0;
	ostringstream locGraphName, locGraphTitle;
	for(int loc_i = 1; loc_i <= locHist->GetNbinsX(); ++loc_i)
	{
		if(locHist->GetBinContent(loc_i) > 0.0)
			++locNumNonZeroBins;
	}

	int locNumFinalXBins = locExcludeZeroesFlag ? locNumNonZeroBins : locHist->GetNbinsX();
	if(locNumFinalXBins == 0)
		return NULL;

	double *locXData = new double[locNumFinalXBins];
	double *locYData = new double[locNumFinalXBins];
	double *locYDataError = new double[locNumFinalXBins];

	unsigned int locIndex = 0;
	for(int loc_i = 1; loc_i <= locHist->GetNbinsX(); ++loc_i)
	{
		if((locHist->GetBinContent(loc_i) > 0.0) || !locExcludeZeroesFlag)
		{
			locXData[locIndex] = locHist->GetXaxis()->GetBinCenter(loc_i);
			locYData[locIndex] = locHist->GetBinContent(loc_i);
			locYDataError[locIndex] = locHist->GetBinError(loc_i);
			++locIndex;
		}
	}

	TGraphErrors *locGraph = new TGraphErrors(locNumFinalXBins, locXData, locYData, NULL, locYDataError);
	locGraphTitle << locHist->GetTitle() << ";" << locHist->GetXaxis()->GetTitle() << ";" << locHist->GetYaxis()->GetTitle();
	locGraph->SetTitle(locGraphTitle.str().c_str());
	locGraphName << locHist->GetName() << "_Graph";
	locGraph->SetName(locGraphName.str().c_str());
	return locGraph;
}

TObjArray* DDataConversion::Convert_HistArrayToGraphArray(TObjArray *locHistArray, bool locExcludeZeroesFlag)
{
	unsigned int locNumHists = locHistArray->GetEntriesFast();
	TObjArray *locGraphErrorsArray = new TObjArray(locNumHists);
	for(unsigned int loc_i = 0; loc_i < locNumHists; ++loc_i)
	{
		TH1D* locHist = (TH1D*)(*locHistArray)[loc_i];
		TGraphErrors* locGraphErrors = Convert_HistToGraph(locHist, locExcludeZeroesFlag);
		locGraphErrorsArray->AddAt(locGraphErrors, loc_i);
	}
	return locGraphErrorsArray;
}

TGraph* DDataConversion::Convert_HistToGraph_NoErrors(TH1D *locHist, bool locExcludeZeroesFlag)
{
	int locNumXBins = locHist->GetNbinsX();
	int locNumNonZeroBins = 0;

	for(int loc_i = 1; loc_i <= locNumXBins; ++loc_i)
	{
		if(locHist->GetBinContent(loc_i) > 0.0)
			locNumNonZeroBins++;
	}

	int locNumFinalXBins = locExcludeZeroesFlag ? locNumNonZeroBins : locNumXBins;

//cout << "numfinal, flag, numnonzero, numx = " << locNumFinalXBins << ", " << locExcludeZeroesFlag << ", " << locNumNonZeroBins << ", " << locNumXBins << endl;
	if(locNumFinalXBins == 0)
		return NULL;

	double *locXData = new double[locNumFinalXBins];
	double *locYData = new double[locNumFinalXBins];

	int locIndex = 0;
	for(int loc_i = 1; loc_i <= locNumXBins; ++loc_i)
	{
		if((locHist->GetBinContent(loc_i) > 0.0) || !locExcludeZeroesFlag)
		{
			locXData[locIndex] = locHist->GetXaxis()->GetBinCenter(loc_i);
			locYData[locIndex] = locHist->GetBinContent(loc_i);
			++locIndex;
		}
	}

	TGraph *locGraph = new TGraph(locNumFinalXBins, locXData, locYData);
	string locGraphTitle = string(locHist->GetTitle()) + string(";") + string(locHist->GetXaxis()->GetTitle()) + string(";") + string(locHist->GetYaxis()->GetTitle());
	locGraph->SetTitle(locGraphTitle.c_str());
	string locGraphName = string(locHist->GetName()) + string("_Graph");
	locGraph->SetName(locGraphName.c_str());
	return locGraph;
}

TObjArray* DDataConversion::Convert_HistArrayToGraphArray_NoErrors(TObjArray *locHistArray, bool locExcludeZeroesFlag)
{
	int locNumHists = locHistArray->GetEntriesFast();
	TObjArray *locGraphArray = new TObjArray(locNumHists);
	for(int loc_i = 0; loc_i < locNumHists; ++loc_i)
	{
		TH1D* locHist = (TH1D*)(*locHistArray)[loc_i];
		TGraph* locGraph = Convert_HistToGraph_NoErrors(locHist, locExcludeZeroesFlag);
		locGraphArray->AddAt(locGraph, loc_i);
	}
	return locGraphArray;
}

TObjArray* DDataConversion::Convert_DistroToMatch_GraphErrors(TObjArray *locDistroArrayToInterpolate, TObjArray *locDistroArrayToMatchTo, double locMinX, double locMaxX)
{
	int locNumDistros = locDistroArrayToInterpolate->GetEntriesFast();
	TObjArray *locInterpolatedDistroArray = new TObjArray(locNumDistros);
	for(int loc_i = 0; loc_i < locNumDistros; ++loc_i)
	{
		TGraphErrors* locGraph_DistroToConvert = (TGraphErrors*)(*locDistroArrayToInterpolate)[loc_i];
		TGraphErrors* locGraph_DistroToMatchTo = (TGraphErrors*)(*locDistroArrayToMatchTo)[loc_i];

		TGraphErrors *locGraph_InterpolatedDistro = Convert_DistroToMatch(locGraph_DistroToConvert, locGraph_DistroToMatchTo, locMinX, locMaxX);
		locInterpolatedDistroArray->AddAt(locGraph_InterpolatedDistro, loc_i);
	}
	return locInterpolatedDistroArray;
}

TGraphErrors* DDataConversion::Convert_DistroToMatch(TGraphErrors *locGraph_DistroToConvert, TGraphErrors *locGraph_DistroToMatchTo, double locMinX, double locMaxX)
{
	//interpolate/extrapolate one distro to match another, such that they each have values for the same x-axis points
		//extrapolates if necessary up to min/max x
	if((locGraph_DistroToConvert == NULL) || (locGraph_DistroToMatchTo == NULL))
		return NULL;

	int locNumPoints_ToMatchTo = locGraph_DistroToMatchTo->GetN();
	double *locXArray_ToMatchTo = locGraph_DistroToMatchTo->GetX();

	int locNumPoints_ToInterpolate = locGraph_DistroToConvert->GetN();
	double *locXArray_ToInterpolate = locGraph_DistroToConvert->GetX();
	double *locYArray_ToInterpolate = locGraph_DistroToConvert->GetY();
	double *locYErrorArray_ToInterpolate = locGraph_DistroToConvert->GetEY();

	double locTolerance = 0.00001;
	if(locNumPoints_ToInterpolate == 1)
	{
		//only one point: return it if exact match, else bail
		double locX = locXArray_ToInterpolate[0];
		if((locX < locMinX) || (locX > locMaxX))
			return NULL;
		for(int loc_i = 0; loc_i < locNumPoints_ToMatchTo; ++loc_i)
		{
			if(fabs(locXArray_ToMatchTo[loc_i] - locX) < locTolerance)
				return locGraph_DistroToConvert;
		}
		return NULL;
	}

	int locNumPoints_Final = 0;
	int locPointOffset = 0;
//cout << "min, max = " << locMinX << ", " << locMaxX << endl;
//cout << "# points to match to: " << locNumPoints_ToMatchTo << endl;
	for(int loc_i = 0; loc_i < locNumPoints_ToMatchTo; ++loc_i)
	{
		if((locXArray_ToMatchTo[loc_i] + locTolerance) < locMinX)
		{
			++locPointOffset;
			continue;
		}
		if((locXArray_ToMatchTo[loc_i] - locTolerance) > locMaxX)
			break;
		++locNumPoints_Final;
	}

	bool locMissedLowEdgeFlag = false;
	int locIndexToCheck = locPointOffset - 1;
	if(locIndexToCheck > 0)
	{
		if(locXArray_ToMatchTo[locPointOffset - 1] < locMinX) //have skipped the low edge!
			locMissedLowEdgeFlag = true;
	}
	int locMissedHighEdgeFlag = 0;
	locIndexToCheck = locPointOffset + locNumPoints_Final - 1 + 1; //index corresponding to point after current max x point
	if(locIndexToCheck < locNumPoints_ToMatchTo)
	{
		if(locXArray_ToMatchTo[locIndexToCheck] > locMaxX) //have skipped the high edge!
			locMissedHighEdgeFlag = 1;
	}

	if(locMissedLowEdgeFlag)
		locNumPoints_Final++;
	if(locMissedHighEdgeFlag == 1)
		locNumPoints_Final++;

//cout << "# points final = " << locNumPoints_Final << endl;
	double *locXArray_Final = new double[locNumPoints_Final];
	double *locYArray_Final = new double[locNumPoints_Final];
	double *locYErrorArray_Final = new double[locNumPoints_Final];

//cout << "numpoints tomatchto, tointerp, final = " << locNumPoints_ToMatchTo << ", " << locNumPoints_ToInterpolate << ", " << locNumPoints_Final << endl;
//cout << "missed low/high edge flags = " << locMissedLowEdgeFlag << ", " << locMissedHighEdgeFlag << endl;
	for(int loc_i = 0; loc_i < locNumPoints_Final; ++loc_i)
	{
		double locX;
		if((loc_i == 0) && locMissedLowEdgeFlag)
			locX = locMinX;
		else if((loc_i == (locNumPoints_Final - 1)) && (locMissedHighEdgeFlag == 1))
			locX = locMaxX;
		else if(locMissedLowEdgeFlag)
			locX = locXArray_ToMatchTo[loc_i + locPointOffset - 1];
		else
			locX = locXArray_ToMatchTo[loc_i + locPointOffset];
		int locInterpolateDistroXLowBin = 0;
		int locInterpolateDistroXHighBin = locNumPoints_ToInterpolate - 1;
		for(int loc_j = 1; loc_j < locNumPoints_ToInterpolate; ++loc_j) //start at 1: high can't be same as low! (0)
		{
//cout << "x, j, xarray = " << locX << ", " << loc_j << ", " << locXArray_ToInterpolate[loc_j] << endl;
			if(locXArray_ToInterpolate[loc_j] > locX) //if true at j = 2: may be a low-side extrapolation: fine!
			{
				locInterpolateDistroXHighBin = loc_j;
				break;
			}
			locInterpolateDistroXLowBin = loc_j; //keep updating until true! (true when break)
		}
		if(locInterpolateDistroXLowBin == locInterpolateDistroXHighBin) //need high-side extrapolation
			locInterpolateDistroXLowBin--;

//cout << "i, low, high bins = " << loc_i << ", " << locInterpolateDistroXLowBin << ", " << locInterpolateDistroXHighBin << endl;

		double locInterpolateDistroXLow = locXArray_ToInterpolate[locInterpolateDistroXLowBin];
		double locInterpolateDistroXHigh = locXArray_ToInterpolate[locInterpolateDistroXHighBin];
		double locInterpolateDistroYLow = locYArray_ToInterpolate[locInterpolateDistroXLowBin];
		double locInterpolateDistroYHigh = locYArray_ToInterpolate[locInterpolateDistroXHighBin];
		double locInterpolateDistroErrorYLow = locYErrorArray_ToInterpolate[locInterpolateDistroXLowBin];
		double locInterpolateDistroErrorYHigh = locYErrorArray_ToInterpolate[locInterpolateDistroXHighBin];

//cout << "i, xlow, xhigh, ylow, yhigh, eylow, eyhigh = " << loc_i << ", " << locInterpolateDistroXLow << ", " << locInterpolateDistroXHigh << ", " << locInterpolateDistroYLow << ", " << locInterpolateDistroYHigh << ", " << locInterpolateDistroErrorYLow << ", " << locInterpolateDistroErrorYHigh << endl;

		double locSlope = (locInterpolateDistroYHigh - locInterpolateDistroYLow)/(locInterpolateDistroXHigh - locInterpolateDistroXLow);
		double locIntercept = locInterpolateDistroYLow - locSlope*locInterpolateDistroXLow;
		double locY = locSlope*locX + locIntercept;

		double locUncertaintyTerm = (locX - locInterpolateDistroXHigh)/(locInterpolateDistroXHigh - locInterpolateDistroXLow);
		double locErrorY = sqrt(locInterpolateDistroErrorYHigh*locInterpolateDistroErrorYHigh*(1.0 + locUncertaintyTerm)*(1.0 + locUncertaintyTerm) + locUncertaintyTerm*locUncertaintyTerm*locInterpolateDistroErrorYLow*locInterpolateDistroErrorYLow);

//cout << "i, x, y, ey = " << loc_i << ", " << locX << ", " << locY << ", " << locErrorY << endl;
		locXArray_Final[loc_i] = locX;
		locYArray_Final[loc_i] = locY;
		locYErrorArray_Final[loc_i] = locErrorY;
	}

	string locGraphTitle = locGraph_DistroToConvert->GetTitle();
	string locGraphTitle_XAxis = locGraph_DistroToConvert->GetXaxis()->GetTitle();
	string locGraphTitle_YAxis = locGraph_DistroToConvert->GetYaxis()->GetTitle();
	string locGraphName = string(locGraph_DistroToConvert->GetName()) + string("_Interpolated");

	TGraphErrors *locGraph_InterpolatedDistro = new TGraphErrors(locNumPoints_Final, locXArray_Final, locYArray_Final, NULL, locYErrorArray_Final);
	locGraph_InterpolatedDistro->SetName(locGraphName.c_str());
	locGraph_InterpolatedDistro->SetTitle(locGraphTitle.c_str());
	locGraph_InterpolatedDistro->GetXaxis()->SetTitle(locGraphTitle_XAxis.c_str());
	locGraph_InterpolatedDistro->GetYaxis()->SetTitle(locGraphTitle_YAxis.c_str());

	return locGraph_InterpolatedDistro;
}

TGraphErrors* DDataConversion::Reverse_PointOrder(TGraphErrors *locInputGraph)
{
	int locNumPoints = locInputGraph->GetN();
	double *locX_Input = locInputGraph->GetX();
	double *locY_Input = locInputGraph->GetY();
	double *locEY_Input = locInputGraph->GetEY();

	double *locX_Output = new double[locNumPoints];
	for(int loc_i = 0; loc_i < locNumPoints; ++loc_i)
		locX_Output[loc_i] = locX_Input[locNumPoints - 1 - loc_i];

	double *locY_Output = new double[locNumPoints];
	for(int loc_i = 0; loc_i < locNumPoints; ++loc_i)
		locY_Output[loc_i] = locY_Input[locNumPoints - 1 - loc_i];
	double *locEY_Output = new double[locNumPoints];
	for(int loc_i = 0; loc_i < locNumPoints; ++loc_i)
		locEY_Output[loc_i] = locEY_Input[locNumPoints - 1 - loc_i];

	TGraphErrors *locOutputGraph = new TGraphErrors(locNumPoints, locX_Output, locY_Output, NULL, locEY_Output);
	string locGraphName = string(locInputGraph->GetName()) + string("_Reversed");
	locOutputGraph->SetName(locGraphName.c_str());
	string locGraphTitle = string(locInputGraph->GetTitle()) + string(";") + string(locInputGraph->GetXaxis()->GetTitle()) + string(";") + string(locInputGraph->GetYaxis()->GetTitle());
	locOutputGraph->SetTitle(locGraphTitle.c_str());
	return locOutputGraph;
}

TObjArray* DDataConversion::Convert_DistroArrayToMatch_Hists(TObjArray *locDistroArrayToInterpolate, TObjArray *locDistroArrayToMatchTo)
{
	int locNumDistros = locDistroArrayToInterpolate->GetEntriesFast();
	TObjArray *locInterpolatedDistroArray = new TObjArray(locNumDistros);
	for(int loc_i = 0; loc_i < locNumDistros; ++loc_i)
	{
		TH1D* locHist_DistroToInterpolate = (TH1D*)(*locDistroArrayToInterpolate)[loc_i];
		TH1D* locHist_DistroToMatchTo = (TH1D*)(*locDistroArrayToMatchTo)[loc_i];
		TH1D* locHist_InterpolatedDistro = NULL;
		if((locHist_DistroToInterpolate != NULL) && (locHist_DistroToMatchTo != NULL))
			locHist_InterpolatedDistro = Convert_DistroToMatch(locHist_DistroToInterpolate, locHist_DistroToMatchTo);
		else
			locHist_InterpolatedDistro = NULL;
		locInterpolatedDistroArray->AddAt(locHist_InterpolatedDistro, loc_i);
	}
	return locInterpolatedDistroArray;
}

TH1D* DDataConversion::Convert_DistroToMatch(TH1D *locHist_DistroToInterpolate, TH1D *locHist_DistroToMatchTo)
{
	//also extrapolates if necessary!!!!!!!
	if((locHist_DistroToInterpolate == NULL) || (locHist_DistroToMatchTo == NULL))
		return NULL;

	string locHistTitle = locHist_DistroToInterpolate->GetTitle();
	string locHistTitle_XAxis = locHist_DistroToInterpolate->GetXaxis()->GetTitle();
	string locHistTitle_YAxis = locHist_DistroToInterpolate->GetYaxis()->GetTitle();
	string locHistName = string(locHist_DistroToInterpolate->GetName()) + string("_Interpolated");

	TH1D* locHist_InterpolatedDistro = (TH1D*)locHist_DistroToMatchTo->Clone(locHistName.c_str());
	locHist_InterpolatedDistro->Reset();
	locHist_InterpolatedDistro->SetTitle(locHistTitle.c_str());
	locHist_InterpolatedDistro->GetXaxis()->SetTitle(locHistTitle_XAxis.c_str());
	locHist_InterpolatedDistro->GetYaxis()->SetTitle(locHistTitle_YAxis.c_str());

	for(int loc_i = 1; loc_i <= locHist_InterpolatedDistro->GetNbinsX(); ++loc_i)
	{
		double locX = locHist_InterpolatedDistro->GetBinCenter(loc_i);
		int locInterpolateDistroXLowBin = 1;
		int locInterpolateDistroXHighBin = 2;
		for(int loc_j = 2; loc_j <= locHist_DistroToInterpolate->GetNbinsX(); ++loc_j) //start at 2: high can't be same as low! (1)
		{
			if(locHist_DistroToInterpolate->GetBinCenter(loc_j) > locX) //if true at j = 2: may be a low-side extrapolation: fine!
			{
				locInterpolateDistroXHighBin = loc_j;
				break;
			}
			locInterpolateDistroXLowBin = loc_j; //keep updating until true! (true when break)
		}
		if(locInterpolateDistroXLowBin == locInterpolateDistroXHighBin) //need high-side extrapolation
			locInterpolateDistroXLowBin--;

//cout << "low, high bins = " << locInterpolateDistroXLowBin << ", " << locInterpolateDistroXHighBin << endl;

		double locInterpolateDistroXLow = locHist_DistroToInterpolate->GetBinCenter(locInterpolateDistroXLowBin);
		double locInterpolateDistroXHigh = locHist_DistroToInterpolate->GetBinCenter(locInterpolateDistroXHighBin);
		double locInterpolateDistroYLow = locHist_DistroToInterpolate->GetBinContent(locInterpolateDistroXLowBin);
		double locInterpolateDistroYHigh = locHist_DistroToInterpolate->GetBinContent(locInterpolateDistroXHighBin);
		double locInterpolateDistroErrorYLow = locHist_DistroToInterpolate->GetBinError(locInterpolateDistroXLowBin);
		double locInterpolateDistroErrorYHigh = locHist_DistroToInterpolate->GetBinError(locInterpolateDistroXHighBin);
//cout << "i, xlow, xhigh, ylow, yhigh = " << loc_i << ", " << locInterpolateDistroXLow << ", " << locInterpolateDistroXHigh << ", " << locInterpolateDistroYLow << ", " << locInterpolateDistroYHigh << endl;

		double locSlope = (locInterpolateDistroYHigh - locInterpolateDistroYLow)/(locInterpolateDistroXHigh - locInterpolateDistroXLow);
		double locIntercept = locInterpolateDistroYLow - locSlope*locInterpolateDistroXLow;
		double locY = locSlope*locX + locIntercept;

		double locUncertaintyTerm = (locX - locInterpolateDistroXHigh)/(locInterpolateDistroXHigh - locInterpolateDistroXLow);
		double locErrorY = sqrt(locInterpolateDistroErrorYHigh*locInterpolateDistroErrorYHigh*(1.0 + locUncertaintyTerm)*(1.0 + locUncertaintyTerm) + locUncertaintyTerm*locUncertaintyTerm*locInterpolateDistroErrorYLow*locInterpolateDistroErrorYLow);

//cout << "i, x, y, ey = " << loc_i << ", " << locX << ", " << locY << ", " << locErrorY << endl;

		locHist_InterpolatedDistro->SetBinContent(loc_i, locY);
		locHist_InterpolatedDistro->SetBinError(loc_i, locErrorY);
	}

	return locHist_InterpolatedDistro;
}

void DDataConversion::Interpolate_DistroToPoint(TH1D* locHist, double locPoint, double &locValue, double &locUncertainty, bool locPrintFlag)
{
	//find the value at a given point in x
	if(locPoint <= locHist->GetBinCenter(1))
	{
		locValue = locHist->GetBinContent(1);
		locUncertainty = locHist->GetBinError(1);
		if(locPrintFlag)
			cout << "Point, value, uncertainty = " << locPoint << ", " << locValue << ", " << locUncertainty << endl; //y = mx + b
		return;
	}

	int locNumBins = locHist->GetNbinsX();
	if(locPoint >= locHist->GetBinCenter(locNumBins))
	{
		locValue = locHist->GetBinContent(locNumBins);
		locUncertainty = locHist->GetBinError(locNumBins);
		if(locPrintFlag)
			cout << "Point, value, uncertainty = " << locPoint << ", " << locValue << ", " << locUncertainty << endl; //y = mx + b
		return;
	}

	for(int loc_i = 2; loc_i <= locNumBins; ++loc_i) //don't start at beginning: will break (would extrapolate, not interpolate!)
	{
		double locCurrentBinCenter = locHist->GetBinCenter(loc_i);
		if(!(fabs(locCurrentBinCenter - locPoint) > 0.0))
		{
			locValue = locHist->GetBinContent(loc_i);
			locUncertainty = locHist->GetBinError(loc_i);
			if(locPrintFlag)
				cout << "Point, value, uncertainty = " << locPoint << ", " << locValue << ", " << locUncertainty << endl; //y = mx + b
			return;
		}

		double locPreviousBinCenter = locHist->GetBinCenter(loc_i - 1);
		if((locCurrentBinCenter > locPoint) && (locPreviousBinCenter < locPoint)) //interpolate!
		{
			double locSlope = (locHist->GetBinContent(loc_i) - locHist->GetBinContent(loc_i - 1))/(locHist->GetBinCenter(loc_i) - locHist->GetBinCenter(loc_i - 1));
			double locIntercept = locHist->GetBinContent(loc_i) - locSlope*locHist->GetBinCenter(loc_i); //b = y - mx
			locValue = locSlope*locPoint + locIntercept;

			double locUncertaintyTerm = (locPoint - locHist->GetBinCenter(loc_i))/(locHist->GetBinCenter(loc_i) - locHist->GetBinCenter(loc_i - 1));
			locUncertainty = sqrt(locHist->GetBinError(loc_i)*locHist->GetBinError(loc_i)*(1.0 + locUncertaintyTerm)*(1.0 + locUncertaintyTerm) + locUncertaintyTerm*locUncertaintyTerm*locHist->GetBinError(loc_i - 1)*locHist->GetBinError(loc_i - 1));
			if(locPrintFlag)
				cout << "Point, value, uncertainty = " << locPoint << ", " << locValue << ", " << locUncertainty << endl; //y = mx + b
			return;
		}
	}
}

TObjArray* DDataConversion::Offset_Points_HistArray(TObjArray *locHistArray, double locOffset)
{
	int locNumHists = locHistArray->GetEntriesFast();
	TObjArray* locOffsetHistArray = new TObjArray(locNumHists);
	for(int loc_i = 0; loc_i < locNumHists; ++loc_i)
	{
		TH1D* locHist = (TH1D*)(*locHistArray)[loc_i];
		TH1D* locOffsetHist = Offset_Points(locHist, locOffset);
		locOffsetHistArray->AddAt(locOffsetHist, loc_i);
	}
	return locOffsetHistArray;
}

TH1D* DDataConversion::Offset_Points(TH1D *locHist, double locOffset)
{
	int locNumBins = locHist->GetNbinsX();
	double *locBinEdgeArray = new double[locNumBins + 1];
	for(int loc_i = 0; loc_i < locNumBins; ++loc_i)
		locBinEdgeArray[loc_i] = locHist->GetXaxis()->GetBinLowEdge(loc_i + 1) + locOffset;
	locBinEdgeArray[locNumBins] = locHist->GetXaxis()->GetBinUpEdge(locNumBins);

	ostringstream locHistName, locHistTitle;
	locHistName << locHist->GetName() << "_Offset";
	locHistTitle << locHist->GetTitle() << ";" << locHist->GetXaxis()->GetTitle() << ";" << locHist->GetYaxis()->GetTitle() << ";" << locHist->GetZaxis()->GetTitle();
	TH1D *locOffsetHist = new TH1D(locHistName.str().c_str(), locHistTitle.str().c_str(), locNumBins, locBinEdgeArray);

	//set values, errors
	for(int loc_i = 1; loc_i <= locNumBins; ++loc_i)
	{
		locOffsetHist->SetBinContent(loc_i, locHist->GetBinContent(loc_i));
		locOffsetHist->SetBinError(loc_i, locHist->GetBinError(loc_i));
	}
	return locOffsetHist;
}

void DDataConversion::Calculate_AverageDeviationPercentageArray_GraphErrors(TObjArray *locDistroArrayToCompare, TObjArray *locDistroArrayToCompareTo, double locRangeMin, double locRangeMax)
{
	double locWeightedMeanNumerator_Total = 0.0, locWeightedMeanDenominator_Total = 0.0;
	for(int loc_i = 0; loc_i < locDistroArrayToCompare->GetEntriesFast(); ++loc_i)
	{
		TGraphErrors* locGraph_DistroToCompare = (TGraphErrors*)(*locDistroArrayToCompare)[loc_i];
		TGraphErrors* locGraph_DistroToCompareTo = (TGraphErrors*)(*locDistroArrayToCompareTo)[loc_i];
		if((locGraph_DistroToCompare != NULL) && (locGraph_DistroToCompareTo != NULL))
			Calculate_AverageDeviationPercentage(locGraph_DistroToCompare, locGraph_DistroToCompareTo, locWeightedMeanNumerator_Total, locWeightedMeanDenominator_Total, locRangeMin, locRangeMax);
	}
}

double DDataConversion::Calculate_AverageDeviationPercentage(TGraphErrors* locGraph_DistroToCompare, TGraphErrors* locGraph_DistroToCompareTo, double& locWeightedMeanNumerator_Total, double& locWeightedMeanDenominator_Total, double locRangeMin, double locRangeMax)
{
	if((locGraph_DistroToCompare == NULL) || (locGraph_DistroToCompareTo == NULL))
		return -100.0;

	int locNumPoints = locGraph_DistroToCompare->GetN();
	if(locGraph_DistroToCompareTo->GetN() != locNumPoints)
		return -100.0;

	double *locXArray = locGraph_DistroToCompare->GetX();

	double *locYArray_ToCompareTo = locGraph_DistroToCompareTo->GetY();
	double *locYErrorArray_ToCompareTo = locGraph_DistroToCompareTo->GetEY();

	double *locYArray_ToCompare = locGraph_DistroToCompare->GetY();
	double *locYErrorArray_ToCompare = locGraph_DistroToCompare->GetEY();

	double locWeightedMeanNumerator = 0.0, locWeightedMeanDenominator = 0.0;
	int locNumPointsCompared = 0;
	for(int loc_i = 0; loc_i < locNumPoints; ++loc_i)
	{
		double locX = locXArray[loc_i];
		if((locX < locRangeMin) || (locX > locRangeMax))
			continue; //don't include!!!
		++locNumPointsCompared;
		double locValue1 = locYArray_ToCompare[loc_i];
		double locValue2 = locYArray_ToCompareTo[loc_i];
		double locError1 = locYErrorArray_ToCompare[loc_i];
		double locError2 = locYErrorArray_ToCompareTo[loc_i];
		double locDifference = (locValue2 - locValue1)/locValue2; //= 1 - Val1/Val2
		double locUncertaintyTerm = 1.0/(locValue2*locValue2);
		double locDifferenceVariance = locError1*locError1*locUncertaintyTerm + locValue1*locValue1*locError2*locError2*locUncertaintyTerm*locUncertaintyTerm;
		locWeightedMeanNumerator += locDifference/locDifferenceVariance;
		locWeightedMeanDenominator += 1.0/locDifferenceVariance;
		locWeightedMeanNumerator_Total += locDifference/locDifferenceVariance;
		locWeightedMeanDenominator_Total += 1.0/locDifferenceVariance;
	}
//cout << "numbins compared, range min/max = " << locNumPointsCompared << ", " << locRangeMin << ", " << locRangeMax << endl;
	double locWeightedMean = (locWeightedMeanDenominator > 0.0) ? 100.0*locWeightedMeanNumerator/locWeightedMeanDenominator : -100.0;
	cout << "weighted mean % difference = " << locWeightedMean << endl;
	double locWeightedMean_Total = (locWeightedMeanDenominator_Total > 0.0) ? 100.0*locWeightedMeanNumerator_Total/locWeightedMeanDenominator_Total : -100.0;
	cout << "total weighted mean % difference = " << locWeightedMean_Total << endl;
	return locWeightedMean; //if > 0: Cross1 is x% < Cross2 (% of Cross2)
}

void DDataConversion::Calculate_AverageDeviationPercentageArray_Hists(TObjArray *locDistroArrayToCompare, TObjArray *locDistroArrayToCompareTo, double locRangeMin, double locRangeMax)
{
	double locWeightedMeanNumerator_Total = 0.0, locWeightedMeanDenominator_Total = 0.0;
	for(int loc_i = 0; loc_i < locDistroArrayToCompare->GetEntriesFast(); ++loc_i)
	{
		TH1D* locHist_DistroToCompare = (TH1D*)(*locDistroArrayToCompare)[loc_i];
		TH1D* locHist_DistroToCompareTo = (TH1D*)(*locDistroArrayToCompareTo)[loc_i];
		if((locHist_DistroToCompare != NULL) && (locHist_DistroToCompareTo != NULL))
			Calculate_AverageDeviationPercentage(locHist_DistroToCompare, locHist_DistroToCompareTo, locWeightedMeanNumerator_Total, locWeightedMeanDenominator_Total, locRangeMin, locRangeMax);
	}
}

double DDataConversion::Calculate_AverageDeviationPercentage(TH1D* locHist_DistroToCompare, TH1D* locHist_DistroToCompareTo, double& locWeightedMeanNumerator_Total, double& locWeightedMeanDenominator_Total, double locRangeMin, double locRangeMax)
{
	if((locHist_DistroToCompare == NULL) || (locHist_DistroToCompareTo == NULL))
		return -100.0;

	int locNumBins = locHist_DistroToCompare->GetNbinsX();
	if(locHist_DistroToCompareTo->GetNbinsX() != locNumBins)
		return -100.0;

	double locWeightedMeanNumerator = 0.0, locWeightedMeanDenominator = 0.0;
	for(int loc_i = 1; loc_i <= locNumBins; ++loc_i)
	{
		double locX = locHist_DistroToCompare->GetBinCenter(loc_i);
		if((locX < locRangeMin) || (locX > locRangeMax))
			continue; //don't include!!!
		double locValue1 = locHist_DistroToCompare->GetBinContent(loc_i);
		double locValue2 = locHist_DistroToCompareTo->GetBinContent(loc_i);
		double locError1 = locHist_DistroToCompare->GetBinError(loc_i);
		double locError2 = locHist_DistroToCompareTo->GetBinError(loc_i);
		double locDifference = (locValue2 - locValue1)/locValue2; //= 1 - Val1/Val2
		double locUncertaintyTerm = 1.0/(locValue2*locValue2);
		double locDifferenceVariance = locError1*locError1*locUncertaintyTerm + locValue1*locValue1*locError2*locError2*locUncertaintyTerm*locUncertaintyTerm;
		locWeightedMeanNumerator += locDifference/locDifferenceVariance;
		locWeightedMeanDenominator += 1.0/locDifferenceVariance;
		locWeightedMeanNumerator_Total += locDifference/locDifferenceVariance;
		locWeightedMeanDenominator_Total += 1.0/locDifferenceVariance;
	}
	double locWeightedMean = (locWeightedMeanDenominator > 0.0) ? 100.0*locWeightedMeanNumerator/locWeightedMeanDenominator : -100.0;
	cout << "weighted mean % difference = " << locWeightedMean << endl;
	double locWeightedMean_Total = (locWeightedMeanDenominator_Total > 0.0) ? 100.0*locWeightedMeanNumerator_Total/locWeightedMeanDenominator_Total : -100.0;
	cout << "Total weighted mean % difference = " << locWeightedMean_Total << endl;
	return locWeightedMean;
}

DCrossSectionData* DDataConversion::Convert_CrossSection_WithGaps(DCrossSectionData* locCrossSectionData_ToConvert, DCrossSectionData* locCrossSectionData_ConvertTo)
{
	vector<double>& locFirstDimValues_ConvertTo = locCrossSectionData_ConvertTo->dFirstDimValues;
	vector<double>& locFirstDimValues_ToConvert = locCrossSectionData_ToConvert->dFirstDimValues;

	TObjArray* locNewCrossSectionArray = new TObjArray();

	int locNumCreated = 0;
	int locLowBinIndex_ConvertTo = -1;
	for(size_t loc_i = 0; loc_i < locFirstDimValues_ToConvert.size(); ++loc_i)
	{
		//search first-dim bins of to-convert cross section
		double locToConvertValue = locFirstDimValues_ToConvert[loc_i];
if(gDebugLevel > 0)
cout << "to convert: " << locToConvertValue << endl;

		//find bins just on either side //resume search at previous low index
		bool locExactMatchFlag = false;
		bool locMatchFoundFlag = false;
		for(size_t loc_j = locLowBinIndex_ConvertTo + 1; loc_j < locFirstDimValues_ConvertTo.size(); ++loc_j) //find the target
		{
			double locConvertToValue = locFirstDimValues_ConvertTo[loc_j];
if(gDebugLevel > 0)
cout << "j, convert_to_E, to-convert: " << loc_j << ", " << locConvertToValue << ", " << locToConvertValue << endl;

			//check if exact match
			if(fabs(locConvertToValue - locToConvertValue) < 0.0001)
			{
if(gDebugLevel > 0)
cout << "exact match" << endl;
				locMatchFoundFlag = true;
				locExactMatchFlag = true;
				locLowBinIndex_ConvertTo = loc_j;
				break;
			}

			if((loc_i == 0) && (locConvertToValue < locToConvertValue))
				continue; //would have to extrapolate below the range: nope

			if(loc_j < (locFirstDimValues_ConvertTo.size() - 1))
			{
				double locNextConvertToValue = locFirstDimValues_ConvertTo[loc_j + 1];
if(gDebugLevel > 0)
	cout << "next convert-to: " << locNextConvertToValue << endl;
				if(fabs(locNextConvertToValue - locToConvertValue) < fabs(locToConvertValue - locConvertToValue))
					continue; //closer to the next one: continue
			}

			//convert-to is now the closest one to to-convert
			//but, is there another to-convert that is closer to convert-to?
			while((loc_i + 1) < locFirstDimValues_ToConvert.size())
			{
if(gDebugLevel > 0)
cout << "next toconvert = " << locFirstDimValues_ToConvert[loc_i + 1] << endl;
				double locNextToConvertValue = locFirstDimValues_ToConvert[loc_i + 1];
				double locNextDistance = fabs(locNextToConvertValue - locConvertToValue);
				double locThisDistance = fabs(locFirstDimValues_ToConvert[loc_i] - locConvertToValue);
if(gDebugLevel > 0)
cout << "this distance, next distance = " << locThisDistance << ", " << locNextDistance << endl;
				if(locThisDistance < locNextDistance)
					break; //nope, on the closest one

				//yes, there is. but, before i skip this one, is the next to-convert closer to the NEXT convert-to instead of this one?
				double locNextConvertToValue = locFirstDimValues_ConvertTo[loc_j + 1];
				double locNextNextDistance = fabs(locNextConvertToValue - locNextToConvertValue);
if(gDebugLevel > 0)
cout << "next distance, next-next distance = " << locNextDistance << ", " << locNextNextDistance << endl;
				if(locNextNextDistance < locNextDistance)
					break; //the next to-convert is closer to the NEXT convert-to, don't use it

if(gDebugLevel > 0)
cout << "skip and move on" << endl;
				++loc_i; //The next one is closer, skip this one and move on. 
			}

			//update, in case it has changed
			locToConvertValue = locFirstDimValues_ToConvert[loc_i];
if(gDebugLevel > 0)
cout << "updated to-convert: " << locToConvertValue << endl;

			//check if exact match
			if(fabs(locConvertToValue - locToConvertValue) < 0.0001)
			{
if(gDebugLevel > 0)
cout << "exact match" << endl;
				locMatchFoundFlag = true;
				locExactMatchFlag = true;
				locLowBinIndex_ConvertTo = loc_j;
				break;
			}

			if(((loc_i + 1) == locFirstDimValues_ToConvert.size()) && (locToConvertValue < locConvertToValue))
				break; //at the end, we're done

			locMatchFoundFlag = true;
			locLowBinIndex_ConvertTo = loc_j;
if(gDebugLevel > 0)
cout << "bins found: " << locLowBinIndex_ConvertTo << ", " << locLowBinIndex_ConvertTo + 1 << endl;
			break;
		}

		if(!locMatchFoundFlag)
			break; //we are completely done

		if(locExactMatchFlag)
		{
			//fill in up to this point
			while(locNewCrossSectionArray->GetEntriesFast() < locLowBinIndex_ConvertTo)
				locNewCrossSectionArray->AddLast(NULL);

			//copy over
			locNewCrossSectionArray->AddLast(locCrossSectionData_ToConvert->dCrossSectionArray->At(loc_i));
			++locNumCreated;
			continue;
		}

		int locTargetIndex = locLowBinIndex_ConvertTo;
		double locTargetEnergy = locFirstDimValues_ConvertTo[locLowBinIndex_ConvertTo];

		//fill in up to this point
		while(locNewCrossSectionArray->GetEntriesFast() < locTargetIndex)
			locNewCrossSectionArray->AddLast(NULL);


		//ok, what is the target energy in relation to the to-convert energy?
		int locLowBinIndex_ToConvert = loc_i;
		if(loc_i == 0)
			locLowBinIndex_ToConvert = 0;
		else if(loc_i == (locFirstDimValues_ToConvert.size() - 1))
			locLowBinIndex_ToConvert = loc_i - 1;
		else if(locTargetEnergy > locToConvertValue)
			locLowBinIndex_ToConvert = loc_i;
		else
			locLowBinIndex_ToConvert = loc_i - 1;

		double locToConvertValue_Low = locFirstDimValues_ToConvert[locLowBinIndex_ToConvert];
		double locToConvertValue_High = locFirstDimValues_ToConvert[locLowBinIndex_ToConvert + 1];

		TGraphErrors* locGraphErrors_Low = (TGraphErrors*)locCrossSectionData_ToConvert->dCrossSectionArray->At(locLowBinIndex_ToConvert);
		TGraphErrors* locGraphErrors_High = (TGraphErrors*)locCrossSectionData_ToConvert->dCrossSectionArray->At(locLowBinIndex_ToConvert + 1);

		TGraphErrors* locGraphErrors_Converted = Interpolate_CrossSection_ToOtherBin(locGraphErrors_Low, locGraphErrors_High, locToConvertValue_Low, locToConvertValue_High, locTargetEnergy);
		locNewCrossSectionArray->AddLast(locGraphErrors_Converted);
		if(locGraphErrors_Converted != NULL)
			++locNumCreated;
	}

	if(locNumCreated == 0)
		return NULL;

	string locGraphName = locCrossSectionData_ToConvert->dGraphName + string("_Converted");
	return (new DCrossSectionData(locCrossSectionData_ToConvert, locGraphName, locNewCrossSectionArray, locFirstDimValues_ConvertTo));
}

DCrossSectionData* DDataConversion::Convert_CrossSection_NearbyOnly(DCrossSectionData* locCrossSectionData_ToConvert, DCrossSectionData* locCrossSectionData_ConvertTo, double locMaxDistance)
{
	vector<double>& locFirstDimValues_ConvertTo = locCrossSectionData_ConvertTo->dFirstDimValues;
	vector<double>& locFirstDimValues_ToConvert = locCrossSectionData_ToConvert->dFirstDimValues;

	TObjArray* locNewCrossSectionArray = new TObjArray();

	int locNumCreated = 0;
	int locLowBinIndex_ConvertTo = -1;
	for(size_t loc_i = 0; loc_i < locFirstDimValues_ToConvert.size(); ++loc_i)
	{
		//search first-dim bins of to-convert cross section
		double locToConvertValue = locFirstDimValues_ToConvert[loc_i];
if(gDebugLevel > 0)
cout << "to convert: " << locToConvertValue << endl;

		//find bins just on either side //resume search at previous low index
		bool locExactMatchFlag = false;
		bool locMatchFoundFlag = false;
		for(size_t loc_j = locLowBinIndex_ConvertTo + 1; loc_j < locFirstDimValues_ConvertTo.size(); ++loc_j) //find the target
		{
			double locConvertToValue = locFirstDimValues_ConvertTo[loc_j];
if(gDebugLevel > 0)
cout << "j, convert_to_E, to-convert: " << loc_j << ", " << locConvertToValue << ", " << locToConvertValue << endl;

			//check if exact match
			if(fabs(locConvertToValue - locToConvertValue) < 0.0001)
			{
if(gDebugLevel > 0)
cout << "exact match" << endl;
				locMatchFoundFlag = true;
				locExactMatchFlag = true;
				locLowBinIndex_ConvertTo = loc_j;
				break;
			}

			if((loc_i == 0) && (locConvertToValue < locToConvertValue))
				continue; //would have to extrapolate below the range: nope

			if(loc_j < (locFirstDimValues_ConvertTo.size() - 1))
			{
				double locNextConvertToValue = locFirstDimValues_ConvertTo[loc_j + 1];
if(gDebugLevel > 0)
	cout << "next convert-to: " << locNextConvertToValue << endl;
				if(fabs(locNextConvertToValue - locToConvertValue) < fabs(locToConvertValue - locConvertToValue))
					continue; //closer to the next one: continue
			}

			//convert-to is now the closest one to to-convert
			//but, is there another to-convert that is closer to convert-to?
			while((loc_i + 1) < locFirstDimValues_ToConvert.size())
			{
if(gDebugLevel > 0)
cout << "next toconvert = " << locFirstDimValues_ToConvert[loc_i + 1] << endl;
				double locNextToConvertValue = locFirstDimValues_ToConvert[loc_i + 1];
				double locNextDistance = fabs(locNextToConvertValue - locConvertToValue);
				double locThisDistance = fabs(locFirstDimValues_ToConvert[loc_i] - locConvertToValue);
if(gDebugLevel > 0)
cout << "this distance, next distance = " << locThisDistance << ", " << locNextDistance << endl;
				if(locThisDistance < locNextDistance)
					break; //nope, on the closest one

				//yes, there is. but, before i skip this one, is the next to-convert closer to the NEXT convert-to instead of this one?
				double locNextConvertToValue = locFirstDimValues_ConvertTo[loc_j + 1];
				double locNextNextDistance = fabs(locNextConvertToValue - locNextToConvertValue);
if(gDebugLevel > 0)
cout << "next distance, next-next distance = " << locNextDistance << ", " << locNextNextDistance << endl;
				if(locNextNextDistance < locNextDistance)
					break; //the next to-convert is closer to the NEXT convert-to, don't use it

if(gDebugLevel > 0)
cout << "skip and move on" << endl;
				++loc_i; //The next one is closer, skip this one and move on. 
			}

			//update, in case it has changed
			locToConvertValue = locFirstDimValues_ToConvert[loc_i];
if(gDebugLevel > 0)
cout << "updated to-convert: " << locToConvertValue << endl;

			//check if exact match
			if(fabs(locConvertToValue - locToConvertValue) < (locMaxDistance + 0.0001))
			{
if(gDebugLevel > 0)
cout << "exact match" << endl;
				locMatchFoundFlag = true;
				locExactMatchFlag = true;
				locLowBinIndex_ConvertTo = loc_j;
				break;
			}
		}

		if(locExactMatchFlag)
		{
			//fill in up to this point
			while(locNewCrossSectionArray->GetEntriesFast() < locLowBinIndex_ConvertTo)
				locNewCrossSectionArray->AddLast(NULL);

			//copy over
			locNewCrossSectionArray->AddLast(locCrossSectionData_ToConvert->dCrossSectionArray->At(loc_i));
			++locNumCreated;
			continue;
		}
	}

	if(locNumCreated == 0)
		return NULL;

	string locGraphName = locCrossSectionData_ToConvert->dGraphName + string("_Converted");
	return (new DCrossSectionData(locCrossSectionData_ToConvert, locGraphName, locNewCrossSectionArray, locFirstDimValues_ConvertTo));
}

DCrossSectionData* DDataConversion::Convert_CrossSection_FillGaps(DCrossSectionData* locCrossSectionData_ToConvert, DCrossSectionData* locCrossSectionData_ConvertTo)
{
	vector<double>& locFirstDimValues_ConvertTo = locCrossSectionData_ConvertTo->dFirstDimValues;
	vector<double>& locFirstDimValues_ToConvert = locCrossSectionData_ToConvert->dFirstDimValues;

	TObjArray* locNewCrossSectionArray = new TObjArray();

	int locNumCreated = 0;
	size_t locBinIndex_Low = 0, locBinIndex_High = 1, locBinIndex_PreviousLow = 99999, locBinIndex_PreviousHigh = 99999;
	for(size_t loc_i = 0; loc_i < locFirstDimValues_ConvertTo.size(); ++loc_i)
	{
		//search first-dim bins of to-convert cross section
		double locFirstDimValue = locFirstDimValues_ConvertTo[loc_i];
if(gDebugLevel > 0)
cout << "convert to: " << locFirstDimValue << endl;
		if((locFirstDimValue < locFirstDimValues_ToConvert.front()) || (locFirstDimValue > locFirstDimValues_ToConvert.back()))
		{
			locNewCrossSectionArray->AddLast(NULL);
			continue; //would have to extrapolate
		}
if(gDebugLevel > 0)
cout << "find bins" << endl;

		//find bins just on either side //resume search at previous low index
		bool locExactMatchFlag = false;
		for(size_t loc_j = locBinIndex_Low; loc_j < locFirstDimValues_ToConvert.size(); ++loc_j)
		{
if(gDebugLevel > 0)
cout << "j, to_convert_E, first dim value: " << loc_j << ", " << locFirstDimValues_ToConvert[loc_j] << ", " << locFirstDimValue << endl;
			if(locFirstDimValues_ToConvert[loc_j] < locFirstDimValue)
				continue;
			locBinIndex_Low = loc_j - 1;
			locBinIndex_High = loc_j;
if(gDebugLevel > 0)
cout << "bins found: " << locBinIndex_Low << ", " << locBinIndex_High << endl;
if(gDebugLevel > 0)
cout << "compare: " << locFirstDimValues_ToConvert[loc_j] << ", " << locFirstDimValue << ", " << endl;
			if(fabs(locFirstDimValues_ToConvert[loc_j] - locFirstDimValue) < 0.001)
				locExactMatchFlag = true;
if(gDebugLevel > 0)
cout << "exact match flag = " << locExactMatchFlag << endl;
			break;
		}

		if(locExactMatchFlag)
		{
if(gDebugLevel > 0)
cout << "exact match" << endl;
			locBinIndex_Low = locBinIndex_High;
			locBinIndex_PreviousLow = locBinIndex_High;
			locBinIndex_PreviousHigh = locBinIndex_High + 1;
			++locNumCreated;
			locNewCrossSectionArray->AddLast(locCrossSectionData_ToConvert->dCrossSectionArray->At(locBinIndex_High));
			continue;
		}

		double locFirstDimValue_Low = locFirstDimValues_ToConvert[locBinIndex_Low];
		double locFirstDimValue_High = locFirstDimValues_ToConvert[locBinIndex_High];

if(gDebugLevel > 0)
cout << "bin index low, high = " << locBinIndex_Low << ", " << locBinIndex_High << endl;
		TGraphErrors* locGraphErrors_Low = (TGraphErrors*)locCrossSectionData_ToConvert->dCrossSectionArray->At(locBinIndex_Low);
		TGraphErrors* locGraphErrors_High = (TGraphErrors*)locCrossSectionData_ToConvert->dCrossSectionArray->At(locBinIndex_High);

		TGraphErrors* locGraphErrors_Converted = Interpolate_CrossSection_ToOtherBin(locGraphErrors_Low, locGraphErrors_High, locFirstDimValue_Low, locFirstDimValue_High, locFirstDimValue);
		locNewCrossSectionArray->AddLast(locGraphErrors_Converted);
		if(locGraphErrors_Converted != NULL)
			++locNumCreated;

		locBinIndex_PreviousLow = locBinIndex_Low;
		locBinIndex_PreviousHigh = locBinIndex_High;
	}

	if(locNumCreated == 0)
		return NULL;

	string locGraphName = locCrossSectionData_ToConvert->dGraphName + string("_Converted");
	return (new DCrossSectionData(locCrossSectionData_ToConvert, locGraphName, locNewCrossSectionArray, locFirstDimValues_ConvertTo));
}

TGraphErrors* DDataConversion::Interpolate_CrossSection_ToOtherBin(TGraphErrors* locGraphErrors_Low, TGraphErrors* locGraphErrors_High, double locFirstDimValue_Low, double locFirstDimValue_High, double locFirstDimValue_Target)
{
if(gDebugLevel > 0)
cout << "energies low/high/target: " << locFirstDimValue_Low << ", " << locFirstDimValue_High << ", " << locFirstDimValue_Target << endl;
	int locNumPoints_Low = locGraphErrors_Low->GetN();
	int locNumPoints_High = locGraphErrors_High->GetN();
	double* locXValues_Low = locGraphErrors_Low->GetX();
	double* locXValues_High = locGraphErrors_High->GetX();
	double locTolerance = 0.000001;

	//find range of x-axis to use (highest low to lowest high)
	double locMinX = (locXValues_Low[0] > locXValues_High[0]) ? locXValues_Low[0] : locXValues_High[0];
	locMinX -= locTolerance;
	double locMaxX = (locXValues_Low[locNumPoints_Low - 1] > locXValues_High[locNumPoints_High - 1]) ? locXValues_High[locNumPoints_High - 1] : locXValues_Low[locNumPoints_Low - 1];
	locMaxX += locTolerance;
if(gDebugLevel > 0)
cout << "min, max X = " << locMinX << ", " << locMaxX << endl;

	//interpolate within first-dim-bin to match X bins
	int locClosestEBinToTargetFlag = 0; //0 is low, 1 is high
	if((locFirstDimValue_High - locFirstDimValue_Target) < (locFirstDimValue_Target - locFirstDimValue_Low))
		locClosestEBinToTargetFlag = 1; //high side closest

	if(locClosestEBinToTargetFlag == 1)
	{
		//high side closest to target, interpolate low-side to match high-side x-axis binning
		locGraphErrors_Low = Convert_DistroToMatch(locGraphErrors_Low, locGraphErrors_High, locMinX, locMaxX); //interpolate
		locGraphErrors_High = Convert_DistroToMatch(locGraphErrors_High, locGraphErrors_Low, locMinX, locMaxX); //truncate to match!
	}
	else
	{
		//low side closest to target, interpolate high-side to match low-side x-axis binning
		locGraphErrors_High = Convert_DistroToMatch(locGraphErrors_High, locGraphErrors_Low, locMinX, locMaxX); //interpolate
		locGraphErrors_Low = Convert_DistroToMatch(locGraphErrors_Low, locGraphErrors_High, locMinX, locMaxX); //truncate to match!
	}
	if((locGraphErrors_Low == NULL) || (locGraphErrors_High == NULL))
		return NULL;

	double* locCrossValues_Low = locGraphErrors_Low->GetY();
	double* locCrossValues_High = locGraphErrors_High->GetY();
	if((locCrossValues_Low == NULL) || (locCrossValues_High == NULL))
		return NULL;

	double* locCrossUncertainties_Low = locGraphErrors_Low->GetEY();
	double* locCrossUncertainties_High = locGraphErrors_High->GetEY();
if(gDebugLevel > 1)
cout << "pointers: " << locCrossValues_Low << ", " << locCrossValues_High << ", " << locCrossUncertainties_Low << ", " << locCrossUncertainties_High << endl;
	double locFirstDimWidth = locFirstDimValue_High - locFirstDimValue_Low;
if(gDebugLevel > 1)
cout << "first-dim's: " << locFirstDimValue_High << ", " << locFirstDimValue_Low << ", " << locFirstDimWidth << endl;

	int locNumPoints_Combined = locGraphErrors_High->GetN(); //same now
	double* locXValues_Combined = locGraphErrors_High->GetX(); //same now
	double* locCrossValues_Combined = new double[locNumPoints_Combined];
	double* locCrossUncertainties_Combined = new double[locNumPoints_Combined];

	//find cross values, uncertainties
if(gDebugLevel > 1)
cout << "num points = " << locNumPoints_Combined << endl;
	for(int loc_i = 0; loc_i < locNumPoints_Combined; ++loc_i)
	{
if(gDebugLevel > 2)
cout << "i, #points = " << loc_i << ", " << locNumPoints_Combined << endl;
		double locCrossUncertaintyLow = locCrossUncertainties_Low[loc_i];
		double locCrossUncertaintyHigh = locCrossUncertainties_High[loc_i];
		double locSlope = (locCrossValues_High[loc_i] - locCrossValues_Low[loc_i])/locFirstDimWidth;
		double locIntercept = locCrossValues_High[loc_i] - locSlope*locFirstDimValue_High;
		locCrossValues_Combined[loc_i] = locSlope*locFirstDimValue_Target + locIntercept;
		double locUncertaintyTerm = (locFirstDimValue_Target - locFirstDimValue_High)/locFirstDimWidth;
if(gDebugLevel > 2)
cout << "crosshigh, crosslow, slope, intercept = " << locCrossValues_High[loc_i] << ", " << locCrossValues_Low[loc_i] << ", " << locSlope << ", " << locIntercept << endl;
		if((locCrossUncertaintyLow > 0.0) && (locCrossUncertaintyHigh > 0.0))
			locCrossUncertainties_Combined[loc_i] = sqrt(locCrossUncertaintyHigh*locCrossUncertaintyHigh*(1.0 + locUncertaintyTerm)*(1.0 + locUncertaintyTerm) + locUncertaintyTerm*locUncertaintyTerm*locCrossUncertaintyLow*locCrossUncertaintyLow);
		else
			locCrossUncertainties_Combined[loc_i] = 0.0;
//if(locCrossValues_Combined[loc_i] < 0.0)
if(gDebugLevel > 2)
cout << "i, cos, value = " << loc_i << ", " << locXValues_Combined[loc_i] << ", " << locCrossValues_Combined[loc_i] << endl;
	}

	TGraphErrors* locGraphErrors = new TGraphErrors(locNumPoints_Combined, locXValues_Combined, locCrossValues_Combined, NULL, locCrossUncertainties_Combined);
	ostringstream locTitleStream;
	locTitleStream << "E_{#gamma} = " << locFirstDimValue_Target << " GeV;" << locGraphErrors_Low->GetXaxis()->GetTitle() << ";" << locGraphErrors_Low->GetYaxis()->GetTitle();
	locGraphErrors->SetTitle(locTitleStream.str().c_str());
	return (locGraphErrors);
}