#include "StandardLabels.C"

void dtheta_vs_theta(void)
{
	TFile *fil = new TFile("rezest.root");
	TH1D *pres_vs_theta = (TH1D*)gROOT->FindObject("pres_vs_theta");
	
	TH1D *dtheta_vs_theta = (TH1D*)pres_vs_theta->Clone("dtheta_vs_theta");
	
	double f = 0.10; // maximum fractional error is allowed to increase due to delta_theta
	
	for(int ibin=1; ibin<=pres_vs_theta->GetNbinsX(); ibin++){
		double theta = pres_vs_theta->GetXaxis()->GetBinCenter(ibin);
		double pres = pres_vs_theta->GetBinContent(ibin);
		double lambda = 3.1415927/2.0 - theta/57.3;
		double lambda_res = pres/tan(lambda)*f*(f+2.0);
		
		lambda_res *= 1000.0; // convert to mrad

		dtheta_vs_theta->SetBinContent(ibin, lambda_res);
	}


	TCanvas *c1 = new TCanvas("c1");
	c1->SetGrid();
	c1->SetTicks();
	
	TH2D *axes = new TH2D("axes", "Tolerance in angular uncertainty", 100, 0.0, 90.0, 100, 0.0, 10.0);
	axes->SetStats(0);
	axes->SetXTitle("#theta (degrees)");
	axes->SetYTitle("#Delta#theta (milliradians)");
	axes->Draw();

	dtheta_vs_theta->GetYaxis()->UnZoom(); // why do I have to do this???
	dtheta_vs_theta->Draw("same");
	
	char str[256];
	sprintf(str, "Tolerance set to increase relative momentum error by %d%%", (int)(100.0*f));
	StandardLabels(axes, str);

	c1->SaveAs("dtheta_vs_theta.png");
	c1->SaveAs("dtheta_vs_theta.pdf");
}