#include <TROOT.h>
#include <TH1.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TColor.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
//-----------------
void tres_avg(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);
	
	// Create histo to hold uncertainty
	TH1D *tavg_err = new TH1D("tavg_err", "BCAL time average uncertainty for full showers", 100, 0.0, 500.0);
	tavg_err->SetXTitle("Uncertainty in t_{avg} for shower (ps)");

	// 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;
	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 = 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; 
			
			tavg_err->Fill(t_tot);
			
			t_tot = 0.0;
			A_tot = 0.0;
			last_event = hit.event;
		}
		
		if(hit.layer<1 || hit.layer>10)continue; // bulletproof
		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;
	}

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	tavg_err->SetFillColor(kAzure);
	tavg_err->SetFillStyle(3000);
	tavg_err->Draw();
	
	StandardLabels1D(tavg_err, "#theta=90^{o}  ;  0#leqE_{#gamma}#leq2GeV","", "hdgeant steps convoluted with electronic pulse shape     threshold=44.7mV");
	
	c1->SaveAs("tres_avg.png");
	c1->SaveAs("tres_avg.pdf");
}