#include "DCalcCrossSection.h"

ClassImp(DCalcCrossSection)
ClassImp(DCrossSectionFactors)

//Constants
const double DCrossSectionFactors::dAvogadrosNumber = 6.02214179E23; // (1/mol) (from http://physics.nist.gov/cgi-bin/cuu/Value?na)

/***************************************************** CROSS SECTION: STATISTICAL UNCERTAINTIES *****************************************************/

TPair* DCalcCrossSection::Calc_CrossSection(DCrossSectionFactors locFactors, TH1D* locACYHist, double locFlux, double locYCF)
{
	//STATISTICAL ERROR ONLY!! (Systematics (target factor, branching ratio, etc) need to be included elsewhere!!)
	double locCrossSectionFactor = locFactors.Calc_CrossSectionFactor();

	//Make cross section histogram
	string locHistName = string(locACYHist->GetName()) + string("_CrossSection");
	TH1D* locCrossSectionHist = (TH1D*)locACYHist->Clone(locHistName.c_str());
	locCrossSectionHist->Reset("M");

	//Set Cross-section-axis title
	string locCrossSectionTitle = Build_CrossSectionTitle(locFactors.dBarnesFactor, locFactors.dDCosThetaFlag);
	locCrossSectionHist->GetYaxis()->SetTitle(locCrossSectionTitle.c_str());

	//Make error-fraction histogram
	locHistName = string(locACYHist->GetName()) + string("_CrossSectionStatErrors");
	TH1D* locStatErrorHist = (TH1D*)locACYHist->Clone(locHistName.c_str());
	locStatErrorHist->Reset("M");

	//Set error-y-axis title
	locStatErrorHist->GetYaxis()->SetTitle("#sigma_{Statistical} Fraction");

	for(int loc_i = 1; loc_i <= locACYHist->GetNbinsX(); ++loc_i)
	{
		pair<double, double> locACY(locACYHist->GetBinContent(loc_i), locACYHist->GetBinError(loc_i));
		double locBinWidth = locACYHist->GetBinWidth(loc_i);

		if(!(locACY.first > 0.0) || !(locFlux > 0.0))
			continue;

		pair<double, double> locCrossSection = Calc_CrossSection(locCrossSectionFactor, locACY, locBinWidth, locFlux, locYCF);

		locCrossSectionHist->SetBinContent(loc_i, locCrossSection.first);
		locCrossSectionHist->SetBinError(loc_i, locCrossSection.second);
		locStatErrorHist->SetBinContent(loc_i, locCrossSection.second/locCrossSection.first);
	}

	locCrossSectionHist->SetEntries(locACYHist->GetEntries());
	locStatErrorHist->SetEntries(locACYHist->GetEntries());

	return (new TPair(locCrossSectionHist, locStatErrorHist));
}

TPair* DCalcCrossSection::Calc_CrossSection(DCrossSectionFactors locFactors, TH1D* locACYHist, TH1D* locFluxHist, double locYCF)
{
	//STATISTICAL ERROR ONLY!! (Systematics (target factor, branching ratio, etc) need to be included elsewhere!!)
	double locCrossSectionFactor = locFactors.Calc_CrossSectionFactor();

	//Make cross section histogram
	string locHistName = string(locACYHist->GetName()) + string("_CrossSection");
	TH1D* locCrossSectionHist = (TH1D*)locACYHist->Clone(locHistName.c_str());
	locCrossSectionHist->Reset("M");

	//Set Cross-section-axis title
	string locCrossSectionTitle = Build_CrossSectionTitle(locFactors.dBarnesFactor, locFactors.dDCosThetaFlag);
	locCrossSectionHist->GetYaxis()->SetTitle(locCrossSectionTitle.c_str());

	//Make error-fraction histogram
	locHistName = string(locACYHist->GetName()) + string("_CrossSectionStatErrors");
	TH1D* locStatErrorHist = (TH1D*)locACYHist->Clone(locHistName.c_str());
	locStatErrorHist->Reset("M");

	//Set error-y-axis title
	locStatErrorHist->GetYaxis()->SetTitle("#sigma_{Statistical} Fraction");

	for(int loc_i = 1; loc_i <= locACYHist->GetNbinsX(); ++loc_i)
	{
		pair<double, double> locACY(locACYHist->GetBinContent(loc_i), locACYHist->GetBinError(loc_i));
		double locFlux = locFluxHist->GetBinContent(loc_i);

		if(!(locACY.first > 0.0) || !(locFlux > 0.0))
			continue;

		pair<double, double> locCrossSection = Calc_CrossSection(locCrossSectionFactor, locACY, 2.0, locFlux, locYCF); //2: integrating over all costheta (-1 to 1)

		locCrossSectionHist->SetBinContent(loc_i, locCrossSection.first);
		locCrossSectionHist->SetBinError(loc_i, locCrossSection.second);
		locStatErrorHist->SetBinContent(loc_i, locCrossSection.second/locCrossSection.first);
	}

	locCrossSectionHist->SetEntries(locACYHist->GetEntries());
	locStatErrorHist->SetEntries(locACYHist->GetEntries());

	return (new TPair(locCrossSectionHist, locStatErrorHist));
}

TPair* DCalcCrossSection::Calc_CrossSection(DCrossSectionFactors locFactors, TH2D* locACYHist, TH1D* locFluxHist, double locYCF)
{
	//ASSUMES: locFluxHist corresponds to locACYHist x-axis!
	//STATISTICAL ERROR ONLY!! (Systematics (target factor, branching ratio, etc) need to be included elsewhere!!)
	double locCrossSectionFactor = locFactors.Calc_CrossSectionFactor();

	//Make cross section histogram
	string locHistName = string(locACYHist->GetName()) + string("_CrossSection");
	TH1D* locCrossSectionHist = (TH1D*)locACYHist->Clone(locHistName.c_str());
	locCrossSectionHist->Reset("M");

	//Set Cross-section title
	string locCrossSectionTitle = Build_CrossSectionTitle(locFactors.dBarnesFactor, locFactors.dDCosThetaFlag);
	locCrossSectionHist->SetTitle(locCrossSectionTitle.c_str());

	//Make error-fraction histogram
	locHistName = string(locACYHist->GetName()) + string("_CrossSectionStatErrors");
	TH1D* locStatErrorHist = (TH1D*)locACYHist->Clone(locHistName.c_str());
	locStatErrorHist->Reset("M");

	//Set error title
	locStatErrorHist->SetTitle("#sigma_{Statistical} Fraction");

	for(int loc_i = 1; loc_i <= locACYHist->GetNbinsX(); ++loc_i)
	{
		double locFlux = locFluxHist->GetBinContent(loc_i);
		for(int loc_j = 1; loc_j <= locACYHist->GetNbinsY(); ++loc_j)
		{
			pair<double, double> locACY(locACYHist->GetBinContent(loc_i, loc_j), locACYHist->GetBinError(loc_i, loc_j));
			double locBinWidth = locACYHist->GetYaxis()->GetBinWidth(loc_j);

			if(!(locACY.first > 0.0) || !(locFlux > 0.0))
				continue;

			pair<double, double> locCrossSection = Calc_CrossSection(locCrossSectionFactor, locACY, locBinWidth, locFlux, locYCF);

			locCrossSectionHist->SetBinContent(loc_i, loc_j, locCrossSection.first);
			locCrossSectionHist->SetBinError(loc_i, loc_j, locCrossSection.second);
			locStatErrorHist->SetBinContent(loc_i, loc_j, locCrossSection.second/locCrossSection.first);
		}
	}

	locCrossSectionHist->SetEntries(locACYHist->GetEntries());
	locStatErrorHist->SetEntries(locACYHist->GetEntries());

	return (new TPair(locCrossSectionHist, locStatErrorHist));
}

pair<double, double> DCalcCrossSection::Calc_CrossSection(DCrossSectionFactors locFactors, pair<double, double> locACY, double locBinWidth, double locFlux, double locYCF)
{
	double locCrossSectionFactor = locFactors.Calc_CrossSectionFactor();
	return Calc_CrossSection(locCrossSectionFactor, locACY, locBinWidth, locFlux, locYCF);
}

//INTERNAL, PRIVATE
pair<double, double> DCalcCrossSection::Calc_CrossSection(double locCrossSectionFactor, pair<double, double> locACY, double locBinWidth, double locFlux, double locYCF)
{
	//STATISTICAL ERROR ONLY!! (Systematics (target factor, branching ratio, etc) need to be included elsewhere!!)

	//locACY: Acceptance-Corrected Yield
	//locYCF: Yield-Correction Factor
	if(!(locACY.first > 0.0) || !(locBinWidth > 0.0) || !(locFlux > 0.0))
		return pair<double, double>(0.0, 0.0);

	//ycf & flux uncertainties are scale uncertainties: incorporate later
	vector<pair<double, double>	> locDataNumerator;
	locDataNumerator.push_back(locACY);
	locDataNumerator.push_back(pair<double, double>(locYCF, 0.0));
	return DStatUtilities::Calc_Ratio(locDataNumerator, pair<double, double>(locFlux, 0.0), locCrossSectionFactor/locBinWidth);
}

/************************************************* CROSS SECTION: SYSTEMATIC & SCALE UNCERTAINTIES **************************************************/

TPair* DCalcCrossSection::Calc_CrossSection_SystematicErrors(TH1* locCrossHist, TH1* locACYSystFractionHist)
{
	//Works for TH1*, TH2*, and TH3*
	//Pair: first is cross section with syst errors, second is final systematic uncertainty fraction hist

	//Make total syst histogram
		//Errors are from ACY only
	string locHistName = string(locCrossHist->GetName()) + string("SystErrors");
	TH1* locTotalSystFractionHist = (TH1*)locACYSystFractionHist->Clone(locHistName.c_str());

	//Make cross section histogram
	locHistName = string(locCrossHist->GetName()) + string("WithSyst");
	TH1* locCrossSectionWithSystHist = (TH1*)locCrossHist->Clone(locHistName.c_str());

	//LOOP OVER HISTOGRAM, SET THE ERROR
	for(int loc_i = 1; loc_i <= locCrossHist->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locCrossHist->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locCrossHist->GetNbinsZ(); ++loc_k)
			{
				double locCrossSection = locCrossHist->GetBinContent(loc_i, loc_j, loc_k);
				if(!(fabs(locCrossSection) > 0.0))
					continue;

				double locTotalSystFraction = locACYSystFractionHist->GetBinContent(loc_i, loc_j, loc_k);
				locCrossSectionWithSystHist->SetBinError(loc_i, loc_j, loc_k, locCrossSection*locTotalSystFraction);
			}
		}
	}

	locCrossSectionWithSystHist->SetEntries(locCrossHist->GetEntries());
	locTotalSystFractionHist->SetEntries(locCrossHist->GetEntries());

	return (new TPair(locCrossSectionWithSystHist, locTotalSystFractionHist));
}

TPair* DCalcCrossSection::Calc_CrossSection_ScaleErrors(DCrossSectionFactors locFactors, TH1* locCrossHist, pair<double, double> locFluxErrorFractionPair, pair<double, double> locYCFErrorFractionPair)
{
	//Works for TH1*
	//Output TPair: first is cross section with scale errors, second is final scale uncertainty fraction hist
	//locFluxErrorFractionPair, locYCFErrorFractionPair: first is stat, second is syst

	//Make total syst histogram
		//Errors are from ACY only
	string locHistName = string(locCrossHist->GetName()) + string("ScaleErrors");
	TH1* locScaleFractionHist = (TH1*)locCrossHist->Clone(locHistName.c_str());
	locScaleFractionHist->Reset("M");

	//Make cross section histogram
	locHistName = string(locCrossHist->GetName()) + string("WithScale");
	TH1* locCrossSectionWithScaleHist = (TH1*)locCrossHist->Clone(locHistName.c_str());

	//Systematic uncertainties add in quadrature
		//Think of each one as an independent dimension in some Hilbert space

	double locBranchingRatioErrorFraction = locFactors.dBranchingRatio.second/locFactors.dBranchingRatio.first;
	double locIndependentVariance = locBranchingRatioErrorFraction*locBranchingRatioErrorFraction; //Branching ratio
	locIndependentVariance += locFactors.dTargetErrorFraction*locFactors.dTargetErrorFraction; //Target
	locIndependentVariance += locFluxErrorFractionPair.first*locFluxErrorFractionPair.first; //Flux stat error
	locIndependentVariance += locFluxErrorFractionPair.second*locFluxErrorFractionPair.second; //Flux syst error
	locIndependentVariance += locYCFErrorFractionPair.first*locYCFErrorFractionPair.first; //YCF stat error
	locIndependentVariance += locYCFErrorFractionPair.second*locYCFErrorFractionPair.second; //YCF syst error
	double locScaleErrorFraction = sqrt(locIndependentVariance);

	//LOOP OVER HISTOGRAM, SQUARE THE ERROR AND SET
	for(int loc_i = 1; loc_i <= locCrossHist->GetNbinsX(); ++loc_i)
	{
		double locCrossSection = locCrossHist->GetBinContent(loc_i);
		if(!(fabs(locCrossSection) > 0.0))
			continue;

		locCrossSectionWithScaleHist->SetBinError(loc_i, locCrossSection*locScaleErrorFraction);
		locScaleFractionHist->SetBinContent(loc_i, locScaleErrorFraction);
	}

	locCrossSectionWithScaleHist->SetEntries(locCrossHist->GetEntries());
	locScaleFractionHist->SetEntries(locCrossHist->GetEntries());

	return (new TPair(locCrossSectionWithScaleHist, locScaleFractionHist));
}

TPair* DCalcCrossSection::Calc_CrossSection_ScaleErrors(DCrossSectionFactors locFactors, TH2* locCrossHist, TH1* locFluxHist, double locFluxSystErrorFraction, pair<double, double> locYCFErrorFractionPair)
{
	//Works for TH2*
	//Output TPair: first is cross section with scale errors, second is final scale uncertainty fraction hist
	//pair<double, double> locYCFErrorPair: first is stat, second is syst

	//Make total syst histogram
		//Errors are from ACY only
	string locHistName = string(locCrossHist->GetName()) + string("ScaleErrors");
	TH1* locScaleFractionHist = (TH1*)locCrossHist->Clone(locHistName.c_str());
	locScaleFractionHist->Reset("M");

	//Make cross section histogram
	locHistName = string(locCrossHist->GetName()) + string("WithScale");
	TH1* locCrossSectionWithScaleHist = (TH1*)locCrossHist->Clone(locHistName.c_str());

	//Systematic uncertainties add in quadrature
		//Think of each one as an independent dimension in some Hilbert space

	double locBranchingRatioErrorFraction = locFactors.dBranchingRatio.second/locFactors.dBranchingRatio.first;
	double locIndependentVariance = locBranchingRatioErrorFraction*locBranchingRatioErrorFraction; //Branching ratio
	locIndependentVariance += locFactors.dTargetErrorFraction*locFactors.dTargetErrorFraction; //Target
	locIndependentVariance += locFluxSystErrorFraction*locFluxSystErrorFraction; //Flux syst error
	locIndependentVariance += locYCFErrorFractionPair.first*locYCFErrorFractionPair.first; //YCF stat error
	locIndependentVariance += locYCFErrorFractionPair.second*locYCFErrorFractionPair.second; //YCF syst error

	//LOOP OVER HISTOGRAM, SQUARE THE ERROR AND SET
	for(int loc_i = 1; loc_i <= locCrossHist->GetNbinsX(); ++loc_i)
	{
		double locFluxStatErrorFraction = locFluxHist->GetBinError(loc_i)/locFluxHist->GetBinContent(loc_i);
		for(int loc_j = 1; loc_j <= locCrossHist->GetNbinsY(); ++loc_j)
		{
			double locCrossSection = locCrossHist->GetBinContent(loc_i, loc_j);
			if(!(fabs(locCrossSection) > 0.0))
				continue;

			double locScaleFraction = sqrt(locIndependentVariance + locFluxStatErrorFraction*locFluxStatErrorFraction);

			locCrossSectionWithScaleHist->SetBinError(loc_i, loc_j, locCrossSection*locScaleFraction);
			locScaleFractionHist->SetBinContent(loc_i, loc_j, locScaleFraction);
		}
	}

	locCrossSectionWithScaleHist->SetEntries(locCrossHist->GetEntries());
	locScaleFractionHist->SetEntries(locCrossHist->GetEntries());

	return (new TPair(locCrossSectionWithScaleHist, locScaleFractionHist));
}

/******************************************************** CROSS SECTION: TOTAL UNCERTAINTIES ********************************************************/

TPair* DCalcCrossSection::Calc_CrossSection_TotalErrors(TH1* locCrossHistWithStat, TH1* locCrossHistWithSyst, TH1* locCrossHistWithScale)
{
	//Works for TH1*, TH2*, and TH3*
	//Pair: first is cross section with total errors, second is total systematic uncertainty fraction hist

	//Make total error histogram
		//Errors are from ACY only
	string locHistName = string(locCrossHistWithStat->GetName()) + string("TotalErrors");
	TH1* locTotalErrorFractionHist = (TH1*)locCrossHistWithStat->Clone(locHistName.c_str());
	locTotalErrorFractionHist->Reset("M");

	//Make cross section histogram
	locHistName = string(locCrossHistWithStat->GetName()) + string("WithTotal");
	TH1* locCrossSectionWithTotalHist = (TH1*)locCrossHistWithStat->Clone(locHistName.c_str());

	//LOOP OVER HISTOGRAM, SQUARE THE ERROR AND SET
	for(int loc_i = 1; loc_i <= locCrossHistWithStat->GetNbinsX(); ++loc_i)
	{
		for(int loc_j = 1; loc_j <= locCrossHistWithStat->GetNbinsY(); ++loc_j)
		{
			for(int loc_k = 1; loc_k <= locCrossHistWithStat->GetNbinsZ(); ++loc_k)
			{
				double locCrossSection = locCrossHistWithStat->GetBinContent(loc_i, loc_j, loc_k);
				if(!(fabs(locCrossSection) > 0.0))
					continue;

				double locStatError = locCrossHistWithStat->GetBinError(loc_i, loc_j, loc_k);
				double locSystError = locCrossHistWithSyst->GetBinError(loc_i, loc_j, loc_k);
				double locScaleError = locCrossHistWithScale->GetBinError(loc_i, loc_j, loc_k);
				double locTotalError = sqrt(locStatError*locStatError + locSystError*locSystError + locScaleError*locScaleError);

				locCrossSectionWithTotalHist->SetBinError(loc_i, loc_j, loc_k, locTotalError);
				locTotalErrorFractionHist->SetBinContent(loc_i, loc_j, loc_k, locTotalError/locCrossSection);
			}
		}
	}

	locCrossSectionWithTotalHist->SetEntries(locCrossHistWithStat->GetEntries());
	locTotalErrorFractionHist->SetEntries(locCrossHistWithStat->GetEntries());

	return (new TPair(locCrossSectionWithTotalHist, locTotalErrorFractionHist));
}

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

string DCalcCrossSection::Build_CrossSectionTitle(double locBarnesFactor, bool locDCosThetaFlag)
{
	string locBarnesString = "??";
	if(!(fabs(locBarnesFactor - 1.0E24) > 0.0))
		locBarnesString = "b";
	else if(!(fabs(locBarnesFactor - 1.0E27) > 0.0))
		locBarnesString = "mb";
	else if(!(fabs(locBarnesFactor - 1.0E30) > 0.0))
		locBarnesString = "#mub";
	else if(!(fabs(locBarnesFactor - 1.0E33) > 0.0))
		locBarnesString = "nb";
	else if(!(fabs(locBarnesFactor - 1.0E36) > 0.0))
		locBarnesString = "pb";
	else if(!(fabs(locBarnesFactor - 1.0E39) > 0.0))
		locBarnesString = "fb";

	string locUnitsString = "";
	if(locDCosThetaFlag && (locBarnesString != ""))
		locUnitsString = string(" (") + locBarnesString + string(")");
	else
		locUnitsString = string(" (") + locBarnesString + string("/sr)");

	string locCrossSectionTitle = locDCosThetaFlag ? "d#sigma/dcos(#theta)" : "d#sigma/d#Omega";
	if(locUnitsString != "")
		locCrossSectionTitle += locUnitsString;

	return locCrossSectionTitle;
}

/*********************************************************************** FLUX ***********************************************************************/

void DCalcCrossSection::Calc_FluxForBinning(TH1D *locInputFluxHist, TH1D* locReBinnedFluxHist)
{
	//locReBinnedFluxHist must exist prior to entry: informs binning
	double locTotalFlux = 0.0;
	for(int loc_k = 1; loc_k <= locReBinnedFluxHist->GetNbinsX(); ++loc_k)
	{
		double locRangeMin = locReBinnedFluxHist->GetXaxis()->GetBinLowEdge(loc_k);
		double locRangeMax = locReBinnedFluxHist->GetXaxis()->GetBinUpEdge(loc_k);

		pair<double, double> locFlux = Calc_FluxInRange(locRangeMin, locRangeMax, locInputFluxHist);

		locReBinnedFluxHist->SetBinContent(loc_k, locFlux.first);
		locReBinnedFluxHist->SetBinError(loc_k, locFlux.second);
		locTotalFlux += locFlux.first;
	}
	locReBinnedFluxHist->SetEntries(locTotalFlux);
}

pair<double, double> DCalcCrossSection::Calc_FluxInRange(double locRangeMin, double locRangeMax, TH1D *locInputFluxHist)
{
	//assumes locRangeMax & locRangeMin aren't in the same bin of the input hist (true unless ridiculous)
	double locFlux = 0.0;
	double locFluxUncertainty = 0.0;

	TAxis* locAxis = locInputFluxHist->GetXaxis();
	int locRangeMinBin = locAxis->FindBin(locRangeMin);
	int locRangeMaxBin = locAxis->FindBin(locRangeMax);

	for(int loc_i = locRangeMinBin; loc_i <= locRangeMaxBin; ++loc_i)
	{
		double locBinFlux = locInputFluxHist->GetBinContent(loc_i);
		double locBinFluxError = locInputFluxHist->GetBinError(loc_i);
		if(!(locBinFlux > 0.0))
			continue;
		if(loc_i == locRangeMinBin)
		{
			double locFluxFraction = (locAxis->GetBinUpEdge(loc_i) - locRangeMin)/locAxis->GetBinWidth(loc_i);
			double locNewFluxContribution = locFluxFraction * locBinFlux;
			double locFluxFractionVariance = locNewFluxContribution*(1.0 - locFluxFraction)/(locBinFlux * locBinFlux);
			double locNewFluxContributionVariance = locFluxFraction * locFluxFraction * locBinFluxError * locBinFluxError + locFluxFractionVariance * locBinFlux * locBinFlux;

			locFlux += locNewFluxContribution;
			locFluxUncertainty += locNewFluxContributionVariance;
		}
		else if(loc_i == locRangeMaxBin)
		{
			double locFluxFraction = (locRangeMax - locAxis->GetBinLowEdge(loc_i))/locAxis->GetBinWidth(loc_i);
			double locNewFluxContribution = locFluxFraction * locBinFlux;
			double locFluxFractionVariance = locNewFluxContribution*(1.0 - locFluxFraction)/(locBinFlux * locBinFlux);
			double locNewFluxContributionVariance = locFluxFraction * locFluxFraction * locBinFluxError * locBinFluxError + locFluxFractionVariance * locBinFlux * locBinFlux;

			locFlux += locNewFluxContribution;
			locFluxUncertainty += locNewFluxContributionVariance;
		}
		else
		{
			locFlux += locBinFlux;
			locFluxUncertainty += locBinFluxError * locBinFluxError;
		}
	}
	locFluxUncertainty = sqrt(locFluxUncertainty);
	return pair<double, double>(locFlux, locFluxUncertainty);
}