#include <string>

#include "TFile.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.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"

using namespace std;

double Calc_InitGaussHeight(TH1* locMassHist, double locGaussMean, double locGaussWidth, double locNumSigmas);
DFitFunctor* Build_PolyBackground_Linear(TH1* locMassHist, double locFitRangeMin, double locFitRangeMax);
DFitFunctor* Build_PolyBackground_Quadratic(TH1* locMassHist, double locFitRangeMin, double locFitRangeMax);
DFitFunctor_DoubleGaussian* Build_DoubleGaussian(double locGaussHeight, double locGaussMean, double locGaussWidth);
DFitControl* Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);

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

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_MMVsBeamE_New("very_final/cascade_mm.root");
Draw_MMVsBeamE_New("very_final/cascade_mm_lambda.root");
Draw_MMVsBeamE_New("very_final/cascade_inclusive.root");
*/
void Draw_MMVsBeamE_New(string locInputFileName, bool locOffKPlusFlag = false)
{
	DPlotDrawer locPlotDrawer;
	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("MM");
	locCanvasInstructions->dStyle->SetPadLeftMargin(0.12);
	locCanvasInstructions->dStyle->SetTitleOffset(1.05, "y");

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH2 *locMissingMass_Center, *locMissingMass_Sidebands;
	if(locOffKPlusFlag)
	{
		locMissingMass_Center = (TH2*)locInputFile->Get("MissMassOffKPluses_BestCombo_RFSignal");
		locMissingMass_Sidebands = (TH2*)locInputFile->Get("MissMassOffKPluses_BestCombo_RFSideband");
	}
	else
	{
		locMissingMass_Center = (TH2*)locInputFile->Get("Hist_MissingMass_PostKinFitCut_KinFit_Center/MissingMassVsBeamE");
		locMissingMass_Sidebands = (TH2*)locInputFile->Get("Hist_MissingMass_PostKinFitCut_KinFit_Sideband/MissingMassVsBeamE");
	}
	TH2* locMissingMass_Subtracted = (TH2*)locMissingMass_Center->Clone("Clone");

	int locRebinX = 2, locRebinY = 10; //cascade_mm_lambda: 2, 10, showprojy with 100 bins  //cascade_mm: hopeless //cascade_inclusive: 2, 15, showprojy with 150 bins (E-6)
	locMissingMass_Center->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Sidebands->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Subtracted->Rebin2D(locRebinX, locRebinY);

	locMissingMass_Subtracted->Add(locMissingMass_Sidebands, -0.5);
	locMissingMass_Subtracted->GetYaxis()->SetRangeUser(1.15, 3.0);
	locPlotDrawer.Draw_Object(locMissingMass_Subtracted);

	double locMinBeamE = 3.0;
	double locMaxBeamE = 8.0;
	int locMinProjBin = locMissingMass_Subtracted->GetXaxis()->FindBin(locMinBeamE);
	int locMaxProjBin = locMissingMass_Subtracted->GetXaxis()->FindBin(locMaxBeamE);
	TH1* loc1DMissingMassHist = locMissingMass_Subtracted->ProjectionY("MassHist_ProjY", locMinProjBin, locMaxProjBin);
	loc1DMissingMassHist->GetYaxis()->SetTitle("Counts / 20 MeV/c^{2}");
	if(locInputFileName == "very_final/cascade_mm.root")
		loc1DMissingMassHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}(X)");
	else if(locInputFileName == "very_final/cascade_mm_lambda.root")
		loc1DMissingMassHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Lambda(X)");
	else if(locInputFileName == "very_final/cascade_inclusive.root")
	{
		loc1DMissingMassHist->Rebin(2);
		loc1DMissingMassHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}(X)");
		loc1DMissingMassHist->GetYaxis()->SetTitle("Counts / 40 MeV/c^{2}");
	}
	locPlotDrawer.Draw_Object(loc1DMissingMassHist, locCanvasInstructions);
//	locMissingMass_Subtracted->Draw("COLZ");
}

/*
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_Dalitz("final/kpi0sigma0.root");
Draw_Dalitz("final/kpisigmaplus.root");
*/
void Draw_Dalitz(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;
	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Dalitz");
	locCanvasInstructions->dStyle->SetOptFit(102);
	locCanvasInstructions->dStyle->SetOptStat(10);
	locCanvasInstructions->dStyle->SetPadLeftMargin(0.10);
	locCanvasInstructions->dPadLargeRightMargin = 0.1;

	TH2* locDalitzHist = (TH2*)locInputFile->Get("Hist_Dalitz_KinFit_PostKinFitCut/Dalitz");
	locDalitzHist->Rebin2D(2, 2);
//	locDalitzHist->GetXaxis()->SetTitle("K^{#plus}#pi^{0} Invariant Mass^{2} (GeV/c^{2})^{2}");
//	locDalitzHist->GetYaxis()->SetTitle("#pi^{0}#Sigma^{0} Invariant Mass^{2} (GeV/c^{2})^{2}");
	locDalitzHist->GetXaxis()->SetTitle("K^{#plus}#pi^{#minus} Invariant Mass^{2} (GeV/c^{2})^{2}");
	locDalitzHist->GetYaxis()->SetTitle("#pi^{#minus}#Sigma^{#plus} Invariant Mass^{2} (GeV/c^{2})^{2}");
	locPlotDrawer.Draw_Object(locDalitzHist, locCanvasInstructions);
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_MMVsConLev("analysis/cascade_mm.root");
Draw_MMVsConLev("analysis/cascade_mm_lambda.root");
Draw_MMVsConLev("analysis/cascade_inclusive.root");
*/
void Draw_MMVsConLev(string locInputFileName)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH2* locMissingMass_Center = (TH2*)locInputFile->Get("Hist_MissingMass_KinFit_Center/MissingMassVsConfidenceLevel_LogX");
	TH2* locMissingMass_Sidebands = (TH2*)locInputFile->Get("Hist_MissingMass_KinFit_Sideband/MissingMassVsConfidenceLevel_LogX");
	TH2* locMissingMass_Subtracted = (TH2*)locMissingMass_Center->Clone("Clone");

	int locRebinX = 9, locRebinY = 10;
	locMissingMass_Center->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Sidebands->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Subtracted->Rebin2D(locRebinX, locRebinY);

	locMissingMass_Subtracted->Add(locMissingMass_Sidebands, -0.5);
	locMissingMass_Center->Draw("COLZ");
	new TCanvas();
	locMissingMass_Sidebands->Draw("COLZ");
	new TCanvas();
	locMissingMass_Subtracted->Draw("COLZ");
}


/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_MMVsBeamE("final/cascade_mm.root");
Draw_MMVsBeamE("final/cascade_mm_lambda.root");
Draw_MMVsBeamE("final/cascade_inclusive_E6.root");
Draw_MMVsBeamE("final/klambdaX.root", true);
*/
void Draw_MMVsBeamE(string locInputFileName, bool locOffKPlusFlag = false)
{
	DPlotDrawer locPlotDrawer;
	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("MM");
	locCanvasInstructions->dStyle->SetPadLeftMargin(0.12);
	locCanvasInstructions->dStyle->SetTitleOffset(1.05, "y");

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH2 *locMissingMass_Center, *locMissingMass_Sidebands;
	if(locOffKPlusFlag)
	{
		locMissingMass_Center = (TH2*)locInputFile->Get("Hist_MissingMass_OffKPlus_KinFit_PostKinFitCut_Center/MissingMassVsBeamE");
		locMissingMass_Sidebands = (TH2*)locInputFile->Get("Hist_MissingMass_OffKPlus_KinFit_PostKinFitCut_Sideband/MissingMassVsBeamE");
	}
	else
	{
		locMissingMass_Center = (TH2*)locInputFile->Get("Hist_MissingMass_PostKinFitCut_KinFit_Center/MissingMassVsBeamE");
		locMissingMass_Sidebands = (TH2*)locInputFile->Get("Hist_MissingMass_PostKinFitCut_KinFit_Sideband/MissingMassVsBeamE");
	}
	TH2* locMissingMass_Subtracted = (TH2*)locMissingMass_Center->Clone("Clone");

	int locRebinX = 2, locRebinY = 10; //cascade_mm_lambda: 2, 10, showprojy with 100 bins  //cascade_mm: hopeless //cascade_inclusive: 2, 15, showprojy with 150 bins (E-6)
	locMissingMass_Center->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Sidebands->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Subtracted->Rebin2D(locRebinX, locRebinY);

	locMissingMass_Subtracted->Add(locMissingMass_Sidebands, -0.5);
	locPlotDrawer.Draw_Object(locMissingMass_Subtracted);

	double locMinBeamE = 3.0;
	double locMaxBeamE = 8.0;
	int locMinProjBin = locMissingMass_Subtracted->GetXaxis()->FindBin(locMinBeamE);
	int locMaxProjBin = locMissingMass_Subtracted->GetXaxis()->FindBin(locMaxBeamE);
	TH1* loc1DMissingMassHist = locMissingMass_Subtracted->ProjectionY("MassHist_RFSignal_PARA", locMinProjBin, locMaxProjBin);
	loc1DMissingMassHist->GetYaxis()->SetTitle("Counts / 20 MeV/c^{2}");
	if(locInputFileName == "final/cascade_mm.root")
		loc1DMissingMassHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}(X)");
	else if(locInputFileName == "final/cascade_mm_lambda.root")
		loc1DMissingMassHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Lambda(X)");
	else if(locInputFileName == "final/cascade_inclusive_E6.root")
	{
		loc1DMissingMassHist->Rebin(2);
		loc1DMissingMassHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}(X)");
		loc1DMissingMassHist->GetYaxis()->SetTitle("Counts / 40 MeV/c^{2}");
	}
	locPlotDrawer.Draw_Object(loc1DMissingMassHist, locCanvasInstructions);
//	locMissingMass_Subtracted->Draw("COLZ");
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_MM("analysis/cascade_mm.root");
Draw_MM("analysis/cascade_mm_lambda.root");
Draw_MM("analysis/cascade_inclusive.root");
Draw_MM("final/klambdaX.root", true);
*/
void Draw_MM(string locInputFileName, bool locOffKPlusFlag = false)
{
	DPlotDrawer locPlotDrawer;
	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("MM");
	locCanvasInstructions->dStyle->SetTitleOffset(1.3, "y");

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1 *locMissingMass_Center, *locMissingMass_Sidebands;
	if(locOffKPlusFlag)
	{
		locMissingMass_Center = (TH1*)locInputFile->Get("Hist_MissingMass_OffKPlus_KinFit_PostKinFitCut_Center/MissingMass");
		locMissingMass_Sidebands = (TH1*)locInputFile->Get("Hist_MissingMass_OffKPlus_KinFit_PostKinFitCut_Sideband/MissingMass");
	}
	else
	{
		locMissingMass_Center = (TH1*)locInputFile->Get("Hist_MissingMass_PostKinFitCut_KinFit_Center/MissingMass");
		locMissingMass_Sidebands = (TH1*)locInputFile->Get("Hist_MissingMass_PostKinFitCut_KinFit_Sideband/MissingMass");
	}
	TH1* locMissingMass_Subtracted = (TH1*)locMissingMass_Center->Clone("Clone");

	int locRebinAmount = 4;
	locMissingMass_Center->Rebin(locRebinAmount);
	locMissingMass_Sidebands->Rebin(locRebinAmount);
	locMissingMass_Subtracted->Rebin(locRebinAmount);

	locMissingMass_Subtracted->Add(locMissingMass_Sidebands, -0.5);

	locPlotDrawer.Draw_Object(locMissingMass_Subtracted, locCanvasInstructions);
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_MMVsConLev("analysis/cascade_mm.root", false);
Draw_MMVsConLev("analysis/cascade_mm_lambda.root", false);
Draw_MMVsConLev("analysis/cascade_inclusive.root", false);
Draw_MMVsConLev("klambdaX.root", true);
*/
void Draw_MMVsConLev(string locInputFileName, bool locOffKPlusFlag = false)
{
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH2 *locMissingMass_Center, *locMissingMass_Sidebands;
	if(locOffKPlusFlag)
	{
		locMissingMass_Center = (TH2*)locInputFile->Get("Hist_MissingMass_OffKPlus_KinFit_Center/MissingMassVsConfidenceLevel_LogX");
		locMissingMass_Sidebands = (TH2*)locInputFile->Get("Hist_MissingMass_OffKPlus_KinFit_Sideband/MissingMassVsConfidenceLevel_LogX");
	}
	else
	{
		locMissingMass_Center = (TH2*)locInputFile->Get("Hist_MissingMass_KinFit_Center/MissingMassVsConfidenceLevel_LogX");
		locMissingMass_Sidebands = (TH2*)locInputFile->Get("Hist_MissingMass_KinFit_Sideband/MissingMassVsConfidenceLevel_LogX");
	}
	TH2* locMissingMass_Subtracted = (TH2*)locMissingMass_Center->Clone("Clone");

	int locRebinX = 9, locRebinY = 4;
	locMissingMass_Center->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Sidebands->Rebin2D(locRebinX, locRebinY);
	locMissingMass_Subtracted->Rebin2D(locRebinX, locRebinY);

	locMissingMass_Subtracted->Add(locMissingMass_Sidebands, -0.5);
	locMissingMass_Center->Draw("COLZ");
	new TCanvas();
	locMissingMass_Sidebands->Draw("COLZ");
	new TCanvas();
	locMissingMass_Subtracted->Draw("COLZ");
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_Cascade_ConLev("final/cascade_vertp4.root");
*/
void Draw_Cascade_ConLev(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	DCanvasInstructions* locInstructions = new DCanvasInstructions("ConLev");
	locInstructions->dPadInstructions[0]->dLogYFlag = true;
	locInstructions->dStyle->SetPadLeftMargin(0.1);

	TH1* locHist = (TH1*)locInputFile->Get("Hist_KinFitResults/ConfidenceLevel");
	locHist->Rebin(10);
	locPlotDrawer.Draw_Object(locHist, locInstructions);
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_Cascade_Vertex("final/cascade_vertp4.root");
*/
void Draw_Cascade_Vertex(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;

	//lifetime
	DCanvasInstructions* locLifetimeInstructions = new DCanvasInstructions("Lifetime");
	locLifetimeInstructions->dStyle->SetOptFit(112);

	//VertexZ
	DCanvasInstructions* locVertexZInstructions = new DCanvasInstructions("VertexZ");
	locVertexZInstructions->dPadInstructions[0]->dLegendNames.push_back("Photoproduction Vertex");
	locVertexZInstructions->dPadInstructions[0]->dLegendNames.push_back("#Xi^{#minus} Decay Vertex");
	locVertexZInstructions->dPadInstructions[0]->dLegendNames.push_back("#Lambda Decay Vertex");
	locVertexZInstructions->dPadInstructions[0]->dHistFillColorFlag = false;
	locVertexZInstructions->dStyle->SetPadLeftMargin(0.08);
	locVertexZInstructions->dStyle->SetHistLineWidth(3);

	//Vertex-Z
	TObjArray* locVertZArray = new TObjArray(3);
	locHist = (TH1*)locInputFile->Get("Hist_ParticleComboKinematics_KinFit/Step0__Photon_Proton_->_K+_K+_Xi-/StepVertexZ");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}");
	locHist->GetXaxis()->SetTitle("Vertex-Z (cm)");
	locHist->Rebin(8);
	locHist->GetXaxis()->SetRangeUser(30.0, 120.0);
	locVertZArray->AddAt(locHist, 0);
	locHist = (TH1*)locInputFile->Get("Hist_ParticleComboKinematics_KinFit/Step1__Xi-_->_Lambda_Pi-/StepVertexZ");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}");
	locHist->GetXaxis()->SetTitle("Vertex-Z (cm)");
	locHist->Rebin(8);
	locHist->GetXaxis()->SetRangeUser(30.0, 120.0);
	locVertZArray->AddAt(locHist, 1);
	locHist = (TH1*)locInputFile->Get("Hist_ParticleComboKinematics_KinFit/Step2__Lambda_->_Proton_Pi-/StepVertexZ");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}");
	locHist->GetXaxis()->SetTitle("Vertex-Z (cm)");
	locHist->Rebin(8);
	locHist->GetXaxis()->SetRangeUser(30.0, 120.0);
	locVertZArray->AddAt(locHist, 2);
	locCanvas = locPlotDrawer.Draw_Superimpose(locVertZArray, locVertexZInstructions);
//	locCanvas->Write();
return;
	//Lambda Lifetime
	locHist = (TH1*)locInputFile->Get("Hist_ParticleComboKinematics_KinFit/Step2__Lambda_->_Proton_Pi-/LambdaLifetime_RestFrame");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}");
	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.05, 0.8); //0.01 -> 1.2
	locCanvas = locPlotDrawer.Draw_Object(locHist, locLifetimeInstructions);

	//Xi- Lifetime
	locHist = (TH1*)locInputFile->Get("Hist_ParticleComboKinematics_KinFit/Step1__Xi-_->_Lambda_Pi-/Xi-Lifetime_RestFrame");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}K^{#plus}#Xi^{#minus}");
	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.05, 0.8); //0.01 -> 1.2
	locCanvas = locPlotDrawer.Draw_Object(locHist, locLifetimeInstructions);
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_Lambda_Vertex("analysis/pass2/klambda.root");
*/
void Draw_Lambda_Vertex(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* locLifetimeInstructions = new DCanvasInstructions("Lifetime");
	locLifetimeInstructions->dStyle->SetOptFit(112);

	//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
	TObjArray* locVertZArray = new TObjArray(2);
	locHist = (TH1*)locInputFile->Get("Hist_ParticleComboKinematics_KinFit/Step0__Photon_Proton_->_K+_Lambda/StepVertexZ");
	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("Hist_ParticleComboKinematics_KinFit/Step1__Lambda_->_Proton_Pi-/StepVertexZ");
	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("Hist_ParticleComboKinematics_KinFit/Step1__Lambda_->_Proton_Pi-/LambdaLifetime_RestFrame");
	locHist->SetTitle("#gammap#rightarrowK^{#plus}#Lambda");
	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.05, 0.8); //0.01 -> 1.2
	locCanvas = locPlotDrawer.Draw_Object(locHist, locLifetimeInstructions);
}

void Fit_Lifetime(TH1* locHist, vector<double>& locInitParams, double locFitRangeMin, double locFitRangeMax)
{
	DFitFunctor_Lifetime* locFitFunctor = new DFitFunctor_Lifetime(locInitParams, kRed);
	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;
}

/*
root -l
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Draw_Cascades.C+
Draw_CascadeInvariantMass("final/cascade_vertp4.root");
*/
void Draw_CascadeInvariantMass(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;
	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dMinPadWidth = 400;
	locCanvasInstructions->dStyle->SetOptFit(102);
	locCanvasInstructions->dStyle->SetTitleOffset(1.1, "y");
	locCanvasInstructions->dStyle->SetPadLeftMargin(0.12);

	TH1* locMassHist = (TH1*)locInputFile->Get("Hist_InvariantMass_XiMinus_PostKinFitCut_KinFit/InvariantMass");
	locMassHist->Rebin(4);
	locMassHist->GetXaxis()->SetRangeUser(1.2, 1.5);
	locMassHist->GetYaxis()->SetTitle("# Combos / 2 MeV/c^{2}");
//	locMassHist->SetTitle("All Beam Energies");

	double locGaussMean = 1.321;
	double locGaussWidth = 0.005;
	double locNumSigmas = 5.0;
	double locGaussHeight = Calc_InitGaussHeight(locMassHist, locGaussMean, locGaussWidth, locNumSigmas);

	double locFitRangeMin = 1.26;
	double locFitRangeMax = 1.44;
	double locCutRangeMin, locCutRangeMax;
	DFitControl* locFitControl = Fit_DoubleGauss(locMassHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);

	DFitControl::DYieldData locYieldData;
	locFitControl->Calc_Yield(locMassHist, false, true, locYieldData);
	cout << "Yield: " << locYieldData.dYield << " +/- " << 100.0*locYieldData.dYieldUncertainty/locYieldData.dYield << "%" << endl;

	locPlotDrawer.Draw_Object(locMassHist, locCanvasInstructions);

	DFitFunctor_DoubleGaussian* locDoubleGaussian = static_cast<DFitFunctor_DoubleGaussian*>(locFitControl->dSignalFitGroup->dFitFunctors[0]);
	cout << "total width = " << locDoubleGaussian->Calc_StandardDeviation().first << " +/- " << locDoubleGaussian->Calc_StandardDeviation().second << endl;
}

DFitControl* Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax)
{
	//Build Functors
	DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = Build_DoubleGaussian(locGaussHeight, locGaussMean, locGaussWidth);
	//DFitFunctor* locFitFunctor_Background = Build_PolyBackground_Linear(locHist, locFitRangeMin, locFitRangeMax);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground_Quadratic(locHist, locFitRangeMin, locFitRangeMax);

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_DoubleGaussian, locFitFunctor_Background);
	locFitControl->dLogLikelihoodFlag = false;

	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);
//	cout << "FIT STATUS = " << locFitStatus << endl;

	//cut range
	locFitFunctor_DoubleGaussian->Calc_CutRange(3.0, locCutRangeMin, locCutRangeMax);

	return locFitControl;
}

DFitFunctor_DoubleGaussian* Build_DoubleGaussian(double locGaussHeight, double locGaussMean, double locGaussWidth)
{
	vector<double> locInitParams(5, 0.0);
	//peak gaussian
	locInitParams[0] = 0.85*locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = 0.9*locGaussWidth;
	//non-gaussian-tail gaussian
	locInitParams[3] = 0.15*locGaussHeight;  locInitParams[4] = 3.0*locGaussWidth;
	DFitFunctor_DoubleGaussian* locDoubleGaussian = new DFitFunctor_DoubleGaussian(locInitParams, kBlue);
	return locDoubleGaussian;
}

double Calc_InitGaussHeight(TH1* locMassHist, double locGaussMean, double locGaussWidth, double locNumSigmas)
{
	int locLowSideBin = locMassHist->GetXaxis()->FindBin(locGaussMean - locNumSigmas*locGaussWidth);
	if(locLowSideBin < 1)
		locLowSideBin = 1;
	int locHighSideBin = locMassHist->GetXaxis()->FindBin(locGaussMean + locNumSigmas*locGaussWidth);
	if(locHighSideBin > locMassHist->GetNbinsX())
		locHighSideBin = locMassHist->GetNbinsX();
	double locAvgBackground = 0.5*(locMassHist->GetBinContent(locLowSideBin) + locMassHist->GetBinContent(locHighSideBin));

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

DFitFunctor* Build_PolyBackground_Linear(TH1* locMassHist, double locFitRangeMin, double locFitRangeMax)
{
	double locSlope = 0.0, locIntercept = 0.0;
	DFitUtilities::Guess_Params_Linear(locMassHist, 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));
}

DFitFunctor* Build_PolyBackground_Quadratic(TH1* locMassHist, double locFitRangeMin, double locFitRangeMax)
{
	double locFlatParam = 0.0, locLinearParam = 0.0, locQuadParam = 0.0;
	double locCriticalPointX = 1.38;
	double locCriticalPointY = locMassHist->GetBinContent(locMassHist->GetXaxis()->FindBin(locCriticalPointX));
	double locCurve_VerticalChange = locMassHist->GetBinContent(locMassHist->GetXaxis()->FindBin(locFitRangeMin)) - locCriticalPointY;
	double locCurve_HorizontalRangeOfVerticalChange = locCriticalPointX - locFitRangeMin;
	DFitUtilities::Calc_ParabolaParams(locCurve_VerticalChange, locCurve_HorizontalRangeOfVerticalChange, locCriticalPointX, locCriticalPointY, locFlatParam, locLinearParam, locQuadParam);

	vector<double> locInitParams(3, 0.0);
	locInitParams[0] = locFlatParam;  locInitParams[1] = locLinearParam;  locInitParams[2] = locQuadParam;
	cout << "params are: " << locInitParams[0] << ", " << locInitParams[1] << ", " << locInitParams[2] << endl;
	return (new DFitFunctor_Polynomial(locInitParams, kMagenta));
}