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

//-----------------
// DynamicRange
//-----------------
void DynamicRange(void)
{
	//gStyle->SetStatY(0.90);
	//gStyle->SetOptStat(2200);
	TColor::CreateColorWheel();
	
	// Set energy bands (thrown) to look at
	double Ehi_max = 2.0;
	double Ehi_min = 1.8;
	double Elo_max = 0.140;
	double Elo_min = 0.040;

	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 *h_fADCup_Ehi = new TH2D("h_fADCup_Ehi", "", 100, 0.025, 3.25, MaxSummedLayers(), 0.5, (double)MaxSummedLayers()+0.5);
	TH2D *h_fADCdn_Ehi = (TH2D*)h_fADCup_Ehi->Clone("h_fADCdn_Ehi");
	TH2D *h_fADCup_Elo = new TH2D("h_fADCup_Elo", "", 20, 0.025, 0.15, MaxSummedLayers(), 0.5, (double)MaxSummedLayers()+0.5);
	TH2D *h_fADCdn_Elo = (TH2D*)h_fADCup_Elo->Clone("h_fADCdn_Elo");
	
	// 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 fADCup_Ehi[10]; // index is summed layer number
	double fADCdn_Ehi[10]; // index is summed layer number
	double fADCup_Elo[10]; // index is summed layer number
	double fADCdn_Elo[10]; // index is summed layer number
	int last_event=1;
	for(int i=1; i<=Nentries; i++){
		tree->GetEntry(i);
		
		if(hit.event!=last_event){
			
			// Fill histos and reset variables
			for(int ilayer=1; ilayer<=MaxSummedLayers(); ilayer++){
				h_fADCup_Ehi->Fill(BCAL_fADC2E(fADCup_Ehi[ilayer]), ilayer);
				h_fADCdn_Ehi->Fill(BCAL_fADC2E(fADCdn_Ehi[ilayer]), ilayer);
				h_fADCup_Elo->Fill(BCAL_fADC2E(fADCup_Elo[ilayer]), ilayer);
				h_fADCdn_Elo->Fill(BCAL_fADC2E(fADCdn_Elo[ilayer]), ilayer);
			
				fADCup_Ehi[ilayer] = 0;
				fADCdn_Ehi[ilayer] = 0;
				fADCup_Elo[ilayer] = 0;
				fADCdn_Elo[ilayer] = 0;
			}

			last_event = hit.event;
		}
		
		if(hit.layer<1 || hit.layer>MaxSummedLayers())continue; //
		
		// Remember max ADC for thrown particles in upper energy band
		if(hit.E_thrown<=Ehi_max && hit.E_thrown>=Ehi_min){
			if(hit.fADC_up > fADCup_Ehi[hit.layer])fADCup_Ehi[hit.layer] = hit.fADC_up;
			if(hit.fADC_dn > fADCdn_Ehi[hit.layer])fADCdn_Ehi[hit.layer] = hit.fADC_dn;
		}

		// Remember max ADC for thrown particles in lower energy band
		if(hit.E_thrown<=Elo_max && hit.E_thrown>=Elo_min){
			if(hit.fADC_up > fADCup_Elo[hit.layer])fADCup_Elo[hit.layer] = hit.fADC_up;
			if(hit.fADC_dn > fADCdn_Elo[hit.layer])fADCdn_Elo[hit.layer] = hit.fADC_dn;
		}
	}
	

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	

	//StandardLabels(h, Scheme(), "", AngleStr("#theta_{proton}="));
	
	//c1->SaveAs("DynamicRange.png");
	//c1->SaveAs("DynamicRange.pdf");
	
}