void pres_vs_p_vs_theta_discrete_eugene(const char *suffix="")
{
	gROOT->Reset();
	gStyle->SetPalette(1);
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	c1->SetFillColor(kWhite);
	
	char fname[256];
	sprintf(fname, "hd_res_charged%s.root", suffix);
	TFile *f = new TFile(fname);
	TH2D *dp_over_p_sigma = (TH2D*)gROOT->FindObject("dp_over_p_sigma");
	
	// Draw empty axes
	TH2D *axes = new TH2D("axes", "Total momentum resolution from M.C.", 100, 0.0, 135, 100, 0.0, 8.0);
	axes->SetStats(0);
	axes->SetXTitle("Polar angle #theta (degrees)");
	axes->SetYTitle("Total Momentum Resolution (%)");
	axes->SetLabelSize(0.06, "Y");
	axes->SetLabelSize(0.055, "X");
	axes->GetYaxis()->SetTitleSize(0.06);
	axes->GetYaxis()->SetTitleOffset(0.60);
	axes->GetXaxis()->SetTitleSize(0.05);
	axes->GetXaxis()->SetTitleOffset(0.90);
	axes->Draw();
	
	// Make legend
	TLegend *leg = new TLegend(0.507, 0.530, 0.873, 0.841);
	leg->SetFillColor(kWhite);
	
	// Loop over momenta to plot
	TAxis *yaxis = dp_over_p_sigma->GetYaxis();
	double mom[]={ 1.0, 2.0, 3.0};
	int Nmom = 3;
	int colors[] = {kRed, kBlue, kMagenta, kGreen};
	int styles[] = {20, 21, 22, 23};
	for(int i=0; i<Nmom; i++){
		int ybin = yaxis->FindBin(mom[i]);
		cout<<"Projecting for momentum="<<mom[i]<<" (ybin="<<ybin<<")"<<endl;
		dp_over_p_sigma->ProjectionX("_px",ybin-2,ybin+2);
		dp_over_p_sigma_px = (TH1D*)gROOT->FindObject("dp_over_p_sigma_px");
		dp_over_p_sigma_px->Scale(1.0/5.0);
		
		char hname[256];
		sprintf(hname, "h%d", i);
		TH1D *h = (TH1D*)dp_over_p_sigma_px->Clone(hname);
		h->Rebin(3);
		h->Scale(1.0/3.0);
		h->SetLineColor(colors[i]);
		h->SetMarkerColor(colors[i]);
		h->SetMarkerStyle(styles[i]);
		
		h->Scale(100.0); // convert to percentage
		
		// Fix bin 9 of 3GeV histo
		if(i==2){
			h->SetBinContent(9, (h->GetBinContent(8)+h->GetBinContent(10))/2.0);
		}

		h->Draw("sameP");
		
		h->Fit("pol4","0Q","",22.0, 135.0);
		TF1 *fun = h->GetFunction("pol4");
		TF1 *fun4 = (TF1*)fun->Clone("fun4");
		fun4->SetRange(28.0,135.0);
		fun4->SetLineColor(colors[i]);
		fun4->Draw("same");

		h->Fit("pol5","0Q","", 0.0, 40.0);
		TF1 *fun = h->GetFunction("pol5");
		fun->SetLineColor(colors[i]);
		fun->SetRange(0.0, 28.0);
		fun->Draw("same");

		char lab[256];
		sprintf(lab, "%3.3f GeV/c", mom[i]);
		leg->AddEntry(h, lab);

		// Fix bins 7 and 9 and 50
		h->SetBinContent(7, fun->Eval(h->GetBinCenter(7)));
		if(i==2)h->SetBinContent(9, fun->Eval(h->GetBinCenter(9)));
		h->SetBinContent(50, fun4->Eval(h->GetBinCenter(50)));
		
	}
	
	leg->Draw();
	
	// Add standard labeling
	StandardLabels2D(axes, suffix);
	
	char fname[256];
	sprintf(fname, "pres_vs_p_vs_theta_discrete_eugene%s.gif", suffix);
	c1->SaveAs(fname);
	sprintf(fname, "pres_vs_p_vs_theta_discrete_eugene%s.pdf", suffix);
	c1->SaveAs(fname);
}





//----------------
// StandardLabels2D
//----------------
void StandardLabels2D(TH2D *axes=NULL, string suffix="")
{
	// 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, "August 13, 2009 DL");
	ConvertFromNDC2D(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(33);
	lab->SetTextColor(kBlue);
	lab->Draw();

	// SVN Revision
	lab = new TLatex(0.7, 0.645, "svn revision: 5431");
	ConvertFromNDC2D(lab, axes);
	lab->SetTextSize(0.02);
	lab->SetTextAlign(31);
	lab->SetTextColor(kBlue);
	lab->Draw();
	
	// Event type
	lab = new TLatex(0.43, 0.660, "Single #pi^{+}");
	ConvertFromNDC2D(lab, axes);
	lab->SetTextSize(0.055);
	lab->SetTextAlign(32);
	lab->Draw();

	// Conditions
	if(suffix.length()>0 && suffix[0]=='_')suffix = suffix.substr(1,suffix.length());
	string cond = "Full Reconstruction ("+suffix+")";
	lab = new TLatex(0.7, 0.59, cond.c_str());
	ConvertFromNDC2D(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(31);
	lab->Draw();
}

//----------------
// ConvertFromNDC
//----------------
void ConvertFromNDC2D(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();
	TAxis *yaxis = h->GetYaxis();

	// These seem to be sensitive to values set by "SetRangeUser" as well as pad margins
	double xmin = xaxis->GetBinCenter(xaxis->GetFirst());
	double xmax = xaxis->GetBinCenter(xaxis->GetLast());

	double ymin = yaxis->GetBinCenter(yaxis->GetFirst());
	double ymax = yaxis->GetBinCenter(yaxis->GetLast());

	double x = obj->GetX();
	double y = obj->GetY();
	
	x = xmin + (xmax-xmin)*(0.5+x/1.15);
	y = ymin + (ymax-ymin)*(0.5+y/1.15);
	
	obj->SetX(x);
	obj->SetY(y);

}