#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 <TGraph.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 npoints=10000;
const int nbin=100;
const double low_edge = 0.;
const double up_edge = 100.;
int npfits;
double x_point[npoints];
TH1D *h1_Mn;
TH1D *h2_Mn;

void myFCN(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *par, Int_t /*iflag */)
{
	double chi2 = 0.;
	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]*par[0]);

	const double currnorm = h1_Mn->Integral()/h2_Mn->Integral();
	h2_Mn->Scale(currnorm);

	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.);

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

	chi2 = fval;

	delete h2_Mn;
}

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

	//Fill Hist

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

	//Fit
	const double initfact = 5.;

	const int nfitpar = 1;
	double par0[nfitpar] = {initfact};

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

	ROOT::Fit::Fitter fitter;
	fitter.Config().SetParamsSettings(nfitpar,par0);
	//fitter.Config().ParSettings(0).Set("f",0.5);
	fitter.Config().ParSettings(0).SetLimits(0., 10.);
	fitter.Config().ParSettings(0).SetStepSize(5);

	ROOT::Math::MinimizerOptions opt;
	// print the default minimizer option values
	opt.SetMinimizerType("Minuit2");
	opt.SetTolerance(1);
	opt.SetMaxFunctionCalls(1000000);
	opt.SetMaxIterations(1000000);
	opt.SetPrintLevel(3);

	ROOT::Fit::FitResult result;

	opt.SetMinimizerAlgorithm("Simplex");
	opt.Print();
	fitter.SetFCN(myFCN);
	fitter.FitFCN();
	result = fitter.Result();
	result.Print(std::cout);

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

	opt.SetMinimizerAlgorithm("Minos");
	opt.Print();
	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;

	// Draw Hist

	TH1D *h2;
	h2 = new TH1D("h2", "h2", 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;
}