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


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

	// Convert from Tesla/cm to Gauss/cm
	dBtot_vs_r_vs_z->Scale(10000.0);

	// Create histogram to hold gradient population
	TH1D *grad = new TH1D("grad", "Inner Tracking Volume", 100, 0.0, 250.0);
	grad->SetXTitle("Max B-field Gradient (Gauss/cm)");
	grad->SetYTitle("Volume fraction");

	// Find bins defining inner volume
	double ZMIN = 17.0;
	double ZMAX = 360.0;
	double RMIN = 0.0;
	double RMAX = 65.0;
	
	int izmin = dBtot_vs_r_vs_z->GetXaxis()->FindBin(ZMIN);
	int izmax = dBtot_vs_r_vs_z->GetXaxis()->FindBin(ZMAX);
	int irmin = dBtot_vs_r_vs_z->GetYaxis()->FindBin(RMIN);
	int irmax = dBtot_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++){

			double r = dBtot_vs_r_vs_z->GetYaxis()->GetBinCenter(jbin);

			double dBtot = dBtot_vs_r_vs_z->GetBinContent(ibin, jbin);
			grad->Fill(dBtot, TMath::TwoPi()*r);
		}
	}

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

	StandardLabels1D(grad,"solenoid_1500_poisson_20100920_01");
	
	// Save
	c1->SaveAs("volume_frac_vs_gradient.pdf");
	c1->SaveAs("volume_frac_vs_gradient.gif");
}

//--------------------
// AddSpecLines
//--------------------
void AddSpecLines(double frac, TH1D *h)
{
	// 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(0.0, 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/2.0, frac+0.005, str);
	lab->SetTextAlign(21);
	lab->SetTextSize(0.025);
	lab->Draw();

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