#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 "PlotFitter/DFitFunctors_BreitWigners.h"

#include "PlotData/DStatUtilities.h"

//#include "../particleType.h"

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 locCriticalPointX, double locChangeX);
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);

double Fit_Sigma0_Pol(TH1* locMassHist);

void Fix_HistNegatives(TH1* locHist)
{
	for(int loc_i = 1; loc_i <= locHist->GetNbinsX(); ++loc_i)
	{
		if(locHist->GetBinContent(loc_i) < 0.0)
			locHist->SetBinContent(loc_i, 0.0);
	}
}
/*
.x ../Load.C
.L ~/Code/Work/gluex/ROOT_Scripts/strange_v2/Fit_Sigma0.C+
Calc_TSlope("~/Scratch/gluex/strange_v2/ksigma0.root");
*/
void Calc_TSlope(string locInputFileName)
{
	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->SetOptStat(10);

	TH2* loc2DHist = (TH2*)locInputFile->Get("BestCombo/Sigma0MassVsT");
	TH1* locTSlopeHist = loc2DHist->ProjectionX("TProjMassBin", loc2DHist->GetYaxis()->FindBin(1.16001), loc2DHist->GetYaxis()->FindBin(1.2299));

	vector<double> locInitParams(2, 0.0);
	locInitParams[0] = 500.0;  locInitParams[1] = 3.0;
	DFitFunctor_Exponential* locExponential = new DFitFunctor_Exponential(locInitParams);

	//Do fit
	DFitControl* locFitControl = new DFitControl(-1.5, -0.6, locExponential);
	locFitControl->dLogLikelihoodFlag = false;

	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

	int locFitStatus = locPlotFitter.Fit_Hist(locTSlopeHist, locFitControl);

	DPlotDrawer locPlotDrawer;
	locPlotDrawer.Draw_Object(locTSlopeHist, locCanvasInstructions);
}

/*
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Fit_Sigma0.C+
Calc_Yields("final/ksigma0_kinfit.root");
*/
void Calc_Yields(string locInputFileName, bool locCoherentPeakFlag = false);
void Calc_Yields(string locInputFileName, bool locCoherentPeakFlag)
{
	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 = NULL;
	if(locCoherentPeakFlag)
		locMassHist = (TH1*)locInputFile->Get("Hist_InvariantMass_Sigma0_KinFit_PostBeamECut/InvariantMass");
	else
		locMassHist = (TH1*)locInputFile->Get("Hist_InvariantMass_Sigma0_KinFit_PostKinFitCut/InvariantMass");
	locMassHist->Rebin(5);
	locMassHist->GetYaxis()->SetTitle("# Combos / MeV/c^{2}");
	locMassHist->SetTitle("All Beam Energies");

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

	double locFitRangeMin = 1.13;
	double locFitRangeMax = 1.28;
	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;
}

/*
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Fit_Sigma0.C+
Subtract_RFSidebands("final/ksigma0_golden_flux.root");
Subtract_RFSidebands("final/ksigma0_fluxratio_all.root");
*/
void Subtract_RFSidebands(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;

	DCanvasInstructions* locCanvasInstructions_Fit = new DCanvasInstructions("Fit");
	locCanvasInstructions_Fit->dStyle->SetOptFit(102);

	//JUST PICK 1 t-bin: -1.5 -> -0.2
	TFile* locInputFile = new TFile(locInputFileName.c_str(), "READ");
	double locMinT = -1.5, locMaxT = -0.2; //-1.5, -0.2
	int locMinTBin, locMaxTBin;
	//mass hist: 400, phi: 360
	int locMassHistRebin = 8, locPhiHistRebin = 30; //mass = 4 for all, 8 for golden //phi = 15 for all, 30 for golden
	double locSigma0MassCutMin = 1.16, locSigma0MassCutMax = 1.23;

	// GET SIGMA0 SIDEBAND SUBTRACTION FACTOR: PARA
	double locSigma0SubtractFactor_PARA = 0.0;
	{
		TH2* locSigma0MassVsT_RFSignal_PARA = (TH2*)locInputFile->Get("BeamAsymmetry/Sigma0MassVsT_RFSignal_PARA");
		locMinTBin = locSigma0MassVsT_RFSignal_PARA->GetXaxis()->FindBin(locMinT + 0.0001);
		locMaxTBin = locSigma0MassVsT_RFSignal_PARA->GetXaxis()->FindBin(locMaxT - 0.0001);
		TH1* locMassHist_RFSignal_PARA = locSigma0MassVsT_RFSignal_PARA->ProjectionY("MassHist_RFSignal_PARA", locMinTBin, locMaxTBin);
		locMassHist_RFSignal_PARA->Rebin(locMassHistRebin);

		TH2* locSigma0MassVsT_RFSideband_PARA = (TH2*)locInputFile->Get("BeamAsymmetry/Sigma0MassVsT_RFSideband_PARA");
		TH1* locMassHist_RFSideband_PARA = locSigma0MassVsT_RFSideband_PARA->ProjectionY("MassHist_RFSideband_PARA", locMinTBin, locMaxTBin);
		locMassHist_RFSideband_PARA->Rebin(locMassHistRebin);

		//Sideband subtract RF
		locMassHist_RFSignal_PARA->Add(locMassHist_RFSideband_PARA, -0.5);
		//Fix_HistNegatives(locMassHist_RFSignal_PARA);

		locSigma0SubtractFactor_PARA = Fit_Sigma0_Pol(locMassHist_RFSignal_PARA);
		locPlotDrawer.Draw_Object(locMassHist_RFSignal_PARA, locCanvasInstructions_Fit);
	}

	//GET PROD PLANE PHI HIST: PARA
	TH1* locPhiHist_PARA = NULL;
	{
		//Prod Plane Phi: PARA, Sigma0 Signal
		TH2* locProdPlanePhiVsT_RFSignal_Sigma0Signal_PARA = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSignal_Sigma0Signal_PARA");
		TH1* locPhiHist_RFSignal_Sigma0Signal_PARA = locProdPlanePhiVsT_RFSignal_Sigma0Signal_PARA->ProjectionY("PhiHist_RFSignal_Sigma0Signal_PARA", locMinTBin, locMaxTBin);
		locPhiHist_RFSignal_Sigma0Signal_PARA->Rebin(locPhiHistRebin);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSignal_Sigma0Signal_PARA, locCanvasInstructions_Fit);

		TH2* locProdPlanePhiVsT_RFSideband_Sigma0Signal_PARA = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSideband_Sigma0Signal_PARA");
		TH1* locPhiHist_RFSideband_Sigma0Signal_PARA = locProdPlanePhiVsT_RFSideband_Sigma0Signal_PARA->ProjectionY("PhiHist_RFSideband_Sigma0Signal_PARA", locMinTBin, locMaxTBin);
		locPhiHist_RFSideband_Sigma0Signal_PARA->Rebin(locPhiHistRebin);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSideband_Sigma0Signal_PARA, locCanvasInstructions_Fit);

		//Sideband subtract RF
		locPhiHist_RFSignal_Sigma0Signal_PARA = (TH1*)locPhiHist_RFSignal_Sigma0Signal_PARA->Clone("Phi_Sigma0Signal_PARA");
		locPhiHist_RFSignal_Sigma0Signal_PARA->Add(locPhiHist_RFSideband_Sigma0Signal_PARA, -0.5);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSignal_Sigma0Signal_PARA, locCanvasInstructions_Fit);
		locPhiHist_PARA = locPhiHist_RFSignal_Sigma0Signal_PARA;

/*
		//Prod Plane Phi: PARA, Sigma0 Sideband
		TH2* locProdPlanePhiVsT_RFSignal_Sigma0Sideband_PARA = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSignal_Sigma0Sideband_PARA");
		TH1* locPhiHist_RFSignal_Sigma0Sideband_PARA = locProdPlanePhiVsT_RFSignal_Sigma0Sideband_PARA->ProjectionY("PhiHist_RFSignal_Sigma0Sideband_PARA", locMinTBin, locMaxTBin);
		locPhiHist_RFSignal_Sigma0Sideband_PARA->Rebin(locPhiHistRebin);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSignal_Sigma0Sideband_PARA, locCanvasInstructions_Fit);

		TH2* locProdPlanePhiVsT_RFSideband_Sigma0Sideband_PARA = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSideband_Sigma0Sideband_PARA");
		TH1* locPhiHist_RFSideband_Sigma0Sideband_PARA = locProdPlanePhiVsT_RFSideband_Sigma0Sideband_PARA->ProjectionY("PhiHist_RFSideband_Sigma0Sideband_PARA", locMinTBin, locMaxTBin);
		locPhiHist_RFSideband_Sigma0Sideband_PARA->Rebin(locPhiHistRebin);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSideband_Sigma0Sideband_PARA, locCanvasInstructions_Fit);

		//Sideband subtract RF
		locPhiHist_RFSignal_Sigma0Sideband_PARA = (TH1*)locPhiHist_RFSignal_Sigma0Sideband_PARA->Clone("Phi_Sigma0Sideband_PARA");
		locPhiHist_RFSignal_Sigma0Sideband_PARA->Add(locPhiHist_RFSideband_Sigma0Sideband_PARA, -0.5);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSignal_Sigma0Sideband_PARA, locCanvasInstructions_Fit);

		//Prod Plane Phi: Sigma0 Sideband Subtraction
		locPhiHist_PARA = (TH1*)locPhiHist_RFSignal_Sigma0Signal_PARA->Clone("ProdPlanePhi_PARA");
		locPhiHist_PARA->Add(locPhiHist_RFSignal_Sigma0Sideband_PARA, locSigma0SubtractFactor_PARA); */
		locPlotDrawer.Draw_Object(locPhiHist_PARA, locCanvasInstructions_Fit);
	}

	// GET SIGMA0 SIDEBAND SUBTRACTION FACTOR: PERP
	double locSigma0SubtractFactor_PERP = 0.0;
	{
		TH2* locSigma0MassVsT_RFSignal_PERP = (TH2*)locInputFile->Get("BeamAsymmetry/Sigma0MassVsT_RFSignal_PERP");
		locMinTBin = locSigma0MassVsT_RFSignal_PERP->GetXaxis()->FindBin(locMinT + 0.0001);
		locMaxTBin = locSigma0MassVsT_RFSignal_PERP->GetXaxis()->FindBin(locMaxT - 0.0001);
		TH1* locMassHist_RFSignal_PERP = locSigma0MassVsT_RFSignal_PERP->ProjectionY("MassHist_RFSignal_PERP", locMinTBin, locMaxTBin);
		locMassHist_RFSignal_PERP->Rebin(locMassHistRebin);

		TH2* locSigma0MassVsT_RFSideband_PERP = (TH2*)locInputFile->Get("BeamAsymmetry/Sigma0MassVsT_RFSideband_PERP");
		TH1* locMassHist_RFSideband_PERP = locSigma0MassVsT_RFSideband_PERP->ProjectionY("MassHist_RFSideband_PERP", locMinTBin, locMaxTBin);
		locMassHist_RFSideband_PERP->Rebin(locMassHistRebin);

		//Sideband subtract RF
		locMassHist_RFSignal_PERP->Add(locMassHist_RFSideband_PERP, -0.5);
		//Fix_HistNegatives(locMassHist_RFSignal_PERP);

		locSigma0SubtractFactor_PERP = Fit_Sigma0_Pol(locMassHist_RFSignal_PERP);
		locPlotDrawer.Draw_Object(locMassHist_RFSignal_PERP, locCanvasInstructions_Fit);
	}

	//GET PROD PLANE PHI HIST: PERP
	TH1* locPhiHist_PERP = NULL;
	{
		//Prod Plane Phi: PERP, Sigma0 Signal
		TH2* locProdPlanePhiVsT_RFSignal_Sigma0Signal_PERP = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSignal_Sigma0Signal_PERP");
		TH1* locPhiHist_RFSignal_Sigma0Signal_PERP = locProdPlanePhiVsT_RFSignal_Sigma0Signal_PERP->ProjectionY("PhiHist_RFSignal_Sigma0Signal_PERP", locMinTBin, locMaxTBin);
		locPhiHist_RFSignal_Sigma0Signal_PERP->Rebin(locPhiHistRebin);

		TH2* locProdPlanePhiVsT_RFSideband_Sigma0Signal_PERP = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSideband_Sigma0Signal_PERP");
		TH1* locPhiHist_RFSideband_Sigma0Signal_PERP = locProdPlanePhiVsT_RFSideband_Sigma0Signal_PERP->ProjectionY("PhiHist_RFSideband_Sigma0Signal_PERP", locMinTBin, locMaxTBin);
		locPhiHist_RFSideband_Sigma0Signal_PERP->Rebin(locPhiHistRebin);

		//Sideband subtract RF
		locPhiHist_RFSignal_Sigma0Signal_PERP = (TH1*)locPhiHist_RFSignal_Sigma0Signal_PERP->Clone("Phi_Sigma0Signal_PERP");
		locPhiHist_RFSignal_Sigma0Signal_PERP->Add(locPhiHist_RFSideband_Sigma0Signal_PERP, -0.5);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSignal_Sigma0Signal_PERP, NULL);
		locPhiHist_PERP = locPhiHist_RFSignal_Sigma0Signal_PERP;
/*
		//Prod Plane Phi: PERP, Sigma0 Sideband
		TH2* locProdPlanePhiVsT_RFSignal_Sigma0Sideband_PERP = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSignal_Sigma0Sideband_PERP");
		TH1* locPhiHist_RFSignal_Sigma0Sideband_PERP = locProdPlanePhiVsT_RFSignal_Sigma0Sideband_PERP->ProjectionY("PhiHist_RFSignal_Sigma0Sideband_PERP", locMinTBin, locMaxTBin);
		locPhiHist_RFSignal_Sigma0Sideband_PERP->Rebin(locPhiHistRebin);

		TH2* locProdPlanePhiVsT_RFSideband_Sigma0Sideband_PERP = (TH2*)locInputFile->Get("BeamAsymmetry/ProdPlanePhiVsT_RFSideband_Sigma0Sideband_PERP");
		TH1* locPhiHist_RFSideband_Sigma0Sideband_PERP = locProdPlanePhiVsT_RFSideband_Sigma0Sideband_PERP->ProjectionY("PhiHist_RFSideband_Sigma0Sideband_PERP", locMinTBin, locMaxTBin);
		locPhiHist_RFSideband_Sigma0Sideband_PERP->Rebin(locPhiHistRebin);

		//Sideband subtract RF
		locPhiHist_RFSignal_Sigma0Sideband_PERP = (TH1*)locPhiHist_RFSignal_Sigma0Sideband_PERP->Clone("Phi_Sigma0Sideband_PERP");
		locPhiHist_RFSignal_Sigma0Sideband_PERP->Add(locPhiHist_RFSideband_Sigma0Sideband_PERP, -0.5);
//		locPlotDrawer.Draw_Object(locPhiHist_RFSignal_Sigma0Sideband_PERP, NULL);

		//Prod Plane Phi: Sigma0 Sideband Subtraction
		locPhiHist_PERP = (TH1*)locPhiHist_RFSignal_Sigma0Signal_PERP->Clone("ProdPlanePhi_PERP");
		locPhiHist_PERP->Add(locPhiHist_RFSignal_Sigma0Sideband_PERP, locSigma0SubtractFactor_PERP); */
		locPlotDrawer.Draw_Object(locPhiHist_PERP, locCanvasInstructions_Fit);
	}
/*
	P_PARA = 0.440 +/- 0.009 +/- 0.013
	P_PERP = 0.382 +/- 0.008 +/- 0.011
*/
	//FIT PARA/PERP RATIO
	{
		//calc ratio
		double locFluxRatio = 1.04; //+/- 0.05 //GOLDEN PERIOD
		//double locFluxRatio = 0.943; //FULL PERIOD
		TH1* locNumeratorHist = (TH1*)locPhiHist_PARA->Clone("Ratio_Num");
		locNumeratorHist->Add(locPhiHist_PERP, -1.0*locFluxRatio); //numerator
		TH1* locDenominatorHist = (TH1*)locPhiHist_PARA->Clone("Ratio_Denom");
		locDenominatorHist->Add(locPhiHist_PERP, locFluxRatio); //denominator
		TH1* locRatioHist = DStatUtilities::Calc_Ratio(locNumeratorHist, locDenominatorHist);

		/*
		//Fit ratio
		vector<double> locInitParams(2, 0.0);
		double locAmplitude = 0.33; //locRatioHist->GetMaximum();
		locInitParams[0] = locAmplitude;  locInitParams[1] = 0.0;
		cout << "init params = " << locAmplitude << ", " << 0.0 << endl;
		DFitFunctor_Cosine* locCosine = new DFitFunctor_Cosine(locInitParams, true, kRed);

		DPlotFitter locPlotFitter;
		DFitControl* locFitControl = new DFitControl(-180.0, 180.0, locCosine, "Cosine");
		int locFitStatus = locPlotFitter.Fit_Hist(locRatioHist, locFitControl);
		//locRatioHist->GetYaxis()->SetTitle("P_{#gamma} #Sigma");
		*/


		//Fit To: (P_PARA + P_PERP)*Sigma*cos(2*phi)/(2 + (P_PERP - P_PARA)*Sigma*cos(2*phi))
		string locFitString = "(0.440 + 0.382)*[0]*cos(2.0*x*3.14159/180.0)/(2.0 + (0.382 - 0.440)*[0]*cos(2.0*x*3.14159/180.0))";
		TF1* locFunc = new TF1("Ratio_Fit", locFitString.c_str(), -180.0, 180.0);
		locFunc->SetParameter(0, 0.9);
		locFunc->SetLineColor(kRed);
		locRatioHist->Fit(locFunc, "QNR+");
		locRatioHist->GetListOfFunctions()->Add(locFunc);

		pair<double, double> locPolarization(0.401, 0.010);
		pair<double, double> locFitAmplitude(0.3779, 0.2212);
		pair<double, double> locBeamAsymmetry = DStatUtilities::Calc_Ratio(locFitAmplitude, locPolarization);
		cout << "beam asymm = " << locBeamAsymmetry.first << ", " << locBeamAsymmetry.second << endl;

		locCanvasInstructions_Fit->dStyle->SetOptFit(112);
		locPlotDrawer.Draw_Object(locRatioHist, locCanvasInstructions_Fit);
	}
}

double Fit_Sigma0_Pol(TH1* locMassHist)
{
	DPlotFitter locPlotFitter;

	//Do fit
	double locGaussMean = 1.1926;
	double locGaussWidth = 0.005;
	double locNumSigmas = 5.0;
	double locGaussHeight = Calc_InitGaussHeight(locMassHist, locGaussMean, locGaussWidth, locNumSigmas);

	double locFitRangeMin = 1.12;
	double locFitRangeMax = 1.28;
	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;

	//1.16, 1.23
	double locSidebandIntegral = locMassHist->Integral(1, locMassHist->GetXaxis()->FindBin(1.1599));
	locSidebandIntegral += locMassHist->Integral(locMassHist->GetXaxis()->FindBin(1.23001), locMassHist->GetNbinsX());
	locSidebandIntegral *= locMassHist->GetBinWidth(1);
cout << "sideband integral = " << locSidebandIntegral << endl;
	DFitFunctor_Polynomial* locPolynomial = (DFitFunctor_Polynomial*)locFitControl->dBackgroundFitGroup->dFitFunctors[0];
	double locPeakBackgroundIntegral = locPolynomial->dFunc->Integral(1.16, 1.23);
	cout << "peak integral = " << locPeakBackgroundIntegral << endl;
	double locSidebandScaleFactor = -1.0*locPeakBackgroundIntegral/locSidebandIntegral;
	cout << "subtract factor = " << locSidebandScaleFactor << endl;
	return locSidebandScaleFactor;
}

/*
.x $ROOT_ANALYSIS_HOME/Load.C
.L scripts_analysis/Fit_Sigma0.C+
Fit_Sigma0("final/ksigma0.root");
TPOL: Divide amplitude by 0.401 + 0.010
THIS FUNCTION IS OLD AND DEPRECATED
*/
void Fit_Sigma0(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);

	TH2* loc2DHist = (TH2*)locInputFile->Get("BeamAsymmetry/Sigma0MassVsT");
	TH1* locMassHist = loc2DHist->ProjectionY("MassProjTBin", loc2DHist->GetXaxis()->FindBin(-1.4999), loc2DHist->GetXaxis()->FindBin(-0.2001));
	locMassHist->Rebin(4);

	int locRebinAmount = 15; //360 //12: 30 //15: 24

	//Sideband phi
	loc2DHist = (TH2*)locInputFile->Get("BestCombo/ProdPlanePhiVsT_Sideband_PARA");
	loc2DHist->Rebin2D(1, locRebinAmount);
	TH1* locProdPlanePhiHist_Sideband_PARA = loc2DHist->ProjectionY("ProdPlanePhi_Sideband_PARA", loc2DHist->GetXaxis()->FindBin(-1.4999), loc2DHist->GetXaxis()->FindBin(-0.2001));
	locProdPlanePhiHist_Sideband_PARA->Sumw2();
	locPlotDrawer.Draw_Object(locProdPlanePhiHist_Sideband_PARA, locCanvasInstructions);

	loc2DHist = (TH2*)locInputFile->Get("BestCombo/ProdPlanePhiVsT_Sideband_PERP");
	loc2DHist->Rebin2D(1, locRebinAmount);
	TH1* locProdPlanePhiHist_Sideband_PERP = loc2DHist->ProjectionY("ProdPlanePhi_Sideband_PERP", loc2DHist->GetXaxis()->FindBin(-1.4999), loc2DHist->GetXaxis()->FindBin(-0.2001));
	locProdPlanePhiHist_Sideband_PERP->Sumw2();
	locPlotDrawer.Draw_Object(locProdPlanePhiHist_Sideband_PERP, locCanvasInstructions);

	//Signal phi
	loc2DHist = (TH2*)locInputFile->Get("BestCombo/ProdPlanePhiVsT_Signal_PARA");
	loc2DHist->Rebin2D(1, locRebinAmount);
	TH1* locProdPlanePhiHist_Signal_PARA = loc2DHist->ProjectionY("ProdPlanePhi_Signal_PARA", loc2DHist->GetXaxis()->FindBin(-1.4999), loc2DHist->GetXaxis()->FindBin(-0.2001));
	locProdPlanePhiHist_Signal_PARA->Sumw2();

	loc2DHist = (TH2*)locInputFile->Get("BestCombo/ProdPlanePhiVsT_Signal_PERP");
	loc2DHist->Rebin2D(1, locRebinAmount);
	TH1* locProdPlanePhiHist_Signal_PERP = loc2DHist->ProjectionY("ProdPlanePhi_Signal_PERP", loc2DHist->GetXaxis()->FindBin(-1.4999), loc2DHist->GetXaxis()->FindBin(-0.2001));
	locProdPlanePhiHist_Signal_PERP->Sumw2();


	//Do fit
	double locGaussMean = 1.1926;
	double locGaussWidth = 0.015;
	double locNumSigmas = 5.0;
	double locGaussHeight = Calc_InitGaussHeight(locMassHist, locGaussMean, locGaussWidth, locNumSigmas);

	double locFitRangeMin = 1.12;
	double locFitRangeMax = 1.27;
	double locCutRangeMin, locCutRangeMax;
	DFitControl* locFitControl = Fit_DoubleGauss(locMassHist, locGaussHeight, locGaussMean, locGaussWidth, locFitRangeMin, locFitRangeMax, locCutRangeMin, locCutRangeMax);
	locPlotDrawer.Draw_Object(locMassHist, locCanvasInstructions);

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

//1.16, 1.23
	double locSidebandIntegral = locMassHist->Integral(1, locMassHist->GetXaxis()->FindBin(1.1599), "width");
	locSidebandIntegral += locMassHist->Integral(locMassHist->GetXaxis()->FindBin(1.23001), locMassHist->GetNbinsX(), "width");

	DFitFunctor_Polynomial* locPolynomial = (DFitFunctor_Polynomial*)locFitControl->dBackgroundFitGroup->dFitFunctors[0];
	double locPeakBackgroundIntegral = locPolynomial->dFunc->Integral(1.16, 1.23);
	double locSidebandScaleFactor = -1.0*locPeakBackgroundIntegral/locSidebandIntegral;
cout << "side integral, peak back interal, factor = " << locSidebandIntegral << ", " << locPeakBackgroundIntegral << ", " << locSidebandScaleFactor << endl;

	//Subtract sidebands
	locProdPlanePhiHist_Signal_PARA->Add(locProdPlanePhiHist_Sideband_PARA, locSidebandScaleFactor);
	locProdPlanePhiHist_Signal_PERP->Add(locProdPlanePhiHist_Sideband_PERP, locSidebandScaleFactor);
	locPlotDrawer.Draw_Object(locProdPlanePhiHist_Signal_PARA, locCanvasInstructions);
	locPlotDrawer.Draw_Object(locProdPlanePhiHist_Signal_PERP, locCanvasInstructions);

	TH1* locNumeratorHist = (TH1*)locProdPlanePhiHist_Signal_PARA->Clone("Ratio");
	locNumeratorHist->Add(locProdPlanePhiHist_Signal_PERP, -1.0); //numerator
	TH1* locDenominatorHist = (TH1*)locProdPlanePhiHist_Signal_PARA->Clone("Ratio");
	locDenominatorHist->Add(locProdPlanePhiHist_Signal_PERP, 1.0); //denominator
	TH1* locRatioHist = DStatUtilities::Calc_Ratio(locNumeratorHist, locDenominatorHist);
//	TH1* locRatioHist = locProdPlanePhiHist_Signal_PARA->GetAsymmetry(locProdPlanePhiHist_Signal_PERP);

	//Fit ratio
	vector<double> locInitParams(2, 0.0);
	double locAmplitude = locRatioHist->GetMaximum();
	locInitParams[0] = locAmplitude;  locInitParams[1] = 0.0;
	DFitFunctor_Cosine* locCosine = new DFitFunctor_Cosine(locInitParams, true, kRed);

	locFitControl = new DFitControl(-180.0, 180.0, locCosine, "Cosine");
	locFitControl->dLogLikelihoodFlag = true;
	int locFitStatus = locPlotFitter.Fit_Hist(locRatioHist, locFitControl);
	locRatioHist->GetYaxis()->SetTitle("P_{#gamma} #Sigma");

	pair<double, double> locPolarization(0.401, 0.010);
	pair<double, double> locFitAmplitude(0.3779, 0.2212);
	pair<double, double> locBeamAsymmetry = DStatUtilities::Calc_Ratio(locFitAmplitude, locPolarization);
	cout << "beam asymm = " << locBeamAsymmetry.first << ", " << locBeamAsymmetry.second << endl;

	locCanvasInstructions->dStyle->SetOptFit(112);
	locPlotDrawer.Draw_Object(locRatioHist, locCanvasInstructions);
}

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);
	double locCriticalPointX = locFitRangeMax;
	double locChangeX = locFitRangeMin;
	DFitFunctor* locFitFunctor_Background = Build_PolyBackground_Quadratic(locHist, locCriticalPointX, locChangeX);

	//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 locCriticalPointX, double locChangeX)
{
	int locCriticalPointBin = locMassHist->GetXaxis()->FindBin(locCriticalPointX);
	int locChangeBin = locMassHist->GetXaxis()->FindBin(locChangeX);

	double locCriticalPointY = locMassHist->GetBinContent(locCriticalPointBin);
	double locCurve_VerticalChange = locMassHist->GetBinContent(locChangeBin) - locCriticalPointY;
	double locCurve_HorizontalRangeOfVerticalChange = locCriticalPointX - locChangeX;

	double locFlatParam = 0.0, locLinearParam = 0.0, locQuadParam = 0.0;
	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));
}