// 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
// profile histogram as a function of fADC value. The profile
// histos are fit to a specific functional form and code is
// generated to capture those fit results into a C++ accessible
// format: timewalk.cc.

void timewalk_calibrate_profile(void)
{

	TFile *f = new TFile("hd_root.root");
	
	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, 15.0, 55.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_dn->SetLineColor(kRed);
	
	double parms_up[10+1][4];
	double parms_dn[10+1][4];

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

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

		TProfile *hup = new TProfile("hup", "", 100, 0.0, 800.0);
		TProfile *hdn = new TProfile("hdn", "", 100, 0.0, 800.0);

		char cut_up[256];
		char cut_dn[256];
		sprintf(cut_up, "layer==%d && tup>0 && tup<100.0", ilayer);
		sprintf(cut_dn, "layer==%d && tdn>0 && tdn<100.0", ilayer);

		tree->Project("hup", "tup:fADC_up", cut_up);
		tree->Project("hdn", "tdn:fADC_dn", cut_dn);
	
		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_up, "same");
		tree->Draw("tdn:fADC_dn", cut_dn, "same");

		fun_up->SetParameter(0, 29.0);
		fun_up->SetParameter(1, 237.0);
		fun_up->SetParameter(2, 1.2);
		fun_up->SetParameter(3, -35);

		fun_dn->SetParameter(0, 21.0);
		fun_dn->SetParameter(1, 0.4);
		fun_dn->SetParameter(2, 0.061);
		fun_dn->SetParameter(3, -1.2);
		
		hup->Fit(fun_up,"", "same", 20.0, 800.0);
		for(ipar=0; ipar<=3; ipar++)parms_up[ilayer][ipar] = fun_up->GetParameter(ipar);

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

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

	// Write parameters out in form of C++ file
	const char *cppfname = "timewalk_profile.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(unsigned 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(unsigned 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_profile.ps\". Convert to PDF with:"<<endl;
	cout<<endl;
	cout<<"  ps2pdf timewalk_profile.ps"<<endl;
	cout<<endl;
}