#include <TROOT.h>
#include <TH1.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TColor.h>
#include <TGraph.h>
#include <TF1.h>
#include <TLatex.h>
#include <TProfile.h>

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


// The file tavg_res.cc is created by the tavg_check_tgraph.C macro
#include "tavg_res.cc"

#include "StandardLabels.C"

//------------------------------------
// 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;
//------------------------------------

//-----------------
// tres_avg_vs_E
//-----------------
void tres_avg_vs_E(void)
{
	gStyle->SetStatY(0.90);
	gStyle->SetOptStat(2200);
	TColor::CreateColorWheel();

	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);
	
	TProfile *h = new TProfile("h","", 400, 0.0, 0.750);
	
	// Loop over all events and fill histogram. We do this
	// rather than TTree::Project since we need to call the 
	// BCAL_tavg_res for each entry
	int Nentries = branch->GetEntries();
	double t_tot = 0.0;
	double A_tot = 0.0;
	double G_tot = 0.0;
	int last_event=1;
	vector<double> tavg_errs;
	vector<double> geomns;
	for(int i=1; i<=Nentries; i++){
		tree->GetEntry(i);
		
		// Fill histogram on event boundaries and reset for next event
		if(hit.event!=last_event){
			t_tot = sqrt(t_tot/pow(A_tot,2.0));
			t_tot *=1000.0;  // convert to ps
			
			// Convert to (t1-t2)/2 since that is proportional to
			// the position and more directly comparable to the
			// time average.
			//t_tot /= 2.0; 
			
			// See 1GeV dir for how conversion factor was obtained
			double E = G_tot/2686.0; // convert to GeV
			
			tavg_errs.push_back(t_tot);
			geomns.push_back(E);
			
			h->Fill(E, t_tot);
			
			t_tot = 0.0;
			A_tot = 0.0;
			G_tot = 0.0;
			last_event = hit.event;
		}
		
		if(hit.layer<1 || hit.layer>4)continue; // mimic 2006 analysis
		if(hit.sector!=2)continue; // mimic 2006 analysis
		if(hit.tup_corrected<-100)continue;
		if(hit.tdn_corrected<-100)continue;
		
		// Accumulate info for this event
		double geomn = hit.geometric_mean;
		double sigma_t = BCAL_tavg_res(geomn, hit.layer);
		double weight = 1.0/(sigma_t*sigma_t);
		t_tot += pow(weight*sigma_t, 2.0);
		A_tot += weight;
		G_tot += geomn;
	}

	// Put data in a TGraph
	TGraph *g = new TGraph(geomns.size(), &geomns[0], &tavg_errs[0]);	

	// Create a fit function
	TF1 *fun = new TF1("fun", "sqrt([0]*[0] + pow([1]/sqrt(x) ,2.0))");
	fun->SetParameter(0, 0.0);
	fun->SetParameter(1, 60.0);
	fun->SetLineColor(kRed);
	fun->SetLineWidth(2);
	
	//g->Fit(fun, "", "", 0.02, 2.4);
	
	fun->SetLineColor(kBlue);
	h->Fit(fun, "", "", 0.02, 0.750);

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TH2D *axes = new TH2D("axes","BCAL Time resolution vs. reconstructed energy", 100, 0.0, 2.5, 100, 0.0, 500.0);
	axes->SetStats(0);
	axes->SetXTitle("Geometric mean (GeV)");
	axes->SetYTitle("Uncertainty in t_{avg} for shower (ps)");
	axes->Draw();
	
	g->SetMarkerStyle(8);
	g->SetMarkerSize(0.25);
;
	g->Draw("Psame");
	h->Draw("same");
	
	char str[256];
	sprintf(str, "#sigma_{t_{avg}} = #frac{%4.1f ps}{#sqrt{E}} #oplus %4.1f ps", fun->GetParameter(1), fun->GetParameter(0));
	TLatex *lab = new TLatex(1.0, 305, str);
	lab->SetTextSize(0.05);
	lab->Draw();
	
	StandardLabels(axes, "#theta=90^{o}  ;  0#leqE_{#gamma}#leq2GeV","", "hdgeant steps convoluted with electronic pulse shape     threshold=44.7mV");
	
	c1->SaveAs("tres_avg_vs_E.png");
	c1->SaveAs("tres_avg_vs_E.pdf");
}