#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 nbin=100;
const double low_edge = 0.;
const double up_edge = 2.;

const int npoints=10000;
double x_point[npoints];

TH1D *h1_Mn;
TH1D *h2_Mn;
TGraph *g1;
TGraph *g2;

void myFCN(Double_t &fval, Double_t *par, Int_t &npfits)
{
	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;

	for (unsigned int i=0; i<nbin; i++)
	{
		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++;
		}
	}

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

	delete h2_Mn;

	delete g1;
	delete g2;
}

//int main()
int chi2vspar0()
{
	//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.);

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

	TH1D *hmin = new TH1D("hmin", "hmin", nbin, low_edge, up_edge);

	const double step = (up_edge-low_edge)/(double)(nbin);
	for(int iter=0; iter<nbin; ++iter)
	{
		Double_t par=step*(1+iter);
		Double_t chi2 = 0.;
		Int_t npfits = 0;
		myFCN(chi2, &par, npfits);
		hmin->Fill(par, chi2);
		cout << iter << " " << par << " " << chi2 << endl;
	}

	TCanvas *cc = new TCanvas();

	hmin->SetMarkerStyle(kPlus);
	hmin->SetLineColor(1);
	hmin->SetMarkerColor(1);
	hmin->Draw("E0");
	hmin->GetXaxis()->SetTitle("");
	hmin->GetXaxis()->SetTitleSize(0.05);
	hmin->GetXaxis()->SetLabelSize(0.05);
	hmin->SetTitle("chi2 vs par0");
	cc->Update();

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

	return 0;
}