#include <TROOT.h>
#include <TTree.h>
#include <TH2.h>
#include <TGraphErrors.h>
#include <TFile.h>
#include <TF1.h>
#include <TCanvas.h>

#include <iostream>
#include <fstream>
using namespace std;


//------------------------------------
// Copied from JEventProcessor_bcal_timing.cc
typedef struct{
	int event;
	int layer;
	int sector;
	int fADC_up;
	int fADC_dn;
	float Etot;
	float geometric_mean;
	float tup;
	float tdn;
	float tup_corrected;
	float tdn_corrected;
}HIT_t;
//------------------------------------

//-----------------
// tavg_calibrate_tgraph
//-----------------
void tavg_calibrate_tgraph(void)
{

	TFile *f = new TFile("hd_root.root");
	f->cd();
	
	TTree *tree = (TTree*)gROOT->FindObject("tree");

	HIT_t hit;
	HIT_t *hitptr = &hit;
	TBranch *branch = tree->GetBranch("T");
	branch->SetAddress(hitptr);
	vector<float> tavgs[10+1];
	vector<float> geomns[10+1];
	int Nentries = branch->GetEntries();
	for(int i=1; i<=Nentries; i++){
		tree->GetEntry(i);
		if(hit.layer<1 || hit.layer>10)continue; // bulletproof
		if(hit.tup_corrected<-100)continue;
		if(hit.tdn_corrected<-100)continue;
		float tavg = (hit.tup_corrected + hit.tdn_corrected)/2.0;
		float geomn = hit.geometric_mean;
		if(fabs(tavg)>15.0)continue;
		tavgs[hit.layer].push_back(tavg);
		geomns[hit.layer].push_back(geomn);
	}
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TH2D *axes = new TH2D("axes","", 100, 0.0, 1750.0, 100, -15.0, 15.0);
	axes->SetStats(0);
	axes->SetXTitle("fADC");
	axes->SetYTitle("t (ns)");

	TF1 *fun = new TF1("fun", "[0]+log(x)*([1]+log(x)*([2]+log(x)*([3]+log(x)*([4]+log(x)*([5]+log(x)*[6])))))", 0.0, 3500.0);
	
	double parms[10+1][7];

	// Open PDF file for plots
	c1->Print("tavg_tgraph.ps[");

	// Loop over layers
	int Nlayers=0;
	for(int ilayer=1; ilayer<=10 ; ilayer++, Nlayers++){

		TGraph *g = new TGraph(tavgs[ilayer].size(), &geomns[ilayer][0], &tavgs[ilayer][0]);

		char title[256];
		sprintf(title, "#Deltat for layer %d", ilayer);
		axes->SetTitle(title);
		
		fun->SetLineColor(kRed);
		fun->SetParameter(0, 0.0);
		fun->SetParameter(1, 0.0);
		fun->SetParameter(2, 0.0);
		fun->SetParameter(3, 0.0);
		fun->SetParameter(4, 0.0);
		fun->SetParameter(5, 0.0);
		fun->SetParameter(6, 0.0);
		g->Fit("fun","","",30.0, 1750.0);

		axes->Draw();
		g->Draw("Psame");

		for(int ipar=0; ipar<=6; ipar++)parms[ilayer][ipar] = fun->GetParameter(ipar);

		cout<<"Number of points:"<<g->GetN()<<endl;

		c1->Update();
		c1->Print("tavg_tgraph.ps");		
	}
	
	// Close PDF file 
	c1->Print("tavg_tgraph.ps]");

	// Write parameters out in form of C++ file
	const char *cppfname = "tavg_tgraph.cc";
	ofstream pout(cppfname);
	pout<<endl;
	pout<<"#include<cmath>"<<endl;
	pout<<"using namespace std;"<<endl;
	pout<<endl;
	pout<<"#define NLAYERS_BCAL "<<Nlayers<<endl;
	pout<<"double BCAL_tavg(double geometric_mean, int layer)"<<endl;
	pout<<"{"<<endl;
	pout<<"	double tw_par[NLAYERS_BCAL][7] ={"<<endl;
	for(int i=1; i<=Nlayers; i++){
		pout<<"		{";
		for(int j=0; j<=6; j++){
			if(j!=0)pout<<", ";
			pout<<parms[i][j];
		}
		pout<<"}";
		if(i<Nlayers)pout<<",";
		pout<<"\t// layer "<<i<<endl;
	}
	pout<<"		};"<<endl;
	pout<<endl;
	pout<<"	layer = layer<1 ? 1:layer>NLAYERS_BCAL ? NLAYERS_BCAL:layer;"<<endl;
	pout<<"	double *p = tw_par[layer-1];"<<endl;
	pout<<"	double x = log(geometric_mean);"<<endl;
	pout<<"	return p[0]+x*(p[1]+x*(p[2]+x*(p[3]+x*(p[4]+x*(p[5]+x*p[6])))));"<<endl;
	pout<<"}"<<endl;
	pout<<endl;
	pout<<endl;
	pout.close();
	cout<<endl;
	cout<<"Wrote timewalk parameters to \""<<cppfname<<"\""<<endl;
	cout<<endl;
	
	cout<<"Plots written to \"tavg_tgraph.ps\". Convert to PDF with:"<<endl;
	cout<<endl;
	cout<<"  ps2pdf tavg_tgraph.ps"<<endl;
	cout<<endl;
	cout<<"attempting to run for you ...."<<endl;
	system("ps2pdf tavg_tgraph.ps");
	cout<<"Done."<<endl;
}