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

#include <TFile.h>
#include <TH2.h>
#include <TROOT.h>

#include <time.h>
#include <stdlib.h>

#include <iostream>
#include <sstream>
#include <fstream>
#include <map>
using namespace std;

map<double, double> FindSteps(double r, TH2D *h, double z_max=550.0);
double GetStandardStepSize(double step_size);
void PrintLatexTable(ostream &ofs, double r_inches, map<double,double> &points);
void PrintDataPoints(double r_inches, map<double,double> &points);

//--------------------
// step_map
//--------------------
void step_map(void)
{
	
	new TFile("bfield.root");
	TH2D *dBtot_vs_r_vs_z = (TH2D*)gROOT->FindObject("dBtot_vs_r_vs_z");
	
	TH2D *step_size_vs_r_vs_z = (TH2D*)dBtot_vs_r_vs_z->Clone("step_size_vs_r_vs_z");

	// Convert from Tesla/cm to Gauss/cm
	step_size_vs_r_vs_z->Scale(10000.0);
	
	// Convert to step sizes
	int Nrbins = step_size_vs_r_vs_z->GetYaxis()->GetNbins();
	int Nzbins = step_size_vs_r_vs_z->GetXaxis()->GetNbins();
	for(int ir=1; ir<=Nrbins; ir++){
		for(int iz=1; iz<=Nzbins; iz++){
			double val = step_size_vs_r_vs_z->GetBinContent(iz, ir);
			double step_size = 200.0/val;
			step_size_vs_r_vs_z->SetBinContent(iz, ir, step_size);
		}
	}
	
	vector <map<double, double> > step_maps;
	step_maps.push_back(FindSteps(1.0 * 2.54, step_size_vs_r_vs_z));
	step_maps.push_back(FindSteps(12.0 * 2.54, step_size_vs_r_vs_z));
	step_maps.push_back(FindSteps(22.0 * 2.54, step_size_vs_r_vs_z));
	step_maps.push_back(FindSteps(32.0 * 2.54, step_size_vs_r_vs_z, 355.0));
	
	ofstream ofs("step_map.tex");
	PrintLatexTable(ofs, 0.0, step_maps[0]);
	PrintLatexTable(ofs, 12.0, step_maps[1]);
	PrintLatexTable(ofs, 22.0, step_maps[2]);
	PrintLatexTable(ofs, 32.0, step_maps[3]);
	ofs.close();
	
	PrintDataPoints(0.0, step_maps[0]);
	PrintDataPoints(12.0, step_maps[1]);
	PrintDataPoints(22.0, step_maps[2]);
	PrintDataPoints(32.0, step_maps[3]);	
}

//--------------------
// FindSteps
//--------------------
map<double, double> FindSteps(double r, TH2D *h, double z_max)
{
	// Get bin in r
	int rbin = h->GetYaxis()->FindBin(r);

	// Keep track of points and step sizes
	map<double, double> points;
	
	// Step downstream
	double z_start = 0.0;
	double z = z_start;
	while (z < z_max){
		
		// Get step size for the point we're at
		int zbin = h->GetXaxis()->FindBin(z);
		double step_size1 = GetStandardStepSize(h->GetBinContent(zbin, rbin));
	
		// Get step_size for point we're stepping to
		zbin = h->GetXaxis()->FindBin(z + step_size1);
		double step_size2 = GetStandardStepSize(h->GetBinContent(zbin, rbin));
		double step_size = (step_size1 > step_size2 ? step_size2:step_size1);
		
		points[z] = step_size;
		
		z += step_size;
	}

	return points;
}

//--------------------
// GetStandardStepSize
//--------------------
double GetStandardStepSize(double step_size)
{
	// Convert the floating point value from the map to
	// one of the standard step sizes
	int Nsizes = 5;
	double step_sizes[] = {1.0, 2.0, 4.0, 8.0, 12.0};

	double s = 1.0; // cm (minimum step size)
	for(int i=0; i<Nsizes-1; i++){
		if(step_size > step_sizes[i]){
			s= step_sizes[i+1];
		}else{
			break;
		}
	}

	return s;
}

//--------------------
// PrintLatexTable
//--------------------
void PrintLatexTable(ostream &ofs, double r_inches, map<double,double> &points)
{
	unsigned int Ncols = 6;
	unsigned int max_rows = 34;

	char r_cm[256];
	sprintf(r_cm, "%3.1f", r_inches*2.54);
	
	char ref[256];
	sprintf(ref,"step_tab%d",(int)r_inches);
	
	stringstream caption;
	caption << "Measurement locations in Z for when probe is placed at r="<<r_cm<<" cm ("<<r_inches<<" inches) from the beam line.";
	caption << "Values are in cm and measured downstream along the beamline ";
	caption << "with zero being even with the upstream face of the most upstream coil in the magnet";
	caption << "(see figure \\ref{GlueX_Angular_Tolerance}).";

	ofs << "%============== "<< r_cm <<" cm ("<<r_inches<<"'') =====================================" << endl;
	ofs << "\\begin{table}[htdp]" << endl;
	ofs << "\\caption{"<<caption.str()<<"}" << endl;
	ofs << "\\begin{center}" << endl;
	ofs << "\\begin{tabular}{|";  for(unsigned int i=0; i<Ncols; i++) ofs << "r|"; ofs << "}" << endl;
	ofs << "\\hline" << endl;
	ofs << "\\multicolumn{"<<Ncols<<"}{|c|}{" << r_cm <<" cm ("<<r_inches<<"'')   -   "<<points.size()<<" points} \\\\" <<endl;
	ofs << "\\hline" << endl;

	
	unsigned int Nvals = 0;
	unsigned int Ntables=1;
	map<double,double>::iterator iter=points.begin();
	for(unsigned int i=0; i<points.size(); i++, iter++, Nvals++){
		double z = iter->first;
		double step = iter->second;
		ofs << z << " \\emph{\\footnotesize(+" << step << "cm)}";
		if(i%Ncols == (Ncols-1)){
			ofs << " \\\\" << endl;
			if((i/Ncols + 1)%max_rows == 0){
				// exceeded maximum number of rows. Close out this table and
				// start another since it won't all fit on on page otherwise.
				ofs << "\\hline" << endl;
				ofs << "\\end{tabular}" << endl;
				ofs << "\\end{center}" << endl;
				ofs << "\\label{"<<ref<<"}" << endl;
				ofs << "\\end{table}%" << endl;
				ofs << endl;
				
				if(Ntables%4 == 0)ofs<<"\\clearpage"<<endl;

				ofs << "\\begin{table}[htdp]" << endl;
				ofs << "\\caption{Continuation of table \\ref{"<<ref<<"}.}" << endl;
				ofs << "\\begin{center}" << endl;
				ofs << "\\begin{tabular}{|";  for(unsigned int ii=0; ii<Ncols; ii++) ofs << "r|"; ofs << "}" << endl;
				ofs << "\\hline" << endl;
				ofs << "\\multicolumn{"<<Ncols<<"}{|c|}{" << r_cm <<" cm ("<<r_inches<<"'')   -   "<<points.size()<<" points} \\\\" <<endl;
				ofs << "\\hline" << endl;

				sprintf(ref,"step_tab%d_%d",(int)r_inches, Ntables++);
			}
		}else{
			ofs << " & ";
		}
	}
	
	// Close off last line
	unsigned int i=Nvals%Ncols;
	if(i != (Ncols-1)){
		for(; i<(Ncols-1); i++){
			ofs << " & ";
		}
		ofs << " \\\\" << endl;
	}

	ofs << "\\hline" << endl;
	ofs << "\\end{tabular}" << endl;
	ofs << "\\end{center}" << endl;
	ofs << "\\label{"<<ref<<"}" << endl;
	ofs << "\\end{table}%" << endl;
	ofs << endl;
}

//--------------------
// PrintDataPoints
//--------------------
void PrintDataPoints(double r_inches, map<double,double> &points)
{
	char fname[256];
	sprintf(fname, "step_map_r%02dinches.txt", (int)r_inches);
	
	ofstream ofs(fname);
	time_t t = time(NULL);

	ofs << "#" <<endl;
	ofs << "# Generated by step_map.C" <<endl;
	ofs << "#" <<endl;
	ofs << "# "<<ctime(&t);
	ofs << "#" <<endl;
	ofs << "# user: "<<getenv("USER")<<endl;
	ofs << "# host: "<<getenv("HOST")<<endl;
	ofs << "#" <<endl;
	ofs << "#   z(cm) r(cm)"<<endl;
	ofs << "# =============="<<endl;
	
	
	double r_cm = r_inches*2.54;
	map<double,double>::iterator iter=points.begin();
	for(unsigned int i=0; i<points.size(); i++, iter++){
		double z = iter->first;
		ofs << z << " " << r_cm << endl;
	}
	
	ofs.close();
}