#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 tdiff_res.cc is created by the tdiff_check_tgraph.C macro
#include "tdiff_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_diff_vs_E4
//-----------------
void tres_diff_vs_E4(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, 2.5);
	TH2D *h = new TH2D("h", "", 50, 0.0, 0.85, 40, -600.0, 600.0);
	
	// Loop over all events and fill histogram. We do this
	// rather than TTree::Project since we need to call the 
	// BCAL_tdiff_res for each entry
	int Nentries = branch->GetEntries();
	double t_tot = 0.0;
	double A_tot = 0.0;
	double E_tot = 0.0;
	int last_event=1;
	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 /= A_tot;
			
			h->Fill(E_tot, t_tot);
			
			t_tot = 0.0;
			A_tot = 0.0;
			E_tot = 0.0;
			last_event = hit.event;
		}
		
		if(hit.layer<1 || hit.layer>10)continue; // match 2006 beam test
		//if(hit.sector!=2)continue; // match 2006 beam test
		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_tdiff_res(geomn, hit.layer);
		double weight = 1.0/(sigma_t*sigma_t);
		double E = geomn/2686.0; // convert to GeV
		//double weight = E;
		t_tot += weight*1000.0*(hit.tup_corrected - hit.tdn_corrected)/2.0;
		A_tot += weight;
		E_tot += E;
	}
	
	h->FitSlicesY();
	TH1D *h_2 = (TH1D*)gROOT->FindObject("h_2");
	h_2->SetLineColor(kRed);

	// Create a fit function
	TF1 *fun = new TF1("fun", "sqrt([0]*[0] + pow([1]/sqrt(x) ,2.0))");
	fun->FixParameter(0, 0.0);
	fun->SetParameter(1, 50.0);
	fun->SetLineColor(kRed);
	fun->SetLineWidth(2);
	
	//g->Fit(fun, "", "", 0.01, 2.4);
	
	fun->SetLineColor(kBlue);
	h_2->Fit(fun, "", "", 0.1, 0.8);

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TH2D *axes = new TH2D("axes","", 100, 0.0, 0.85, 100, -600.0, 600.0);
	axes->SetStats(0);
	axes->SetXTitle("Geometric mean (GeV)");
	axes->SetYTitle("Energy weighted #Deltat/2 (ps)");
	axes->Draw();
	
	//g->SetMarkerStyle(8);
	//g->SetMarkerSize(0.25);
;
	//g->Draw("Psame");
	h->Draw("same");
	h_2->Draw("Psame");
	
	char str[256];
	sprintf(str, "#sigma_{#Deltat/2} = #frac{%4.1f ps}{#sqrt{E}} #oplus %4.1f ps", fun->GetParameter(1), fun->GetParameter(0));
	TLatex *lab = new TLatex(0.4, 250, str);
	lab->SetTextSize(0.05);
	lab->Draw();

	StandardLabels(axes, "#theta=90^{o} (center of module) ;  150#leqE_{#gamma}#leq650MeV","", "2006 Beam Test conditions     threshold=44.7mV", "NIM Article Method  weighting: 1/#sigma_{t}^{2}");
	
	c1->SaveAs("tres_diff_vs_E4.png");
	c1->SaveAs("tres_diff_vs_E4.pdf");
}