#include "TFile.h"
#include "TList.h"
#include "TH1.h"
#include "TH2.h"
#include "TKey.h"
#include "TDirectory.h"
#include "TObjArray.h"
#include "TClass.h"

#include "PlotDrawer/DCanvasInstructions.h"
#include "PlotDrawer/DPadInstructions.h"
#include "PlotDrawer/DPlotDrawer.h"
#include "PlotFitter/DPlotFitter.h"
#include "PlotFitter/DFitControl.h"
#include "PlotFitter/DFitFunctors_Miscellaneous.h"
#include "PlotFitter/DFitFunctors_Gaussians.h"

#include <map>
#include <set>
#include <vector>
#include <string>

bool gFinalFlag = false;

void Draw_PreLambda(string locInputFileName);
void Draw_Lambda_FromTree(string locInputFileName);
void Draw_Lambda(string locInputFileName);
void Draw_Phi(string locInputFileName);
void Draw_Sigma0(string locInputFileName);
void Draw_SigmaPlus(string locInputFileName);
void Draw_K0(string locInputFileName);

void Draw_KStarPlus(string locInputFileName);
void Draw_SigmaStar0(string locInputFileName);
void Draw_KStarZero(string locInputFileName);
void Draw_SigmaStarPlus(string locInputFileName);

void Draw_Analyses(string locInputFileName, bool locGroupPlotsFlag = false);
void Draw_Analysis(TDirectory* locReactionDir, TDirectory* locOutputDir, int locMassRebinAmount, bool locGroupPlotsFlag);

void Fit_Lifetime(TH1* locHist, vector<double>& locInitParams, double locFitRangeMin, double locFitRangeMax);

double Calc_InitPeakHeight(TH1* locHist, double locGaussMean, double locPeakEdgeDistance);

DFitFunctor_DoubleGaussian* Build_DoubleGaussian(double locGaussHeight, double locGaussMean, double locGaussWidth, double loc2ndGaussWidth);
DFitFunctor_Voigtian* Build_Voigtian(double locPeakHeight, double locPeakMean, double locGaussWidth, double locBWGamma);
DFitFunctor* Build_PolyBackground(TH1* locHist, double locFitRangeMin, double locFitRangeMax, int locPolyBackPower);
DFitFunctor_Gaussian* Build_Gaussian(double locGaussHeight, double locGaussMean, double locGaussWidth);

bool Fit_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax);
bool Fit_Voigtian(TH1* locHist, double locPeakHeight, double locPeakMean, double locGaussWidth, double locBWGamma, double locFitRangeMin, double locFitRangeMax);
bool Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double loc2ndGaussWidth, double locFitRangeMin, double locFitRangeMax, pair<double, double>& locTotalSigma);

void Draw_All(void)
{
	Draw_Lambda_FromTree("klambda_lifetime_hists.root");
	Draw_Lambda("hd_root.root");
	Draw_Phi("hd_root.root");
	Draw_Sigma0("hd_root.root");
	Draw_SigmaPlus("kpisigmaplus_hists.root");
	Draw_K0("k0pilambda_hists_allfit.root");

	Draw_KStarPlus("kpi0lambda_hists_kstar.root");
	Draw_SigmaStar0("kpi0lambda_hists_sigmastar.root");
	Draw_KStarZero("hd_root.root");
	Draw_SigmaStarPlus("hd_root.root");
}

void Draw_Analyses(string locInputFileName, bool locGroupPlotsFlag)
{
	gROOT->SetBatch(kTRUE);

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	string locLocalInputFileName = locInputFileName.substr(locInputFileName.rfind("/") + 1);
	string locOutputFileName = locLocalInputFileName.substr(0, locLocalInputFileName.rfind(".")) + string(".canvii.root");
	TFile* locOutputFile = new TFile(locOutputFileName.c_str(), "RECREATE");

	// manually define rebinning
	map<string, int> locMassRebinMap; //key is reaction name

	// loop over all keys (reactions) in this directory
	TList* locListOfReactionKeys = locInputFile->GetListOfKeys();
	for(int loc_i = 0; loc_i < locListOfReactionKeys->GetEntries(); ++loc_i)
	{
		//get object
		TKey* locReactionKey = (TKey*)locListOfReactionKeys->At(loc_i);
		cout << "Key " << loc_i << " of " << locListOfReactionKeys->GetEntries() << endl;
		string locReactionKeyName = locReactionKey->GetName();
		TObject* locReactionObject = locInputFile->Get(locReactionKeyName.c_str());

		//check if is reaction directory
		if(!locReactionObject->IsA()->InheritsFrom(TDirectory::Class()))
			continue;
		if(locReactionKeyName == "Independent")
			continue;

		cout << "Reaction Name: " << locReactionKeyName << endl;
/*
		if(locReactionKeyName == "lambdaX")
			continue; //too much background
		if((locReactionKeyName == "omega_missp") || (locReactionKeyName == "omegaX"))
			continue; //too much background
		if((locReactionKeyName == "omega") || (locReactionKeyName == "omega_p4fit"))
			continue; //play with cuts in TTree (beware, conlev not flat!)
		if((locReactionKeyName == "k0sigma") || (locReactionKeyName == "k0sigma_p4fit") || (locReactionKeyName == "k0sigmaX"))
			continue; //WAY too much background: final state is p, pi+, pi-, pi0!!
		if((locReactionKeyName == "ksigma0_p4fit") || (locReactionKeyName == "ksigma0X"))
			continue; //too much background
		if((locReactionKeyName == "k0pi0sigmaplus_p4fit") || (locReactionKeyName == "k0pi0sigmaplusX"))
			continue; //WAY too much background
		if((locReactionKeyName == "k0pisigma0") || (locReactionKeyName == "k0pisigma0_p4fit") || (locReactionKeyName == "k0pisigma0X"))
			continue; //WAY too much background
*/
		TDirectory* locReactionDir = (TDirectory*)locReactionObject;
		TDirectory* locOutputDir = locOutputFile->mkdir(locReactionKeyName.c_str(), locReactionKeyName.c_str());

//		Draw_Analysis(locReactionDir, locOutputDir, locMassRebinMap[locReactionKeyName], locGroupPlotsFlag);
		Draw_Analysis(locReactionDir, locOutputDir, 5, locGroupPlotsFlag);
	} //end reaction loop

	locOutputFile->Close();
}

void Draw_Analysis(TDirectory* locReactionDir, TDirectory* locOutputDir, int locMassRebinAmount, bool locGroupPlotsFlag)
{
	map<string, vector<pair<string, TH1*> > > locMassHistMap; //first key is reaction name, second is base-folder-name, pair string is legend name
	TH1* locEventsHist = NULL; //key is reaction name, value is #-events-survived hist
	TH1* locCombos1DHist = NULL; //key is reaction name, value is #-combos-survived hist
	TH2* locCombos2DHist = NULL; //key is reaction name, value is #-combos-survived hist
	TH1* locKinFitConfidenceHist = NULL; //key is reaction name
	TH1* locBeamRFDeltaTHist = NULL; //key is reaction name
	TH1* locNumBeamHist = NULL; //key is reaction name

	// Loop over action keys
	TList* locListOfActionKeys = locReactionDir->GetListOfKeys();
	set<string> locMassBaseNames;
	for(int loc_j = 0; loc_j < locListOfActionKeys->GetEntries(); ++loc_j)
	{
		//get object
		TKey* locActionKey = (TKey*)locListOfActionKeys->At(loc_j);
		string locActionKeyName = locActionKey->GetName();
		TObject* locActionObject = locReactionDir->Get(locActionKeyName.c_str());

		if(locActionKeyName == "NumEventsSurvivedAction")
			locEventsHist = (TH1*)locActionObject;
		else if(locActionKeyName == "NumCombosSurvivedAction1D")
			locCombos1DHist = (TH1*)locActionObject;
		else if(locActionKeyName == "NumCombosSurvivedAction")
			locCombos2DHist = (TH2*)locActionObject;

		//check if is directory
		if(!locActionObject->IsA()->InheritsFrom(TDirectory::Class()))
			continue;

		TDirectory* locActionDir = (TDirectory*)locActionObject;

		if(locActionKeyName == "Hist_ComboConstruction")
		{
			locBeamRFDeltaTHist = (TH1*)locActionDir->Get("BeamPhotonRFDeltaT");
			locNumBeamHist = (TH1*)locActionDir->Get("NumSurvivingBeamParticles");
		}
		else if(locActionKeyName == "Hist_KinFitResults")
			locKinFitConfidenceHist = (TH1*)locActionDir->Get("ConfidenceLevel");
		else if(locActionKeyName.substr() == "Hist_KinFitResults")
			locKinFitConfidenceHist = (TH1*)locActionDir->Get("ConfidenceLevel");
		else if(locActionKeyName.find("Mass") != string::npos)
		{
			//get hist
			TH1* locHist = NULL;
			if(locActionKeyName.find("Invariant") != string::npos)
				locHist = (TH1*)locActionDir->Get("InvariantMass");
			else if(locActionKeyName.find("Squared") != string::npos)
				locHist = (TH1*)locActionDir->Get("MissingMassSquared");
			else if(locActionKeyName.find("2DInvariant") != string::npos)
				continue;
			else
				locHist = (TH1*)locActionDir->Get("MissingMass");
			if(locHist == NULL)
				continue;

			//search to see if this matches a previous mass action
			set<string>::iterator locIterator = locMassBaseNames.begin();
			string locMatchName = "";
			string locLegendName = "Initial";
			for(; locIterator != locMassBaseNames.end(); ++locIterator)
			{
				string locPreviousMassAction = *locIterator;
				if(locActionKeyName.substr(0, locPreviousMassAction.size()) != locPreviousMassAction)
					continue; //not a match
				locMatchName = locPreviousMassAction;
				locLegendName = locActionKeyName.substr(locPreviousMassAction.size() + 1);
				break;
			}

			//insert hist
			pair<string, TH1*> locHistPair(locLegendName, locHist);
			if(locMatchName != "")
				locMassHistMap[locMatchName].push_back(locHistPair);
			else //new
			{
				locMassHistMap[locActionKeyName].push_back(locHistPair);
				locMassBaseNames.insert(locActionKeyName);
			}
		} //end mass search
	} //end action loop
	locOutputDir->cd();

	//DRAW
	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;
	TObjArray* locCanvasArray = NULL;

	//log-y
	DCanvasInstructions* locLogYInstructions = new DCanvasInstructions("LogY");
	locLogYInstructions->dPadInstructions[0]->dLogYFlag = true;

	//# combos at the end
	string locHistName = string(locCombos2DHist->GetName()) + string("_py");
	TH1* locNumCombosAtEnd = locCombos2DHist->ProjectionY(locHistName.c_str(), locCombos2DHist->GetNbinsX(), locCombos2DHist->GetNbinsX());
	locNumCombosAtEnd->GetXaxis()->SetRangeUser(1, 10);
	locNumCombosAtEnd->GetXaxis()->SetTitle("# Particle Combos Survived Final Action");

	//survival
	TObjArray* locSurvivalArray = new TObjArray(3);
	locSurvivalArray->AddAt(locEventsHist, 0);
	locSurvivalArray->AddAt(locCombos1DHist, 1);
	locSurvivalArray->AddAt(locNumCombosAtEnd, 2);

	if(locGroupPlotsFlag)
		locCanvasArray = locPlotDrawer.Draw_Array(locSurvivalArray, locLogYInstructions);
	else
	{
		for(int loc_j = 0; loc_j < locSurvivalArray->GetEntriesFast(); ++loc_j)
		{
			locCanvas = locPlotDrawer.Draw_Object(locSurvivalArray->At(loc_j), locLogYInstructions);
			locCanvas->Write();
		}
	}

	//kinfit
	if(locKinFitConfidenceHist != NULL)
	{
		locCanvas = locPlotDrawer.Draw_Object(locKinFitConfidenceHist, locLogYInstructions);
		locCanvas->Write();
	}

	//masses
	TObjArray* locResultsArray = new TObjArray();
	map<string, vector<pair<string, TH1*> > >::iterator locMassIterator = locMassHistMap.begin();
	DCanvasInstructions* locMassInstructions = new DCanvasInstructions(locReactionDir->GetName());
	for(; locMassIterator != locMassHistMap.end(); ++locMassIterator)
	{
		vector<pair<string, TH1*> >& locMassHists = locMassIterator->second;

		//SET LEGEND NAMES, SCALE FLAG (optional)
		TObjArray* locMassArray = new TObjArray(locMassHists.size());
		vector<string> locLegendNames;
		for(size_t loc_i = 0; loc_i < locMassHists.size(); ++loc_i)
		{
			locLegendNames.push_back(locMassHists[loc_i].first);
			locMassHists[loc_i].second->Rebin(locMassRebinAmount);
			locMassArray->AddAt(locMassHists[loc_i].second, loc_i);
		}

		int locPadIndex = locGroupPlotsFlag ? locResultsArray->GetEntriesFast() : 0;
		locMassInstructions->Get_PadInstructions(locPadIndex)->dLegendNames = locLegendNames;

		//SUPERIMPOSE HISTOGRAMS:
		if(!locGroupPlotsFlag)
		{
			locCanvas = locPlotDrawer.Draw_Superimpose(locMassArray, locMassInstructions);
			locCanvas->Write();
		}
		else
			locResultsArray->AddLast(locMassArray);
	}
	if(locGroupPlotsFlag)
		locCanvasArray = locPlotDrawer.Draw_2DArray(locResultsArray, locMassInstructions);
}

void Fit_Lifetime(TH1* locHist, vector<double>& locInitParams, double locFitRangeMin, double locFitRangeMax)
{
	DFitFunctor_Lifetime* locFitFunctor = new DFitFunctor_Lifetime(locInitParams);
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor);

	DPlotFitter locPlotFitter;
	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);
	cout << "hist name, fit status = " << locHist->GetName() << ", " << locFitStatus << endl;
}

void Draw_Phi(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//log-y
	DCanvasInstructions* locLogYInstructions = new DCanvasInstructions("LogY");
	locLogYInstructions->dPadInstructions[0]->dLogYFlag = true;

	//Confidence Level
	locHist = (TH1*)locInputFile->Get("phi/Hist_KinFitResults/ConfidenceLevel");
	locHist->GetYaxis()->SetTitle("");
	locHist->Rebin(4);
	locCanvas = locPlotDrawer.Draw_Object(locHist, locLogYInstructions);
//	locCanvas->Write();

	//Mass
	locHist = (TH1*)locInputFile->Get("phi/Hist_InvariantMass_phiMeson_TightKinFitCutKinFit/InvariantMass");
	locHist->Rebin(8);
	locHist->SetTitle("#gammap#rightarrow#phip");
	locHist->GetXaxis()->SetRangeUser(0.96, 1.08);
	locHist->GetYaxis()->SetTitle("");
	if(!gFinalFlag)
	{
		double locPeakMean = 1.019;
		double locGaussWidth = 0.002;
		double locPeakHeight = Calc_InitPeakHeight(locHist, locPeakMean, 5.0*locGaussWidth);
		double locBWGamma = 0.0043;
		double locFitRangeMin = 1.0;
		double locFitRangeMax = 1.055;
		Fit_Voigtian(locHist, locPeakHeight, locPeakMean, locGaussWidth, locBWGamma, locFitRangeMin, locFitRangeMax);
	}
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_Sigma0(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//log-y
	DCanvasInstructions* locLogYInstructions = new DCanvasInstructions("LogY");
	locLogYInstructions->dPadInstructions[0]->dLogYFlag = true;

	//Confidence Level
	locHist = (TH1*)locInputFile->Get("ksigma0/Hist_KinFitResults/ConfidenceLevel");
	locHist->GetYaxis()->SetTitle("");
	locHist->Rebin(4);
	locCanvas = locPlotDrawer.Draw_Object(locHist, locLogYInstructions);
//	locCanvas->Write();

	//Mass
	locHist = (TH1*)locInputFile->Get("ksigma0/Hist_InvariantMass_Sigma0_TightKinFitCutKinFit/InvariantMass");
	locHist->Rebin(8);
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Sigma^{0}");
	locHist->GetYaxis()->SetTitle("");
	if(!gFinalFlag)
	{
		double locGaussMean = 1.192;
		double locGaussWidth = 0.002;
		double locGaussHeight = Calc_InitPeakHeight(locHist, locGaussMean, 0.03);
		double locFitRangeMin = 1.16;
		double locFitRangeMax = 1.25;
		pair<double, double> locTotalSigma;
		Fit_DoubleGauss(locHist, locGaussHeight, locGaussMean, locGaussWidth, 2.0*locGaussWidth, locFitRangeMin, locFitRangeMax, locTotalSigma);
		cout << "total sigma, error = " << locTotalSigma.first << ", " << locTotalSigma.second << endl;
	}
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_K0(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//Mass
	locHist = (TH1*)locInputFile->Get("KZeroMass_KinFit");
	locHist->SetTitle("#gammap#rightarrowK^{0}#pi^{#plus}#Lambda");
	locHist->Rebin(10);
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_SigmaPlus(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//Mass
	locHist = (TH1*)locInputFile->Get("SigmaPlus_KinFit");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#pi^{#minus}#Sigma^{+}");
	locHist->Rebin(8);
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_KStarZero(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//Mass
	locHist = (TH1*)locInputFile->Get("k2pilambda/Hist_InvariantMass_KStar0_KinFitCutKinFit/InvariantMass");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#pi^{#minus}#pi^{#plus}#Lambda");
	locHist->Rebin(15);
	locHist->GetYaxis()->SetTitle("");
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_SigmaStarPlus(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//Mass
	locHist = (TH1*)locInputFile->Get("k2pilambda_p4fit/Hist_InvariantMass_SigmaStarPlus_TightKinFitCutKinFit/InvariantMass");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#pi^{#minus}#pi^{#plus}#Lambda");
	locHist->Rebin(15);
	locHist->GetYaxis()->SetTitle("");
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_KStarPlus(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	DPlotDrawer locPlotDrawer;

	//Mass
	TH1* locHist = (TH1*)locInputFile->Get("KStarMass_OverlapCut_KinFit");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#pi^{0}#Lambda");
	locHist->Rebin(12);
	TCanvas* locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_SigmaStar0(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	DPlotDrawer locPlotDrawer;

	//Mass
	TH1* locHist = (TH1*)locInputFile->Get("SigmaStarMass_OverlapCut_KinFit");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#pi^{0}#Lambda");
	locHist->Rebin(15);
	TCanvas* locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_PreLambda(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//log-y, wide-screen
	DCanvasInstructions* locLogYInstructions = new DCanvasInstructions("LogY");
	locLogYInstructions->dPadInstructions[0]->dLogYFlag = true;

	//Missing Mass Sq
	locHist = (TH1*)locInputFile->Get("MissingMassSquared");
	locHist->Rebin(2); //or 4
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	locHist->GetYaxis()->SetTitle("");
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();

	//Lambda Mass
	locHist = (TH1*)locInputFile->Get("LambdaMass");
	locHist->Rebin(2); //or 4
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	locHist->GetXaxis()->SetRangeUser(1.08, 1.15);
	locHist->GetYaxis()->SetTitle("");
	locCanvas = locPlotDrawer.Draw_Object(locHist);
}

void Draw_Lambda(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//log-y, wide-screen
	DCanvasInstructions* locLogYInstructions = new DCanvasInstructions("LogY");
	locLogYInstructions->dPadInstructions[0]->dLogYFlag = true;

	//Confidence Level
	locHist = (TH1*)locInputFile->Get("klambda_nomassfit/Hist_KinFitResults/ConfidenceLevel");
	locHist->GetYaxis()->SetTitle("");
	locHist->Rebin(4);
	locCanvas = locPlotDrawer.Draw_Object(locHist, locLogYInstructions);
//	locCanvas->Write();

	//Lambda Mass
	locHist = (TH1*)locInputFile->Get("klambda_nomassfit/Hist_InvariantMass_Lambda_TightKinFitCutKinFit/InvariantMass");
	locHist->Rebin(2); //or 4
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	locHist->GetXaxis()->SetRangeUser(1.08, 1.15);
	locHist->GetYaxis()->SetTitle("");
	if(!gFinalFlag)
	{
		double locGaussMean = 1.116;
		double locGaussWidth = 0.002;
		double locGaussHeight = Calc_InitPeakHeight(locHist, locGaussMean, 4.0*locGaussWidth);
		double locFitRangeMin = 1.095;
		double locFitRangeMax = 1.135;
		pair<double, double> locTotalSigma;
		Fit_DoubleGauss(locHist, locGaussHeight, locGaussMean, locGaussWidth, 2.0*locGaussWidth, locFitRangeMin, locFitRangeMax, locTotalSigma);
		cout << "total sigma, error = " << locTotalSigma.first << ", " << locTotalSigma.second << endl;
	}
	locCanvas = locPlotDrawer.Draw_Object(locHist);
//	locCanvas->Write();
}

void Draw_Lambda_FromTree(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	vector<double> locInitParams; // [0]*e^([1]*x);
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;

	DPlotDrawer locPlotDrawer;
	TCanvas* locCanvas = NULL;

	//log-y, wide-screen
	DCanvasInstructions* locLogYInstructions = new DCanvasInstructions("LogY");
	locLogYInstructions->dPadInstructions[0]->dLogYFlag = true;

	//VertexZ
	DCanvasInstructions* locVertexZInstructions = new DCanvasInstructions("VertexZ");
	locVertexZInstructions->dPadInstructions[0]->dLegendNames.push_back("#Lambda Production Vertex");
	locVertexZInstructions->dPadInstructions[0]->dLegendNames.push_back("#Lambda Decay Vertex");
	locVertexZInstructions->dPadInstructions[0]->dHistFillColorFlag = false;
	locVertexZInstructions->dStyle->SetHistLineWidth(2);

	//Vertex-Z
	//_NoPIDHit
	TObjArray* locVertZArray = new TObjArray(2);
	locHist = (TH1*)locInputFile->Get("ProductionVertexZ");
	locHist->Rebin(2);
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	locHist->GetXaxis()->SetTitle("Vertex-Z (cm)");
	locHist->GetXaxis()->SetRangeUser(20.0, 120.0);
	locVertZArray->AddAt(locHist, 0);
	locHist = (TH1*)locInputFile->Get("LambdaDecayVertexZ");
	locHist->Rebin(2);
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	locHist->GetXaxis()->SetTitle("Vertex-Z (cm)");
	locHist->GetXaxis()->SetRangeUser(20.0, 120.0);
	locVertZArray->AddAt(locHist, 1);
	locCanvas = locPlotDrawer.Draw_Superimpose(locVertZArray, locVertexZInstructions);
//	locCanvas->Write();

	//Lifetime
	locHist = (TH1*)locInputFile->Get("LambdaLifetime");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	locHist->Rebin(4);
	locHist->GetXaxis()->SetRangeUser(0.0, 2.0);
	locInitParams.clear();
	locInitParams.resize(2);
	locInitParams[0] = locHist->GetMaximum();  locInitParams[1] = 0.26;
	Fit_Lifetime(locHist, locInitParams, 0.01, 1.2);
	locCanvas = locPlotDrawer.Draw_Object(locHist);
}

double Calc_InitPeakHeight(TH1* locHist, double locGaussMean, double locPeakEdgeDistance)
{
	int locLowSideBin = locHist->GetXaxis()->FindBin(locGaussMean - locPeakEdgeDistance);
	int locHighSideBin = locHist->GetXaxis()->FindBin(locGaussMean + locPeakEdgeDistance);
	double locAvgBackground = 0.5*(locHist->GetBinContent(locLowSideBin) + locHist->GetBinContent(locHighSideBin));

	int locPeakBin = locHist->GetXaxis()->FindBin(locGaussMean);
	return locHist->GetBinContent(locPeakBin) - locAvgBackground;
}

bool Fit_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax)
{
	//Build Functors
	DFitFunctor_Gaussian* locFitFunctor_Gaussian = Build_Gaussian(locGaussHeight, locGaussMean, locGaussWidth);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground(locHist, locFitRangeMin, locFitRangeMax, 2);

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_Gaussian, locFitFunctor_Background);
	DPlotFitter locPlotFitter;
	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);

	DFitControl::DYieldData locYieldData;
	locFitControl->Calc_Yield(locHist, false, true, locYieldData);
cout << "yield = " << locYieldData.dYield << endl;

	return true;
}

bool Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double loc2ndGaussWidth, double locFitRangeMin, double locFitRangeMax, pair<double, double>& locTotalSigma)
{
	//Build Functors
	DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = Build_DoubleGaussian(locGaussHeight, locGaussMean, locGaussWidth, loc2ndGaussWidth);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground(locHist, locFitRangeMin, locFitRangeMax, 1);

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_DoubleGaussian, locFitFunctor_Background);
	DPlotFitter locPlotFitter;
	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);

	DFitControl::DYieldData locYieldData;
	locFitControl->Calc_Yield(locHist, false, true, locYieldData);
cout << "yield = " << locYieldData.dYield << endl;

	locTotalSigma = locFitFunctor_DoubleGaussian->Calc_StandardDeviation();
	return true;
}

bool Fit_Voigtian(TH1* locHist, double locPeakHeight, double locPeakMean, double locGaussWidth, double locBWGamma, double locFitRangeMin, double locFitRangeMax)
{
	//Build Functors
	DFitFunctor_Voigtian* locFitFunctor_Voigtian = Build_Voigtian(locPeakHeight, locPeakMean, locGaussWidth, locBWGamma);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground(locHist, locFitRangeMin, locFitRangeMax, 1);

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_Voigtian, locFitFunctor_Background);
	DPlotFitter locPlotFitter;
	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);

	DFitControl::DYieldData locYieldData;
	locFitControl->Calc_Yield(locHist, false, true, locYieldData);
cout << "yield = " << locYieldData.dYield << endl;

	return true;
}

DFitFunctor_Voigtian* Build_Voigtian(double locPeakHeight, double locPeakMean, double locGaussWidth, double locBWGamma)
{
	vector<double> locInitParams(4, 0.0);
	locInitParams[0] = locPeakHeight;  locInitParams[1] = locPeakMean;
	locInitParams[2] = locGaussWidth;  locInitParams[3] = locBWGamma;
	DFitFunctor_Voigtian* locVoigtian = new DFitFunctor_Voigtian(locInitParams, kBlue);
	locVoigtian->dParamLimits[2] = pair<double, double>(0.0005, 0.5);
	locVoigtian->dParamLimits[3] = pair<double, double>(0.0005, 0.5);
	return locVoigtian;
}

DFitFunctor_DoubleGaussian* Build_DoubleGaussian(double locGaussHeight, double locGaussMean, double locGaussWidth, double loc2ndGaussWidth)
{
	vector<double> locInitParams(5, 0.0);
	//peak gaussian
	locInitParams[0] = 0.925*locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = 0.9*locGaussWidth;
	//non-gaussian-tail gaussian
	locInitParams[3] = 0.075*locGaussHeight;  locInitParams[4] = loc2ndGaussWidth;
	DFitFunctor_DoubleGaussian* locDoubleGaussian = new DFitFunctor_DoubleGaussian(locInitParams, kBlue);
	locDoubleGaussian->dParamLimits[2] = pair<double, double>(0.0005, 0.5);
	locDoubleGaussian->dParamLimits[4] = pair<double, double>(0.0005, 0.5);
	return locDoubleGaussian;
}

DFitFunctor_Gaussian* Build_Gaussian(double locGaussHeight, double locGaussMean, double locGaussWidth)
{
	vector<double> locInitParams(3, 0.0);
	//peak gaussian
	locInitParams[0] = locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = locGaussWidth;
	DFitFunctor_Gaussian* locGaussian = new DFitFunctor_Gaussian(locInitParams, kBlue);
	locGaussian->dParamLimits[2] = pair<double, double>(0.0005, 0.5);
	return locGaussian;
}

DFitFunctor* Build_PolyBackground(TH1* locHist, double locFitRangeMin, double locFitRangeMax, int locPolyBackPower)
{
	if(locPolyBackPower == 0)
	{
		double locIntercept = 0.0;
		DFitUtilities::Guess_Params_Flat(locHist, locFitRangeMin, locFitRangeMax, locIntercept);
		vector<double> locInitParams(1, locIntercept);
		return (new DFitFunctor_Polynomial(locInitParams, kMagenta));
	}
	else if(locPolyBackPower == 1)
	{
		double locSlope = 0.0, locIntercept = 0.0;
		DFitUtilities::Guess_Params_Linear(locHist, locFitRangeMin, locFitRangeMax, locIntercept, locSlope, true);

		vector<double> locInitParams;
		locInitParams.resize(2);
		locInitParams[0] = locIntercept;  locInitParams[1] = locSlope;

		return (new DFitFunctor_LinearAboveZero(locInitParams, locFitRangeMin, locFitRangeMax, kMagenta));
	}

	double locSlope = 0.0, locIntercept = 0.0;
	DFitUtilities::Guess_Params_Linear(locHist, locFitRangeMin, locFitRangeMax, locIntercept, locSlope, true);

	vector<double> locInitParams(locPolyBackPower + 1, 0.0);
	locInitParams[0] = locIntercept;  locInitParams[1] = locSlope;

	return (new DFitFunctor_Polynomial(locInitParams, kMagenta));
}