//
// CPPsim
//
// GEANT4 based simulation of Charged Pion Polarizability 
// experiment in Hall-D

#include <unistd.h>
#include <signal.h>
#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
#include <thread>
using namespace std;

#include "G4RunManager.hh"
#include "G4MTRunManager.hh"
#include "G4UImanager.hh"

#include "G4HadronicProcessStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4TransportationManager.hh"

#include "QGSP_FTFP_BERT.hh"

#include "CPPUserActionInitialization.hh"
#include "CPPDetectorConstruction.hh"

//#undef G4VIS_USE

#ifdef G4VIS_USE
#include "G4VisExecutive.hh"
#endif

#include "G4UIExecutive.hh"

#include <DANA/DApplication.h>
#include <particleType.h>


// Global variables
DApplication *dapp = NULL;
string gGDMLfilename = "";
unsigned int MAXEVENTS = -1;
string CONFIG_FILE = "control.in";
string INPUT_FILENAME = "";
string OUTPUT_FILENAME = "CPPsim.hddm";
string BFIELD_MAP = "default";
vector<double> KINE_VALS;
vector<double> SCAP_VALS;
vector<double> BEAM_VALS;
double BGRATE;
vector<double> BGGATE_VALS;
G4ThreeVector  BEAM_START(0.0, 0.0, -1000.0);
vector<string> GEANT4_MACROS;
bool INTERACTIVE_MODE = false;
bool BATCH_MODE = false;
bool NO_AUTO_VIS_MAC = false;
bool MODE_MULS = true;
bool MODE_DCAY = true;
bool MODE_HADR = true;
int MODE_ELOSS = 1;
int NTHREADS = 1;
int RNDM = 123;
int RUN_NUMBER = 2;
int HDDM_USE_COMPRESSION = 2; // 0=none, 1=bz2,  2=z
bool HDDM_USE_INTEGRITY_CHECKS = true;
bool GEOMOPT = true;
double FMWPC_DEAD_RADIUS = 0.0;
double FMWPC_DEAD_WIDTH  = 12.0;

// Cobrems values come in through BEAM_VALS. defaults are in GlueXPrimaryGeneratorAction.cc
bool   ONLY_COHERENT = false;
bool   ONLY_INCOHERENT = false;


// Function declarations
void Usage(bool exit_when_done=false);
void ParseCommandLineArguments(int argc,char **argv);
void AdjustPhysicsProcesses(G4UImanager* UImanager);
void ReadControl_in(void);
void PrintConfigInfo(void);

//..............................
// SIGINThandler
// Signal handler to allow us to ctl-C out of a run cleanly
//..............................
void SIGINThandler(int x)
{
	G4cerr << endl << "SIGINT received. Aborting run ..." << endl;
	G4RunManager::GetRunManager()->AbortRun();
}

// COMMENT OUT THE FOLLOWING LINE TO COMPILE WITH MULTI-THREAD SUPPORT
//#undef G4MULTITHREADED

//--------------------
// main
//--------------------
int main(int argc,char **argv)
{
	// Parse all command line arguments
	ParseCommandLineArguments(argc, argv);

	// Create DApplication so we can get at the calibration database and magnetic field maps
	dapp = new DApplication(argc, argv);

	// Read the control.in file
	ReadControl_in();
	
	// Confirm existence of all specified macros
	bool all_macros_exist = true;
	for(unsigned int i=0; i<GEANT4_MACROS.size(); i++){
		ifstream ifs(GEANT4_MACROS[i].c_str());
		if(!ifs.is_open()){
			if(all_macros_exist) Usage(); // print usage statement on first bad macro
			G4cout << "XXXX    Macro file \"" << GEANT4_MACROS[i] << "\" does not exist !   XXXX" << endl;
			all_macros_exist = false;
		}
	}
	if(!all_macros_exist) exit(-1);

	G4UIExecutive* ui = new G4UIExecutive(argc, argv);
	
	// Create G4RunManager and use it to initialize GEANT4
#ifdef G4MULTITHREADED
	G4MTRunManager* runManager = new G4MTRunManager;
	runManager->SetNumberOfThreads(NTHREADS);
#else
	G4RunManager* runManager = new G4RunManager;
#endif
	runManager->SetGeometryToBeOptimized(GEOMOPT);
	runManager->SetUserInitialization(new CPPDetectorConstruction(gGDMLfilename));
	runManager->SetUserInitialization(new QGSP_FTFP_BERT(0));
	runManager->SetUserInitialization(new CPPUserActionInitialization);

#ifdef G4VIS_USE
	G4VisManager* visManager = new G4VisExecutive;
	visManager->Initialize();
	if(!NO_AUTO_VIS_MAC){
		ifstream ifs("vis.mac");
		if(ifs.is_open()){
			ifs.close();
			if(!BATCH_MODE) GEANT4_MACROS.insert(GEANT4_MACROS.begin(), "vis.mac");
		}
	}
#endif

	// Suppress printing the hadronic interactions list.
	// (this trick obtained in response to bug report 1759
	// http://bugzilla-geant4.kek.jp/show_bug.cgi?id=1759)
	G4HadronicProcessStore::Instance()->SetVerbose(0);

	G4UImanager* UImanager = G4UImanager::GetUIpointer();

	runManager->Initialize();
	
	// Set sighandler to catch sigints. Do this after the
	// magnetic field map is read in since that uses the
	// HDGeometry library and installs JANA sig handlers
	signal(SIGINT,SIGINThandler);

	// Set up user interface and graphical user interface

	// Adjust physics processes
	AdjustPhysicsProcesses(UImanager);
	
	// Print selected config. info
	PrintConfigInfo();
	
	// Run all macros
	for(unsigned int i=0; i<GEANT4_MACROS.size(); i++){
		string &macro = GEANT4_MACROS[i];
		G4cout << endl << "---- Running macro \"" << macro << "\" ... ----" << endl;
		UImanager->ApplyCommand(("/control/execute " + macro).c_str());
	}
	
	// Optionally start interactive mode
	if(INTERACTIVE_MODE){
		G4cout << "Starting interactive session ..." << endl;
		ui->SessionStart();
	}else if(MAXEVENTS>0){
		char cmd[256];
		sprintf(cmd,"/run/beamOn %d", MAXEVENTS);
		G4cout << "===== Starting event processing =====" << endl;
		G4cout << cmd << endl;
		UImanager->ApplyCommand(cmd);
	}

	// Tell user were done	
	G4cout << endl << "------- Done -------" << endl;

	delete ui;

#ifdef G4VIS_USE
	delete visManager;
#endif
	
	delete runManager;
	
	return 0;
}

//--------------------
// Usage
//--------------------
void Usage(bool exit_when_done)
{
	G4cout << endl;
	G4cout << " Charged Pion Polarizability in Hall-D simulation." << endl;
	G4cout << endl;
	G4cout << "Usage:" << endl;
	G4cout << "     CPPsim [options]" << endl;
	G4cout << endl;
	G4cout << "  options:" << endl;
	G4cout << "     -h, --help          print help" << endl;
	G4cout << "     -c, --config file   set config file name (def. control.in)" << endl;
	G4cout << "     -i, --interactive   start in interactive mode" << endl;
	G4cout << "     -b, --batch         execute in batch mode (see below)" << endl;	
	G4cout << "     -m, --macro macro   execute the specified macro file" << endl;
	G4cout << "     -nv, --novis        do not automatically execute vis.mac" << endl;
	G4cout << endl;
	G4cout << "This will automatically look for the file \"control.in\" in the" << endl;
	G4cout << "current working directory for configuration options. This file" << endl;
	G4cout << "is of the same format as is used for GlueX's hdgeant. It implements" << endl;
	G4cout << "a few additional flags and does not implement others. See the example" << endl;
	G4cout << "distributed with the CPPsim source code for full details." << endl;
	G4cout << endl;
	G4cout << "If a file named \"vis.mac\" is found in the current working directory" << endl;
	G4cout << "then it is executed automatically unless batch mode is specified" << endl;
	G4cout << "or the -nv or --novis option is given." << endl;
	G4cout << "Batch mode and interactive mode are mutually exclusive and the program" << endl;
	G4cout << "will quit with an error if both are specified." << endl;
	G4cout << endl;
	G4cout << "The -m option can be used to specify a macro to be automatically executed." << endl;
	G4cout << "This option can be specified more than once and all macros will be executed" << endl;
	G4cout << "in the order given on the command line. In addition, a macro can be specified" << endl;
	G4cout << "in the control.in file. If one or more are specified there, then they will" << endl;
	G4cout << "be executed after any macros given on the command line." << endl;
	G4cout << endl;
	
	if(exit_when_done) exit(0);
}

//--------------------
// ParseCommandLineArguments
//--------------------
void ParseCommandLineArguments(int argc,char **argv)
{
	for(int i=1; i<argc; i++){
		string arg(argv[i]);
		string next = (i+1)<argc ? argv[i+1]:"";
		bool next_is_switch = true;
		if(next.size()>0) next_is_switch = (next[0]=='-');
		
		bool switch_missing_arg = false;
		if(arg=="-h" || arg=="--help"       ) Usage(true);
		if(arg=="-i" || arg=="--interactive") INTERACTIVE_MODE = true;
		if(arg=="-b" || arg=="--batch"      ) BATCH_MODE = true;
		if(arg=="-nv"|| arg=="--novis"      ) NO_AUTO_VIS_MAC = true;
		if(arg=="-m" || arg=="--macro"      ){
			if(next_is_switch || next.empty()){
				switch_missing_arg = true;
			}else{
				GEANT4_MACROS.push_back(next);
			}
		}
		if(arg=="-c" || arg=="--config"      ){
			if(next_is_switch || next.empty()){
				switch_missing_arg = true;
			}else{
				CONFIG_FILE = next;
			}
		}
		
		if(switch_missing_arg){
			G4cout << "Argument \"" << arg << "\" requires an argument!" << endl;
			G4cout << "(run with --help for help)" << endl;
			exit(-1);
		}
	}
	
	if(BATCH_MODE && INTERACTIVE_MODE){
		Usage();
		G4cout << "XXXX   Both BATCH and INTERACTIVE modes cannot be specified!   XXXX" << endl << endl;
		exit(-1);
	}

}

//--------------------
// AdjustPhysicsProcesses
//--------------------
void AdjustPhysicsProcesses(G4UImanager* UImanager)
{
	/// This adjusts which physics processes are turn on or off and
	/// and configurations of them based on user values found in
	/// control.in. The relevant settings are stored in global variables.

	G4cout << endl;

	// Multiple Scattering
	if(!MODE_MULS){
		string cmd = "/process/inactivate msc all";
		G4cout << "Turning off multiple scattering (MULS)" << endl;
		G4cout << cmd << endl << endl;
		UImanager->ApplyCommand(cmd.c_str());
	}
	
	// Energy Loss
	vector<string> cmds;
	switch(MODE_ELOSS){
		case 0:
			G4cout << "Turning off energy loss (LOSS=0)" << endl;
			cmds.push_back("/process/inactivate eIoni");
			cmds.push_back("/process/inactivate muIoni");
			cmds.push_back("/process/inactivate ionIoni");
			cmds.push_back("/process/inactivate hIoni");
			break;
		case 1:
		case 2:
			G4cout << "Default energy loss (LOSS=1,2)" << endl;
			break;
		case 4:
			G4cout << "Turning off energy loss fluctuations (LOSS=4)" << endl;
			cmds.push_back("/process/eLoss/fluct false");
			break;
		default:
			G4cout << "Unsuppoerted LOSS value ("<<MODE_ELOSS<<") shuould be 0,1,2, or 4" << endl;
			exit(-1);
	}
	for(uint32_t i=0; i<cmds.size(); i++){
		G4cout << cmds[i] << endl;
		UImanager->ApplyCommand(cmds[i].c_str());
	}
	G4cout << endl;
	
	// Decay
	if(!MODE_DCAY){
		string cmd = "/process/inactivate Decay all";
		G4cout << "Turning off decays (DCAY)" << endl;
		G4cout << cmd << endl << endl;
		UImanager->ApplyCommand(cmd.c_str());
	}
	
	// HADR
	cmds.clear();
	if(!MODE_HADR){
		// This is a subset of the possible reactions. A more complete review is needed.
		G4cout << "Tunining off (selected) Hadronic Interactions" << endl;
		cmds.push_back("/process/inactivate hadElastic");
		cmds.push_back("/process/inactivate neutronInelastic");
		cmds.push_back("/process/inactivate protonInelastic");
		cmds.push_back("/process/inactivate pi+Inelastic");
		cmds.push_back("/process/inactivate pi-Inelastic");
		cmds.push_back("/process/inactivate kaon+Inelastic");
		cmds.push_back("/process/inactivate kaon-Inelastic");
		cmds.push_back("/process/inactivate dInelastic");
		cmds.push_back("/process/inactivate ionInelastic");
	}
	for(uint32_t i=0; i<cmds.size(); i++){
		G4cout << cmds[i] << endl;
		UImanager->ApplyCommand(cmds[i].c_str());
	}
	
	// Random number sequence seeds
	char ccmd[256];
	sprintf(ccmd, "/random/setSequence %d", RNDM);
	G4cout << ccmd << endl;
	UImanager->ApplyCommand(ccmd);
	
	G4cout << endl;

}

//--------------------
// ReadControl_in
//--------------------
void ReadControl_in(void)
{
	ifstream ifs(CONFIG_FILE);
	if(!ifs.is_open()){
		G4cout << endl << "No \"control.in\" file found!" << endl << endl;
		exit(-1);
	}
	G4cout << "Reading configuration from " << CONFIG_FILE << " ... " << endl;
	
	// Loop over file
	int iline=0;
	while(ifs.good()){

		// Read in line from file
		char line[512];
		ifs.getline(line, 512);
		iline++;
		if(ifs.gcount()<1) continue;

		// Tokenize line
		stringstream ss(line);
		vector<string> tokens;
		while(ss.good()){
			string s;
			ss >> s;
			if(!s.empty()) tokens.push_back(s);
		}

		// Skip empty or comment lines
		if(tokens.size() <2)continue; // we need the flag and at least one value for line to be useful
		string flag = tokens[0];
		for(unsigned int i=0; i<flag.size(); i++) flag[i] = toupper(flag[i]);
		if(flag[0] == 'C') continue;
		if(flag[0] == '*') continue;
		
		// Copy values into vals vectors, converting to appropriate forms
		vector<string> vals;
		vector<int> ivals;
		vector<double> fvals;
		for(unsigned int i=1; i<tokens.size(); i++){

			// Remove any single quotes from strings
			string val = "";
			for(unsigned int j=0; j<tokens[i].size(); j++){
				if(tokens[i][j] != '\'') val += tokens[i][j];
			}
			vals.push_back(val);
			
			// Convert ints and floats
			ivals.push_back(atoi(tokens[i].c_str()));
			fvals.push_back(atof(tokens[i].c_str()));
		}
		
		// Look for lines of interest
		if(flag == "TRIG"            ) MAXEVENTS = ivals[0];
		if(flag == "INFILE"          ) INPUT_FILENAME = vals[0];
		if(flag == "OUTFILE"         ) OUTPUT_FILENAME = vals[0];
		if(flag == "GEOMFILE"        ) gGDMLfilename = vals[0];
		if(flag == "GEOMOPT"         ) GEOMOPT = ivals[0];
		if(flag == "GEANT4_MACRO"    ) GEANT4_MACROS.push_back(vals[0]);
		if(flag == "BFIELDMAP"       ) BFIELD_MAP = vals[0];
		if(flag == "RNDM"            ) RNDM = ivals[0];
		if(flag == "MULS"            ) MODE_MULS = ivals[0];
		if(flag == "LOSS"            ) MODE_ELOSS = ivals[0];
		if(flag == "DCAY"            ) MODE_DCAY = ivals[0];
		if(flag == "HADR"            ) MODE_HADR = ivals[0];
		if(flag == "NTHREADS"        ) NTHREADS = vals[0]=="max" ? 0:ivals[0];
		if(flag == "RUNG"            ) RUN_NUMBER = ivals[0];
		if(flag == "BGRATE"          ) BGRATE = fvals[0];
		if(flag == "DEAD_RADIUS"     ) FMWPC_DEAD_RADIUS = fvals[0];
		if(flag == "DEAD_WIDTH"      ) FMWPC_DEAD_WIDTH = fvals[0];
		if(flag == "HDDM_USE_COMPRESSION"      ) HDDM_USE_COMPRESSION = ivals[0];
		if(flag == "HDDM_USE_INTEGRITY_CHECKS" ) HDDM_USE_INTEGRITY_CHECKS = ivals[0];
		if(flag == "KINE"){
			if(fvals.size() != 7){
				G4cout << "KINE flag requires 7 arguments!!" << endl;
				exit(-1);
			}
			KINE_VALS = fvals;
		}
		if(flag == "SCAP"){
			if(fvals.size() != 3){
				G4cout << "SCAP flag requires 3 arguments!!" << endl;
				exit(-1);
			}
			SCAP_VALS = fvals;
		}
		if(flag == "BEAM"){
			if(fvals.size() != 7){
				G4cout << "BEAM flag requires 7 arguments!!" << endl;
				exit(-1);
			}
			BEAM_VALS = fvals;
		}
		if(flag == "BGGATE"){
			if(fvals.size() != 2){
				G4cout << "BGGATE flag requires 2 arguments!!" << endl;
				exit(-1);
			}
			BGGATE_VALS = fvals;
		}
		
	}
	
	ifs.close();
	
	if( (!KINE_VALS.empty()) && SCAP_VALS.empty() ){
		SCAP_VALS.push_back(0.0);
		SCAP_VALS.push_back(0.0);
		SCAP_VALS.push_back(65.0);
	}
	
	// If geometry files was not set by user, then look for default
	if(gGDMLfilename == ""){
		G4cout << "GDML geometry file not specified." << G4endl;
		G4cout << "  looking for ./cpproot.gdml ... "; G4cout.flush();
		ifstream ifs("./cpproot.gdml");
		if(ifs.is_open()){
			ifs.close();
			G4cout << "yes" << G4endl;
			gGDMLfilename = "./cpproot.gdml";
		}else {
			G4cout << "no" << G4endl;

			const char *HDDS_HOME  = getenv("HDDS_HOME");
			const char *BMS_OSNAME = getenv("BMS_OSNAME");
			if(!HDDS_HOME ) G4cout << "HDDS_HOME  environment variable not set!" << G4endl;
			if(!BMS_OSNAME) G4cout << "BMS_OSNAME environment variable not set!" << G4endl;
			if(HDDS_HOME!=NULL && BMS_OSNAME!=NULL){
		
				string tmpname = string(HDDS_HOME) + "/" + BMS_OSNAME + "/src/cpproot.gdml";
				G4cout << "  looking for " + tmpname + " ... "; G4cout.flush();
				ifstream ifs(tmpname.c_str());
				if(ifs.is_open()){
					ifs.close();
					G4cout << "yes" << G4endl;
					gGDMLfilename = tmpname;
				}else{
					G4cout << "no" << G4endl;
				}
			}
		}
	}
	if(gGDMLfilename == ""){
		G4cout << "--- Unable to identify GDML geometry file! ---" << G4endl;
		G4cout << "--- Please specify it using \"GEOMFILE\" in control.in ---" << G4endl;
		exit(-1);
	}
}

//--------------------
// PrintConfigInfo
//--------------------
void PrintConfigInfo(void)
{	
#ifdef G4MULTITHREADED
	bool mt_enabled = true;
	if(NTHREADS==0) NTHREADS = sysconf(_SC_NPROCESSORS_ONLN);
#else
	bool mt_enabled = false;
	if(NTHREADS==0) NTHREADS = 1;
	if(NTHREADS != 1) {
		G4cout << "---------------------------- WARNING ! ----------------------------" << endl;
		G4cout << "You have specified to run with " << NTHREADS << " processing threads" << endl;
		G4cout << "but this executable is not capable of running with multiple" << endl;
		G4cout << "threads. The number of threads is therefore reduced to 1" << endl;
		G4cout << "--------------------------------------------------------------------" << endl;
		NTHREADS = 1;
	}
#endif
	
	
	G4cout << "         MAXEVENTS : " << MAXEVENTS << endl;
	if(INPUT_FILENAME!=""){
		G4cout << "    INPUT_FILENAME : " << (INPUT_FILENAME=="" ? "<none>":INPUT_FILENAME) << endl;
	}else if( !KINE_VALS.empty() ){
		G4cout << "    Particle gun ------------------------------" << endl;
		G4cout << "                       type: " << (int)KINE_VALS[0] << " (" << ParticleType((Particle_t)(KINE_VALS[0]-100)) << ")" << endl;
		G4cout << "                          x: " << SCAP_VALS[0] << endl;
		G4cout << "                          y: " << SCAP_VALS[1] << endl;
		G4cout << "                          z: " << SCAP_VALS[2] << endl;
		G4cout << "                   momentum: " << KINE_VALS[1] << endl;
		G4cout << "                      theta: " << KINE_VALS[2] << endl;
		G4cout << "                        phi: " << KINE_VALS[3] << endl;
		G4cout << "             delta_momentum: " << KINE_VALS[4] << endl;
		G4cout << "                delta_theta: " << KINE_VALS[5] << endl;
		G4cout << "                  delta_phi: " << KINE_VALS[6] << endl;
		G4cout << "    -------------------------------------------" << endl;
	}
	G4cout << "   OUTPUT_FILENAME : " << OUTPUT_FILENAME << endl;
	G4cout << "     gGDMLfilename : " << gGDMLfilename << endl;
	G4cout << "          NTHREADS : " << NTHREADS << (mt_enabled ? " (MT enabled)":"") << endl;
	if(FMWPC_DEAD_RADIUS>0.0) G4cout << " FMWPC Dead Region : " << FMWPC_DEAD_RADIUS << " cm radius" << endl;
	if( FMWPC_DEAD_WIDTH>0.0) G4cout << " FMWPC Dead Region : " << FMWPC_DEAD_WIDTH << "cm x " << FMWPC_DEAD_WIDTH <<"cm square" << endl;
}