#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"

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

/*
.x ../Load.C
.L ~/Code/Work/gluex/ROOT_Scripts/strange_v2/Fit_Lambda.C+
Plot_Misc("~/Scratch/gluex/strange_v2/klambda.root");
*/
void Plot_Misc(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* locPIDHist = (TH2*)locInputFile->Get("Hist_ParticleID/Step0__Photon_Proton_->_K+_Lambda/K+/DeltaTVsP_TOF");
	locPIDHist->GetYaxis()->SetRangeUser(-0.75, 0.75);
	locPlotDrawer.Draw_Object(locPIDHist, locCanvasInstructions);

	TH1* locMassHist = (TH1*)locInputFile->Get("Hist_MissingMass_OffKPlus_PostLambdaCut/MissingMass");
	locPlotDrawer.Draw_Object(locMassHist, locCanvasInstructions);
}

/*
.x ../Load.C
.L ~/Code/Work/gluex/ROOT_Scripts/strange_v2/Fit_Lambda.C+
Fit_Lambda("~/Scratch/gluex/strange_v2/klambda.root");
*/
void Fit_Lambda(string locInputFileName)
{
	DPlotDrawer locPlotDrawer;
	DPlotFitter locPlotFitter;
	locPlotFitter.dQuietFitFlag = true;

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

	DCanvasInstructions* locCanvasInstructions = new DCanvasInstructions("Heya");
	locCanvasInstructions->Set_TargetCanvasAspectRatio(16.0/9.0);
	locCanvasInstructions->dTargetPadAspectRatio = 16.0/9.0;
	locCanvasInstructions->dMinPadWidth = 400;
	locCanvasInstructions->dStyle->SetOptFit(102);
	locCanvasInstructions->dStyle->SetOptStat(10);

	TH1* locMassHist = (TH1*)locInputFile->Get("Hist_InvariantMass_Lambda_KinFit_PostKinFitCut/InvariantMass");

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

	double locFitRangeMin = 1.095;
	double locFitRangeMax = 1.14;
	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;

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

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

	//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] = 1.5*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));
}