#include "StandardLabels.C"

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

	// Get tree from file
	// TFile *f = new TFile("bfield_20090814_01.root");
                    TFile *f = new TFile("bfield_20100701_03.root");
	TTree *Bfield = (TTree*)gROOT->FindObject("Bfield");

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

                    // Design concept 2
	TFile *f1 = new TFile("concept2_v1.root");
	TTree *Bfield1 = (TTree*)gROOT->FindObject("Bfield");
	
	// Make Canvas
	TCanvas *c1 = new TCanvas("c1", "boundaries canvas");
	c1->SetTicks();
	c1->SetGrid();

	// Make empty axes
	TH2D *axes = new TH2D("axes", "Maximum #DeltaB_{r}/B_{z} vs. z", 100, -85.0, 650.0, 100, -0.85, 0.85);
	axes->SetStats(0);
	axes->SetXTitle("Z (cm)");
	axes->SetYTitle("|B_{r}(r)-B_{r}(0)|/B_{z}(0)");
                    axes->GetXaxis()->SetRangeUser(0,360);
                    axes->GetYaxis()->SetRangeUser(0,1);
	axes->Draw();

	// Make a histo to project into
	double zmin = -85.0;
	double zmax = 650.0;
	int Nbins = (zmax-zmin)/2.54;
	zmax = zmin + (double)(Nbins)*2.54;
	TH1D *Br_vs_z_r0 = new TH1D("Br_vs_z_r0", "", Nbins, zmin, zmax);
	TH1D *Bzabs_vs_z_r0 = new TH1D("Bzabs_vs_z_r0", "", Nbins, zmin, zmax);
	
	Bfield->Project("Br_vs_z_r0", "z","abs(Br)*(abs(r-2.54)<(0.1))");
	Bfield->Project("Bzabs_vs_z_r0", "z","abs(Bz)*(abs(r-2.54)<(0.1))");

	TH1D *Br_vs_z_r01 = new TH1D("Br_vs_z_r01", "", Nbins, zmin, zmax);
	TH1D *Bzabs_vs_z_r01 = new TH1D("Bzabs_vs_z_r01", "", Nbins, zmin, zmax);
	
	Bfield1->Project("Br_vs_z_r01", "z","abs(Br)*(abs(r-2.54)<(0.1))");
	Bfield1->Project("Bzabs_vs_z_r01", "z","abs(Bz)*(abs(r-2.54)<(0.1))");

	// Make legend
	TLegend *leg = new TLegend(0.60, 0.75, 0.80, 0.85);
	leg->SetFillColor(kWhite);

	// Make histos for r=2.54-55 cm in 1 inch steps
	int Ncolors=9;
	int color[]={kRed,kYellow,kOrange,kSpring, kGreen, kCyan,kBlue,kMagenta,kBlack};
	TH1D *Br_vs_z[100];
	TH1D *Br_vs_z_upper;
	double r_upper = 0.0;
	int Nhists=0;
	// for(double r=2.0*2.54; r<=69.0; r+=3.0*2.54){	
                   for(double r=5.08; r<=69.0; r+=60.92){
		char hname[256];
		sprintf(hname, "Br_vs_z_r%d", (int)(r/2.54+0.5));
		TH1D *h = Br_vs_z[Nhists++] = (TH1D*)Br_vs_z_r0->Clone(hname);
		h->Reset();
		char weight[256];
		sprintf(weight, "(abs(Br)+0.3*r/65)*(abs(r-%3.2f)<(0.1))", r);   
		// sprintf(weight, "(abs(Br)+0.2*r/65)*(abs(r-%3.2f)<(0.1))", r);   // original specification
		Bfield->Project(hname, "z",weight);
		h->Add(Br_vs_z_r0,-1.0);
		h->Divide(Bzabs_vs_z_r0);
		// h->SetLineColor(color[Nhists%Ncolors]);
		h->SetLineColor(4);
		
		if (r>30) h->Draw("same");
                                        // Bzabs_vs_z_r0->Draw("");

		char lab[256];
		sprintf(lab, "r=%3.1f cm", r);
		leg->AddEntry(h, lab);
		
		if(r<=67.0){
			Br_vs_z_upper = h;
			r_upper = r;
		}
                                         if (r > 65) {

                                        //           plot data from identified file

	 	              sprintf(hname, "concept2");
		             TH1D *h1 = (TH1D*)Br_vs_z_r0->Clone(hname);
		             h1->Reset();
		            sprintf(weight, "(abs(Br))*(abs(r-%3.2f)<(0.1))", r);
		            Bfield1->Project(hname, "z",weight);
		            h1->Add(Br_vs_z_r01,-1.0);
		            h1->Divide(Bzabs_vs_z_r01);
                                                    //h1->Scale(-1.0);
		            h1->SetLineColor(2);
		            h1->Draw("same");
                                       }
	}

	leg->Draw();

	// Add standard Labels
	// StandardLabels(axes,"","","solenoid_1500_poisson_20090814_01");
	StandardLabels(axes,"","","solenoid_poisson_concept2");
	
	c1->SaveAs("dBr_rel_vs_z2.png");
	c1->SaveAs("dBr_rel_vs_z2.pdf");

	// Print table of values at selected Z
	ofstream ofile("dBr_rel_vs_z2.txt");
	ofile<<"# For r="<<r<<" cm"<<endl;
	ofile<<"# z   Br_nom   dBr_max"<<endl;
	for(int bin=4; bin<=Nbins; bin+=10){
		char line[1024];
		double z = Br_vs_z_r0->GetXaxis()->GetBinCenter(bin);
		double Br_nominal = Br_vs_z_r0->GetBinContent(bin);
		double dBr_max = Br_vs_z_upper->GetBinContent(bin);
		sprintf(line, "%3.0f\t%3.2f\t%3.2f", z, Br_nominal, dBr_max);
		ofile<<line<<endl;
	}
	ofile.close();
}