void photon_Eres(void)
{
	gROOT->Reset();
	TColor::CreateColorWheel();

	TFile *f = new TFile("hd_res_photon.root");
	TH2D *dE_over_E_vs_p_vs_theta = (TH2D*)gROOT->FindObject("dE_over_E_vs_p_vs_theta");

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGridy();
	
	TH2D *axes = new TH2D("axes","Reconstructed photon energy resolution vs E_{#gamma}", 100, 0.0, 8.0, 100, 0.0, 0.2);
	axes->SetStats(0);
	axes->SetYTitle("#sigma(#deltaE/E)");
	axes->SetXTitle("E_{#gamma} (GeV)");
	axes->Draw();
	
	TLegend *leg = new TLegend(0.3, 0.5, 0.80, 0.80);

	int colors[] = {kRed, kBlue, kGreen, kMagenta};
	double thetas[] = {3.0, 8.0, 15.0, 30.0};
	int Nthetas=4;
	for(int i=0; i<Nthetas; i++){
		char hname[256];
		sprintf(hname, "dE_over_E_vs_p_%d", (int)thetas[i]);
		int bin = dE_over_E_vs_p_vs_theta->GetXaxis()->FindBin(thetas[i]);
		dE_over_E_vs_p_vs_theta->ProjectionY(hname, bin, bin);
		
		TH1D *h = (TH1D*)gROOT->FindObject(hname);

		h->SetXTitle("E_{#gamma} (GeV)");
		h->SetYTitle("Relative energy resolution");
		h->SetTitle("Reconstructed photon energy resolution vs E_{#gamma}");
		h->SetStats(0);
		h->SetLineColor(colors[i]);
		h->SetLineWidth(2.0);

		char fname[256];
		sprintf(fname, "fun%d", i);
		TF1 *fun = new TF1(fname,"[0]/sqrt(x)+[1]");
		fun->SetLineColor(colors[i]);
		fun->SetParameter(0, 0.05);
		fun->SetParameter(1, 0.02);

		h->Fit(fun,"0");
		fun->SetRange(0.01, 8.0);
		h->Draw("same");
		fun->Draw("same");
		
		char str[256];
		sprintf(str,"#theta=%02d^{o}     %3.1f%%/#sqrt{E} + %3.1f%%", thetas[i], fun->GetParameter(0)*100.0, fun->GetParameter(1)*100.0);
		leg->AddEntry(fun, str);
	}
	
	leg->Draw();
	
	AddStandardLabels(axes);
	
	c1->SaveAs("photon_Eres.gif");
	c1->SaveAs("photon_Eres.pdf");

}

//----------------
// AddStandardLabels
//----------------
void AddStandardLabels(TH2D *axes=NULL)
{
	// This will draw a label or two on the
	// current plot using the NDC coordinates.
	// It is put here to make sure all plots have
	// a consistent labeling.
	
	// Date, Author
	TLatex *lab = new TLatex(0.7, 0.7, "Feb. 15, 2009 DL");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(33);
	lab->Draw();

	// SVN Revision
	lab = new TLatex(0.7, 0.645, "svn revision: 4865");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.02);
	lab->SetTextAlign(31);
	lab->Draw();
	
	// Event type
	lab = new TLatex(0.45, 0.66, "Single #gamma");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(31);
	lab->Draw();

	// Condition
	lab = new TLatex(0.7, 0.61, "E_{recon} #leq #pm20% of E_{thrown}");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(31);
	lab->Draw();

}

//----------------
// ConvertFromNDC
//----------------
void ConvertFromNDC(TLatex *obj, TH2D *h=NULL)
{
	// Bugs in ROOT make it hard to plot labels consistently.
	// For 1D plots, the histogram axes define the coordinate
	// system. For 2D plots, we seem to be forced to use the
	// NDC. There does not seem to be an obvious way to tell
	// which we're using so we pass the information in in the
	// form of the "axes" histogram. If it is not NULL, then
	// we use it to define the limits. Otherwise, we do nothing.

	if(h==NULL)return;

	TAxis *xaxis = h->GetXaxis();
	int Nbinsx = xaxis->GetNbins();
	double xmin = xaxis->GetBinLowEdge(1);
	double xmax = xmin + xaxis->GetBinLowEdge(Nbinsx);

	TAxis *yaxis = h->GetYaxis();
	int Nbinsy = yaxis->GetNbins();
	double ymin = yaxis->GetBinLowEdge(1);
	double ymax = ymin + yaxis->GetBinLowEdge(Nbinsy);
	
	
	double x = obj->GetX();
	double y = obj->GetY();
	cout<<" in:  x="<<x<<"  y="<<y<<endl;
	
	x = xmin + (xmax-xmin)*(0.5+x/1.15);
	y = ymin + (ymax-ymin)*(0.5+y/1.15);
	
	obj->SetX(x);
	obj->SetY(y);

	cout<<"out:  x="<<x<<"  y="<<y<<endl;
}