TGraphErrors* GetPoints(TH2D *h);

void tres_single(void)
{

	TFile *f = new TFile("hd_root.root");
	
	TTree *tree = (TTree*)gROOT->FindObject("tree");

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	// Open PDF file for plots
	c1->Print("tres_single.ps[");

	// Loop over layers
	int Nlayers=0;
	for(int ilayer=1; ilayer<=10 ; ilayer++, Nlayers++){
	
		// Resolution is determined by fitting slices in fADC
		// values where the slices are determined by having a
		// minimum of 300 events. We project onto a 2D histo
		// first, and then project enough x-bins onto the y-axis
		// to have enough events for a fit.	
	
		TH2D *hup = new TH2D("hup", "", 800, 0.0, 800.0, 100, -5.0, 5.0);
		TH2D *hdn = new TH2D("hdn", "", 800, 0.0, 800.0, 100, -5.0, 5.0);

		char cut[256];
		sprintf(cut, "layer==%d", ilayer);
		tree->Project("hup", "tup_corrected:fADC_up", cut);
		tree->Project("hdn", "tdn_corrected:fADC_dn", cut);
		
		TGraphErrors *gup = GetPoints(hup);
		TGraphErrors *gdn = GetPoints(hdn);

		gup->SetLineColor(kCyan);
		gdn->SetLineColor(kMagenta);
		gup->SetLineWidth(2);
		gdn->SetLineWidth(2);
		gup->SetMarkerColor(gup->GetLineColor());
		gdn->SetMarkerColor(gdn->GetLineColor());
		gup->SetMarkerStyle(22);
		gdn->SetMarkerStyle(21);
		
		gup->GetFunction("fun")->SetLineColor(kBlue);
		gdn->GetFunction("fun")->SetLineColor(kRed);

		TH2D *axes = new TH2D("axes","", 100, 0.0, 800.0, 100, 0.0, ilayer<6 ? 1.0:20);
		axes->SetStats(0);
		axes->SetXTitle("fADC");
		axes->SetYTitle("t (ns)");

		char title[256];
		sprintf(title, "Single ended timing resolution for layer %d", ilayer);
		axes->SetTitle(title); 

		axes->Draw();
		
		gup->Draw("Psame");
		gdn->Draw("Psame");
		
		c1->Update();
		c1.Print("tres_single.ps");		
	}
	
	// Close PDF file 
	c1->Print("tres_single.ps]");

	
	cout<<"Plots written to \"tres_single.ps\". Convert to PDF with:"<<endl;
	cout<<endl;
	cout<<"  ps2pdf tres_single.ps"<<endl;
	cout<<endl;
}

//------------------
// GetPoints
//------------------
TGraphErrors* GetPoints(TH2D *h)
{
	int min_entries = 300;

	// Project onto the x-axis so we can find the limits for
	// for the projection on the y-axis based on number of events
	TH1D *h_px = h->ProjectionX("_px", 1);

	// Loop over (uneven) bins
	int Npoints = 0;
	vector<double> x;
	vector<double> y;
	vector<double> yerr;
	int Nbins = h->GetNbinsX();
	int Nevents_total = h_px->Integral();
	for(int start_bin=1; start_bin<=Nbins; ){

		// Find limits to give us min_entries entries
		int end_bin = start_bin;
		while(h_px->Integral(start_bin, end_bin)<min_entries && end_bin<Nbins){
			++end_bin;
		}
		
		// If only a few events are left, go ahead and include them here.
		if(end_bin<Nbins && h_px->Integral(end_bin, Nbins)<min_entries)end_bin=Nbins;

		// Project onto y-axis
		cout<<"projecting bins "<<start_bin<<" to "<<end_bin<<"  Nentries="<<h_px->Integral(start_bin, end_bin)<<endl;
		TH1D *h_py = h->ProjectionY("_py",start_bin, end_bin);
		
		// Fit to Gaussian
		h_py->Fit("gaus","0");
		//c1->Update();

		// Get fit results
		double sigma = h_py->GetFunction("gaus")->GetParameter(2);
		double sigma_err = h_py->GetFunction("gaus")->GetParError(2);
		double fADC = (h->GetXaxis()->GetBinCenter(start_bin) + h->GetXaxis()->GetBinCenter(end_bin))/2.0;

		x.push_back(fADC);
		y.push_back(sigma);
		yerr.push_back(sigma_err);
		Npoints++;

		// Increment for next iteration
		start_bin=end_bin+1;
	}
	
	TGraphErrors *g = new TGraphErrors(Npoints, &x[0], &y[0], 0, &yerr[0]);
	
	// Fit the graph. If there are less than 8 points, just fit to 
	// a flat line.
	TF1 *fun;
	if(Npoints<8){
		fun  = new TF1("fun", "[0]", 0.0, 800.0);
	}else{
		fun = new TF1("fun", "[0] + [1]/(pow(x,[2])+[3])", 0.0, 800.0);
		fun->SetParameter(0, 0.037);
		fun->SetParameter(1, 1000.0);
		fun->SetParameter(2, 1.0);
		fun->SetParameter(3, -3.0);
	}
	
	g->Fit(fun);
	
	return g;
}