#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

#include <TROOT.h>
#include <TApplication.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <THStack.h>
#include <TPaveStats.h>
#include <TGaxis.h>
#include <TGraphErrors.h>
#include <TMath.h>
#include <TFitResult.h>
#include <TVirtualFitter.h>

#include <TRandom3.h>
#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"

using namespace std;

const int nmax_calls=1000;
const int npoints=10000;
const int nbin=100;
const double low_edge = 0.;
const double up_edge = 100.;
int icall;
int npfits;
double x_point[npoints];
TH1D *h1_Mn;
TH1D *h2_Mn;
TGraphErrors *g1;
TGraphErrors *g2;
TGraphErrors *gchi2;
TGraphErrors *gscatter;
TCanvas *c_debug;

void myFCN(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *par, Int_t /*iflag */)
{
	npfits = 0;

	// make MC hist
	h2_Mn = new TH1D("h2_Mn", "h2_Mn", nbin, low_edge, up_edge);
	h2_Mn->Sumw2();

	for (unsigned int i=0; i<npoints; i++)
		h2_Mn->Fill(x_point[i]);

	g1 = new TGraphErrors(h1_Mn);
	g2 = new TGraphErrors(h2_Mn);

	/* Multiply by factor */
	for (int i=0;i<g2->GetN();i++) g2->GetX()[i] = par[0]*g2->GetX()[i];

	/* Normalize */
	h2_Mn->Reset();
	for (unsigned int i=0; i<npoints; i++)
	{
	    const double x = h2_Mn->GetBinCenter(i+1);
	    const double c = g2->Eval(x, 0, "s");
		h2_Mn->SetBinContent(i+1, c);
	}
	const double currnorm = h1_Mn->Integral()/h2_Mn->Integral();
	for (int i=0;i<g2->GetN();i++) g2->GetY()[i] = g2->GetY()[i]/currnorm;

	c_debug->cd(1);
	g2->SetLineColor(kRed);
	g2->SetMarkerColor(kRed);
	if(g1->GetMaximum() > g2->GetMaximum())
	{
		g1->SetTitle("Graph from Histograms");
		g1->Draw("AP");
		g2->Draw("P same");
		g1->SetMaximum(1000.);
	}
	else
	{
		g2->SetTitle("Graph from Histograms");
		g2->Draw("AP");
		g1->Draw("P same");
		g2->SetMaximum(1000.);
	}

	for (unsigned int i=0; i<nbin; i++)
	{
//		const double h2cont = h2_Mn->GetBinContent(i+1);
//		const double h2err2 = TMath::Power(h2_Mn->GetBinError(i+1), 2.);
		const double h1cont = h1_Mn->GetBinContent(i+1);
		const double h1err2 = TMath::Power(h1_Mn->GetBinError(i+1), 2.);

	    const double x = h1_Mn->GetBinCenter(i+1);

	    const double h2cont = g2->Eval(x, 0, "s");
	    const double h2err2 = h2cont;

		if (h1cont>0)
		{
			fval += TMath::Power(h2cont-h1cont, 2.)/(h1err2 + h2err2);
			npfits++;
		}
	}

	gchi2->SetPoint(icall, icall, fval);
	gscatter->SetPoint(icall, par[0], fval);

	c_debug->cd(2);
	gchi2->SetMarkerStyle(kFullCircle);
	gchi2->SetMaximum(10000.);
	gchi2->Draw("AP");
	c_debug->cd(3);
	gscatter->SetMarkerStyle(kFullCircle);
	gscatter->SetMaximum(10000.);
	gscatter->GetXaxis()->SetLimits(0., 2.);
	gscatter->Draw("AP");
	c_debug->WaitPrimitive();
	c_debug->Update();

	++icall;

	//cout << par[0] << endl;

	delete h2_Mn;
	delete g2;

	delete g1;
}

//int main()
int x_scale_fit()
{
	//TApplication theApp("App", 0, 0);
	gStyle->SetOptStat("emuro");

	//Fill Hist

	TH1D *h1;
	TRandom3 r1(12345);
	h1 = new TH1D("h1st", "h1st", nbin, low_edge, up_edge);
	h1->Sumw2();
	for (unsigned int i=0; i<npoints; i++)
		h1->Fill(r1.Gaus(25., 15.));
	TRandom3 r2(12345);
	for (unsigned int i=0; i<npoints; i++)
		x_point[i] = r2.Gaus(50., 30.);

	c_debug = new TCanvas();
	c_debug->Divide(3);

	//Fit
	const double initfact = 5.;

	/* Copy Local to Global for Fcn */
	h1_Mn = (TH1D *) h1->Clone("h1_clone");

	double chi2min, fitprob, fitfact, errfitfact, chi2, nvpar;

	gchi2 = new TGraphErrors(nmax_calls);
	gchi2->SetTitle("Chi2 vs calls");
	gscatter = new TGraphErrors(nmax_calls);
	gscatter->SetTitle("Chi2 vs par0");

	for(int iter=0; iter<1; ++iter)
	{
		icall=0;

		const int nfitpar = 1;
		double par0[nfitpar] = {iter==0? initfact : fitfact};

		ROOT::Fit::Fitter fitter;
		fitter.Config().SetParamsSettings(nfitpar,par0);
		fitter.Config().ParSettings(0).SetLimits(0., 10.);
		//fitter.Config().ParSettings(0).SetStepSize(5.-iter*0.05);
		fitter.Config().ParSettings(0).SetStepSize(1.);

		ROOT::Math::MinimizerOptions opt;
		opt.SetMinimizerType("Minuit2");
		opt.SetTolerance(1);
		opt.SetMaxFunctionCalls(1000);
		//opt.SetMaxIterations(1000000);
		opt.SetPrintLevel(3);

//		opt.SetMinimizerAlgorithm("Simplex");
//		opt.Print();
//		fitter.Config().SetMinimizerOptions(opt);
//		fitter.SetFCN(myFCN);
//		fitter.FitFCN();
//		ROOT::Fit::FitResult result;
//		result = fitter.Result();
//		result.Print(std::cout);

		opt.SetMinimizerAlgorithm("Migrad");
		opt.Print();
		fitter.Config().SetMinimizerOptions(opt);
		fitter.SetFCN(myFCN);
		fitter.FitFCN();
		ROOT::Fit::FitResult result;
		result = fitter.Result();
		result.Print(std::cout);

		chi2min = result.Chi2();
		fitprob = result.Prob();
		fitfact = result.Parameter(0);
		errfitfact = result.ParError(0);

		chi2=chi2min;
		nvpar = result.NFreeParameters()+npfits;

		cout << "Fit converged toward value " <<  fitfact << " ± " << errfitfact << " with probability " << fitprob << endl;
		cout << chi2 << " " << npfits-nvpar << " " << chi2/(npfits-nvpar) << endl;
	}

	//	fitter.Config().ParSettings(0).SetStepSize(1.);
	//	opt.SetMinimizerAlgorithm("Simplex");
	//	//opt.Print();
	//	fitter.SetFCN(myFCN);
	//	fitter.FitFCN();
	//	result = fitter.Result();
	//	//result.Print(std::cout);
	//
	//	chi2min = result.Chi2();
	//	fitprob = result.Prob();
	//	fitfact = result.Parameter(0);
	//	errfitfact = result.ParError(0);
	//
	//	chi2=chi2min;
	//	nvpar = result.NFreeParameters()+npfits;
	//
	//	cout << "Fit converged toward value " <<  fitfact << " ± " << errfitfact << " with probability " << fitprob << endl;
	//	cout << chi2 << " " << npfits-nvpar << " " << chi2/(npfits-nvpar) << endl;

	//	opt.SetMinimizerAlgorithm("Migrad");
	//	//opt.Print();
	//	fitter.Config().SetMinimizerOptions(opt);
	//	fitter.SetFCN(myFCN);
	//	fitter.FitFCN();
	//	result = fitter.Result();
	//	//result.Print(std::cout);
	//
	//	opt.SetMinimizerAlgorithm("Minos");
	//	opt.Print();
	//	fitter.Config().SetMinimizerOptions(opt);
	//	fitter.SetFCN(myFCN);
	//	fitter.FitFCN();
	//	result = fitter.Result();
	//	result.Print(std::cout);

	//	const double chi2min = result.Chi2();
	//	const double fitprob = result.Prob();
	//	const double fitfact = result.Parameter(0);
	//	const double errfitfact = result.ParError(0);
	//
	//	const double chi2=chi2min;
	//	const double nvpar = result.NFreeParameters()+npfits;
	//
	//	cout << "Fit converged toward value " <<  fitfact << " ± " << errfitfact << " with probability " << fitprob << endl;
	//	cout << chi2 << " " << npfits-nvpar << " " << chi2/(npfits-nvpar) << endl;

	c_debug->Close();

	// Draw Hist

	TH1D *h2 = new TH1D("h2nd", "h2nd", nbin, low_edge, up_edge);
	for (unsigned int i=0; i<npoints; i++)
		h2->Fill(x_point[i]*fitfact);
	const double fitnorm = h1->Integral()/h2->Integral();

	TCanvas *cc = new TCanvas();

	h2->Scale(fitnorm);
	h1->SetMarkerStyle(kPlus);
	h1->SetLineColor(4);
	h1->SetMarkerColor(4);
	h1->Draw("E0");
	h1->GetXaxis()->SetTitle("");
	h1->GetXaxis()->SetTitleSize(0.05);
	h1->GetXaxis()->SetLabelSize(0.05);
	h1->SetTitle("Comparison");
	h2->SetLineColor(2);
	h2->Draw("sames");
	cc->Update();

	cc->Update();
	cc->WaitPrimitive();

	return 0;
}