#include <string>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
using namespace std;

#include <TFile.h>
#include <TTree.h>
#include <TH3D.h>
#include <TH2D.h>

#define _DBG_ cerr<<__FILE__<<":"<<__LINE__<<" "
#define _DGB__ _DBG_<<endl

// Globals
string FILENAME;
string PLANE="000";
int ANGLE=0;
int YEAR, MONTH, DAY, HOUR, MINUTE, SECONDS;
double SCAN_SEP = 2.54; // scan separation distance in cm (0.0254 m = 1 inch) (value of 2.57 copied from Michah's macro)
double SCAN_DIST = 118.0;  // total scan distance (cm)
double DATAPOINT_SEP = SCAN_DIST/867.0;  // approximate distance between measurments (needs to be replaced with a calculated value)
int MAX_X_STEPS = 867;

// Function declarations
void MakeMeanAndRMSHistos(TTree *scan);
void ParseCommandLineArgs(int narg, char *argv[]);
void Usage(void);

//---------------
// main
//---------------
int main(int narg, char *argv[])
{

	// Parse the command line arguments
	ParseCommandLineArgs(narg, argv);
	
	// Open the input file
	ifstream ifile(FILENAME.c_str());
	if(!ifile.is_open()){
		_DBG_<<"Unable to open file \""<<FILENAME<<"\"!"<<endl;
		return -5;
	}
	
	// Open ROOT file
	char rfilename[256];
	sprintf(rfilename, "plane%s %03ddeg %04d_%02d_%02d %02d_%02d_%02d.root", PLANE.c_str(), ANGLE, YEAR, MONTH, DAY, HOUR, MINUTE, SECONDS);
	TFile *rfile = new TFile(rfilename, "RECREATE", "Made by sname_mktree");
	if(!rfile->IsOpen()){
		_DBG_<<"Unable to create ROOT file \""<<rfilename<<"\"!"<<endl;
		return -6;
	}
	
	// Create TTrees
	double vals[10];
	int ivals[10];
	TTree *scan = new TTree("scan", "flatness data from ascii file");
	scan->Branch("S", vals, "x/D:y:depth");
	scan->Branch("P", ivals, "pass/I:n");
	
	// Inform user what's going on
	cout<<endl;
	cout<<"----- Converting FDC snake scan data to ROOT file ----"<<endl;
	cout<<" Input file: "<<FILENAME<<endl;
	cout<<"Output file: "<<rfilename<<endl;
	cout<<endl;
	cout<<"      PLANE: "<<PLANE<<endl;
	cout<<"      ANGLE: "<<ANGLE<<endl;
	cout<<"       YEAR: "<<YEAR<<endl;
	cout<<"      MONTH: "<<MONTH<<endl;
	cout<<"        DAY: "<<DAY<<endl;
	cout<<"       HOUR: "<<HOUR<<endl;
	cout<<"     MINUTE: "<<MINUTE<<endl;
	cout<<"    SECONDS: "<<SECONDS<<endl;
	cout<<endl;
	cout<<" Scanning input ..."<<endl;

	// Loop over lines in file
	int pass = 0;
	int Nmeasurements_this_pass=0;
	while(!ifile.eof()){
	
		// Read in first token as string. This usually is a number, but between
		// passes, it will be the string "horizontal" from a line that says 
		// "horizontal pass X complete".
		string str;
		ifile>>str;
		
		// If the string is "horizontal" then process the end of pass
		if(str=="horizontal"){
			ifile>>str; // read in "pass" string;
			ifile>>pass;
			ifile>>str; // read in "complete" string
			cout<<"finished pass "<<pass<<" ("<<Nmeasurements_this_pass<<" measurements)           \r";
			Nmeasurements_this_pass = 0;
			continue;
		}
		
		// Must be a depth value. Add to tree.
		double x = DATAPOINT_SEP*(double)Nmeasurements_this_pass;
		double y = SCAN_SEP*(double)pass;
		
		// Every other pass is in the backwards direction
		// The 8cm shift is empirical taken from eyeballing the depth-x plot of the
		// first 2 passes in "plane999 999deg 2008_11_06 08_51_09.txt"
		if((pass%2)==0)x = (double)MAX_X_STEPS*DATAPOINT_SEP - x - 8.0;
		
		// Linear correction for slope of laser support structure
		double depth = atof(str.c_str());
		float slope = 500.0E-3/97.0; // slope is 500 microns per 97 cm
		float x_center = 56.0; // approximate
		depth += slope*(x-56.0);		
		
		// Fill tree
		vals[0] = x;
		vals[1] = y;
		vals[2] = depth;
		ivals[0] = pass;
		ivals[1] = Nmeasurements_this_pass;
		scan->Fill();
		
		Nmeasurements_this_pass++;
	}
	cout<<endl;
	cout<<endl;
	
	// Make the mean and RMS histo grams
	MakeMeanAndRMSHistos(scan);

	// Close ROOT file
	rfile->Write();
	delete rfile;

	cout<<"Complete"<<endl;

	return 0;
}

//----------------
// MakeMeanAndRMSHistos
//----------------
void MakeMeanAndRMSHistos(TTree *scan)
{
	cout<<"Generating mean and rms histograms ..."<<endl;

	// Create a histograms to project data onto
	TH3D *depth = new TH3D("depth","", 120/3, -10.0, 110.0, 44, -10.0, 102.0, 2000, 0.0, 8.0);
	TH2D *depth2D = new TH2D("depth2D","", 120/3, -10.0, 110.0, 44, -10.0, 102.0);
	scan->Project("depth", "depth:y:x");
	
	// Create temporary histo to hold z-projections
	TAxis *zaxis = depth->GetZaxis();
	TH1D *h = new TH1D("h","", zaxis->GetNbins(), zaxis->GetXmin(), zaxis->GetXmax());
	
	// Loop over all bins and find mean and RMS for each
	TH2D *depth_mean = (TH2D*)depth2D->Clone("depth_mean");
	TH2D *depth_rms = (TH2D*)depth2D->Clone("depth_rms");
	for(int ibin=1; ibin<=depth->GetXaxis()->GetNbins(); ibin++){
		for(int jbin=1; jbin<=depth->GetYaxis()->GetNbins(); jbin++){
			
			h->Reset();
			for(int kbin=1; kbin<=h->GetNbinsX(); kbin++){
				h->SetBinContent(kbin, depth->GetBinContent(ibin, jbin, kbin));
			}

			double mean = h->GetMean();
			double rms = h->GetRMS();
			
			depth_mean->SetBinContent(ibin, jbin, mean);
			depth_rms->SetBinContent(ibin, jbin, 1000.0*rms);
		}
	}
	
	delete h;

	// Set labels
	depth_mean->SetStats(0);
	depth_mean->SetXTitle("x (cm)");
	depth_mean->SetYTitle("y (cm)");
	depth_mean->SetTitle("Mean of depth measurements per ~1 in^{2}");

	depth_rms->SetStats(0);
	depth_rms->SetXTitle("x (cm)");
	depth_rms->SetYTitle("y (cm)");
	depth_rms->SetTitle("RMS of depth measurements per ~1 in^{2}");
}

//---------------
// ParseCommandLineArgs
//---------------
void ParseCommandLineArgs(int narg, char *argv[])
{
	// Initialize FILENAME global
	FILENAME = "";

	// Loop over arguments
	for(int i=1; i<narg; i++){
		string arg(argv[i]);
		string next((i+1)<narg ? argv[i+1]:"");
		bool used_next = false;
		
		if(arg=="-h" || arg=="--help"){
			Usage();
			exit(0);
		}else{
			if(FILENAME!=""){
				Usage();
				_DBG_<<"More than one filename defined! (first=\""<<FILENAME<<"\", second=\""<<arg<<"\")"<<endl;
				_DBG_<<"Exiting ..."<<endl;
				exit(-2);
			}else{
				FILENAME=arg;
			}
		}
		
		if(used_next && i==(narg-1)){
			Usage();
			_DBG_<<"Value needed for \""<<arg<<"\"!"<<endl;
			exit(-1); 
		}
	}

	// Check that a file was specified
	if(FILENAME==""){
		Usage();
		_DBG_<<"No file specified!"<<endl;
		_DBG_<<"You MUST specify a file to convert!"<<endl;
		exit(-3);
	}
	
	// Parse filename to get info about it
	size_t pos = FILENAME.find("plane");
	if(pos==string::npos){_DBG_<<"Unable to parse filename \""<<FILENAME<<"\"!"<<endl; exit(-4);}
	
	stringstream ss_PLANE (FILENAME.substr(pos+5, 3));
	stringstream ss_ANGLE (FILENAME.substr(pos+5+4, 3));
	stringstream ss_YEAR  (FILENAME.substr(pos+5+4+7, 4));
	stringstream ss_MONTH (FILENAME.substr(pos+5+4+7+5, 2));
	stringstream ss_DAY   (FILENAME.substr(pos+5+4+7+5+3, 2));
	stringstream ss_HOUR  (FILENAME.substr(pos+5+4+7+5+3+3, 2));
	stringstream ss_MINUTE(FILENAME.substr(pos+5+4+7+5+3+3+3, 2));
	stringstream ss_SECONDS(FILENAME.substr(pos+5+4+7+5+3+3+3+3, 2));
	ss_PLANE>>PLANE;
	ss_ANGLE>>ANGLE;
	ss_YEAR>>YEAR;
	ss_MONTH>>MONTH;
	ss_DAY>>DAY;
	ss_HOUR>>HOUR;
	ss_MINUTE>>MINUTE;
	ss_SECONDS>>SECONDS;
}

//---------------
// Usage
//---------------
void Usage(void)
{
	cout<<endl;
	cout<<"Usage:"<<endl;
	cout<<"    snake_mktree [-h|--help] filename.txt"<<endl;
	cout<<endl;
	cout<<" Generate a ROOT file with a TTree containing the data read from"<<endl;
	cout<<"an FDC cathode scan data file. The program should be run with a "<<endl;
	cout<<"single argument indicating the name of the source data file. The"<<endl;
	cout<<"input file can exist anywhere, but the output file will be placed"<<endl;
	cout<<"in the current working directory. Note that the name of the input"<<endl;
	cout<<"file is important. It should be of the form:"<<endl;
	cout<<endl;
	cout<<"    planeZZZ AAAdeg YYYY_MM_DD hh_mm_ss.txt"<<endl;
	cout<<endl;
	cout<<" where:"<<endl;
	cout<<"      ZZZ is the 3 digit plane number (leading zeros if neccessary)"<<endl;
	cout<<"      AAA is the 3 digit angle of the plane in degrees (rounded)"<<endl;
	cout<<"      YYYY, MM, DD are the year month and date of the scan"<<endl;
	cout<<"      hh, mm, ss are the hour minute and second of the scan"<<endl; 
	cout<<endl;
	cout<<" Note that this program is run automatically after every scan so"<<endl;
	cout<<"usually, running it by hand is not needed."<<endl;
	cout<<endl;
	cout<<"Please refer questions/comments to David Lawrence:"<<endl;
	cout<<"      davidl@jlab.org     x5567"<<endl;
	cout<<endl;
}