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



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

//-----------------
// Eres_vs_E
//-----------------
void Eres_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);
	
	TH2D *Ediff = new TH2D("Ediff","", 100, 0.0, 0.750, 100, -0.5, 0.5);
	
	// 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 E_tot = 0.0;
	double G_tot = 0.0;
	int last_event=0;
	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 && last_event!=0){

			E_tot /= 1000.0; // convert to GeV
			
			// See 1GeV dir for how conversion factor was obtained
			//double E = G_tot/2686.0; // convert to GeV
			double E = G_tot/3200.0; // convert to GeV (empirical from eyeballing Ediff)
			
			Ediff->Fill(E, E-E_tot);
			
			E_tot = 0.0;
			G_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
		E_tot += hit.Etot;
		G_tot += hit.geometric_mean;
	}
	
	Ediff->FitSlicesY();
	TH1D *Ediff_2 = (TH1D*)gROOT->FindObject("Ediff_2");
	TH1D *Ediff_over_E = (TH1D*)Ediff_2->Clone("Ediff_over_E");
	
	for(int ibin=2; ibin<=Ediff_over_E->GetNbinsX(); ibin++){
		double E = Ediff_over_E->GetBinCenter(ibin);
		double diff = Ediff_over_E->GetBinContent(ibin);
		Ediff_over_E->SetBinContent(ibin, diff/E);
	}	

	// Create a fit function
	TF1 *fun = new TF1("fun", "sqrt([0]*[0] + pow([1]/sqrt(x) ,2.0))");
	fun->SetParameter(0, 0.013);
	fun->SetParameter(1, 0.042);
	fun->SetLineColor(kRed);
	fun->SetLineWidth(2);
	
	Ediff_over_E->Fit(fun, "", "", 0.02, 1.8);

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TH2D *axes = new TH2D("axes","BCAL Energy resolution vs. reconstructed energy", 100, 0.0, 2.5, 100, 0.0, 0.25);
	axes->SetStats(0);
	axes->SetXTitle("Geometric mean (GeV)");
	axes->SetYTitle("#sigma_{E}/E");
	axes->Draw();

	Ediff_over_E->Draw("Psame");

	char str[256];
	sprintf(str, "#frac{#sigma_{#Delta E}}{E_{recon}} = #frac{%4.1f %%}{#sqrt{E}} #oplus %4.1f %%", fun->GetParameter(1)*100.0, fun->GetParameter(0)*100.0);
	TLatex *lab = new TLatex(1.0, 0.175, str);
	lab->SetTextAlign(22);
	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("Eres_vs_E.png");
	c1->SaveAs("Eres_vs_E.pdf");
}