#include "StandardLabels.C"
#include "GlueX_boundaries.C"


void volume_frac_vs_theta_outer(void)
{
	gROOT->Reset();
	//gStyle->SetPadRightMargin(0.15);
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();
	
	TFile *f = new TFile("bfield.root");
	TH2D *cos_theta_vs_r_vs_z = (TH2D*)gROOT->FindObject("cos_theta_vs_r_vs_z");

	// Create histogram to hold theta population
	TH1D *theta = new TH1D("theta", "Outer Tracking Volume", 100, 0.0, 1000.0);
	theta->SetXTitle("#theta angle of B-field direction (mrad)");
	theta->SetYTitle("Volume fraction");

	// Find bins defining outer volume
	double ZMIN = 360.0;
	double ZMAX = 617.0;
	double RMIN = 0.0;
	double RMAX = 126.0;
	
	int izmin = cos_theta_vs_r_vs_z->GetXaxis()->FindBin(ZMIN);
	int izmax = cos_theta_vs_r_vs_z->GetXaxis()->FindBin(ZMAX);
	int irmin = cos_theta_vs_r_vs_z->GetYaxis()->FindBin(RMIN);
	int irmax = cos_theta_vs_r_vs_z->GetYaxis()->FindBin(RMAX);
	
	// Loop over bins in region of interest
	for(int ibin=izmin; ibin<=izmax; ibin++){
		for(int jbin=irmin; jbin<=irmax; jbin++){

			// Exclude bins outside bore of magnet
			double r = cos_theta_vs_r_vs_z->GetYaxis()->GetBinCenter(jbin);
			double z = cos_theta_vs_r_vs_z->GetXaxis()->GetBinCenter(ibin);
			double r_max = 65.0 + (r-407.0)*(126.0-65.0)/(617.0-360.0);
			if(z<407.0)r_max = 65.0;
			if(r>r_max)continue;

			double cos_theta = cos_theta_vs_r_vs_z->GetBinContent(ibin, jbin);
			
			double val = fabs(acos(cos_theta) - TMath::Pi())*1000.0; 
			theta->Fill(val, TMath::TwoPi()*r);
		}
	}

	// Integrate the theta histogram
	for(int ibin=2; ibin<=theta->GetNbinsX(); ibin++){
		double sum = theta->GetBinContent(ibin) + theta->GetBinContent(ibin-1);
		theta->SetBinContent(ibin, sum);
	}
	
	// Normalize to get fraction
	theta->Scale(1.0/theta->GetBinContent(theta->GetNbinsX()));
	
	double xmin = 0.0;
	double xmax = 1000.0;
	
	// Draw
	theta->SetStats(0);
	theta->SetLineColor(kBlue);
	theta->SetLineWidth(4);
	theta->GetXaxis()->SetRangeUser(xmin, xmax);
	theta->Draw("C");
	
	AddSpecLines(0.50, theta, xmin, xmax);
	AddSpecLines(0.70, theta, xmin, xmax);
	AddSpecLines(0.85, theta, xmin, xmax);
	AddSpecLines(0.98, theta, xmin, xmax);

	StandardLabels1D(theta,"solenoid_1500_poisson_20090814_01");
	
	// Save
	c1->SaveAs("volume_frac_vs_theta_outer.pdf");
	c1->SaveAs("volume_frac_vs_theta_outer.gif");
}

//--------------------
// AddSpecLines
//--------------------
void AddSpecLines(double frac, TH1D *h, double xmin, double xmax)
{
	// For the given fraction (y-value), look up the x-value in
	// the given histogram and draw horizontal and vertical lines
	// with labels.
	
	// Find x-value by looking for bin whose contents are closest to
	// value in frac.
	int ibin_best = 1;
	for(int ibin=2; ibin<h->GetNbinsX(); ibin++){
		double current = h->GetBinContent(ibin);
		double last = h->GetBinContent(ibin_best);
		if(fabs(current-frac) < fabs(last-frac))ibin_best = ibin;
	}
	double x = h->GetBinCenter(ibin_best);
	
	// Draw horizontal line
	TLine *lin = new TLine(xmin, frac, x, frac);
	lin->SetLineStyle(2);
	lin->SetLineColor(kMagenta);
	lin->SetLineWidth(3);
	lin->Draw();
	
	// Draw vertical line
	TLine *lin = new TLine(x, 0.0, x, frac);
	lin->SetLineStyle(2);
	lin->SetLineColor(kMagenta);
	lin->SetLineWidth(3);
	lin->Draw();
	
	// Draw Horizontal label
	char str[256];
	sprintf(str, "%d%%", 100.0*frac);
	TLatex *lab = new TLatex((x-xmin)/2.0 + xmin, frac+0.005, str);
	lab->SetTextAlign(21);
	lab->SetTextSize(0.025);
	lab->Draw();

	// Draw Vertical label
	char str[256];
	sprintf(str, "%d mrad", x);
	TLatex *lab = new TLatex(x+(xmax-xmin)/100.0, frac/2, str);
	lab->SetTextAngle(270);
	lab->SetTextAlign(21);
	lab->SetTextSize(0.025);
	lab->Draw();
}