//
//
// This macro is used to generate ASCII files containing a parameterized 
// magnetic field map in a format suitable for use with the JANA
// JCalibrationFile class. It takes as input the root files 
// BfieldParameters_Bz.root and BfieldParameters_Bx.root which themselves
// are produced by the ParameterizeBField.C macro using bfield.root
// file as input. (The bfield.root file is generated with the bfield2root
// utility from a full map already in the calibDB).
//

//------------------------------
// ParameterizeBField_codegen
//------------------------------
void ParameterizeBField_codegen(const char *based_on="_1500_poisson_20090814_01")
{
	gROOT->Reset();

	// Create master params file
	stringstream fname;
	fname<<"solenoid"<<based_on<<"_params";
	ofstream fout_master(fname.str().c_str());

	// Write header info to file
	time_t now = time(NULL);
	fout_master<<"#"<<endl;
	fout_master<<"# Solenoidal field parameterization"<<endl;
	fout_master<<"#"<<endl;
	fout_master<<"# This file auto-generated by ParameterizeBField_codegen.C macro."<<endl;
	fout_master<<"# "<<ctime(&now);
	fout_master<<"# Produced by "<<getenv("USER")<<" on "<<getenv("HOST")<<endl;
	fout_master<<"#"<<endl;
	fout_master<<"#      input based on: "<<based_on<<endl;
	fout_master<<"#"<<endl;		
	fout_master<<"# This is the master parameters list. The parameterization of the"<<endl;		
	fout_master<<"# field is broken up into several sections in z. Each section has"<<endl;
	fout_master<<"# its own set of parameters kept in its own namepath."<<endl;
	fout_master<<"#"<<endl;
	fout_master<<"# Each line below represents a section. Sections are either representing"<<endl;
	fout_master<<"# Bx or Bz. Given in the table are the limits over with the section"<<endl;
	fout_master<<"# covers."<<endl;
	fout_master<<"#"<<endl;
	fout_master<<"#"<<endl;
	fout_master<<"#% Bi  sec  zmin  zmax  zmid  znorm  rmin  rmax  rmid  rnorm  namepath"<<endl;

	// Write indivdual params files, adding entries to master file
	WriteFile("BfieldParameters_Bz.root", "Bz", fout_master, fname);
	WriteFile("BfieldParameters_Bx.root", "Bx", fout_master, fname);
	
	// Close master params file
	fout_master.close();
}

//------------------------------
// WriteFile
//------------------------------
void WriteFile(const char *input_fname, string element, ofstream &fout_master, stringstream &fname)
{
#define MAX_SECTIONS 100

	// Open input ROOT file
	TFile *f = new TFile(input_fname);
	
	cout<<"------------------- Ignore errors in this section --------------------"<<endl;

	// Get zmin, zmax, rmin, and rmax for each section
	double zmin[MAX_SECTIONS], zmax[MAX_SECTIONS];
	double rmin[MAX_SECTIONS], rmax[MAX_SECTIONS];
	TH1D *z_bounds_hist = (TH1D*)gROOT->FindObject("z_bounds_hist");
	UInt_t Nsec = z_bounds_hist->GetNbinsX()-1;
	for(UInt_t sec=1; sec<=Nsec; sec++){
		zmin[sec-1] = z_bounds_hist->GetBinContent(z_bounds_hist->FindBin(sec));
		zmax[sec-1] = z_bounds_hist->GetBinContent(z_bounds_hist->FindBin(sec+1));

		stringstream hname;
		hname<<"sec"<<(sec)<<"_p0";	// sec1_p0, sec2_p0, ...
		TH1D *sec_p0 = (TH1D*)gROOT->FindObject(hname.str().c_str());
		if(!sec_p0)break;
		rmin[sec-1] = sec_p0->GetXaxis()->GetXmin();
		rmax[sec-1] = sec_p0->GetXaxis()->GetXmax();
	}

	// Get order of Chebyshev polynomials used for 1st level of parameterization
	UInt_t order1;
	for(order1=0; order1<20; order1++){
		stringstream hname;
		hname<<"sec1"<<"_p"<<order1;
		TH1D *sec1_p = (TH1D*)gROOT->FindObject(hname.str().c_str());
		if(!sec1_p){order1--; break;}
	}

	// Get order of Chebyshev polynomials used for 2nd level of parameterization
	TH1D *sec1_pp0 = (TH1D*)gROOT->FindObject("sec1_pp0");
	UInt_t order2 = sec1_pp0->GetNbinsX()-1;

	cout<<"----------------------------------------------------------------------"<<endl;

	// Information messages
	time_t now = time(NULL);
	cout<<endl;
	cout<<"Writing parameters for: "<<element<<endl;
	cout<<"       sections = "<<Nsec<<endl;
	cout<<"1st level order = "<<order1<<endl;
	cout<<"2nd level order = "<<order2<<endl;
	cout<<endl;

	// Parameters
	for(UInt_t sec=1; sec<=Nsec; sec++){
		
		// Open output file
		stringstream my_fname;
		my_fname<<fname.str()<<"_"<<element<<"_sec"<<sec;
		ofstream fout(my_fname.str().c_str());
		cout<<"Writing parameters to "<<fname.str()<<endl;
		
		// Add line to master file
		fout_master<<element;
		fout_master<<"\t"<<sec;
		fout_master<<"\t"<<zmin[sec-1];
		fout_master<<"\t"<<zmax[sec-1];
		fout_master<<"\t"<<(zmax[sec-1]+zmin[sec-1])/2.0;
		fout_master<<"\t"<<(zmax[sec-1]-zmin[sec-1])/2.0;
		fout_master<<"\t"<<rmin[sec-1];
		fout_master<<"\t"<<rmax[sec-1];
		fout_master<<"\t"<<(rmax[sec-1]+rmin[sec-1])/2.0;
		fout_master<<"\t"<<(rmax[sec-1]-rmin[sec-1])/2.0;
		fout_master<<"\t"<<"Magnets/Solenoid/"<<my_fname.str();
		fout_master<<endl;
		
		// Write header info to file
		fout<<"#"<<endl;
		fout<<"# Solenoidal field parameterization"<<endl;
		fout<<"#"<<endl;
		fout<<"# This file auto-generated by ParameterizeBField_codegen.C macro."<<endl;
		fout<<"# "<<ctime(&now);
		fout<<"# Produced by "<<getenv("USER")<<" on "<<getenv("HOST")<<endl;
		fout<<"#"<<endl;
		fout<<"#      input file: "<<f->GetName()<<endl;
		fout<<"#         section: "<<sec<<"  (out of "<<Nsec<<")"<<endl;
		fout<<"# 1st level order: "<<order1<<endl;
		fout<<"# 2nd level order: "<<order2<<endl;
		fout<<"#"<<endl;		
		fout<<"# In the following table, values are coefficients of Chebyshev"<<endl;		
		fout<<"# polynomials. There are 2 levels of parameterization corresponding"<<endl;		
		fout<<"# to the 2 dimensions in which the field map is parameterized (z and r)."<<endl;		
		fout<<"#"<<endl;
		fout<<"# The first level fits the field as a function of z in lab coordinates."<<endl;
		fout<<"# The second level parametrizes those coefficients as a function of r"<<endl;
		fout<<"# also in lab coordinates."<<endl;
		fout<<"#"<<endl;
		fout<<"#"<<endl;
		fout<<"#"<<endl;

		int chars_per_val = 14;
		string header((order2+1)*chars_per_val, ' ');
		for(UInt_t j=0; j<=order2; j++){
			stringstream name;
			name<<"p"<<j;
			string &s = name.str();
			header.replace((j+1)*chars_per_val-s.length()-2, s.length(), s);
		}
		fout<<"#"<<header<<endl;

		// Loop
		for(UInt_t i=0; i<=order1; i++){
			stringstream hname;
			hname<<"sec"<<sec<<"_pp"<<i;
			TH1D *h = (TH1D*)gROOT->FindObject(hname.str().c_str());
			
			string line((order2+1)*chars_per_val, ' ');
			for(UInt_t j=0; j<=order2; j++){
				stringstream val;
				val<<h->GetBinContent(j);
				string &s = val.str();
				if(s.length()>=chars_per_val){
					cerr<<"ERROR: too many characters in val! ("<<__FILE__<<":"<<__LINE__<<")"<<endl;
					return -1;
				}
				line.replace((j+1)*chars_per_val-s.length(), s.length(), s);
			}
			fout<<line<<endl;
		}
		
		fout.close();
	}
}