#include <iostream>
#include <sstream>
#include <string>
#include <fstream>

//masses
void Fit_Pi0Mass(string locInputFileName);
void Fit_OmegaMass(string locInputFileName);
void Draw_MissingMassSquared(string locInputFileName);

//other
void Draw_PID_Timing(string locInputFileName);
void Draw_VertexFit(string locInputFileName);
void Draw_KinFitCutResults(string locInputFileName);

//Utilities
void Draw_HorizontalLine(double locY, int locColor, int locWidth = 3);
void Draw_VerticalLine(double locX, int locColor, int locWidth = 3);

//Fit Utils
double Calc_InitGaussHeight(TH1* locHist, double locGaussMean, double locGaussWidth, double locNumSigmas);
DFitFunctor* Build_PolyBackground(TH1* locHist, double locFitRangeMin, double locFitRangeMax);

//Fit Double Gauss
void Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);
DFitFunctor_DoubleGaussian* Build_DoubleGaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth);

//Fit Voigtian
void Fit_Voigtian(TH1* locHist, double locHeight, double locMean, double locGaussWidth, double locBWGamma, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);
DFitFunctor_Voigtian* Build_Voigtian(TH1* locHist, double locHeight, double locMean, double locGaussWidth, double locBWGamma);

/*
.x Load.C
.L omega_plots.C
//Draw_KinFitCutResults("~/Scratch/gluex/hd_root_3185_p4fitcuts.root");
*/
void Draw_KinFitCutResults(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;
	TObjArray* locHistArray = NULL;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;

	//First do con-lev
	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_KinFitResults/ConfidenceLevel");
	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);

	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("Pre-Cut");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("P4 KinFit Converged");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("P4 KinFit CL Cut > 5#sigma");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("P4 KinFit CL Cut > 3#sigma");
//	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("P4 KinFit CL Cut > 1%");
//	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("P4 KinFit CL Cut > 5%");

	//Missing Mass
	locHistArray = new TObjArray();
	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_MissingMassSquared_PreCut/MissingMassSquared");
	locHist->Rebin(2);
	locHist->GetXaxis()->SetRangeUser(-0.03, 0.03);
	locHist->GetXaxis()->SetTitle("#gammap#rightarrow#pi^{#plus}#pi^{#minus}#gamma#gamma Missing Mass Squared (GeV/c^{2})^{2}");
	locHist->GetYaxis()->SetTitle("# Combos / 100 (MeV/c^{2})^{2}");
	locHistArray->AddLast(locHist);

	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_MissingMassSquared_Converged/MissingMassSquared");
	locHist->Rebin(2);
	locHist->GetXaxis()->SetRangeUser(-0.03, 0.03);
	locHist->GetXaxis()->SetTitle("#gammap#rightarrow#pi^{#plus}#pi^{#minus}#gamma#gamma Missing Mass Squared (GeV/c^{2})^{2}");
	locHist->GetYaxis()->SetTitle("# Combos / 100 (MeV/c^{2})^{2}");
	locHistArray->AddLast(locHist);

	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_MissingMassSquared_Cut_5sigma/MissingMassSquared");
	locHist->Rebin(2);
	locHist->GetXaxis()->SetRangeUser(-0.03, 0.03);
	locHist->GetXaxis()->SetTitle("#gammap#rightarrow#pi^{#plus}#pi^{#minus}#gamma#gamma Missing Mass Squared (GeV/c^{2})^{2}");
	locHist->GetYaxis()->SetTitle("# Combos / 100 (MeV/c^{2})^{2}");
	locHistArray->AddLast(locHist);

	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_MissingMassSquared_Cut_3sigma/MissingMassSquared");
	locHist->Rebin(2);
	locHist->GetXaxis()->SetRangeUser(-0.03, 0.03);
	locHist->GetXaxis()->SetTitle("#gammap#rightarrow#pi^{#plus}#pi^{#minus}#gamma#gamma Missing Mass Squared (GeV/c^{2})^{2}");
	locHist->GetYaxis()->SetTitle("# Combos / 100 (MeV/c^{2})^{2}");
	locHistArray->AddLast(locHist);

	locCanvas = locPlotDrawer.Draw_Superimpose(locHistArray, locCanvasInstructions);

	//Omega
	locHistArray = new TObjArray();
	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_InvariantMass_Omega_PreCut/InvariantMass");
	locHist->Rebin(5);
	locHist->GetXaxis()->SetTitle("#pi^{#plus}#pi^{#minus}#gamma#gamma Invariant Mass (GeV/c^{2})");
	locHist->GetYaxis()->SetTitle("# Combos / 5 (MeV/c^{2})");
	locHistArray->AddLast(locHist);

	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_InvariantMass_Omega_Converged/InvariantMass");
	locHist->Rebin(5);
	locHist->GetXaxis()->SetTitle("#pi^{#plus}#pi^{#minus}#gamma#gamma Invariant Mass (GeV/c^{2})");
	locHist->GetYaxis()->SetTitle("# Combos / 5 (MeV/c^{2})");
	locHistArray->AddLast(locHist);

	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_InvariantMass_Omega_Cut_5sigma/InvariantMass");
	locHist->Rebin(5);
	locHist->GetXaxis()->SetTitle("#pi^{#plus}#pi^{#minus}#gamma#gamma Invariant Mass (GeV/c^{2})");
	locHist->GetYaxis()->SetTitle("# Combos / 5 (MeV/c^{2})");
	locHistArray->AddLast(locHist);

	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_InvariantMass_Omega_Cut_3sigma/InvariantMass");
	locHist->Rebin(5);
	locHist->GetXaxis()->SetTitle("#pi^{#plus}#pi^{#minus}#gamma#gamma Invariant Mass (GeV/c^{2})");
	locHist->GetYaxis()->SetTitle("# Combos / 5 (MeV/c^{2})");
	locHistArray->AddLast(locHist);

	locCanvas = locPlotDrawer.Draw_Superimpose(locHistArray, locCanvasInstructions);
}

/*
.x Load.C
.L omega_plots.C
//Draw_VertexFit("~/Scratch/gluex/hd_root_3185_widepi0.root");
*/
void Draw_VertexFit(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;

	//Confidence Level
	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Cut_OneVertexKinFit/ConfidenceLevel");
	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);

	//Vertex-z
	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Cut_OneVertexKinFit/VertexZ");
	locHist->GetXaxis()->SetRangeUser(30.0, 110.0);
	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);

	//Vertex 2D
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Cut_OneVertexKinFit/VertexYVsX");
	loc2DHist->GetXaxis()->SetRangeUser(-4.0, 4.0);
	loc2DHist->GetYaxis()->SetRangeUser(-4.0, 4.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
}

/*
.x Load.C
.L omega_plots.C
//Draw_PID_Timing("~/Scratch/gluex/hd_root_3185_widepi0.root");
*/
void Draw_PID_Timing(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dStyle->SetOptFit(0);

	//Beta-vs-p functions
	TF1* locBetaVsPFunc_Proton = new TF1("BetaVsPFunc_Proton", "x/sqrt(x*x + [0]*[0])", 0.0, 12.0);
	locBetaVsPFunc_Proton->SetParameter(0, 0.938272046);
	locBetaVsPFunc_Proton->SetLineWidth(4);
	locBetaVsPFunc_Proton->SetLineColor(kMagenta);
	locBetaVsPFunc_Proton->SetNpx(1000);

	TF1* locBetaVsPFunc_Kaon = new TF1("BetaVsPFunc_Kaon", "x/sqrt(x*x + [0]*[0])", 0.0, 12.0);
	locBetaVsPFunc_Kaon->SetParameter(0, 0.493677);
	locBetaVsPFunc_Kaon->SetLineWidth(4);
	locBetaVsPFunc_Kaon->SetLineColor(kMagenta);
	locBetaVsPFunc_Kaon->SetNpx(1000);

	TF1* locBetaVsPFunc_Pion = new TF1("BetaVsPFunc_Pion", "x/sqrt(x*x + [0]*[0])", 0.0, 12.0);
	locBetaVsPFunc_Pion->SetParameter(0, 0.13957018);
	locBetaVsPFunc_Pion->SetLineWidth(4);
	locBetaVsPFunc_Pion->SetLineColor(kMagenta);
	locBetaVsPFunc_Pion->SetNpx(1000);

	//Beta vs p q+: TOF
	loc2DHist = (TH2*)locInputFile->Get("Independent/Hist_DetectorPID/TOF/BetaVsP_q+");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	loc2DHist->GetYaxis()->SetRangeUser(0.0, 1.1);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Proton);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Kaon);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Pion);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);

	//Beta vs p q+: BCAL
	loc2DHist = (TH2*)locInputFile->Get("Independent/Hist_DetectorPID/BCAL/BetaVsP_q+");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	loc2DHist->GetYaxis()->SetRangeUser(0.0, 1.1);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Proton);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Pion);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);

	//Beta vs p q+: FCAL
	loc2DHist = (TH2*)locInputFile->Get("Independent/Hist_DetectorPID/FCAL/BetaVsP_q+");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	loc2DHist->GetYaxis()->SetRangeUser(0.0, 1.1);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Proton);
	loc2DHist->GetListOfFunctions()->Add(locBetaVsPFunc_Pion);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);

	//delta-t vs shower-E q0: BCAL
	loc2DHist = (TH2*)locInputFile->Get("p3pi_pmiss_any/Hist_PID/Photon/BCAL/DeltaTVsShowerE_Photon");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 2.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);

	//delta-t vs shower-E q0: FCAL
	loc2DHist = (TH2*)locInputFile->Get("p3pi_pmiss_any/Hist_PID/Photon/FCAL/DeltaTVsShowerE_Photon");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);

/*
	//Proton TOF Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Proton/TOF/DeltaTVsP");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-1.0, kMagenta, 5);
	Draw_HorizontalLine(1.0, kMagenta, 5);

	//Proton BCAL Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Proton/BCAL/DeltaTVsP");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-1.0, kMagenta, 5);
	Draw_HorizontalLine(1.0, kMagenta, 5);

	//Proton FCAL Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Proton/FCAL/DeltaTVsP");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-2.0, kMagenta, 5);
	Draw_HorizontalLine(2.0, kMagenta, 5);

	//Pi+ TOF Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Pi+/TOF/DeltaTVsP");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-1.0, kMagenta, 5);
	Draw_HorizontalLine(1.0, kMagenta, 5);

	//Pi+ BCAL Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Pi+/BCAL/DeltaTVsP");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-1.0, kMagenta, 5);
	Draw_HorizontalLine(1.0, kMagenta, 5);

	//Pi+ FCAL Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Pi+/FCAL/DeltaTVsP");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-2.0, kMagenta, 5);
	Draw_HorizontalLine(2.0, kMagenta, 5);

	//Photon BCAL Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Photon/BCAL/DeltaTVsShowerE");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 2.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-1.0, kMagenta, 5);
	Draw_HorizontalLine(1.0, kMagenta, 5);

	//Photon FCAL Delta-t
	loc2DHist = (TH2*)locInputFile->Get("p3pi_preco_any/Hist_PID/Photon/FCAL/DeltaTVsShowerE");
	loc2DHist->Rebin2D(1, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 2.0);
	loc2DHist->GetYaxis()->SetRangeUser(-5.0, 5.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	Draw_HorizontalLine(-2.0, kMagenta, 5);
	Draw_HorizontalLine(2.0, kMagenta, 5);
*/
}

/*
.x Load.C
.L omega_plots.C
//Draw_MissingMassSquared("~/Code/Work/gluex/ROOT_Scripts/p3pi_hists_measuredcuts.root");
//Draw_MissingMassSquared("~/Scratch/gluex/hd_root_3185_widepi0.root");
*/
void Draw_MissingMassSquared(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;

//	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_MissingMassSquared/MissingMassSquared");
	locHist = (TH1*)locInputFile->Get("MissingMassSquared");
//	locHist->Rebin(2);
	locHist->GetXaxis()->SetRangeUser(-0.03, 0.03);
	locHist->GetXaxis()->SetTitle("#gammap#rightarrow#pi^{#plus}#pi^{#minus}#gamma#gamma Missing Mass Squared (GeV/c^{2})^{2}");
	locHist->GetYaxis()->SetTitle("# Combos / 200 (MeV/c^{2})^{2}");

	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);
	Draw_VerticalLine(-0.007, kRed);
	Draw_VerticalLine(0.005, kRed);
}

/*
.x Load.C
.L omega_plots.C
//Fit_Pi0Mass("~/Code/Work/gluex/ROOT_Scripts/p3pi_hists_kinfitcuts.root");
//Fit_Pi0Mass("~/Code/Work/gluex/ROOT_Scripts/p3pi_hists_measuredcuts.root");
////Fit_Pi0Mass("~/Scratch/gluex/hd_root_3185_widepi0.root");
*/
void Fit_Pi0Mass(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;

//	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_InvariantMass_Pi0/InvariantMass");
	locHist = (TH1*)locInputFile->Get("Pi0Mass_Measured");
//	locHist = (TH1*)locInputFile->Get("Pi0Mass_KinFit");
	locHist->Rebin(4);
	locHist->GetYaxis()->SetTitle("# Combos / 2 MeV/c^{2}");

	double locGaussMean = 0.135;
	double locGaussWidth = 0.01;
	double locGaussHeight = Calc_InitGaussHeight(locHist, locGaussMean, locGaussWidth, 5.0);
	double locFitRangeMin = 0.05; //measured
	double locFitRangeMax = 0.22; //measured
//	double locFitRangeMin = 0.08; //kinfit
//	double locFitRangeMax = 0.19; //kinfit
	double locCutRangeMin = 0.0, locCutRangeMax = 0.0;
	Fit_DoubleGauss(locHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);

	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);
	Draw_VerticalLine(locCutRangeMin, kRed);
	Draw_VerticalLine(locCutRangeMax, kRed);
}

/*
.x Load.C
.L omega_plots.C
//Fit_OmegaMass("~/Code/Work/gluex/ROOT_Scripts/p3pi_hists_kinfitcuts.root");
//Fit_OmegaMass("~/Code/Work/gluex/ROOT_Scripts/p3pi_hists_measuredcuts.root");
////Fit_OmegaMass("~/Scratch/gluex/hd_root_3185_widepi0.root");
*/
void Fit_OmegaMass(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;

//	locHist = (TH1*)locInputFile->Get("p3pi_preco_any/Hist_InvariantMass_Omega_Kinfit/InvariantMass");
	locHist = (TH1*)locInputFile->Get("OmegaMass_Measured");
//	locHist = (TH1*)locInputFile->Get("OmegaMass_KinFit");
	locHist->Rebin(6); //measured
//	locHist->Rebin(4); //kinfit
	locHist->GetYaxis()->SetTitle("# Combos / 6 MeV/c^{2}"); //measured
//	locHist->GetYaxis()->SetTitle("# Combos / 4 MeV/c^{2}"); //kinfit
	locHist->GetXaxis()->SetTitle("#pi^{#plus}#pi^{#minus}#gamma#gamma Invariant Mass (GeV/c^{2})");

	double locMean = 0.79;
	double locGaussWidth = 0.033; //measured
//	double locGaussWidth = 0.01; //kinfit
	double locHeight = Calc_InitGaussHeight(locHist, locMean, locGaussWidth, 5.0);
	double locBWGamma = 0.00849;
	double locFitRangeMin = 0.62; //measured
	double locFitRangeMax = 0.94; //measured
//	double locFitRangeMin = 0.68; //kinfit
//	double locFitRangeMax = 0.88; //kinfit
	double locCutRangeMin = 0.0, locCutRangeMax = 0.0;
	Fit_Voigtian(locHist, locHeight, locMean, locGaussWidth, locBWGamma, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);

	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);
	Draw_VerticalLine(0.72, kRed);
	Draw_VerticalLine(0.84, kRed);
}

void Draw_HorizontalLine(double locY, int locColor, int locWidth)
{
	TLine* locLine = new TLine(gPad->GetUxmin(), locY, gPad->GetUxmax(), locY);
	locLine->SetLineColor(locColor);
	locLine->SetLineWidth(locWidth);
	locLine->Draw("SAME");
}

void Draw_VerticalLine(double locX, int locColor, int locWidth)
{
	TLine* locLine = new TLine(locX, gPad->GetUymin(), locX, gPad->GetUymax());
	locLine->SetLineColor(locColor);
	locLine->SetLineWidth(locWidth);
	locLine->Draw("SAME");
}

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

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

void 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(locHist, locGaussHeight, locGaussMean, locGaussWidth);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground(locHist, locFitRangeMin, locFitRangeMax);

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

	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = false;

	int locFitStatus = locPlotFitter.Fit_Hist(locHist, locFitControl);
	cout << "FIT STATUS = " << locFitStatus << endl;
	cout << "TOTAL WIDTH = " << locFitFunctor_DoubleGaussian->Calc_StandardDeviation().first << endl;

	//cut range
	locFitFunctor_DoubleGaussian->Calc_CutRange(3.0, locCutRangeMin, locCutRangeMax);
	cout << "CUT RANGE: " << locCutRangeMin << ", " << locCutRangeMax << endl;
}

void Fit_Voigtian(TH1* locHist, double locHeight, double locMean, double locGaussWidth, double locBWGamma, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax)
{
	//Build Functors
	DFitFunctor_Voigtian* locFitFunctor_Voigtian = Build_Voigtian(locHist, locHeight, locMean, locGaussWidth, locBWGamma);
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground(locHist, locFitRangeMin, locFitRangeMax);

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

	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = false;

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

	//cut range
	locFitFunctor_Voigtian->Calc_CutRange(3.0, locCutRangeMin, locCutRangeMax);
	cout << "CUT RANGE: " << locCutRangeMin << ", " << locCutRangeMax << endl;

	//yield
	DFitControl::DYieldData locYieldData;
	locFitControl->Calc_Yield(locHist, false, 0.68, 0.88, locYieldData); //measured
//	locFitControl->Calc_Yield(locHist, false, 0.72, 0.84, locYieldData); //kinfit

	cout << "yield, sig/back = " << locYieldData.dYield << ", " << locYieldData.dYield/locYieldData.dBackground << endl;
}

DFitFunctor* Build_PolyBackground(TH1* locHist, double locFitRangeMin, double locFitRangeMax)
{
	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));
}

DFitFunctor_DoubleGaussian* Build_DoubleGaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth)
{
	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] = 1.5*locGaussWidth;
	DFitFunctor_DoubleGaussian* locDoubleGaussian = new DFitFunctor_DoubleGaussian(locInitParams, kBlue);
	locDoubleGaussian->dParamLimits[2] = pair<double, double>(0.005, 0.5);
	locDoubleGaussian->dParamLimits[4] = pair<double, double>(0.005, 0.5);
	return locDoubleGaussian;
}

DFitFunctor_Voigtian* Build_Voigtian(TH1* locHist, double locHeight, double locMean, double locGaussWidth, double locBWGamma)
{
	vector<double> locInitParams(4, 0.0);
	locInitParams[0] = locHeight;  locInitParams[1] = locMean;
	locInitParams[2] = locGaussWidth;  locInitParams[3] = locBWGamma;

	DFitFunctor_Voigtian* locVoigtian = new DFitFunctor_Voigtian(locInitParams, kBlue);
	locVoigtian->dParamLimits[2] = pair<double, double>(0.005, 0.5);
	locVoigtian->Fix_Parameter(3);
	return locVoigtian;
}