#include <TROOT.h>
#include <TTree.h>
#include <TH2.h>
#include <TGraphErrors.h>
#include <TFile.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TLatex.h>

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

#include "tdiff_tgraph.cc"

#include "StandardLabels.C"
#include "IsOuter.h"

TGraphErrors* GetPoints(TH2D *h);

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

//-----------------
// tdiff_fit_res
//-----------------
void tdiff_fit_res(void)
{

	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);
	
	// Make histograms to hold tdiff vs. geometric mean
	// for each layer.
	TH2D *tdiff_vs_geomn[10+1];
	for(int ilayer=1; ilayer<=10 ; ilayer++){
		char hname[256];
		sprintf(hname, "tdiff_vs_geomn_lay%02d",ilayer);
		tdiff_vs_geomn[ilayer] = new TH2D(hname, "", 2000, 0.0, 2000.0, 100, -5.0, 5.0);
	}

	int Nentries = branch->GetEntries();
	for(int i=1; i<=Nentries; i++){
		tree->GetEntry(i);
		if(hit.layer<1 || hit.layer>10)continue; // bulletproof
		if(hit.tup_corrected<-100)continue;
		if(hit.tdn_corrected<-100)continue;
		float tdiff = hit.tup_corrected - hit.tdn_corrected;
		float geomn = hit.geometric_mean;
		float tdiff_corrected = tdiff - BCAL_tdiff(geomn, hit.layer);
		tdiff_corrected /= 2.0;
		
		tdiff_vs_geomn[hit.layer]->Fill(geomn, tdiff_corrected);
	}
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TH2D *axes = new TH2D("axes","", 100, 0.0, 200.0, 100, 0.0, 2000.0);
	axes->SetStats(0);
	axes->SetXTitle("fADC (MeV)");
	axes->SetYTitle("#Deltat/2 (ps)");

	// Open PDF file for plots
	c1->Print("tdiff_fit_res.ps[");

	// Loop over layers
	int Nlayers=0;
	TF1 *sigma_vs_geomn[10+1];
	for(int ilayer=1; ilayer<=10 ; ilayer++, Nlayers++){

		char title[256];
		sprintf(title, "#Deltat/2 for layer %d", ilayer);
		axes->SetTitle(title);
		
		// Get TGraph and uncertainty function
		TGraphErrors *gres = GetPoints(tdiff_vs_geomn[ilayer]);
		if(gres->GetN()<4)continue;
		TF1 *tmp_fun = gres->GetFunction("fun");
		TF1 *fun = (TF1*)tmp_fun->Clone("fun"); // make copy so we can refit to antoher function below
		
		// Draw axes
		axes->Draw();

		// Plot fADC histo
		TH1D *tmp_px = tdiff_vs_geomn[ilayer]->ProjectionX("tmp_px");
		double xmax = tmp_px->GetXaxis()->GetXmax();
		tmp_px->GetXaxis()->SetLimits(0.0, xmax/2686.0*1000.0);
		tmp_px->Rebin(4);
		tmp_px->Scale(600.0/tmp_px->GetMaximum());
		tmp_px->SetFillStyle(3001);
		tmp_px->SetFillColor(kGray);
		tmp_px->SetLineColor(kBlack);
		tmp_px->Draw("same");

		// Draw data points
		gres->Draw("PEsame");
		
		// Fit using 1/sqrt(E)
		TF1 *fun2 = new TF1("fun2", "[0]/sqrt(x/1000.0)", 0.0, 800.0);
		fun2->SetParameter(0, 0.0);
		fun2->SetLineColor(kBlue);
		gres->Fit(fun2,"", "0", 10.0, 180.0);
		fun2->Draw("same");
		
		// Draw uncertainty function
		fun->Draw("same");
		
		if(fun){
			char funname[256];
			sprintf(funname, "sigma_vs_geomn_lay%02d", ilayer);
			sigma_vs_geomn[ilayer] = (TF1*)fun->Clone(funname);
			
			char str[256];
			sprintf(str, "#sigma_{#Deltat/2} = %4.2fps + #frac{%4.1gps}{E^{%4.2f} + %4.2f}", fun->GetParameter(0), fun->GetParameter(1), fun->GetParameter(2), fun->GetParameter(3));
			TLatex *lab = new TLatex(100.0, 1100.0, str);
			lab->SetTextAlign(22);
			lab->SetTextSize(0.05);
			lab->SetTextColor(kRed);
			lab->Draw();
			
		}else{
			sigma_vs_geomn[ilayer] = NULL;
		}
		
		// Label 1/sqrt(E) function
		char str[256];
		sprintf(str, "#sigma_{#Deltat/2} = #frac{%4.1gps}{#sqrt{E}}", fun2->GetParameter(0));
		TLatex *lab = new TLatex(100.0, 1500.0, str);
		lab->SetTextAlign(22);
		lab->SetTextSize(0.05);
		lab->SetTextColor(kBlue);
		lab->Draw();
		
		StandardLabels(axes);

		c1->Update();
		c1->Print("tdiff_fit_res.ps");
		
		char fname[256];
		sprintf(fname, "tdiff_fit_res_layer%02d.png", ilayer);
		c1->SaveAs(fname);
		
		delete tmp_px;
	}
	
	// Close PDF file 
	c1->Print("tdiff_fit_res.ps]");
	
	cout<<"Plots written to \"tdiff_fit_res.ps\". Convert to PDF with:"<<endl;
	cout<<endl;
	cout<<"  ps2pdf tdiff_fit_res.ps"<<endl;
	cout<<endl;
	cout<<"attempting to run for you ...."<<endl;
	system("ps2pdf tdiff_fit_res.ps");
	cout<<endl;

}


//------------------
// GetPoints
//------------------
TGraphErrors* GetPoints(TH2D *h)
{
	int min_entries = 300;

	// Project onto the x-axis so we can find the limits for
	// for the projection on the y-axis based on number of events
	TH1D *h_px = h->ProjectionX("_px", 1);

	// Loop over (uneven) bins
	int Npoints = 0;
	vector<double> x;
	vector<double> y;
	vector<double> yerr;
	int Nbins = h->GetNbinsX();
	for(int start_bin=1; start_bin<=Nbins; ){

		// Find limits to give us min_entries entries
		int end_bin = start_bin;
		while(h_px->Integral(start_bin, end_bin)<min_entries && end_bin<Nbins){
			++end_bin;
		}
		
		// If only a few events are left, go ahead and include them here.
		if(end_bin<Nbins && h_px->Integral(end_bin, Nbins)<min_entries)end_bin=Nbins;

		// Project onto y-axis
		cout<<"projecting bins "<<start_bin<<" to "<<end_bin<<"  Nentries="<<h_px->Integral(start_bin, end_bin)<<endl;
		TH1D *h_py = h->ProjectionY("_py",start_bin, end_bin);

		if(h_py->Integral()>=20){
		
			// Fit to Gaussian
			h_py->Fit("gaus","0");
			//c1->Update();

			// Get fit results
			double sigma = h_py->GetFunction("gaus")->GetParameter(2);
			double sigma_err = h_py->GetFunction("gaus")->GetParError(2);

			// Find average fADC value for hits in this bin range
			double fADC = 0.0;
			double norm = 0.0;
			for(int bin=start_bin; bin<=end_bin; bin++){
				double weight = h_px->GetBinContent(bin);
				fADC += h_px->GetBinCenter(bin)*weight;
				norm += weight;
			}
			fADC /= norm;
			
			// Add floor term to error obtained from Gaussian fit
			double epsilon = 0.020; // ns
			sigma_err = sqrt(sigma_err*sigma_err + epsilon*epsilon);

			x.push_back(fADC/2686*1000.0);
			y.push_back(sigma*1000.0);
			yerr.push_back(sigma_err);
			Npoints++;
		}

		// Increment for next iteration
		start_bin=end_bin+1;
	}
	
	TGraphErrors *g = new TGraphErrors(Npoints, &x[0], &y[0], 0, &yerr[0]);
	g->SetMarkerColor(kBlack);
	g->SetMarkerStyle(22);
	g->SetMarkerSize(1.5);
	g->SetLineColor(kMagenta);
	
	// Fit the graph. If there are less than 8 points, just fit to 
	// a flat line.
	TF1 *fun;
	if(Npoints<8){
		fun  = new TF1("fun", "[0]", 0.0, 800.0);
	}else{
		fun = new TF1("fun", "[0] + [1]/(pow(x/1000.0,[2])+[3])", 0.0, 800.0);
		fun->SetParameter(0, 0.0);
		fun->SetParameter(1, 1.0);
		fun->SetParameter(2, 0.5);
		fun->SetParameter(3, 0.0);
	}
	
	fun->SetLineColor(kRed);
	g->Fit(fun,"", "", 10.0, 180.0);
	
	return g;
}