// This will fit histograms from the tree in hd_root.root to
// obtain a set of timewalk corrections functions.
//
// This is done by projecting the uncorrected times onto a
// histogram in bins of fADC value. The maximum is then used to
// fill a histo and fit it to a specific functional form and code is
// generated to capture those fit results into a C++ accessible
// format: timewalk_max.cc.

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

void timewalk_calibrate_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, 0.0, 60.0);
	axes->SetStats(0);
	axes->SetXTitle("fADC");
	axes->SetYTitle("t (ns)");

	TF1 *fun_up = new TF1("fun_up", "[0] + [1]/(pow(x,[2])+[3])", 0.0, 800.0);
	TF1 *fun_dn = new TF1("fun_dn", "[0] + [1]/(pow(x,[2])+[3])", 0.0, 800.0);
	fun_up->SetLineWidth(2);
	fun_dn->SetLineWidth(2);

	fun_up->SetLineColor(kBlue);
	fun_up->SetParameter(0, 0.037);
	fun_up->SetParameter(1, 1000.0);
	fun_up->SetParameter(2, 1.0);
	fun_up->SetParameter(3, -3.0);

	fun_dn->SetLineColor(kRed);
	fun_dn->SetParameter(0, 0.020);
	fun_dn->SetParameter(1, 500.0);
	fun_dn->SetParameter(2, 1.0);
	fun_dn->SetParameter(3, -3.0);
	
	fun_up->SetParLimits(1, 0.0, 10.0);
	fun_dn->SetParLimits(1, 0.0, 10.0);
	fun_up->SetParLimits(2, 0.0, 5.0);
	fun_dn->SetParLimits(2, 0.0, 5.0);
	
	double parms_up[10+1][4];
	double parms_dn[10+1][4];

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

	// Loop over layers
	int Nlayers=0;
	for(int ilayer=1; ilayer<=10 ; 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, 600, 0.0, 60.0);
		TH2D *hdn_2D = new TH2D("hdn_2D","", 100, 0.0, 800.0, 600, 0.0, 60.0);
		
		char cut[256];
		sprintf(cut, "layer==%d", ilayer);
		tree->Project("hup_2D", "tup:fADC_up", cut);
		tree->Project("hdn_2D", "tdn: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 for layer %d", ilayer);
		axes->SetTitle(title); 

		axes->Draw();
		tree->Draw("tup:fADC_up", cut, "same");
		tree->Draw("tdn:fADC_dn", cut, "same");
		
		fun_up->SetParameter(0, 35.0);
		fun_up->SetParameter(1, 0.05);
		fun_up->SetParameter(2, 0.03);
		fun_up->SetParameter(3, -1.1);

		fun_dn->SetParameter(0, 17.0);
		fun_dn->SetParameter(1, 0.1);
		fun_dn->SetParameter(2, 0.02);
		fun_dn->SetParameter(3, -1.0);

		hup->Fit(fun_up,"", "same", 25.0, 800.0);
		for(int ipar=0; ipar<=3; ipar++)parms_up[ilayer][ipar] = fun_up->GetParameter(ipar);

		hdn->Fit(fun_dn,"", "same", 25.0, 800.0);
		for(int ipar=0; ipar<=3; ipar++)parms_dn[ilayer][ipar] = fun_dn->GetParameter(ipar);

		c1->Update();
		c1->Print("timewalk.ps");		
	}
	
	// Close PDF file 
	c1->Print("timewalk.ps]");

	// Write parameters out in form of C++ file
	const char *cppfname = "timewalk_max.cc";
	ofstream pout(cppfname);
	pout<<endl;
	pout<<"#include<cmath>"<<endl;
	pout<<"using namespace std;"<<endl;
	pout<<endl;
	pout<<"#define NLAYERS_BCAL "<<Nlayers<<endl;
	pout<<"double BCAL_TimeWalkUp(double fADC, int layer)"<<endl;
	pout<<"{"<<endl;
	pout<<"	double tw_par[NLAYERS_BCAL][4] ={"<<endl;
	for(int i=1; i<=Nlayers; i++){
		pout<<"		{";
		for(int j=0; j<=3; j++){
			if(j!=0)pout<<", ";
			pout<<parms_up[i][j];
		}
		pout<<"}";
		if(i<Nlayers)pout<<",";
		pout<<"\t// layer "<<i<<endl;
	}
	pout<<"		};"<<endl;
	pout<<endl;
	pout<<"	layer = layer<1 ? 1:layer>NLAYERS_BCAL ? NLAYERS_BCAL:layer;"<<endl;
	pout<<"	double *p = tw_par[layer-1];"<<endl;
	pout<<"	return p[0]+p[1]/(pow(fADC,p[2])+p[3]);"<<endl;
	pout<<"}"<<endl;
	pout<<endl;
	pout<<endl;
	pout<<"double BCAL_TimeWalkDn(double fADC, int layer)"<<endl;
	pout<<"{"<<endl;
	pout<<"	double tw_par[NLAYERS_BCAL][4] ={"<<endl;
	for(int i=1; i<=Nlayers; i++){
		pout<<"		{";
		for(int j=0; j<=3; j++){
			if(j!=0)pout<<", ";
			pout<<parms_dn[i][j];
		}
		pout<<"}";
		if(i<Nlayers)pout<<",";
		pout<<"\t// layer "<<i<<endl;
	}
	pout<<"		};"<<endl;
	pout<<endl;
	pout<<"	layer = layer<1 ? 1:layer>NLAYERS_BCAL ? NLAYERS_BCAL:layer;"<<endl;
	pout<<"	double *p = tw_par[layer-1];"<<endl;
	pout<<"	return p[0]+p[1]/(pow(fADC,p[2])+p[3]);"<<endl;
	pout<<"}"<<endl;
	pout.close();
	cout<<endl;
	cout<<"Wrote timewalk parameters to \""<<cppfname<<"\""<<endl;
	cout<<endl;
	
	cout<<"Plots written to \"timewalk.ps\". Convert to PDF with:"<<endl;
	cout<<endl;
	cout<<"  ps2pdf timewalk.ps"<<endl;
	cout<<endl;
	cout<<"attempting to run for you ...."<<endl;
	system("ps2pdf timewalk.ps");
	cout<<"Done."<<endl;
}