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

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

#include "Geom.h"
#include "Angle.h"
#include "StandardLabels.C"

void timewalk_check_max(void)
{

	TFile *f = new TFile("hd_root.root");
	f->cd(); // avoid compiler warnings when running with "+"
	
	TTree *tree = (TTree*)gROOT->FindObject("tree");

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TH2D *axes = new TH2D("axes","", 100, 0.0, 800.0, 100, -5.0, 5.0);
	axes->SetStats(0);
	axes->SetXTitle("fADC");
	axes->SetYTitle("t (ns)");
	
	// Open PDF file for plots
	c1->Print("timewalk_check_max.ps[");

	// Loop over layers
	int Nlayers=0;
	//for(int ilayer=1; ilayer<=MaxSummedLayers() ; ilayer++, Nlayers++){
	for(int ilayer=1; ilayer<=MaxSummedLayers() ; ilayer++, Nlayers++){

		// Create histos to hold maxima
		TH1D *hup = new TH1D("hup", "", 100, 0.0, 800.0);
		TH1D *hdn = new TH1D("hdn", "", 100, 0.0, 800.0);
		
		// First, project onto a 2D histo and then project slices
		// of that onto a temporary histo. This speeds things up
		// considerably compared to projecting each slice individually
		TH2D *hup_2D = new TH2D("hup_2D","", 100, 0.0, 800.0, 300, -15.0, 15.0);
		TH2D *hdn_2D = new TH2D("hdn_2D","", 100, 0.0, 800.0, 300, -15.0, 15.0);
		
		char cut[256];
		sprintf(cut, "layer==%d", ilayer);
		tree->Project("hup_2D", "tup_corrected:fADC_up", cut);
		tree->Project("hdn_2D", "tdn_corrected:fADC_dn", cut);
		
		// Loop over bins in fADC
		for(int bin=1; bin<=hup->GetNbinsX(); bin++){
			
			TH1D *hup_2D_py = hup_2D->ProjectionY("_py",bin,bin);
			TH1D *hdn_2D_py = hdn_2D->ProjectionY("_py",bin,bin);
			
			double hup_max = hup_2D_py->GetBinCenter(hup_2D_py->GetMaximumBin());
			double hup_integral = hup_2D_py->Integral();
			double hup_rms = hup_2D_py->GetRMS();
			double hup_max_err = hup_integral>0.0 ? hup_rms/sqrt(hup_integral):0.0;
			hup->SetBinContent(bin, hup_max);
			hup->SetBinError(bin, hup_max_err);

			double hdn_max = hdn_2D_py->GetBinCenter(hdn_2D_py->GetMaximumBin());
			double hdn_integral = hdn_2D_py->Integral();
			double hdn_rms = hdn_2D_py->GetRMS();
			double hdn_max_err = hdn_integral>0.0 ? hdn_rms/sqrt(hdn_integral):0.0;
			hdn->SetBinContent(bin, hdn_max);
			hdn->SetBinError(bin, hdn_max_err);
		}
	
		hup->SetLineColor(kCyan);
		hdn->SetLineColor(kMagenta);
		hup->SetLineWidth(2);
		hdn->SetLineWidth(2);

		char title[256];
		sprintf(title, "Timewalk Check for layer %d", ilayer);
		axes->SetTitle(title); 

		axes->Draw();
		tree->Draw("tup_corrected:fADC_up", cut, "same");
		tree->Draw("tdn_corrected:fADC_dn", cut, "same");
		
		hup->Draw("same");
		hdn->Draw("same");
		
		StandardLabels(axes, Scheme(), AngleStr("#theta_{#gamma}="));

		c1->Update();
		c1->Print("timewalk_check_max.ps");		
		
		char fname[256];
		sprintf(fname, "timewalk_check_max_layer%0d.png", ilayer);
		c1->SaveAs(fname);
	}
	
	// Close PDF file 
	c1->Print("timewalk_check_max.ps]");

	
	cout<<"Plots written to \"timewalk_check_max.ps\". Convert to PDF with:"<<endl;
	cout<<endl;
	cout<<"  ps2pdf timewalk_check_max.ps"<<endl;
	cout<<endl;
	cout<<"attempting to run for you ...."<<endl;
	system("ps2pdf timewalk_check_max.ps");
	cout<<"Done."<<endl;
}