// This will fit histograms from the tree in hd_root.root to
// obtain a set of timewalk corrections functions.
//
// This is done by doing an unbinned fit to the data points
// using TTree::UnbinnedFit(...). 
// Code is generated to capture the fit results into a C++ 
// accessibleformat: timewalk_unbinned.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_unbinned(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);
	
	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<=1 ; ilayer++, Nlayers++){
		
		char cut_up[256];
		char cut_dn[256];
		sprintf(cut_up, "layer==%d && tup>10", ilayer);	
		sprintf(cut_dn, "layer==%d && tdn>10", ilayer);	

		tree->Fit("fun_up","tup:fADC_up", cut_up);
		for(int ipar=0; ipar<=3; ipar++)parms_up[ilayer][ipar] = fun_up->GetParameter(ipar);

		tree->Fit("fun_dn","tdn:fADC_dn", cut_dn);
		for(int ipar=0; ipar<=3; ipar++)parms_dn[ilayer][ipar] = fun_dn->GetParameter(ipar);

		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->Draw("same");
		fun_dn->Draw("same");

		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_unbinned.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;
}