void timewalk(void)
{
	gROOT->Reset();
	
	// Get pointer to bcal TTree
	TTree *bcal = (TTree*)gROOT->FindObject("bcal");


       	char tn7_corr[256];
	FitTimeWalk(bcal, "n7", "tphoton",  1500.0,  tn7_corr);
	
	char ts7_corr[256];
	FitTimeWalk(bcal, "s7", "tphoton", 1500.0,  ts7_corr);

	char tn8_corr[256];
	FitTimeWalk(bcal, "n8", "tphoton",  2200.0, tn8_corr);

	char ts8_corr[256];
	FitTimeWalk(bcal, "s8", "tphoton", 1800.0, ts8_corr);

	char tn9_corr[256];
	FitTimeWalk(bcal, "n9", "tphoton", 900.0, tn9_corr);

	char ts9_corr[256];
	FitTimeWalk(bcal, "s9", "tphoton",  1000.0, ts9_corr);
	
	cout<<tn7_corr<<endl;
	cout<<ts7_corr<<endl;
	cout<<tn8_corr<<endl;
	cout<<ts8_corr<<endl;
	cout<<tn9_corr<<endl;
	cout<<ts9_corr<<endl;
}

void FitTimeWalk(TTree *bcal, char *chan, char *ref, double chan_max, char *t_corr)
{
	gROOT->Reset();

	// Create canvas with 2 pads, top for no timewalk correction, bottom corrected
	TCanvas *c1 = new TCanvas("c1","",800,800);
	c1->Divide(1,2);
	
	// Make histogram to hold amplitude vs. time plot.
	TH2F *a_vs_t = new TH2F("a_vs_t", "#Deltat vs. ADC", 400, 0.0, chan_max, 100, -10.0, 10.0);
	a_vs_t->SetStats(0);
	char lab[256];
	sprintf(lab, "ADC %s", chan);
	a_vs_t->SetXTitle(lab);
	sprintf(lab, "t%s + %s - 531 (ns)", chan, ref);
	a_vs_t->SetYTitle(lab);

	// Define timewalk function
	TF1 *tw = new TF1("tw","[0]*exp([1]*x)+[2]+[3]*x");
	tw->SetParName(0, "Expo. Amplitude");
	tw->SetParName(1, "Expo. Slope");
	tw->SetParName(2, "Constant");
	tw->SetParName(3, "Linear Slope");
	tw->SetLineColor(kBlue);
        tw->SetParameters(3.0,-0.005,0.5,-0.0006);

	// Here, we project the time difference between the channel
	// we are correcting ("chan") and a reference channel ("ref")
	// onto the a_vs_t histogram we made above.
	// The reference channel is cut to have a large ADC value
	// (using "ref_cut") so that it's timewalk is small and
	// does not contribute much to the shape.
	char str[256], cut[256];
	sprintf(str, "t%s+%s-531:%s", chan, ref, chan);
	sprintf(cut, "Nphotons==1 && abs(t%s+%s-531) <10",chan, ref);
	bcal->Project("a_vs_t", str, cut, "", 10000000);
	
	// Fit slices of the 2-D histo with gaussians for each
	// bin in the x-axis (ADC-axis). The parameters are placed
	// in 1-D histograms. The one of most interest is a_vs_t_1
	// which holds the means of the gaussians as a function of
	// ADC value. Fit this with the timewalk function defined above.
	a_vs_t->FitSlicesY();
	TH1D *a_vs_t_1 = (TH1D*)gROOT->FindObject("a_vs_t_1");
	a_vs_t_1->Fit("tw","0","", 50.0, chan_max);

	// Draw the uncorrected data that is in the 2-D histo. a_vs_t
	// and overlay the fitted timewalk function.
	c1->cd(1);
	c1->GetPad(1)->SetTickx();
	c1->GetPad(1)->SetTicky();
	a_vs_t->Draw();
	tw->Draw("same");

	// Make a pretty label of the timewalk correction function
	double par[4];
	tw->GetParameters(par);
	sprintf(str,"tw_{%s} =  %f*e^{%f*%s} + %f + (%f*%s)", chan , par[0], par[1], chan, par[2], par[3], chan);
	TLatex *l = new TLatex(0.1*chan_max, 8.0, str);
	l->SetTextSize(0.055);
	l->Draw();
	c1->Update();

	// Use the parameters from the fit to build
	// an expression for the time difference that includes
	// the timewalk correction for "chan"
	// t_corr = (tn8 - (p0*exp(p1*n8) + p2 + p3*n8))
	sprintf(t_corr,"(t%s-(%f*exp(%f*%s)+%f+(%f*%s)))", chan, par[0], par[1], chan, par[2], par[3], chan);
	// str = t_corr-tn7:n8
	sprintf(str, "%s+%s-531:%s>>a_vs_t_corr",t_corr, ref, chan);
	// cut = Nphotons==1 && n7>600 && abs(t_corr-tn7)<10
	sprintf(cut, "Nphotons==1 && abs(%s+%s-531)<10", t_corr, ref);
	TH2F *a_vs_t_corr = new TH2F("a_vs_t_corr", "#Deltat vs. ADC with timewalk correction", 400, 0.0, chan_max, 100, -10.0, 10.0);
	a_vs_t_corr->SetStats(0);
	sprintf(lab, "ADC %s", chan);
	a_vs_t_corr->SetXTitle(lab);
	sprintf(lab, "t%s_{tw corrected} + %s - 531 (ns)", chan, ref);
	a_vs_t_corr->SetYTitle(lab);

	// Draw the corrected data in a_vs_t_corr in the bottom pad
	c1->cd(2);
	c1->GetPad(2)->SetTickx();
	c1->GetPad(2)->SetTicky();
	bcal->Draw(str, cut);
	c1->Update();
	
	sprintf(str, "bcal_timewalk_%s.gif", chan);
	c1->SaveAs(str);
}