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

#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TF1.h"
#include "TLine.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"

#include "function_forum.h"





template<class myaddr>
TBranch *MatchBranchAddress (
		TTree *tree,
		const char *branchname,
		myaddr *address )
{
	TBranch *result = tree -> GetBranch ( branchname );
	result -> SetAddress ( address );
	return result;
}

int  nrechit;
std::vector<int>   *hit_xcoor;
std::vector<int>   *hit_ycoor;
std::vector<float> *hit_energy;

TBranch *branch_nrechit;
TBranch *branch_hit_xcoor;
TBranch *branch_hit_ycoor;
TBranch *branch_hit_energy;

void InitBranch ( TTree *tree )
{
	hit_xcoor  = new std::vector<int>( );
	hit_ycoor  = new std::vector<int>( );
	hit_energy = new std::vector<float>( );

	branch_nrechit = MatchBranchAddress (
			tree, "nrechit", &nrechit );
	branch_hit_xcoor = MatchBranchAddress (
			tree, "hit_xcoor", &hit_xcoor );
	branch_hit_ycoor = MatchBranchAddress (
			tree, "hit_ycoor", &hit_ycoor );
	branch_hit_energy = MatchBranchAddress (
			tree, "hit_energy", &hit_energy );
}





int  spectrum_fitter_forum ( int flux )
{
	ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit", "Simplex");
	time_t   time_start = time(0);
	gStyle -> SetOptStat ( 0 );
	gStyle -> SetOptFit ( 0 );

	double peak = 20;
	double width = 3;
	double pedsigma = 6;
	double physigma = 2;



	std::string input_dir = "Output/Tree/";
	std::string input_path = Form ( "%streehit_%04d.root", input_dir.data( ), flux );
	TFile *file_input = new TFile ( input_path.data( ), "read" );
	TTree *tree_input = ( TTree* )file_input -> Get( "tree_rechit" );
	InitBranch ( tree_input );

	TH1D *hist_rechit = new TH1D ( "hist_rechit", "", 200, 0, 200 );
	hist_rechit -> SetLineColor ( kBlack );
	hist_rechit -> SetLineWidth ( 2 );

	long nentries = tree_input -> GetEntriesFast( );
	TRandom3 *rand_gen = new TRandom3( 0 );
	for ( long en=0; en<nentries; en++ )
	{
		branch_nrechit    -> GetEntry( en );
		branch_hit_xcoor  -> GetEntry( en );
		branch_hit_ycoor  -> GetEntry( en );
		branch_hit_energy -> GetEntry( en );

		double energy = 0;
		for ( int i=0; i<nrechit; i++ )
		{
			int xcoor    = hit_xcoor->at(i);
			int ycoor    = hit_ycoor->at(i);
			if ( xcoor!=0 || ycoor!=0 )   continue;
			energy += hit_energy->at(i);
			energy += rand_gen->Gaus(0, physigma);
		}
		if ( energy == 0 )   continue;
		energy += rand_gen -> Gaus( 0, pedsigma );
		hist_rechit -> Fill ( energy );
	}


	TF1 *function_fit = new TF1 ( "function_fit", pdfsum_multilandau, 0, 90, 6 );
	//TF1 *function_fit = new TF1 ( "function_fit", pdfsum_doublelandau, 0, 90, 6 );
	function_fit -> SetLineColor ( kOrange+7 );
	function_fit -> SetLineWidth ( 2 );
	function_fit -> SetNpx       ( 200 );
	function_fit -> SetParameters ( 25.0, 3.5, 600000.0, flux*1.1, 6.5, 3.0 );
	function_fit -> SetParLimits ( 4, 5.0, 7.0 );
	function_fit -> SetParLimits ( 5, 1.0, 4.0 );
	double fbound1 = 25;
	double fbound2 = 45;
	auto result = hist_rechit -> Fit ( function_fit, "MV0BS", "", fbound1, fbound2 );
	double fpeak  = function_fit->GetParameter( 0 );
	double fwidth = function_fit->GetParameter( 1 );
	double fevent = function_fit->GetParameter( 2 );
	double fflux  = function_fit->GetParameter( 3 );
	double fsig1  = function_fit->GetParameter( 4 );
	double fsig2  = function_fit->GetParameter( 5 );
	double chindf = function_fit->GetChisquare( )/function_fit->GetNDF( );
	printf ( "   mean:  %f\n", fpeak );
	printf ( "   width: %f\n", fwidth );
	printf ( "   area:  %f\n", fevent );
	printf ( "   flux:  %f\n", fflux );
	printf ( "   sig1:  %f\n", fsig1 );
	printf ( "   sig2:  %f\n", fsig2 );

	TLine *vmark1 = new TLine (
			fbound1, 0,
			fbound1, function_fit->Eval(fbound1) );
	vmark1 -> SetLineWidth ( 2 );
	vmark1 -> SetLineStyle ( 2 );
	vmark1 -> SetLineColor ( kTeal+4 );

	TLine *vmark2 = new TLine (
			fbound2, 0,
			fbound2, function_fit->Eval(fbound2) );
	vmark2 -> SetLineWidth ( 2 );
	vmark2 -> SetLineStyle ( 2 );
	vmark2 -> SetLineColor ( kTeal+4 );


	TCanvas *canvas = new TCanvas ( "canvas", "", 800, 600 );
	canvas -> cd( );
	TPad *pad = new TPad ( "pad", "", 0.0, 0.0, 1.0, 1.0 );
	pad -> SetLeftMargin   ( 0.12 );
	pad -> SetRightMargin  ( 0.04 );
	pad -> SetTopMargin    ( 0.08 );
	pad -> SetBottomMargin ( 0.12 );
	pad -> Draw( );
	pad -> cd( );

	hist_rechit  -> Draw ( "hist" );
	clone_rechit -> Draw ( "hist same" );
	function_fit -> Draw ( "same" );
	vmark1 -> Draw ( "same" );
	vmark2 -> Draw ( "same" );

	TLegend *legend = new TLegend ( 0.60, 0.75, 0.95, 0.90 );
	legend -> SetTextSize ( 0.04 );
	legend -> SetTextFont ( 42 );
	legend -> AddEntry ( hist_rechit, "Toy spectrum", "l" );
	legend -> AddEntry ( clone_rechit, "Smoothed spectrum", "l" );
	legend -> AddEntry ( function_fit, "Fit function", "l" );
	legend -> Draw ( "same" );

	TLatex *texresult = new TLatex( );
	texresult -> SetNDC( );
	texresult -> SetTextFont ( 42 );
	texresult -> SetTextSize ( 0.04 );
	texresult -> DrawLatex ( 0.6, 0.70,
			Form ( "peak: %.2f (#color[%d]{%.2f})", fpeak, kOrange+10, 20.0 ) );
	texresult -> DrawLatex ( 0.6, 0.65,
			Form ( "width: %.2f (#color[%d]{%.2f})", fwidth, kOrange+10, 3.0 ) );
	texresult -> DrawLatex ( 0.6, 0.60,
			Form ( "events: %.0f (#color[%d]{%.0f})", fevent, kOrange+10, 500000.0 ) );
	texresult -> DrawLatex ( 0.6, 0.55,
			Form ( "flux: %.2f (#color[%d]{%.1f})", fflux, kOrange+10, float(flux) ) );
	texresult -> DrawLatex ( 0.6, 0.50,
			Form ( "ped sig: %.2f (#color[%d]{%.2f})", fsig1, kOrange+10, 6.0 ) );
	texresult -> DrawLatex ( 0.6, 0.45,
			Form ( "phy sig: %.2f (#color[%d]{%.2f})", fsig2, kOrange+10, 2.0 ) );
	texresult -> DrawLatex ( 0.6, 0.40,
			Form ( "#chi^{2}/ndf: %.4f", chindf ) );

	std::string output_dir = "Output/Plot/";
	system ( Form( "mkdir -p %s", output_dir.data( ) ) );

	std::string plot_path = Form ( "%splotvanilla_%04d.png", output_dir.data( ), flux );
	canvas -> SaveAs ( plot_path.data( ) );

	time_t time_end  = time(0);
	float dt_total = (time_end - time_start) / 3600.0;
	result -> Print( );
	printf ("   Job's done in [ %.3f ] hours\n", dt_total);
	return 0;
}