#include <iostream>
#include <iomanip>
#include <fstream>
#include <exception>
#include <string>
#include <sstream>
#include <set>
#include <vector>
#include <iterator>
#include <stdio.h>
#include <sys/time.h>
#include <time.h>
using namespace std;

#define debug 0


#include "GPUInitWires.cu"
#include "GPUFitTrack.cu"
#include "GPUMagneticFieldStepper.cu"
#include "GPUResidualFinder.cu"
#include "GPUErrorCheck.h"


#define MINWIREHITS 6 //Minimum wire hits needed per candidate

//----------------
// main
//----------------
int main(int narg, char *argv[])
{
	int succeed = 0, fail = 0;
	
	InitCUDA();
	GPUInitWires();

	// Create list of mass hypotheses and setup GPU for tracking
	vector<float> mass_hypotheses;
	mass_hypotheses.push_back(0.13965);
	//mass_hypotheses.push_back(0.542);
	//mass_hypotheses.push_back(0.938);
	GPUInitTracking(mass_hypotheses);
	
	// Open output file for fit results
	string resultsfname="final_state_params.txt";
	ofstream resultsfile(resultsfname.c_str());
	
	// Open candidate definitions file
 	string fname="gluex_candidates.txt";
	ifstream file(fname.c_str());
	if(!file.is_open()){
		cout << "Unable to open file \"" << fname << "\"!" << endl;
		exit(-1);
	}
	cout << "Reading " << fname << " ..." << endl;
	
	int c1 = 0;
	int c2 = 0;
	int c3 = 0;
	int c4 = 0;
	int c5 = 0;
	
struct timeval start;
gettimeofday(&start, NULL);

	time_t start_t = time(NULL);
	
	// Read in candidate definitions
	int lineno=0;
	int Nevents = 0;
	while(!file.eof()){
		char line[50000]; // Increased from 256
		file.getline(line, 50000);
		lineno++;  // keep track of line number in case we need to report error
		if(line[0] == '#') continue;    // skip comment lines
		
		// Read in candidate info
		stringstream ss(line);
		vector<float> values; // Store values for line
		while (!ss.fail()) {
			float f;
			ss >> f;
			if (!ss.fail())
				values.push_back(f);
		}
		if(values.empty()){
		 _DBG_ << "empty line in candidates file. Stopping." << endl;
		 continue; // Avoids errors on last line
		}
			if(values.size() < (7 + MINWIREHITS))
			{
				cout << "Error. Candidate values not present on line " << lineno << "." << endl;
				exit(-1);
			}
			
			cout << "Event: " << ++Nevents << "  \r"; cout.flush();

			
			// Extract candidate info
			float q, x, y, z, px, py, pz;
      	  //ss >> q >> x >> y >> z >> px >> py >> pz;
			q = values[0], x = values[1], y = values[2], z = values[3], px = values[4], py = values[5], pz = values[6];
			float3 pos = make_float3(x, y, z);
			float3 mom = make_float3(px, py, pz);
			float p = length(mom);
			float q_over_p = q/p;
			mom *= (1.0/p);

			vector<int> wires_hit; // indices of wires hit on track
			for(int i = 7; i < values.size(); i++) // all other values besides 0-5 should be wires
			{
				float f = values[i];
				wires_hit.push_back(f);
			}

			if(wires_hit.size() > MAX_WIRES_HIT)
			{
				cout << "FATAL: The number of wires found in the wire definitions file (" << wires_hit.size() << ")" << endl;
				cout << "       is greater than number allocated at compile time (" <<  MAX_WIRES_HIT << ")." << endl;
				cout << "       You must edit GPUMagneticFieldStepper.cu and recompile in order to use" <<endl;
				cout << "       the wire definitions in the file \"" << fname << "\"." << endl;
				
				resultsfile << "0 0 0 0 0 0 0 0" << endl; // placeholder indicating fit failed
				
				continue;
			}

#if debug			
cout << "Wires found " << wires_hit.size() << endl;
cout << line << endl;	
#endif	

			vector<fit_result_t> results;
			
			
			int OK = GPUFitTrack(pos, mom, q_over_p, wires_hit, results);
			if(OK == 0){
				succeed++;
				float mass = 0.1395; // temporary
				resultsfile << q_over_p << " "; // q_over_p
				resultsfile << pos.x << " " << pos.y << " " << pos.z << " "; // pos: x y z
				resultsfile << mom.x << " " << mom.y << " " << mom.z << " "; // mom: x y z
				resultsfile << mass << endl; // mass
			}else{
				fail++;
				if(OK == 1) c1++;
				if(OK == 2) c2++;
				if(OK == 3) c3++;
				if(OK == 4) c4++;
				if(OK == 5) c5++;
				resultsfile << "0 0 0 0 0 0 0 0" << endl; // placeholder indicating fit failed
			}
#if debug			
if(Nevents > 0) break;
#endif	
	}
	
struct timeval end;
gettimeofday(&end, NULL);
//printf("%ld\n", ((end.tv_sec*1000000 + end.tv_usec) - (start.tv_sec*1000000 + start.tv_usec)));
long newT = ((end.tv_sec*1000000 + end.tv_usec) - (start.tv_sec*1000000 + start.tv_usec));
cout << "Precision timer: " << (float)((float)Nevents/(float)newT*1000000) << endl;	
	
	time_t end_t = time(NULL);
	double events_per_second = (double)Nevents/(double)(end_t - start_t);
	cout << endl;
	cout << Nevents << " processed at " <<events_per_second << " events/sec" <<endl;
	cout << succeed << " events succeeded, " << fail << " events failed ";
	cout << (float)((float)fail/(float)(succeed+fail))*100 << " percent failed " << endl;
	
	cout << "c1 " << c1 << " c2 " << c2 << " c3 " << c3 << " c4 " << c4 << " c5 " << c5 << endl;
	file.close();
	resultsfile.close();
	cout << "No errors observed." << endl;
 	return 0;
}