#include <stdio.h>
#include <iostream>

#include "TMath.h"
#include "TRandom.h"
#include "TestRooFit_Model.h"

int colour_data = kRed;

int style_mark_data = 20;

//double x_min = -30;
//double x_max = 20;
//int    nbin = 2000;

using namespace RooFit;


void Characterize_hist (TH1F *hist, int colour, int style_mark, double x_min, double x_max, int nbin)
{
	hist -> SetTitle       ("");
	hist -> SetLineColor   (colour);
	hist -> SetMarkerStyle (style_mark);
	hist -> SetMarkerSize  (1);
	hist -> SetMarkerColor (colour);
	hist -> GetXaxis() -> SetTitle       ("ADC counts");
	hist -> GetXaxis() -> SetTitleSize   (0.03);
	hist -> GetXaxis() -> SetTitleOffset (1.2);
	hist -> GetXaxis() -> SetLabelSize   (0.025);
	hist -> GetXaxis() -> SetLabelOffset (0.015);
	hist -> GetXaxis() -> SetRangeUser   (x_min, x_max);
	hist -> GetXaxis() -> SetNdivisions  (nbin, 0);
	hist -> GetYaxis() -> SetTitle       ("Entries");
	hist -> GetYaxis() -> SetTitleSize   (0.03);
	hist -> GetYaxis() -> SetTitleOffset (0.850);
	hist -> GetYaxis() -> SetLabelSize   (0.025);
	hist -> GetYaxis() -> SetLabelOffset (0.005);
}


using namespace RooFit;


void xAna_TestRoofit_Gaussian (double mean = 0, double sigma1 = 1, double sigma2 = 2, double range = 20)
{
	system ("mkdir -p /afs/cern.ch/work/l/lcaophuc/Code_Test/Output/TestRooFit/");
	
	double x_min = mean - range;
	double x_max = mean + range;
	int    nbin  = (int)((x_max-x_min)/0.2);
	
	double  randGaus1, randGaus2, randGausR;
	TRandom *gaus_val1 = new TRandom3 ();
	TRandom *gaus_val2 = new TRandom3 ();
	TRandom *gaus_valR = new TRandom3 ();
	
	TH1F *hist_GausConv   = new TH1F ("hist_GausConv", "", nbin, x_min, x_max);
	TH1F *hist_GausResult = new TH1F ("hist_GausResult", "", nbin, x_min, x_max);
	
	for (long i=0; i<1000000; i++)
	{
		if (i%500000 == 0)   printf ("%9ldth signal particle\n", i+1);
		
		randGaus1 = gaus_val1 -> Gaus (mean, sigma1);
		randGaus2 = gaus_val1 -> Gaus (mean, sigma2);
		
		double sigmaR = sqrt(sigma1*sigma1 + sigma2*sigma2);
		
		randGausR = gaus_valR -> Gaus (mean, sigmaR);
		
		hist_GausConv   -> Fill (randGaus1 + randGaus2);
		hist_GausResult -> Fill (randGausR);
		//hist_Gaus -> Fill (randGaus1);
		//hist_Gaus -> Fill (randGaus2);
	}
	
	//hist_Gaus -> Scale (1/100000);
	
	
	RooRealVar var_x ("var_x", "variable x", x_min, x_max);
	RooDataHist *histroo_GausConv   = new RooDataHist ("histroo_GausConv", "", RooArgList(var_x), hist_GausConv);
	RooDataHist *histroo_GausResult = new RooDataHist ("histroo_GausResult", "", RooArgList(var_x), hist_GausResult);
	
	RooAbsPdf *model_Gaus = Model_Gaussian (&var_x, mean, sigma1, sigma2, 1000000, 1000000);
	
	
	Characterize_hist (hist_GausConv, colour_data, style_mark_data, x_min, x_max, nbin);
	Characterize_hist (hist_GausResult, colour_data, style_mark_data, x_min, x_max, nbin);
	
	
	TF1 *function = model_Gaus -> asTF (RooArgList(var_x));
	double xmax = function -> GetMaximumX (x_min, x_max);
	printf ("The maximum is at: %.14f\n", xmax);
	
	//model_Gaus -> fitTo (*histroo_Gaus, Range(x_min, x_max));
	//RooFitResult *result = model_Gaus -> fitTo (*histroo_Gaus, Range(x_min, x_max));
	
	RooPlot *frame = var_x . frame (-10, 10);
	
	histroo_GausConv   -> plotOn (frame);
	histroo_GausResult -> plotOn (frame);
	model_Gaus         -> plotOn (frame);
	
	frame -> SetYTitle("Events");
	frame -> SetXTitle("ADC counts");
	
	TCanvas *canvas = new TCanvas ("canvas", "", 800, 600);
	canvas -> cd ();
	
	TPad *pad1 = new TPad (Form("pad1"), "", 0.0, 0.0, 1.0, 1.0);
	pad1 -> SetTopMargin    (0.075);
	pad1 -> SetBottomMargin (0.120);
	pad1 -> SetLeftMargin   (0.100);
	pad1 -> SetRightMargin  (0.050);
	pad1 -> SetGrid         (1, 1);
	pad1 -> Draw ("same");
	pad1 -> cd();
	
	frame -> Draw ();
	
	canvas -> SaveAs (Form("/afs/cern.ch/work/l/lcaophuc/Code_Test/Output/TestRooFit/Fit_Gaussian_%.0f.png", range));
}