#include "DStatUtilities.h"

ClassImp(DStatUtilities)

//Control Variables
double DStatUtilities::dMinNumTrialsForAcceptance = 50; // Do not trust acceptance value with less than this many trials

/************************************************************* POINTS *************************************************************/

double DStatUtilities::Calc_BinomialSuccessUncertainty(double locNumSuccesses, double locNumTrials)
{
	return sqrt(locNumSuccesses*(1.0 - locNumSuccesses/locNumTrials));
}

pair<double, double> DStatUtilities::Calc_BinomialRate(double locNumSuccesses, double locNumTrials)
{
	//e.g. acceptance
	double locRate = locNumSuccesses/locNumTrials;
	double locRateError = Calc_BinomialSuccessUncertainty(locNumSuccesses, locNumTrials)/locNumTrials;
	return pair<double, double>(locRate, locRateError);
}

pair<double, double> DStatUtilities::Calc_Sum(pair<double, double> locTerm1, pair<double, double> locTerm2, double locFactor1, double locFactor2)
{
	//Pairs: value, uncertainty
	double locSum = locFactor1*locTerm1.first + locFactor2*locTerm2.first;
	double locSumError = locFactor1*locFactor1*locTerm1.second*locTerm1.second;
	locSumError += locFactor2*locFactor2*locTerm2.second*locTerm2.second;
	locSumError = sqrt(locSumError);
	return pair<double, double>(locSum, locSumError);
}

pair<double, double> DStatUtilities::Calc_Sum(const vector<pair<double, double> >& locTerms)
{
	double locSum = 0.0;
	double locSumError = 0.0;
	for(size_t loc_i = 0; loc_i < locTerms.size(); ++loc_i)
	{
		locSum += locTerms[loc_i].first;
		locSumError += locTerms[loc_i].second*locTerms[loc_i].second;
	}
	locSumError = sqrt(locSumError);
	return pair<double, double>(locSum, locSumError);
}

pair<double, double> DStatUtilities::Calc_Product(pair<double, double> locTerm, double locFactor)
{
	double locProduct = locFactor*locTerm.first;
	double locProductError = locFactor*locTerm.second;
	return pair<double, double>(locProduct, locProductError);
}

pair<double, double> DStatUtilities::Calc_Product(pair<double, double> locTerm1, pair<double, double> locTerm2, double locFactor)
{
	//Pairs: value, uncertainty
/*
	//P = F*V1*V2
	//Variance on P is: 
    = (F*V2, F*V1) (Sig_V1^2,        0) ( F*V2 )
                   (0,        Sig_V2^2) ( F*V1 )

    = (Sig_V1^2*F*V2, Sig_V2^2*F*V1) ( F*V2 )
                                     ( F*V1 )

	 = Sig_V1^2*F^2*V2^2 + Sig_V2^2*F^2*V1^2

	 = F^2*[Sig_V1^2*V2^2 + Sig_V2^2*V1^2]

	 = P^2*[Sig_V1^2/V1^2 + Sig_V2^2/V2^2]

	So Error on P is sqrt of above
*/
	double locProduct = locFactor*locTerm1.first*locTerm2.first;
	double locProductError = locTerm1.second*locTerm1.second/(locTerm1.first*locTerm1.first);
	locProductError += locTerm2.second*locTerm2.second/(locTerm2.first*locTerm2.first);
	locProductError = locProduct*sqrt(locProductError);
	return pair<double, double>(locProduct, locProductError);
}

pair<double, double> DStatUtilities::Calc_Product(const vector<pair<double, double> >& locTerms, double locFactor)
{
	//Pairs: value, uncertainty
	double locProduct = locFactor;
	double locProductError = 0.0;
	for(size_t loc_i = 0; loc_i < locTerms.size(); ++loc_i)
	{
		locProduct *= locTerms[loc_i].first;
		locProductError += locTerms[loc_i].second*locTerms[loc_i].second/(locTerms[loc_i].first*locTerms[loc_i].first);
	}
	locProductError = locProduct*sqrt(locProductError);
	return pair<double, double>(locProduct, locProductError);
}

pair<double, double> DStatUtilities::Calc_Ratio(pair<double, double> locDenominator, double locFactor)
{
	//Pairs: value, uncertainty
/*
	//R = F/D
	//dR/dD = -F/D^2 = -R/D
	//sigma_R = R*sigma_D/D
*/
	double locRatio = locFactor/locDenominator.first;
	double locRatioError = locRatio*locDenominator.second/locDenominator.first;
	return pair<double, double>(locRatio, locRatioError);
}

pair<double, double> DStatUtilities::Calc_Ratio(pair<double, double> locNumerator, pair<double, double> locDenominator, double locFactor)
{
	//Pairs: value, uncertainty
/*
	//R = N/D
	//Variance on R is: 
    = (R/N, -R/D)) (Sig_N^2,       0) (  R/N )
                   (0,       Sig_D^2) ( -R/D )

    = (Sig_N^2*R/N, -Sig_D^2*R/D) (  R/N )
                                  ( -R/D )

	 = Sig_N^2*R^2/(N^2) + Sig_D^2*R^2/(D^2)

	 = R^2*[Sig_N^2/(N^2) + Sig_D^2/(D^2)]

	So Error on R is sqrt of above
*/

	double locRatio = locFactor*locNumerator.first/locDenominator.first;
	double locRatioError = locNumerator.second*locNumerator.second/(locNumerator.first*locNumerator.first);
	locRatioError += locDenominator.second*locDenominator.second/(locDenominator.first*locDenominator.first);
	locRatioError = locRatio*sqrt(locRatioError);
	return pair<double, double>(locRatio, locRatioError);
}

pair<double, double> DStatUtilities::Calc_Ratio(pair<double, double> locNumeratorTerm, const vector<pair<double, double> >& locDenominatorTerms, double locFactor)
{
	//Pairs: value, uncertainty
	vector<pair<double, double> > locNumeratorTerms(1, locNumeratorTerm);
	return Calc_Ratio(locNumeratorTerms, locDenominatorTerms, locFactor);
}

pair<double, double> DStatUtilities::Calc_Ratio(const vector<pair<double, double> >& locNumeratorTerms, pair<double, double> locDenominatorTerm, double locFactor)
{
	//Pairs: value, uncertainty
	vector<pair<double, double> > locDenominatorTerms(1, locDenominatorTerm);
	return Calc_Ratio(locNumeratorTerms, locDenominatorTerms, locFactor);
}

pair<double, double> DStatUtilities::Calc_Ratio(const vector<pair<double, double> >& locNumeratorTerms, const vector<pair<double, double> >& locDenominatorTerms, double locFactor)
{
	//Pairs: value, uncertainty
	double locRatio = locFactor;
	double locRatioError = 0.0;
	for(size_t loc_i = 0; loc_i < locNumeratorTerms.size(); ++loc_i)
	{
		locRatio *= locNumeratorTerms[loc_i].first;
		locRatioError += locNumeratorTerms[loc_i].second*locNumeratorTerms[loc_i].second/(locNumeratorTerms[loc_i].first*locNumeratorTerms[loc_i].first);
	}
	for(size_t loc_i = 0; loc_i < locDenominatorTerms.size(); ++loc_i)
	{
		locRatio /= locDenominatorTerms[loc_i].first;
		locRatioError += locDenominatorTerms[loc_i].second*locDenominatorTerms[loc_i].second/(locDenominatorTerms[loc_i].first*locDenominatorTerms[loc_i].first);
	}

	locRatioError = locRatio*sqrt(locRatioError);
	return pair<double, double>(locRatio, locRatioError);
}

double DStatUtilities::Calc_Average(const vector<double>& locData)
{
	double locAverage = 0.0;
	for(size_t loc_k = 0; loc_k < locData.size(); ++loc_k)
		locAverage += locData[loc_k];
	return locAverage/(double(locData.size()));
}

double DStatUtilities::Calc_SampleSigma(const vector<double>& locData)
{
	return Calc_SampleSigma(Calc_Average(locData), locData);
}

double DStatUtilities::Calc_SampleSigma(double locAverage, const vector<double>& locData)
{
	double locSigma = 0.0;
	for(size_t loc_k = 0; loc_k < locData.size(); ++loc_k)
	{
		double locDelta = locData[loc_k] - locAverage;
		locSigma += locDelta*locDelta;
	}
	return sqrt(locSigma/(locData.size() - 1)); //unbiased
}

pair<double, double> DStatUtilities::Calc_WeightedAverage(const vector<pair<double, double> >& locData)
{
	//Data pair: value, uncertainty
	//returned 2nd double is the uncertainty ON the weighted average, not the sigma of the distribution!!!
	//https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Dealing_with_variance
	if(locData.empty())
		return pair<double, double>(0.0, 0.0);
	if(locData.size() == 1)
		return locData[0];

	double locNumerator = 0.0, locDenominator = 0.0;
	for(size_t loc_k = 0; loc_k < locData.size(); ++loc_k)
	{
		double locWeight = 1.0/(locData[loc_k].second*locData[loc_k].second);
		locNumerator += locData[loc_k].first*locWeight;
		locDenominator += locWeight;
	}
	double locWeightedAverage = locNumerator/locDenominator;
	double locWeightedAverageError = sqrt(1.0/locDenominator); //sigma_error = (d_average/d_data)*sigma_data
	return pair<double, double>(locWeightedAverage, locWeightedAverageError);
}

double DStatUtilities::Calc_WeightedSampleSigma(const vector<pair<double, double> >& locData)
{
	//Data pair: value, uncertainty
	pair<double, double> locWeightedAverage = Calc_WeightedAverage(locData);
	return Calc_WeightedSampleSigma(locWeightedAverage.first, locData);
}

double DStatUtilities::Calc_WeightedSampleSigma(double locWeightedAverage, const vector<pair<double, double> >& locData)
{
	//https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_variance
	//Data pair: value, uncertainty
	if(locData.size() < 2)
		return 0.0;

	double locNumerator = 0.0, locDenominator = 0.0;
	for(size_t loc_i = 0; loc_i < locData.size(); ++loc_i)
	{
		double locDeltaResult = locWeightedAverage - locData[loc_i].first;
		double locWeight = 1.0/(locData[loc_i].second*locData[loc_i].second);
		locNumerator += locDeltaResult*locDeltaResult*locWeight;
		locDenominator += locWeight;
	}

	return sqrt(locNumerator/locDenominator);
}

pair<double, double> DStatUtilities::Calc_SystematicUncertainty(pair<double, double> locTerm1, pair<double, double> locTerm2)
{
	pair<double, double> locSystematicUncertainty = DStatUtilities::Calc_Sum(locTerm1, locTerm2, 1.0, -1.0);
	locSystematicUncertainty.first = fabs(locSystematicUncertainty.first);
	return locSystematicUncertainty;
}

/*********************************************************** HISTOGRAMS ***********************************************************/

TH1* DStatUtilities::Calc_Acceptance(TH1* locSuccessesHist, TH1* locTrialsHist, pair<string, string> locStringReplacement)
{
	TH1* locAcceptanceHist = CloneAndReset_Hist(locSuccessesHist, "_Acceptance", true, locStringReplacement);

	for(int loc_i = 1; loc_i <= locAcceptanceHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locAcceptanceHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locAcceptanceHist->GetNbinsZ(); ++loc_k)
			{
				double locNumSuccesses = locSuccessesHist->GetBinContent(loc_i, loc_j, loc_k);
				double locNumTrials = locTrialsHist->GetBinContent(loc_i, loc_j, loc_k);

				if(!(fabs(locNumSuccesses) > 0.0) || !(fabs(locNumTrials) > dMinNumTrialsForAcceptance))
				{
					locAcceptanceHist->SetBinContent(loc_i, loc_j, loc_k, 0.0);
					locAcceptanceHist->SetBinError(loc_i, loc_j, loc_k, 0.0);
					continue;
				}

				pair<double, double> locAcceptance = Calc_BinomialRate(locNumSuccesses, locNumTrials);
				locAcceptanceHist->SetBinContent(loc_i, loc_j, loc_k, locAcceptance.first);
				locAcceptanceHist->SetBinError(loc_i, loc_j, loc_k, locAcceptance.second);
			}
		}
	}

	locAcceptanceHist->SetEntries(locTrialsHist->GetEntries());
//	locAcceptanceHist->Divide(locSuccessesHist, locTrialsHist, 1.0, 1.0, "B"); //Goes bonkers with error bars
	return locAcceptanceHist;
}

TH1* DStatUtilities::Calc_Acceptance_SeparateTotals(TH1* locSuccessesHist, TH1* locFailuresHist, pair<string, string> locStringReplacement)
{
	TH1* locTrialsHist = CloneAndReset_Hist(locSuccessesHist, "_Trials");
	locTrialsHist->Add(locSuccessesHist, locFailuresHist);
	locTrialsHist->SetEntries(locSuccessesHist->GetEntries() + locFailuresHist->GetEntries());
	return Calc_Acceptance(locSuccessesHist, locTrialsHist, locStringReplacement);
}

TH1* DStatUtilities::Calc_Ratio(TH1* locNumeratorHist, TH1* locDenominatorHist, pair<string, string> locStringReplacement)
{
	TH1* locRatioHist = CloneAndReset_Hist(locNumeratorHist, "_Ratio", true, locStringReplacement);

	for(int loc_i = 1; loc_i <= locRatioHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locRatioHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locRatioHist->GetNbinsZ(); ++loc_k)
			{
				pair<double, double> locNumerator, locDenominator;
				locNumerator.first = locNumeratorHist->GetBinContent(loc_i, loc_j, loc_k);
				locNumerator.second = locNumeratorHist->GetBinError(loc_i, loc_j, loc_k);
				locDenominator.first = locDenominatorHist->GetBinContent(loc_i, loc_j, loc_k);
				locDenominator.second = locDenominatorHist->GetBinError(loc_i, loc_j, loc_k);
				if(!(fabs(locNumerator.first) > 0.0) || !(fabs(locDenominator.first) > 0.0))
				{
					locRatioHist->SetBinContent(loc_i, loc_j, loc_k, 0.0);
					locRatioHist->SetBinError(loc_i, loc_j, loc_k, 0.0);
					continue;
				}

				pair<double, double> locRatio = Calc_Ratio(locNumerator, locDenominator);
				locRatioHist->SetBinContent(loc_i, loc_j, loc_k, locRatio.first);
				locRatioHist->SetBinError(loc_i, loc_j, loc_k, locRatio.second);
			}
		}
	}

	locRatioHist->SetEntries(locNumeratorHist->GetEntries());
//	locRatioHist->Divide(locNumeratorHist, locDenominatorHist, 1.0, 1.0); //Is broken: In bin limit comparison, compares double for equality
	return locRatioHist;
}

TH1* DStatUtilities::Compare(TH1* locCheckHist, TH1* locCompareToHist, pair<string, string> locStringReplacement)
{
	//computes (Check - CompareTo)/CompareTo
	TH1* locComparisonHist = CloneAndReset_Hist(locCheckHist, "_Comparison", true, locStringReplacement);

	for(int loc_i = 1; loc_i <= locCheckHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locCheckHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locCheckHist->GetNbinsZ(); ++loc_k)
			{
				double locCheckValue = locCheckHist->GetBinContent(loc_i, loc_j, loc_k);
				double locCompareToValue = locCompareToHist->GetBinContent(loc_i, loc_j, loc_k);
				if(!(fabs(locCompareToValue) > 0.0) || !(fabs(locCheckValue) > 0.0))
				{
					if(locCheckHist->IsA()->InheritsFrom(TH2::Class()))
						locComparisonHist->SetBinContent(loc_i, loc_j, loc_k, -1.0/0.0); //-inf: So that empty bin is not colored at 0-value
					continue;
				}

				pair<double, double> locCheckPair(locCheckValue, locCheckHist->GetBinError(loc_i, loc_j, loc_k));
				pair<double, double> locCompareToPair(locCompareToValue, locCompareToHist->GetBinError(loc_i, loc_j, loc_k));
				pair<double, double> locRatio = Calc_Ratio(locCheckPair, locCompareToPair);

				locComparisonHist->SetBinContent(loc_i, loc_j, loc_k, locRatio.first - 1.0);
				locComparisonHist->SetBinError(loc_i, loc_j, loc_k, locRatio.second);
			}
		}
	}

	return locComparisonHist;
}

TPair* DStatUtilities::Calc_SystematicUncertainty(TH1* locResultHist, TH1* locAlternativeHist)
{
	TH1* locUncertaintyHist = CloneAndReset_Hist(locResultHist, "_SystUncert");
	TH1* locUncertaintyFractionHist = CloneAndReset_Hist(locResultHist, "_SystUncertFraction");

	//replace title if 2d/3d
	if(locUncertaintyHist->IsA()->InheritsFrom(TH2::Class()) || locUncertaintyHist->IsA()->InheritsFrom(TH3::Class()))
	{
		locUncertaintyHist->SetTitle("#sigma_{Systematic}");
		locUncertaintyFractionHist->SetTitle("#sigma_{Systematic} Fraction");
	}
	else //replace y-axis title if 1d
	{
		locUncertaintyHist->GetYaxis()->SetTitle("#sigma_{Systematic}");
		locUncertaintyFractionHist->GetYaxis()->SetTitle("#sigma_{Systematic} Fraction");
	}

	for(int loc_i = 1; loc_i <= locResultHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locResultHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locResultHist->GetNbinsZ(); ++loc_k)
			{
				double locResultValue = locResultHist->GetBinContent(loc_i, loc_j, loc_k);
				double locAlternativeValue = locAlternativeHist->GetBinContent(loc_i, loc_j, loc_k);

				if(!(fabs(locResultValue) > 0.0) || !(fabs(locAlternativeValue) > 0.0))
					continue; //empty bin, don't compare

				pair<double, double> locResultPair(locResultValue, locResultHist->GetBinError(loc_i, loc_j, loc_k));
				pair<double, double> locAlternativePair(locAlternativeValue, locAlternativeHist->GetBinError(loc_i, loc_j, loc_k));

				pair<double, double> locUncertaintyPair = Calc_SystematicUncertainty(locResultPair, locAlternativePair);
				pair<double, double> locUncertaintyFractionPair = Calc_Ratio(locUncertaintyPair, locResultPair);

				locUncertaintyHist->SetBinContent(loc_i, loc_j, loc_k, locUncertaintyPair.first);
				locUncertaintyHist->SetBinError(loc_i, loc_j, loc_k, locUncertaintyPair.second);

				locUncertaintyFractionHist->SetBinContent(loc_i, loc_j, loc_k, locUncertaintyFractionPair.first);
				locUncertaintyFractionHist->SetBinError(loc_i, loc_j, loc_k, locUncertaintyFractionPair.second);
			}
		}
	}

	locUncertaintyHist->SetEntries(locResultHist->GetEntries());
	locUncertaintyFractionHist->SetEntries(locResultHist->GetEntries());
	return (new TPair(locUncertaintyHist, locUncertaintyFractionHist));
}

/******************************************************* GROUP OF HISTOGRAMS ******************************************************/

TH1* DStatUtilities::Calc_TotalSystematicUncertainty(TObjArray* locUncertaintyHists)
{
	//input array can be value or fraction, it works for both
	TH1* locFirstHist = (TH1*)locUncertaintyHists->At(0);
	TH1* locTotalUncertaintyHist = CloneAndReset_Hist(locFirstHist, "_Total");

	for(int loc_i = 1; loc_i <= locTotalUncertaintyHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locTotalUncertaintyHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locTotalUncertaintyHist->GetNbinsZ(); ++loc_k)
			{
				double locTotalUncertainty = 0.0;
				double locTotalUncertaintyError = 0.0;
				for(int loc_l = 0; loc_l < locUncertaintyHists->GetEntriesFast(); ++loc_l)
				{
					TH1* locHist = (TH1*)locUncertaintyHists->At(loc_l);
					pair<double, double> locStudyError(locHist->GetBinContent(loc_i, loc_j, loc_k), locHist->GetBinError(loc_i, loc_j, loc_k));
					if(!(locStudyError.first > 0.0))
						continue;
					locTotalUncertainty += locStudyError.first*locStudyError.first;
					locTotalUncertaintyError += locStudyError.first*locStudyError.first*locStudyError.second*locStudyError.second;
				}

				if(!(locTotalUncertainty > 0.0))
					continue;
/*
sigma_total = sqrt(var1 + var2 + ...)
d_sigma_total/dvar_i = sigma_i/sigma_total
sigma_total_var = Sum_i [ var_sigma_i * var_i / var_total) ]
sigma_total_var = Sum_i [ var_sigma_i * var_i ] / var_total
so sigma is sqrt of above
*/
				locTotalUncertainty = sqrt(locTotalUncertainty);
				locTotalUncertaintyError = sqrt(locTotalUncertaintyError)/locTotalUncertainty;

				locTotalUncertaintyHist->SetBinContent(loc_i, loc_j, loc_k, locTotalUncertainty);
				locTotalUncertaintyHist->SetBinError(loc_i, loc_j, loc_k, locTotalUncertaintyError);
			}
		}
	}

	locTotalUncertaintyHist->SetEntries(locFirstHist->GetEntries());
	return locTotalUncertaintyHist;
}

TH1* DStatUtilities::Calc_WeightedAverage(TObjArray* locHistArray)
{
	TH1* locFirstHist = (TH1*)locHistArray->At(0);
	TH1* locWeightedAverageHist = CloneAndReset_Hist(locFirstHist, "_WeightedAverage");

	for(int loc_i = 1; loc_i <= locWeightedAverageHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locWeightedAverageHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locWeightedAverageHist->GetNbinsZ(); ++loc_k)
			{
				vector<pair<double, double> > locDataPoints;
				for(int loc_l = 0; loc_l < locHistArray->GetEntriesFast(); ++loc_l)
				{
					TH1* locHist = (TH1*)locHistArray->At(loc_l);
					pair<double, double> locDataPoint(locHist->GetBinContent(loc_i, loc_j, loc_k), locHist->GetBinError(loc_i, loc_j, loc_k));
					if(fabs(locDataPoint.first) > 0.0)
						locDataPoints.push_back(locDataPoint);
				}

				pair<double, double> locWeightedAverage = Calc_WeightedAverage(locDataPoints);
				locWeightedAverageHist->SetBinContent(loc_i, loc_j, loc_k, locWeightedAverage.first);
				locWeightedAverageHist->SetBinError(loc_i, loc_j, loc_k, locWeightedAverage.second);
			}
		}
	}

	locWeightedAverageHist->SetEntries(locFirstHist->GetEntries());
	return locWeightedAverageHist;
}

TH1* DStatUtilities::Calc_WeightedSampleSigmaFraction(TObjArray* locHistArray)
{
	TH1* locFirstHist = (TH1*)locHistArray->At(0);
	TH1* locWeightedSampleSigmaHist = CloneAndReset_Hist(locFirstHist, "_WeightedSampleSigma");

	for(int loc_i = 1; loc_i <= locWeightedSampleSigmaHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locWeightedSampleSigmaHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locWeightedSampleSigmaHist->GetNbinsZ(); ++loc_k)
			{
				vector<pair<double, double> > locDataPoints;
				for(int loc_l = 0; loc_l < locHistArray->GetEntriesFast(); ++loc_l)
				{
					TH1* locHist = (TH1*)locHistArray->At(loc_l);
					pair<double, double> locDataPoint(locHist->GetBinContent(loc_i, loc_j, loc_k), locHist->GetBinError(loc_i, loc_j, loc_k));
					if(fabs(locDataPoint.first) > 0.0)
						locDataPoints.push_back(locDataPoint);
				}

				pair<double, double> locWeightedAverage = Calc_WeightedAverage(locDataPoints);
				if(!(locWeightedAverage.first > 0.0))
					continue;
				double locWeightedSampleSigma = Calc_WeightedSampleSigma(locWeightedAverage.first, locDataPoints);
				locWeightedSampleSigmaHist->SetBinContent(loc_i, loc_j, loc_k, locWeightedSampleSigma/locWeightedAverage.first);
			}
		}
	}

	locWeightedSampleSigmaHist->SetEntries(locFirstHist->GetEntries());
	return locWeightedSampleSigmaHist;
}

/*********************************************************** STRUCTURES ***********************************************************/

TObjArray* DStatUtilities::Perform_Action(const TObjArray* locHistArray1, const TObjArray* locHistArray2, string locAction, pair<string, string> locStringReplacement)
{
	TObjArray* locResultHistArray = new TObjArray();
	for(int loc_i = 0; loc_i < locHistArray1->GetEntriesFast(); ++loc_i)
	{
		TObject* locObject = locHistArray1->At(loc_i);
		if(string(locObject->ClassName()) == "TObjArray")
		{
			TObjArray* locNestedHistArray1 = (TObjArray*)locHistArray1->At(loc_i);
			TObjArray* locNestedHistArray2 = (TObjArray*)locHistArray2->At(loc_i);
			locResultHistArray->AddLast(Perform_Action(locNestedHistArray1, locNestedHistArray2, locAction, locStringReplacement));
		}
		else
		{
			TH1* locHist1 = (TH1*)locHistArray1->At(loc_i);
			TH1* locHist2 = (TH1*)locHistArray2->At(loc_i);
			if(locAction == "Calc_Acceptance")
				locResultHistArray->AddLast(Calc_Acceptance(locHist1, locHist2, locStringReplacement));
			else if(locAction == "Calc_Acceptance_SeparateTotals")
				locResultHistArray->AddLast(Calc_Acceptance_SeparateTotals(locHist1, locHist2, locStringReplacement));
			else if(locAction == "Calc_Ratio")
				locResultHistArray->AddLast(Calc_Ratio(locHist1, locHist2, locStringReplacement));
			else if(locAction == "Compare")
				locResultHistArray->AddLast(Compare(locHist1, locHist2, locStringReplacement));
			else if(locAction == "Calc_SystematicUncertainty") //Returns a TPair!!
				locResultHistArray->AddLast(Calc_SystematicUncertainty(locHist1, locHist2));
		}
	}
	return locResultHistArray;
}

void DStatUtilities::Perform_Action(TDirectory* locDir1, TDirectory* locDir2, TDirectory* locOutputDir, string locAction, pair<string, string> locStringReplacement)
{
	//converted from code by Michael Staib
	cout << "current path = " << locDir1->GetPath() << endl;

	// loop over all keys in this directory
	TList* locListOfKeys = locDir1->GetListOfKeys();
	TKey* locPreviousKey = NULL;
	for(int loc_i = 0; loc_i < locListOfKeys->GetEntries(); ++loc_i)
	{
		TKey* locKey = (TKey*)locListOfKeys->At(loc_i);
		if(locKey == NULL)
			continue;
		string locKeyName = locKey->GetName();

		//keep only the highest cycle number for each key
		if(locPreviousKey != NULL)
		{
			string locPreviousKeyName = locPreviousKey->GetName();
			if(locPreviousKeyName == locKeyName)
				continue;
		}

		TObject* locObject1 = locDir1->Get(locKeyName.c_str());
		TObject* locObject2 = locDir2->Get(locKeyName.c_str());

		if(locObject1->IsA()->InheritsFrom(TDirectory::Class()))
		{
			string locDirName = locObject1->GetName();
			TDirectory* locNewOutputDir = locOutputDir->mkdir(locDirName.c_str(), locDirName.c_str());
			TDirectory* locNewDir1 = locDir1->GetDirectory(locDirName.c_str());
			TDirectory* locNewDir2 = locDir2->GetDirectory(locDirName.c_str());
			Perform_Action(locNewDir1, locNewDir2, locNewOutputDir, locAction, locStringReplacement);
		}
		else if(locObject1->IsA()->InheritsFrom(TH1::Class()))
		{
			// create a new subdir of same name and title in the target file
			locOutputDir->cd();
			if(locAction == "Calc_Acceptance")
				Calc_Acceptance((TH1*)locObject1, (TH1*)locObject2, locStringReplacement);
			else if(locAction == "Calc_Acceptance_SeparateTotals")
				Calc_Acceptance_SeparateTotals((TH1*)locObject1, (TH1*)locObject2, locStringReplacement);
			else if(locAction == "Calc_Ratio")
				Calc_Ratio((TH1*)locObject1, (TH1*)locObject2, locStringReplacement);
			else if(locAction == "Compare")
				Compare((TH1*)locObject1, (TH1*)locObject2, locStringReplacement);
			else if(locAction == "Calc_SystematicUncertainty")
				Calc_SystematicUncertainty((TH1*)locObject1, (TH1*)locObject2);
		}

		locPreviousKey = locKey;
	}
}

/*********************************************************** UTILITIES ************************************************************/

//THIS BELONGS ELSEWHERE
TH1* DStatUtilities::CloneAndReset_Hist(TH1* locReferenceHist, string locNameSuffix, bool locCreateDoubleFlag, pair<string, string> locStringReplacement)
{
	//if locCreateDoubleFlag, created hist will be TH1D, TH2D, or TH3D
		//e.g. calcing ratios of TH1I's!!
	string locHistName = string(locReferenceHist->GetName()) + locNameSuffix;

	TH1* locNewHist = NULL;
	if(!locCreateDoubleFlag)
	{
		locNewHist = (TH1*)locReferenceHist->Clone(locHistName.c_str());
		locNewHist->Reset("M"); //also resets min/max
	}
	else
	{
		string locHistTitle = string(locReferenceHist->GetTitle()) + string(";") + string(locReferenceHist->GetXaxis()->GetTitle()) + string(";");
		locHistTitle += string(locReferenceHist->GetYaxis()->GetTitle()) + string(";") + string(locReferenceHist->GetZaxis()->GetTitle());

		bool locIs3DFlag = locReferenceHist->IsA()->InheritsFrom(TH3::Class());
		bool locIs2DFlag = ((!locIs3DFlag) && locReferenceHist->IsA()->InheritsFrom(TH2::Class()));
		bool locIs1DFlag = ((!locIs3DFlag) && (!locIs2DFlag));

		//get bin edges!!
		Double_t *locXBinEdges = NULL, *locYBinEdges = NULL, *locZBinEdges = NULL;

		//X
		locXBinEdges = new Double_t[locReferenceHist->GetNbinsX() + 1];
		for(Int_t loc_i = 0; loc_i < locReferenceHist->GetNbinsX(); ++loc_i)
			locXBinEdges[loc_i] = locReferenceHist->GetXaxis()->GetBinLowEdge(loc_i + 1);
		locXBinEdges[locReferenceHist->GetNbinsX()] = locReferenceHist->GetXaxis()->GetBinUpEdge(locReferenceHist->GetNbinsX());

		//Y
		if(locIs2DFlag || locIs3DFlag)
		{
			locYBinEdges = new Double_t[locReferenceHist->GetNbinsY() + 1];
			for(Int_t loc_i = 0; loc_i < locReferenceHist->GetNbinsY(); ++loc_i)
				locYBinEdges[loc_i] = locReferenceHist->GetYaxis()->GetBinLowEdge(loc_i + 1);
			locYBinEdges[locReferenceHist->GetNbinsY()] = locReferenceHist->GetYaxis()->GetBinUpEdge(locReferenceHist->GetNbinsY());
		}

		//Z
		if(locIs3DFlag)
		{
			locZBinEdges = new Double_t[locReferenceHist->GetNbinsZ() + 1];
			for(Int_t loc_i = 0; loc_i < locReferenceHist->GetNbinsZ(); ++loc_i)
				locZBinEdges[loc_i] = locReferenceHist->GetZaxis()->GetBinLowEdge(loc_i + 1);
			locZBinEdges[locReferenceHist->GetNbinsZ()] = locReferenceHist->GetZaxis()->GetBinUpEdge(locReferenceHist->GetNbinsZ());
		}

		if(locIs1DFlag)
			locNewHist = new TH1D(locHistName.c_str(), locHistTitle.c_str(), locReferenceHist->GetNbinsX(), locXBinEdges);
		else if(locIs2DFlag)
			locNewHist = new TH2D(locHistName.c_str(), locHistTitle.c_str(), locReferenceHist->GetNbinsX(), locXBinEdges, locReferenceHist->GetNbinsY(), locYBinEdges);
		else if(locIs3DFlag)
		{
			locNewHist = new TH3D(locHistName.c_str(), locHistTitle.c_str(), locReferenceHist->GetNbinsX(), locXBinEdges, \
			locReferenceHist->GetNbinsY(), locYBinEdges, locReferenceHist->GetNbinsZ(), locZBinEdges);
		}
	}

	Replace_TitleSubString(locNewHist, locStringReplacement);
	return locNewHist;
}

//THIS BELONGS ELSEWHERE
void DStatUtilities::Replace_TitleSubString(TH1* locHist, pair<string, string> locStringReplacement)
{
	if(locStringReplacement.first == "")
		return;

	//replace title if 2d/3d
	if(locHist->IsA()->InheritsFrom(TH2::Class()) || locHist->IsA()->InheritsFrom(TH3::Class()))
	{
		TString locTString = locHist->GetTitle();
		locTString.ReplaceAll(locStringReplacement.first.c_str(), locStringReplacement.second.c_str());
		locHist->SetTitle((const char*)locTString);
	}
	else //replace y-axis title if 1d
	{
		TString locTString = locHist->GetYaxis()->GetTitle();
		locTString.ReplaceAll(locStringReplacement.first.c_str(), locStringReplacement.second.c_str());
		locHist->GetYaxis()->SetTitle((const char*)locTString);
	}
}