//	BLTrackFile.cc
/*
This source file is part of G4beamline, http://g4beamline.muonsinc.com
Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved.

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

http://www.gnu.org/copyleft/gpl.html
*/

#include <string.h>
#include "BLTrackFile.hh"
#include "BLAsciiFile.hh"

BLTrackFile::BLTrackFile()
{
	status = BLTF_DUMMY;
	mode = "";
	unit = mm;
	file = 0;
}

BLTrackFile::BLTrackFile(G4String filename, G4String comment, G4String _mode)
{
	status = BLTF_ERROR;
	mode = _mode;
	unit = mm;
	if(mode == "w")
		file = BLAsciiFile::fopen(filename);
	else
		file = fopen(filename.c_str(),"r");
	if(!file || (mode != "r" && mode != "w")) {
		fprintf(stderr,"BLTrackFile: cannot open '%s' in mode '%s'",
					filename.c_str(),mode.c_str());
		return;
	}
	status = BLTF_OK;
	if(mode == "w") {
		fprintf(file,"#BLTrackFile %s\n",comment.c_str());
		fprintf(file,"#x y z Px Py Pz t PDGid EventID TrackID ParentID Weight\n");
		if(unit == mm)
		    fprintf(file,"#mm mm mm MeV/c MeV/c MeV/c ns - - - - -\n");
		else
		    fprintf(file,"#cm cm cm MeV/c MeV/c MeV/c ns - - - - -\n");
	} else if(mode == "r") {
		char line[512];
		if(!fgets(line,sizeof(line),file)) {
			status = BLTF_ENDOFFILE;
			return;
		}
		G4String s(line);
		if(s.substr(0,12) != "#BLTrackFile") {
			G4Exception("BLTrackFIle input","Possible wrong format",
				JustWarning, "Assumed OK");
			rewind(file);
		}
	}
}

BLTrackFile::~BLTrackFile()
{
	if(file) BLAsciiFile::fclose(file);
	file = 0;
	status = BLTF_ERROR;
}

BLTrackFileStatus BLTrackFile::write(G4ThreeVector& pos, G4double time,
			G4ThreeVector& momentum, int PDGid, int eventId,
			int trackId, int parentId, G4double weight)
{
	if(status == BLTF_DUMMY) return BLTF_OK;
	if(status != BLTF_OK || mode != "w") return BLTF_ERROR;

	// left justified to avoid pesky initial spaces in column 1
	fprintf(file,"%.6g %.6g %.6g ",pos[0]/unit,pos[1]/unit,pos[2]/unit);
	fprintf(file,"%.6g %.6g %.6g ",momentum[0]/MeV,momentum[1]/MeV,
						momentum[2]/MeV);
	fprintf(file,"%.6g %d %d %d %d %.6g\n",time/ns,PDGid,eventId,
						trackId,parentId,weight);
	return status;
}

BLTrackFileStatus BLTrackFile::read(G4ThreeVector& pos, G4double& time,
			G4ThreeVector& momentum, int& PDGid, int& eventId,
			int& trackId, int& parentId)
{
	G4double weight;	// ignored
	return read(pos,time,momentum,PDGid,eventId,trackId,parentId,weight);
}

BLTrackFileStatus BLTrackFile::read(G4ThreeVector& pos, G4double& time,
			G4ThreeVector& momentum, int& PDGid, int& eventId,
			int& trackId, int& parentId, G4double &weight)
{
	if(status == BLTF_DUMMY) return BLTF_ENDOFFILE;
	if(status != BLTF_OK || mode != "r") return BLTF_ERROR;

	char line[1024];
	do {
		if(!fgets(line,sizeof(line),file)) {
			status = BLTF_ENDOFFILE;
			return status;
		}
		if(line[0] == '#' && strstr(line,"cm cm cm"))
			unit = cm;
	} while(line[0] == '#' || line[0] == '\n' || line[0] == '\r' || 
							line[0] == '\0');

	// Tom Roberts 2010-04-05
	// Read into doubles, and then convert. This permits event numbers 1.0e7
	// and 10000001 to be read as accurately as possible.
	double data[12];

	//Ajit Kurup 2008-06-02
	//Read in PDGid, eventId,trackId and parentId as integers
	//since the float->int conversion doesn't work properly for
	//large numbers

	data[11] = 1.0;	// in case weight is missing from file
	if(sscanf(line,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
				data+0,data+1,data+2,data+3,data+4,data+5,
				data+6,data+7,data+8,data+9,data+10,
				data+11) < 11) {
		printf("BLTrackFile invalid input line (file abandoned): '%s'\n",
				line);
		fflush(stdout);
		status = BLTF_ERROR;
		return status;
	}

	pos[0] = data[0]*unit;
	pos[1] = data[1]*unit;
	pos[2] = data[2]*unit;
	momentum[0] = data[3]*MeV;
	momentum[1] = data[4]*MeV;
	momentum[2] = data[5]*MeV;
	time = data[6]*ns;
	PDGid = (int)data[7];
	eventId = (int)data[8];
	trackId = (int)data[9];
	parentId = (int)data[10];
	weight = data[11];

	return status;
}