#include "StandardLabels.C"

void Bz_vs_z_beamline2(void)
{
	gROOT->Reset();
	TColor::CreateColorWheel();

	// Make Canvas
	TCanvas *c1 = new TCanvas("c1", "boundaries canvas",200,10,700,500);
	c1->SetTicks();
	c1->SetGrid();

	// Make empty axes
	// TH2D *axes = new TH2D("axes", "B_{z} vs. z along magnet axis", 100, -85, 650, 100, 0.0, 3.0);
	TH2D *axes = new TH2D("axes", "B_{z} vs. z along magnet axis", 100, -85, 550, 100, 0.0, 3.0);
	axes->SetStats(0);
	axes->SetXTitle("Z (cm)");
	axes->SetYTitle("B_{z} (T)");

	// Get trees from files
                    /* LASS field
	TFile *f = new TFile("bfield_20090814_01.root");
	TTree *Bfield_20090814_01 = (TTree*)gROOT->FindObject("Bfield");*/

                    // MIT conceptual design field
	TFile *f2 = new TFile("bfield_MIT20111104_01.root");
	TTree *Bfield_MIT20111104_01 = (TTree*)gROOT->FindObject("Bfield");

                    /* symmetric field
	TFile *f = new TFile("bfield_20100701_03.root");
	TTree *Bfield_20100701_03 = (TTree*)gROOT->FindObject("Bfield");*/

                    // Design concept 2
	TFile *f = new TFile("concept2_v1.root");
	TTree *Bfield_concept2 = (TTree*)gROOT->FindObject("Bfield");

	// Make a histo to project into for LASS field
	double zmin = -85.0;
	double zmax = 660.0;
	int Nbins = (zmax-zmin)/2.54;
	zmax = zmin + (double)(Nbins)*2.54;

	TH1D *Bz_vs_z_MIT20111104_01 = new TH1D("Bz_vs_z_MIT20111104_01", "", Nbins, zmin, zmax);
	Bz_vs_z_MIT20111104_01->SetMarkerStyle(20);
	Bz_vs_z_MIT20111104_01->SetMarkerColor(kRed);
	Bz_vs_z_MIT20111104_01->SetLineColor(Bz_vs_z_MIT20111104_01->GetMarkerColor());
	Bz_vs_z_MIT20111104_01->SetFillColor(kBlue);
	Bz_vs_z_MIT20111104_01->SetLineWidth(3);
	Bz_vs_z_MIT20111104_01->SetFillStyle(3001);
	Bfield_MIT20111104_01->Project("Bz_vs_z_MIT20111104_01", "z","abs(Bz)*(r<=2.55)");

	TH1D *Bz_vs_z_concept2 = new TH1D("Bz_vs_z_concept2", "", Nbins, zmin, zmax);
	Bz_vs_z_concept2->SetMarkerStyle(20);
	Bz_vs_z_concept2->SetMarkerColor(kRed);
	Bz_vs_z_concept2->SetLineColor(Bz_vs_z_concept2->GetMarkerColor());
	Bz_vs_z_concept2->SetFillColor(kBlue);
	Bz_vs_z_concept2->SetLineWidth(3);
	Bz_vs_z_concept2->SetFillStyle(3001);
	Bfield_concept2->Project("Bz_vs_z_concept2", "z","1.0*abs(Bz)*(r<=2.55)");  // scale by factor

	// Do similar for LASS field less 1kG
	TH1D *Bz_vs_z_lower = (TH1D*)Bz_vs_z_MIT20111104_01->Clone("Bz_vs_z_lower");
	Bz_vs_z_lower->SetMarkerColor(kBlue);
	Bz_vs_z_lower->SetLineColor(Bz_vs_z_lower->GetMarkerColor());
	//Bz_vs_z_lower->SetFillStyle(3000);
	Bz_vs_z_lower->SetFillColor(255);
	Bz_vs_z_lower->Reset();
	Bfield_MIT20111104_01->Project("Bz_vs_z_lower", "z","(abs(Bz)-0.100)*(r<=2.55)");

	// And for upper limit field map
	TH1D *Bz_vs_z_upper = (TH1D*)Bz_vs_z_MIT20111104_01->Clone("Bz_vs_z_upper");
	Bz_vs_z_upper->SetMarkerColor(kBlue);
	Bz_vs_z_upper->SetLineColor(Bz_vs_z_upper->GetMarkerColor());
	Bz_vs_z_upper->SetFillColor(kBlue);
	Bfield_MIT20111104_01->Project("Bz_vs_z_upper", "z","(abs(Bz)+0.100)*(r<=2.55)");
	
	// Draw everything
	axes->Draw();
	Bz_vs_z_upper->Draw("same");
	Bz_vs_z_concept2->Draw("same");
	// Bz_vs_z_MIT20111104_01->Draw("same"); concept 2
	Bz_vs_z_lower->Draw("same");

	
	// Bottom histo uses white fill color that eats away part of top axis. Redraw it
	TLine *lin = new TLine(axes->GetXaxis()->GetXmin(), 0.0, axes->GetXaxis()->GetXmax(), 0.0);
	lin->Draw();
	TLine *lin = new TLine(axes->GetXaxis()->GetXmin(), 0.0, axes->GetXaxis()->GetXmin(), -3.0);
	lin->Draw();

	// Add standard Labels
	StandardLabels(axes,"","","solenoid_poisson_concept2");
	
	c1->SaveAs("Bz_vs_z_beamline2.eps"); // causes crashes ???
	c1->SaveAs("Bz_vs_z_beamline2.png");

	// Print table of values at selected Z
	ofstream ofile("Bz_vs_z_beamline2.txt");
	ofile<<"# z   Bz_min   Bz_nom   Bz_max"<<endl;
	for(int bin=4; bin<=Nbins; bin+=10){
		char line[1024];
		double z = Bz_vs_z_concept2->GetXaxis()->GetBinCenter(bin);
		double Bz_min = Bz_vs_z_lower->GetBinContent(bin);
		double Bz_nominal = Bz_vs_z_concept2->GetBinContent(bin);
		double Bz_max = Bz_vs_z_upper->GetBinContent(bin);
		sprintf(line, "%3.0f\t%3.2f\t%3.2f\t%3.2f", z, Bz_min, Bz_nominal, Bz_max);
		ofile<<line<<endl;
	}
	ofile.close();
}