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

#include "PlotDrawer/DCanvasInstructions.h"
#include "PlotDrawer/DPadInstructions.h"
#include "PlotDrawer/DPlotDrawer.h"

#include "PlotFitter/DPlotFitter.h"
#include "PlotFitter/DFitControl.h"
#include "PlotFitter/DFitUtilities.h"
#include "PlotFitter/DFitFunctors_Miscellaneous.h"
#include "PlotFitter/DFitFunctors_Gaussians.h"

#include "PlotData/DStatUtilities.h"

#include "particleType.h"

//Resolution
void Fit_DeltaP(string locInputFileName);
void Fit_DeltaTheta(string locInputFileName);
void Fit_DeltaPhiVsP(string locInputFileName);
void Fit_DeltaPhiVsTheta(string locInputFileName);

//Missing masses
void Draw_MissingMasses(string locInputFileName);

//Efficiency
void Draw_TrackEfficiencies(string locInputFileName);
void Draw_TrackEfficiencyRatio(string locInputFileName_Data, string locInputFileName_MC);
void Draw_MatchingEfficiencies(string locInputFileName);

void Fit_Slices_DoubleGaussians(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locMeanGraph, TGraphErrors*& locSigmaGraph, double locMinimumCounts, double locInitFitRangeMin, double locInitFitRangeMax, double locInitGaussWidth);

DFitControl* Fit_DoubleGauss(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);
DFitControl* Fit_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth, double locFitRangeMin, double locFitRangeMax, double& locCutRangeMin, double& locCutRangeMax);

double Calc_InitGaussHeight(TH1* locHist, double locGaussMean, double locGaussWidth, double locNumSigmas);

DFitFunctor* Build_PolyBackground(TH1* locHist, double locFitRangeMin, double locFitRangeMax);
DFitFunctor_DoubleGaussian* Build_DoubleGaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth);
DFitFunctor_Gaussian* Build_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth);

void Draw_HorizontalLine(double locY, int locColor, int locWidth = 3);
void Draw_VerticalLine(double locX, int locColor, int locWidth = 3);
double* Convert_VectorToArray(const vector<double>& locVector);


/*
.x Load.C
.L effres_plots.C
//Draw_MissingMasses("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Draw_MissingMasses(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;

	//Missing Proton
	locHist = (TH1*)locInputFile->Get("TrackEff_MissingProton/Hist_MissingMass_PostPi0/MissingMass");
	locHist->GetXaxis()->SetTitle("#gammap#rightarrow#pi^{#plus}#pi^{#minus} Missing Mass (GeV/c^{2})");
	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);
	Draw_VerticalLine(0.88, kRed);
	Draw_VerticalLine(1.0, kRed);

	//Missing PiPlus
	locHist = (TH1*)locInputFile->Get("TrackEff_MissingPiPlus/Hist_MissingMassSquared_PostPi0/MissingMassSquared");
	locHist->GetXaxis()->SetTitle("#gammap#rightarrowp#pi^{#minus} Missing Mass Squared (GeV/c^{2})^{2}");
	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);
	Draw_VerticalLine(-0.03, kRed);
	Draw_VerticalLine(0.07, kRed);

	//Missing PiMinus
	locHist = (TH1*)locInputFile->Get("TrackEff_MissingPiMinus/Hist_MissingMassSquared_PostPi0/MissingMassSquared");
	locHist->GetXaxis()->SetTitle("#gammap#rightarrowp#pi^{#plus} Missing Mass Squared (GeV/c^{2})^{2}");
	locCanvas = locPlotDrawer.Draw_Object(locHist, locCanvasInstructions);
	Draw_VerticalLine(-0.03, kRed);
	Draw_VerticalLine(0.07, kRed);
}

/*
.x Load.C
.L effres_plots.C
//Draw_MatchingEfficiencies("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Draw_MatchingEfficiencies(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;
	TObjArray* locCanvasArray = NULL;
	TObjArray* locHistogramSlices = NULL;
	TH2* loc2DHist_Found = NULL;
	TH2* loc2DHist_Missing = NULL;
	TH2* locEfficiencyHist = NULL;
	string locHistName, locNewHistName, locTitle;

	double locEffMax = 1.0;
	double locEffMin = 0.0;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dStyle->SetPalette(55);
	locCanvasInstructions->dStyle->SetNumberContours(20);

	map<Particle_t, string> locDirNameMap;
	locDirNameMap[Proton] = "TrackEff_MissingProton";
	locDirNameMap[PiPlus] = "TrackEff_MissingPiPlus";
	locDirNameMap[PiMinus] = "TrackEff_MissingPiMinus";

	map<Particle_t, string>::iterator locIterator = locDirNameMap.begin();
	for(; locIterator != locDirNameMap.end(); ++locIterator)
	{
		Particle_t locPID = locIterator->first;
		string locDirBaseName = locIterator->second + string("/Custom_TrackingEfficiency/Matching_TimeBased/");

		//SC
		locHistName = locDirBaseName + string("SC/PhiVsTheta_HasHit");
		loc2DHist_Found = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("PhiVsTheta_HasHit_") + string(ParticleType(locPID));
		loc2DHist_Found->SetName(locNewHistName.c_str());
		loc2DHist_Found->Rebin2D(3, 4);

		locHistName = locDirBaseName + string("SC/PhiVsTheta_NoHit");
		loc2DHist_Missing = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("PhiVsTheta_NoHit_") + string(ParticleType(locPID));
		loc2DHist_Missing->SetName(locNewHistName.c_str());
		loc2DHist_Missing->Rebin2D(3, 4);

		locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
//		locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
//		locEfficiencyHist->GetYaxis()->SetRangeUser(0.0, 4.0);
		locTitle = string(ParticleName_ROOT(locPID)) + string(" SC Matching Efficiency");
		locEfficiencyHist->SetTitle(locTitle.c_str());
		locEfficiencyHist->SetMinimum(locEffMin);
		locEfficiencyHist->SetMaximum(locEffMax);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
		locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);

		//BCAL
		locHistName = locDirBaseName + string("BCAL/TrackBCALPhiVsZ_HasHit");
		loc2DHist_Found = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("TrackBCALPhiVsZ_HasHit_") + string(ParticleType(locPID));
		loc2DHist_Found->SetName(locNewHistName.c_str());
		loc2DHist_Found->Rebin2D(6, 8);

		locHistName = locDirBaseName + string("BCAL/TrackBCALPhiVsZ_NoHit");
		loc2DHist_Missing = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("TrackBCALPhiVsZ_NoHit_") + string(ParticleType(locPID));
		loc2DHist_Missing->SetName(locNewHistName.c_str());
		loc2DHist_Missing->Rebin2D(6, 8);

		locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
//		locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
//		locEfficiencyHist->GetYaxis()->SetRangeUser(0.0, 4.0);
		locTitle = string(ParticleName_ROOT(locPID)) + string(" BCAL Matching Efficiency");
		locEfficiencyHist->SetTitle(locTitle.c_str());
		locEfficiencyHist->SetMinimum(locEffMin);
		locEfficiencyHist->SetMaximum(locEffMax);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
		locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);

		//TOF
		locHistName = locDirBaseName + string("TOFPoint/TrackTOFYVsX_HasHit");
		loc2DHist_Found = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("TrackTOFYVsX_HasHit_") + string(ParticleType(locPID));
		loc2DHist_Found->SetName(locNewHistName.c_str());
		loc2DHist_Found->Rebin2D(4, 4);

		locHistName = locDirBaseName + string("TOFPoint/TrackTOFYVsX_NoHit");
		loc2DHist_Missing = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("TrackTOFYVsX_NoHit_") + string(ParticleType(locPID));
		loc2DHist_Missing->SetName(locNewHistName.c_str());
		loc2DHist_Missing->Rebin2D(4, 4);

		locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
//		locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
//		locEfficiencyHist->GetYaxis()->SetRangeUser(0.0, 4.0);
		locTitle = string(ParticleName_ROOT(locPID)) + string(" TOF Matching Efficiency");
		locEfficiencyHist->SetTitle(locTitle.c_str());
		locEfficiencyHist->SetMinimum(locEffMin);
		locEfficiencyHist->SetMaximum(locEffMax);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
		locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);

		//FCAL
		locHistName = locDirBaseName + string("FCAL/TrackFCALYVsX_HasHit");
		loc2DHist_Found = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("TrackFCALYVsX_HasHit_") + string(ParticleType(locPID));
		loc2DHist_Found->SetName(locNewHistName.c_str());
		loc2DHist_Found->Rebin2D(4, 4);

		locHistName = locDirBaseName + string("FCAL/TrackFCALYVsX_NoHit");
		loc2DHist_Missing = (TH2*)locInputFile->Get(locHistName.c_str());
		locNewHistName = string("TrackFCALYVsX_NoHit_") + string(ParticleType(locPID));
		loc2DHist_Missing->SetName(locNewHistName.c_str());
		loc2DHist_Missing->Rebin2D(4, 4);

		locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
//		locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
//		locEfficiencyHist->GetYaxis()->SetRangeUser(0.0, 4.0);
		locTitle = string(ParticleName_ROOT(locPID)) + string(" FCAL Matching Efficiency");
		locEfficiencyHist->SetTitle(locTitle.c_str());
		locEfficiencyHist->SetMinimum(locEffMin);
		locEfficiencyHist->SetMaximum(locEffMax);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//		locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
		locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);
	}
}

/*
.x Load.C
.L effres_plots.C
//Draw_TrackEfficiencies("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Draw_TrackEfficiencies(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;
	TObjArray* locCanvasArray = NULL;
	TObjArray* locHistogramSlices = NULL;
	TH2* loc2DHist_Found = NULL;
	TH2* loc2DHist_Missing = NULL;
	TH2* locEfficiencyHist = NULL;

	double locEffMax = 1.0;
	double locEffMin = 0.0;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dStyle->SetPalette(55);
	locCanvasInstructions->dStyle->SetNumberContours(20);

	//Proton: Found
	loc2DHist_Found = (TH2*)locInputFile->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_Proton");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);
	loc2DHist_Found->GetYaxis()->SetRangeUser(0.0, 4.0);

	//Proton: Missing
	loc2DHist_Missing = (TH2*)locInputFile->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_Proton");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);
	loc2DHist_Missing->GetYaxis()->SetRangeUser(0.0, 4.0);

	//Proton: Efficiency
	locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist->GetYaxis()->SetRangeUser(0.0, 4.0);
	locEfficiencyHist->SetTitle("p Reconstruction Efficiency");
	locEfficiencyHist->SetMinimum(locEffMin);
	locEfficiencyHist->SetMaximum(locEffMax);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);

	//PiPlus: Found
	loc2DHist_Found = (TH2*)locInputFile->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_PiPlus");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiPlus: Missing
	loc2DHist_Missing = (TH2*)locInputFile->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_PiPlus");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiPlus: Efficiency
	locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist->SetTitle("#pi^{#plus} Reconstruction Efficiency");
	locEfficiencyHist->SetMinimum(locEffMin);
	locEfficiencyHist->SetMaximum(locEffMax);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);

	//PiMinus: Found
	loc2DHist_Found = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_PiMinus");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiMinus: Missing
	loc2DHist_Missing = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_PiMinus");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiMinus: Efficiency
	locEfficiencyHist = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist->SetTitle("#pi^{#minus} Reconstruction Efficiency");
	locEfficiencyHist->SetMinimum(locEffMin);
	locEfficiencyHist->SetMaximum(locEffMax);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist, locCanvasInstructions);
}

void Draw_TrackEfficiencyRatio(string locInputFileName_Data, string locInputFileName_MC)
{
	DPlotDrawer locPlotDrawer;

	TFile* locInputFile_Data = new TFile(locInputFileName_Data.c_str(), "READ");
	TFile* locInputFile_MC = new TFile(locInputFileName_MC.c_str(), "READ");

	TH1* locHist = NULL;
	TH2* loc2DHist = NULL;
	TCanvas* locCanvas = NULL;
	TObjArray* locCanvasArray = NULL;
	TObjArray* locHistogramSlices = NULL;
	TH2* loc2DHist_Found = NULL;
	TH2* loc2DHist_Missing = NULL;

	double locEffMax = 1.0;
	double locEffMin = 0.0;

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dStyle->SetPalette(55);
	locCanvasInstructions->dStyle->SetNumberContours(20);

	//Proton: Found
	loc2DHist_Found = (TH2*)locInputFile_Data->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_Proton");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);
	loc2DHist_Found->GetYaxis()->SetRangeUser(0.0, 4.0);

	//Proton: Missing
	loc2DHist_Missing = (TH2*)locInputFile_Data->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_Proton");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);
	loc2DHist_Missing->GetYaxis()->SetRangeUser(0.0, 4.0);

	//Proton: Efficiency
	TH2* locEfficiencyHist_Proton_Data = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist_Proton_Data->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist_Proton_Data->GetYaxis()->SetRangeUser(0.0, 4.0);
	locEfficiencyHist_Proton_Data->SetName("p_eff_data");
	locEfficiencyHist_Proton_Data->SetTitle("p Reconstruction Efficiency");
	locEfficiencyHist_Proton_Data->SetMinimum(locEffMin);
	locEfficiencyHist_Proton_Data->SetMaximum(locEffMax);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist_Proton_Data, locCanvasInstructions);

	//PiPlus: Found
	loc2DHist_Found = (TH2*)locInputFile_Data->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_PiPlus");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiPlus: Missing
	loc2DHist_Missing = (TH2*)locInputFile_Data->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_PiPlus");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiPlus: Efficiency
	TH2* locEfficiencyHist_PiPlus_Data = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist_PiPlus_Data->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist_PiPlus_Data->SetName("pip_eff_data");
	locEfficiencyHist_PiPlus_Data->SetTitle("#pi^{#plus} Reconstruction Efficiency");
	locEfficiencyHist_PiPlus_Data->SetMinimum(locEffMin);
	locEfficiencyHist_PiPlus_Data->SetMaximum(locEffMax);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist_PiPlus_Data, locCanvasInstructions);

	//PiMinus: Found
	loc2DHist_Found = (TH2*)locInputFile_Data->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_PiMinus");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiMinus: Missing
	loc2DHist_Missing = (TH2*)locInputFile_Data->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_PiMinus");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiMinus: Efficiency
	TH2* locEfficiencyHist_PiMinus_Data = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist_PiMinus_Data->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist_PiMinus_Data->SetName("pim_eff_data");
	locEfficiencyHist_PiMinus_Data->SetTitle("#pi^{#minus} Reconstruction Efficiency");
	locEfficiencyHist_PiMinus_Data->SetMinimum(locEffMin);
	locEfficiencyHist_PiMinus_Data->SetMaximum(locEffMax);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist_PiMinus_Data, locCanvasInstructions);

	//MC

	//Proton: Found
	loc2DHist_Found = (TH2*)locInputFile_MC->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_Proton");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);
	loc2DHist_Found->GetYaxis()->SetRangeUser(0.0, 4.0);

	//Proton: Missing
	loc2DHist_Missing = (TH2*)locInputFile_MC->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_Proton");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);
	loc2DHist_Missing->GetYaxis()->SetRangeUser(0.0, 4.0);

	//Proton: Efficiency
	TH2* locEfficiencyHist_Proton_MC = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist_Proton_MC->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist_Proton_MC->GetYaxis()->SetRangeUser(0.0, 4.0);
	locEfficiencyHist_Proton_MC->SetName("p_eff_mc");
	locEfficiencyHist_Proton_MC->SetTitle("p Reconstruction Efficiency");
	locEfficiencyHist_Proton_MC->SetMinimum(locEffMin);
	locEfficiencyHist_Proton_MC->SetMaximum(locEffMax);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist_Proton_MC, locCanvasInstructions);

	//PiPlus: Found
	loc2DHist_Found = (TH2*)locInputFile_MC->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_PiPlus");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiPlus: Missing
	loc2DHist_Missing = (TH2*)locInputFile_MC->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_PiPlus");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiPlus: Efficiency
	TH2* locEfficiencyHist_PiPlus_MC = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist_PiPlus_MC->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist_PiPlus_MC->SetName("pip_eff_mc");
	locEfficiencyHist_PiPlus_MC->SetTitle("#pi^{#plus} Reconstruction Efficiency");
	locEfficiencyHist_PiPlus_MC->SetMinimum(locEffMin);
	locEfficiencyHist_PiPlus_MC->SetMaximum(locEffMax);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist_PiPlus_MC, locCanvasInstructions);

	//PiMinus: Found
	loc2DHist_Found = (TH2*)locInputFile_MC->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackFound_VertexZBin0");
	loc2DHist_Found->SetName("PVsTheta_TrackFound_PiMinus");
	loc2DHist_Found->Rebin2D(4, 5);
	loc2DHist_Found->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiMinus: Missing
	loc2DHist_Missing = (TH2*)locInputFile_MC->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Efficiency_TimeBased/PVsTheta_TrackMissing_VertexZBin0");
	loc2DHist_Missing->SetName("PVsTheta_TrackMissing_PiMinus");
	loc2DHist_Missing->Rebin2D(4, 5);
	loc2DHist_Missing->GetXaxis()->SetRangeUser(0.0, 160.0);

	//PiMinus: Efficiency
	TH2* locEfficiencyHist_PiMinus_MC = (TH2*)DStatUtilities::Calc_Acceptance_SeparateTotals(loc2DHist_Found, loc2DHist_Missing);
	locEfficiencyHist_PiMinus_MC->GetXaxis()->SetRangeUser(0.0, 160.0);
	locEfficiencyHist_PiMinus_MC->SetName("pim_eff_mc");
	locEfficiencyHist_PiMinus_MC->SetTitle("#pi^{#minus} Reconstruction Efficiency");
	locEfficiencyHist_PiMinus_MC->SetMinimum(locEffMin);
	locEfficiencyHist_PiMinus_MC->SetMaximum(locEffMax);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Found, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(loc2DHist_Missing, locCanvasInstructions);
//	locCanvas = locPlotDrawer.Draw_Object(locEfficiencyHist_PiMinus_MC, locCanvasInstructions);


	//computes (MC - Data)/Data
	TH2* locRatioHist_Proton = (TH2*)DStatUtilities::Compare((TH1*)locEfficiencyHist_Proton_MC, (TH1*)locEfficiencyHist_Proton_Data);
	locRatioHist_Proton->SetMaximum(1.0);
	locRatioHist_Proton->SetTitle("p Reconstruction Ratio: (MC - Data)/Data");
	locCanvas = locPlotDrawer.Draw_Object(locRatioHist_Proton, locCanvasInstructions);

	TH2* locRatioHist_PiPlus = (TH2*)DStatUtilities::Compare((TH1*)locEfficiencyHist_PiPlus_MC, (TH1*)locEfficiencyHist_PiPlus_Data);
	locRatioHist_PiPlus->SetMaximum(1.0);
	locRatioHist_PiPlus->SetTitle("#pi^{#plus} Reconstruction Ratio: (MC - Data)/Data");
	locCanvas = locPlotDrawer.Draw_Object(locRatioHist_PiPlus, locCanvasInstructions);

	TH2* locRatioHist_PiMinus = (TH2*)DStatUtilities::Compare((TH1*)locEfficiencyHist_PiMinus_MC, (TH1*)locEfficiencyHist_PiMinus_Data);
	locRatioHist_PiMinus->SetMaximum(1.0);
	locRatioHist_PiMinus->SetTitle("#pi^{#minus} Reconstruction Ratio: (MC - Data)/Data");
	locCanvas = locPlotDrawer.Draw_Object(locRatioHist_PiMinus, locCanvasInstructions);
}

/*
.x Load.C
.L effres_plots.C
//Draw_DeltaPVsTheta("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Draw_DeltaPVsTheta(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

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

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
//	locCanvasInstructions->dPadInstructions[0]->dDefaultDrawOptionMap["TH2"] = "LEGO2";

	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPOverPVsTheta");
	loc2DHist->Rebin2D(4, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 140.0);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
}

/*
.x Load.C
.L effres_plots.C
//Fit_DeltaPhiVsTheta("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Fit_DeltaPhiVsTheta(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

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

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

	TGraphErrors *locMeanGraph, *locSigmaGraph;
	double locMinimumCounts = 500;
	double locInitFitRangeMin = -20.0, locInitFitRangeMax = 20.0;
	double locInitGaussWidth = 4.0;

	TObjArray* locMeanArray = new TObjArray();
	TObjArray* locSigmaArray = new TObjArray();

	//Proton
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPhiVsTheta");
	loc2DHist->SetName("DeltaPhiVsTheta_Proton");
	loc2DHist->Rebin2D(4, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 80.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiPlus
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPhiVsTheta");
	loc2DHist->SetName("DeltaPhiVsTheta_PiPlus");
	loc2DHist->Rebin2D(8, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 160.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiMinus
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPhiVsTheta");
	loc2DHist->SetName("DeltaPhiVsTheta_PiMinus");
	loc2DHist->Rebin2D(8, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 160.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//All
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("p");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#plus}");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#minus}");
	locCanvas = locPlotDrawer.Draw_Superimpose(locMeanArray, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Superimpose(locSigmaArray, locCanvasInstructions);
}

/*
.x Load.C
.L effres_plots.C
//Fit_DeltaPhiVsP("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Fit_DeltaPhiVsP(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

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

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

	TGraphErrors *locMeanGraph, *locSigmaGraph;
	double locMinimumCounts = 500;
	double locInitFitRangeMin = -20.0, locInitFitRangeMax = 20.0;
	double locInitGaussWidth = 8.0;

	TObjArray* locMeanArray = new TObjArray();
	TObjArray* locSigmaArray = new TObjArray();

	//Proton
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPhiVsP");
	loc2DHist->SetName("DeltaPhiVsP_Proton");
	loc2DHist->Rebin2D(2, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiPlus
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPhiVsP");
	loc2DHist->SetName("DeltaPhiVsP_PiPlus");
	loc2DHist->Rebin2D(2, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiMinus
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPhiVsP");
	loc2DHist->SetName("DeltaPhiVsP_PiMinus");
	loc2DHist->Rebin2D(4, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//All
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("p");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#plus}");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#minus}");
	locCanvas = locPlotDrawer.Draw_Superimpose(locMeanArray, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Superimpose(locSigmaArray, locCanvasInstructions);
}

/*
.x Load.C
.L effres_plots.C
//Fit_DeltaTheta("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Fit_DeltaTheta(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

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

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

	TGraphErrors *locMeanGraph, *locSigmaGraph;
	double locMinimumCounts = 500;
	double locInitFitRangeMin = -10.0, locInitFitRangeMax = 10.0;
	double locInitGaussWidth = 0.5;

	TObjArray* locMeanArray = new TObjArray();
	TObjArray* locSigmaArray = new TObjArray();

	//Proton
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaThetaVsP");
	loc2DHist->SetName("DeltaThetaVsP_Proton");
	loc2DHist->Rebin2D(2, 1);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiPlus
	locInitGaussWidth = 0.5;
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaThetaVsP");
	loc2DHist->SetName("DeltaThetaVsP_PiPlus");
	loc2DHist->Rebin2D(4, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiMinus
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaThetaVsP");
	loc2DHist->SetName("DeltaThetaVsP_PiMinus");
	loc2DHist->Rebin2D(4, 2);
	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//All
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("p");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#plus}");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#minus}");
	locCanvas = locPlotDrawer.Draw_Superimpose(locMeanArray, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Superimpose(locSigmaArray, locCanvasInstructions);
}

/*
.x Load.C
.L effres_plots.C
//Fit_DeltaP("~/Scratch/gluex/hd_root_trackeff.root");
*/
void Fit_DeltaP(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

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

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dMinPadWidth = 400;
//	locCanvasInstructions->dPadInstructions[0]

	TGraphErrors *locMeanGraph, *locSigmaGraph;
	double locMinimumCounts = 500;
	double locInitFitRangeMin = -0.2, locInitFitRangeMax = 0.2;
	double locInitGaussWidth = 0.005;

	TObjArray* locMeanArray = new TObjArray();
	TObjArray* locSigmaArray = new TObjArray();

	//Proton
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingProton/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPOverPVsP");
	loc2DHist->SetName("DeltaPOverPVsP_Proton");
	loc2DHist->Rebin2D(4, 2);
//	loc2DHist->GetXaxis()->SetRangeUser(0.0, 3.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiPlus
	locInitGaussWidth = 0.02;
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiPlus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPOverPVsP");
	loc2DHist->SetName("DeltaPOverPVsP_PiPlus");
	loc2DHist->Rebin2D(4, 2);
//	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//PiMinus
	loc2DHist = (TH2*)locInputFile->Get("TrackEff_MissingPiMinus/Custom_TrackingEfficiency/Resolution_TimeBased/DeltaPOverPVsP");
	loc2DHist->SetName("DeltaPOverPVsP_PiMinus");
	loc2DHist->Rebin2D(4, 2);
//	loc2DHist->GetXaxis()->SetRangeUser(0.0, 5.0);
	Fit_Slices_DoubleGaussians(loc2DHist, locHistogramSlices, locMeanGraph, locSigmaGraph, locMinimumCounts, locInitFitRangeMin, locInitFitRangeMax, locInitGaussWidth);
	locMeanGraph->SetTitle("");
	locSigmaGraph->SetTitle("");
	locMeanArray->AddLast(locMeanGraph);
	locSigmaArray->AddLast(locSigmaGraph);
	locCanvas = locPlotDrawer.Draw_Object(loc2DHist, locCanvasInstructions);
	locCanvasArray = locPlotDrawer.Draw_Array(locHistogramSlices, locCanvasInstructions);

	//All
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("p");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#plus}");
	locCanvasInstructions->dPadInstructions[0]->dLegendNames.push_back("#pi^{#minus}");
	locCanvas = locPlotDrawer.Draw_Superimpose(locMeanArray, locCanvasInstructions);
	locCanvas = locPlotDrawer.Draw_Superimpose(locSigmaArray, locCanvasInstructions);
}

void Fit_Slices_DoubleGaussians(TH2* locInput2DHist, TObjArray*& locHistogramSlices, TGraphErrors*& locMeanGraph, TGraphErrors*& locSigmaGraph, double locMinimumCounts, double locInitFitRangeMin, double locInitFitRangeMax, double locInitGaussWidth)
{
	vector<double> locPs, locMeans, locMeanUncertainties, locSigmas, locSigmaUncertainties;

	locHistogramSlices = new TObjArray();
	for(int loc_i = 1; loc_i <= locInput2DHist->GetNbinsX(); ++loc_i)
	{
		ostringstream locHistName;
		locHistName << locInput2DHist->GetName() << "_XBin" << loc_i;
		TH1* loc1DHist = locInput2DHist->ProjectionY(locHistName.str().c_str(), loc_i, loc_i);

		double locIntegral = loc1DHist->Integral(1, loc1DHist->GetNbinsX());
		if(locIntegral < locMinimumCounts)
			continue;
//		if(locIntegral < 1500.0) //use for delta-p/p
//			loc1DHist->Rebin(4);
//		else
//			loc1DHist->Rebin(2);
		if(locIntegral < 1500.0) //use for delta-phi
			loc1DHist->Rebin(2);

		ostringstream locHistTitle;
		locHistTitle << locInput2DHist->GetTitle() << ", " << locInput2DHist->GetXaxis()->GetBinLowEdge(loc_i) << " #leq ";
		locHistTitle << locInput2DHist->GetXaxis()->GetTitle() << " < " << locInput2DHist->GetXaxis()->GetBinUpEdge(loc_i);
		loc1DHist->SetTitle(locHistTitle.str().c_str());
		locHistogramSlices->AddLast(loc1DHist);

		//get init height guess
		double locGaussHeight = Calc_InitGaussHeight(loc1DHist, 0.0, locInitGaussWidth, 5.0);

		//get second guess for width, range min/max
		double locCutRangeMin, locCutRangeMax;
		DFitControl* locFitControl = Fit_Gaussian(loc1DHist, locGaussHeight, 0.0, locInitGaussWidth, locInitFitRangeMin, locInitFitRangeMax, locCutRangeMin, locCutRangeMax);
		DFitFunctor_Gaussian* locFitFunctor_Gaussian = (DFitFunctor_Gaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];

		//set second guess for width, range min/max
		if(loc_i != 1)
			locInitGaussWidth = locFitFunctor_Gaussian->dFunc->GetParameter(2);
		double locFitRangeMin = -5.0*locInitGaussWidth;
		double locFitRangeMax = 5.0*locInitGaussWidth;

		//reset the fit
		locFitControl->Remove_FitFunctions(loc1DHist);
		delete locFitControl;

		//fit
		if(fabs(locFitRangeMax - locFitRangeMin) <= loc1DHist->GetBinWidth(1))
			continue;

		locFitControl = Fit_DoubleGauss(loc1DHist, locGaussHeight, 0.0, locInitGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
		DFitFunctor_DoubleGaussian* locFitFunctor_DoubleGaussian = (DFitFunctor_DoubleGaussian*)locFitControl->dSignalFitGroup->dFitFunctors[0];

		//get fit params
		locPs.push_back(locInput2DHist->GetXaxis()->GetBinCenter(loc_i));
		locMeans.push_back(locFitFunctor_DoubleGaussian->dFunc->GetParameter(1));
		locMeanUncertainties.push_back(locFitFunctor_DoubleGaussian->dFunc->GetParError(1));
		pair<double, double> locSigma = locFitFunctor_DoubleGaussian->Calc_StandardDeviation();
		locSigmas.push_back(locSigma.first);
		locSigmaUncertainties.push_back(locSigma.second);
	}

	unsigned int locArraySize = locMeans.size();
	double* locXArray = Convert_VectorToArray(locPs);
	double* locMeanArray = Convert_VectorToArray(locMeans);
	double* locMeanUncertaintyArray = Convert_VectorToArray(locMeanUncertainties);
	double* locSigmaArray = Convert_VectorToArray(locSigmas);
	double* locSigmaUncertaintyArray = Convert_VectorToArray(locSigmaUncertainties);

	locMeanGraph = new TGraphErrors(locMeans.size(), locXArray, locMeanArray, NULL, locMeanUncertaintyArray);
	locMeanGraph->SetName((string(locInput2DHist->GetName()) + string("_Means")).c_str());
	locMeanGraph->SetTitle(locInput2DHist->GetTitle());
	locMeanGraph->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locMeanGraph->GetYaxis()->SetTitle((string("Mean ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());

	locSigmaGraph = new TGraphErrors(locMeans.size(), locXArray, locSigmaArray, NULL, locSigmaUncertaintyArray);
	locSigmaGraph->SetName((string(locInput2DHist->GetName()) + string("_Sigmas")).c_str());
	locSigmaGraph->SetTitle(locInput2DHist->GetTitle());
	locSigmaGraph->GetXaxis()->SetTitle(locInput2DHist->GetXaxis()->GetTitle());
	locSigmaGraph->GetYaxis()->SetTitle((string("#sigma ") + string(locInput2DHist->GetYaxis()->GetTitle())).c_str());
}

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(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 = 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;
}

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

	//Do fit
	DFitControl* locFitControl = new DFitControl(locFitRangeMin, locFitRangeMax, locFitFunctor_Gaussian, 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_Gaussian->Calc_CutRange(3.0, locCutRangeMin, locCutRangeMax);

	return locFitControl;
}

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;
}

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.85*locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = 0.9*locGaussWidth;
	//non-gaussian-tail gaussian
	locInitParams[3] = 0.15*locGaussHeight;  locInitParams[4] = 1.5*locGaussWidth;
	DFitFunctor_DoubleGaussian* locDoubleGaussian = new DFitFunctor_DoubleGaussian(locInitParams, kBlue);
	return locDoubleGaussian;
}

DFitFunctor_Gaussian* Build_Gaussian(TH1* locHist, double locGaussHeight, double locGaussMean, double locGaussWidth)
{
	vector<double> locInitParams(3, 0.0);
	locInitParams[0] = locGaussHeight;  locInitParams[1] = locGaussMean;  locInitParams[2] = locGaussWidth;

	DFitFunctor_Gaussian* locGaussian = new DFitFunctor_Gaussian(locInitParams, kBlue);
	return locGaussian;
}

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* Convert_VectorToArray(const vector<double>& locVector)
{
	double* locArray = new double[locVector.size()];
	for(size_t loc_i = 0; loc_i < locVector.size(); ++loc_i)
		locArray[loc_i] = locVector[loc_i];
	return locArray;
}