#include "DFitUtilities.h"

ClassImp(DFitUtilities)

using namespace std;

double DFitUtilities::Calc_EntriesInRange(const TH1* locHist, double locRangeMin, double locRangeMax, double& locAverageBinWidth)
{
	//get edge bins
	unsigned int locRangeMinBin = locHist->GetXaxis()->FindBin(locRangeMin);
	unsigned int locRangeMaxBin = locHist->GetXaxis()->FindBin(locRangeMax);

	//compute avg bin width
	double locLowBinLowEdge = locHist->GetXaxis()->GetBinLowEdge(locRangeMinBin);
	double locHighBinHighEdge = locHist->GetXaxis()->GetBinUpEdge(locRangeMaxBin);
	locAverageBinWidth = (locHighBinHighEdge - locLowBinLowEdge)/(double(locRangeMaxBin - locRangeMinBin + 1));

	//compute #events
	double locNumEvents = locHist->Integral(locRangeMinBin, locRangeMaxBin);

	//subtract partial counts in first & last bins
	//first bin
	double locEdgeBinCounts = locHist->GetBinContent(locRangeMinBin);
	double locExcludeFraction = (locRangeMin - locLowBinLowEdge)/locHist->GetXaxis()->GetBinWidth(locRangeMinBin);
	locNumEvents -= locEdgeBinCounts*locExcludeFraction;
	//last bin
	locEdgeBinCounts = locHist->GetBinContent(locRangeMaxBin);
	locExcludeFraction = (locHighBinHighEdge - locRangeMax)/locHist->GetXaxis()->GetBinWidth(locRangeMaxBin);
	locNumEvents -= locEdgeBinCounts*locExcludeFraction;

	return locNumEvents;
}

void DFitUtilities::Calc_CutRange_SignalFraction(TF1* locFunc, double locSignalFraction, double& locCutRangeMin, double& locCutRangeMax, double locMean, double locSigma)
{
	double locCutDistanceFromCenter = 3.0*locSignalFraction*locSigma; //init guess at input % of 3 sigma

	double locRangeMin, locRangeMax;
	locFunc->GetRange(locRangeMin, locRangeMax);
	double locTotalRange = locRangeMax - locRangeMin;

	double locTotalIntegralRangeMin = locRangeMin - locTotalRange;
	double locTotalIntegralRangeMax = locRangeMax + locTotalRange;
	if((3.0*locSigma + locMean) > locTotalIntegralRangeMax)
	{
		locTotalIntegralRangeMax = locMean + 4.0*locSigma;
		locTotalIntegralRangeMin = locMean - 4.0*locSigma;
	}

	locFunc->SetRange(locTotalIntegralRangeMin, locTotalIntegralRangeMax); //temp, not sure if needed
	double locTotalSignalIntegral = locFunc->Integral(locTotalIntegralRangeMin, locTotalIntegralRangeMax);

	//the cut distance must be greater than this initial value, because this is the desired value for only one gaussian, and we have 2

	double locDistanceStep = 0.1*locCutDistanceFromCenter;
	if(locDistanceStep < 0.01)
		locDistanceStep = 0.01;
	if(locCutDistanceFromCenter < 0.005)
		locCutDistanceFromCenter = 0.005;
	bool locPreviouslyTooLowFlag = true;
	double locTolerance = 1.0E-12;
	double locIntegralFraction = 0.0;

//ostringstream locStream;
//locStream << "values = " << locRangeMin << ", " << locRangeMax << ", " << locTotalIntegralRangeMin << ", " << locTotalIntegralRangeMax << ", " << locMean << ", " << locCutDistanceFromCenter << ", " << locTotalSignalIntegral << ", " << locSigma1 << ", " << locSigma2 << ", " << locHeight1 << ", " << locHeight2 << ", " << locDistanceStep;
//gProofServ->SendAsynMessage(TString(locStream.str()));

	//vary locCutDistanceFromCenter until either the integral is within the %-tolerance of the desired amount, or until the max # of iterations
		//step size is decreased 
	unsigned int locLoopCounter = 0;
	do
	{
		double locIntegral = locFunc->Integral(locMean - locCutDistanceFromCenter, locMean + locCutDistanceFromCenter);
		locIntegralFraction = locIntegral/locTotalSignalIntegral;
		++locLoopCounter;
		if(locLoopCounter == 200)
		{
			 locCutDistanceFromCenter = 9.9E9;
			 break;
		}

/*
if(locCounter < 100)
{
locStream.str("");
locStream << "frac, dist, step = " << locIntegralFraction << ", " << locCutDistanceFromCenter << ", " << locDistanceStep;
gProofServ->SendAsynMessage(TString(locStream.str()));
}
++locCounter;
*/
		//integral faction is too large, but was previously too low: have just crossed the value: decrease step size
		if(locIntegralFraction > locSignalFraction)
		{
			if(locPreviouslyTooLowFlag)
			{
				locDistanceStep /= 2.0;
				locPreviouslyTooLowFlag = false;
			}
			locCutDistanceFromCenter -= locDistanceStep;
		}
		else if(locIntegralFraction < locSignalFraction)
		{
			if(!locPreviouslyTooLowFlag)
			{
				locDistanceStep /= 2.0;
				locPreviouslyTooLowFlag = true;
			}
			locCutDistanceFromCenter += locDistanceStep;
		}
		else
			break;
		locIntegral = locFunc->Integral(locMean - locCutDistanceFromCenter, locMean + locCutDistanceFromCenter);
	}
	while(fabs(1.0 - locIntegralFraction/locSignalFraction) > locTolerance);

	locCutRangeMin = locMean - locCutDistanceFromCenter;
	locCutRangeMax = locMean + locCutDistanceFromCenter;
	locFunc->SetRange(locRangeMin, locRangeMax); //restore
}

void DFitUtilities::Fit_Slices_Gaussian(TH2 *loc2DHist, TH1*& locMeanHist, TH1*& locSigmaHist, TH1*& locChiSqHist)
{
	string locHistName, locHistTitle;
	TH1 *loc1DHist;

	loc2DHist->FitSlicesY();

	locHistName = loc2DHist->GetName();
	locHistName += "_0";
	loc1DHist = (TH1*)gDirectory->Get(locHistName.c_str());
	delete loc1DHist;

	locHistName = loc2DHist->GetName();
	locHistName += "_1";
	locMeanHist = (TH1*)gDirectory->Get(locHistName.c_str());
	locHistTitle = loc2DHist->GetTitle();
	locHistTitle += " Mean";
	locMeanHist->SetTitle(locHistTitle.c_str());
	locMeanHist->GetXaxis()->SetTitle(loc2DHist->GetXaxis()->GetTitle());
	locMeanHist->GetYaxis()->SetTitle(loc2DHist->GetYaxis()->GetTitle());
	locMeanHist->SetEntries(loc2DHist->GetEntries());

	locHistName = loc2DHist->GetName();
	locHistName += "_2";
	locSigmaHist = (TH1*)gDirectory->Get(locHistName.c_str());
	locHistTitle = loc2DHist->GetTitle();
	locHistTitle += " Sigma";
	locSigmaHist->SetTitle(locHistTitle.c_str());
	locSigmaHist->GetXaxis()->SetTitle(loc2DHist->GetXaxis()->GetTitle());
	locSigmaHist->GetYaxis()->SetTitle(loc2DHist->GetYaxis()->GetTitle());
	locSigmaHist->SetEntries(loc2DHist->GetEntries());

	locHistName = loc2DHist->GetName();
	locHistName += "_chi2";
	locChiSqHist = (TH1*)gDirectory->Get(locHistName.c_str());
	locHistTitle = loc2DHist->GetTitle();
	locHistTitle += " ChiSq / NDof";
	locChiSqHist->SetTitle(locHistTitle.c_str());
	locChiSqHist->SetEntries(loc2DHist->GetEntries());
}