#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 <TGraphErrors.h>
#include <TF1.h>
#include <TLatex.h>
#include <TProfile.h>
#include <TFitResult.h>

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


// The file tdiff_res.cc is created by the tdiff_check_tgraph.C macro
#include "tdiff_res.cc"
#include "tdiff.cc"
#include "E_calibrate.cc"


#include "Angle.h"
#include "Geom.h"

#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;
	float theta_thrown;
	float E_thrown;
}HIT_t;
//------------------------------------

//-----------------
// tdiff_res
//-----------------
void tdiff_res(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);
	
	TH1D *h = new TH1D("h", "", 200, -5000.0, 5000.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){

			if(A_tot>0.0){
cout<<"t_tot="<<t_tot<<" A_tot="<<A_tot<<"  ratio:"<<t_tot/A_tot<<endl;			
				t_tot /= A_tot;
			
				E_tot = BCAL_fADC2E(E_tot); // convert to GeV

				h->Fill(t_tot);
			}
			t_tot = 0.0;
			A_tot = 0.0;
			E_tot = 0.0;
			last_event = hit.event;
		}
		
		if(hit.layer<1 || hit.layer>MaxSummedLayers())continue; //
		if(hit.tup_corrected<-100)continue;
		if(hit.tdn_corrected<-100)continue;
		
		// Accumulate info for this event
		double geomn = hit.geometric_mean;
		double E = geomn; // convert to GeV just before filling histo above
		double weight = E;
		//double tdiff_corrected = (hit.tup_corrected - hit.tdn_corrected)/2.0 - BCAL_tdiff(geomn, hit.layer);
		//double tdiff_corrected = (hit.tup_corrected - hit.tdn_corrected)/2.0;
		double tdiff_corrected = (hit.tup - hit.tdn)/2.0;
		if(!IsOuter(hit.layer)){
			t_tot += weight*1000.0*tdiff_corrected;
			A_tot += weight;
		}
		E_tot += E;
	}
	
	// Create a fit function
	TF1 *fun = new TF1("fun", "gaus(0)+gaus(3)");
	fun->SetParameter(0, 300.0);
	fun->SetParameter(1,   0.0);
	fun->SetParameter(2, 100.0);
	fun->SetParameter(3, 50.0);
	fun->SetParameter(4,   0.0);
	fun->SetParameter(5, 300.0);
	fun->SetLineColor(kRed);
	fun->SetLineWidth(2);
	TFitResultPtr fitresultptr = h->Fit(fun, "ES");

double A1 = fun->GetParameter(0);
double mean1 = fun->GetParameter(1);
double sigma1 = fun->GetParameter(2);
double A2 = fun->GetParameter(3);
double mean2 = fun->GetParameter(4);
double sigma2 = fun->GetParameter(5);

double mean = (A1*sigma1*mean1 + A2*sigma2*mean2)/(A1*sigma1 + A2*sigma2);
cout<<"mean="<<mean<<endl;

	if(fitresultptr==0){
		TFitResult *fitresult = fitresultptr.Get();
		if(fitresult){
			//fitparm_covariance = fitresult->CovMatrix(1,0);
			//cout<<"Covariance(0,0) = "<<fitresult->CovMatrix(0,0)<<endl;
			//cout<<"Covariance(1,0) = "<<fitresult->CovMatrix(1,0)<<endl;
			//cout<<"Covariance(0,1) = "<<fitresult->CovMatrix(0,1)<<endl;
			//cout<<"Covariance(1,1) = "<<fitresult->CovMatrix(1,1)<<endl;
		}else{
			cout<<"NULL TFitResult !!!"<<endl;
		}
	}else{
		cout<<"******* Bad Fit !!! **********"<<endl;
	}

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	h->Draw();
	
	//char str[256];
	//sprintf(str, "#sigma_{#Deltat/2} = #frac{(%4.1f#pm%4.1f) ps}{#sqrt{E}} #oplus (%4.1f#pm%4.1f) ps"
	//	, fun->GetParameter(1), fun->GetParError(1)
	//	, fun->GetParameter(0), fun->GetParError(0));
	//TLatex *lab = new TLatex(0.4, 250, str);
	//lab->SetTextSize(0.05);
	//lab->Draw();

	StandardLabels1D(h, Scheme(), "", AngleStr("#theta_{#gamma}="));
	
	c1->SaveAs("tdiff_res.png");
	c1->SaveAs("tdiff_res.pdf");
	
	// Write parameter to a file so they can be read in and used
	// to automatically generate the LaTeX table later
	//ofstream ofs("tdiff_parms.dat");
	//ofs<<Scheme()<<endl;
	//ofs<<Angle()<<endl;
	//ofs<<fun->GetParameter(1)<<" "<<fun->GetParError(1)<<endl;
	//ofs<<fun->GetParameter(0)<<" "<<fun->GetParError(0)<<endl;
	//ofs<<fitparm_covariance<<endl;
	//ofs.close();
}