#include "StandardLabels.C"


void dtheta_vs_p_vs_theta_discrete(const char *suffix="")
{
	gROOT->Reset();
	gStyle->SetPalette(1);
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	char fname[256];
	sprintf(fname, "hd_res_charged%s.root", suffix);
	TFile *f = new TFile(fname);
	TH2D *dtheta_sigma = (TH2D*)gROOT->FindObject("dtheta_sigma");
	
	// Draw empty axes
	TH2D *axes = new TH2D("axes", "#theta angle resolution", 100, 0.0, 165, 100, 0.0, 100.0);
	axes->SetStats(0);
	axes->SetXTitle("Polar angle #theta (degrees)");
	axes->SetYTitle("Reconstructed #theta resolution (mrad)");
	axes->Draw();
	
	// Make legend
	TLegend *leg = new TLegend(0.425, 0.718, 0.626, 0.871);
	leg->SetFillColor(kWhite);
	
	// Loop over momenta to plot
	TAxis *yaxis = dtheta_sigma->GetYaxis();
	double mom[]={0.350, 0.5, 1.0, 2.0};
	int Nmom = 4;
	int colors[] = {kBlue, kRed, kGreen, kMagenta};
	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;
		dtheta_sigma->ProjectionX("_px",ybin,ybin);
		dtheta_sigma_px = (TH1D*)gROOT->FindObject("dtheta_sigma_px");

		char hname[256];
		sprintf(hname, "h%d", i);
		TH1D *h = (TH1D*)dtheta_sigma_px->Clone(hname);
		h->SetLineColor(colors[i]);
		h->SetMarkerColor(colors[i]);
		h->SetMarkerStyle(styles[i]);
		
		h->Draw("sameP");
		
		char lab[256];
		sprintf(lab, "%3.3f GeV/c", mom[i]);
		leg->AddEntry(h, lab);
	}
	
	leg->Draw();
	
	// Add standard labeling
	StandardLabels2D(axes,suffix);
	
	char fname[256];
	sprintf(fname, "dtheta_vs_p_vs_theta_discrete%s.gif", suffix);
	c1->SaveAs(fname);
	sprintf(fname, "dtheta_vs_p_vs_theta_discrete%s.pdf", suffix);
	c1->SaveAs(fname);
}