// VertexGen_ES.cc
// Contact person:Nuno Barros <nfbarros@hep.upenn.edu>
// See VertexGen_ES.hh for more details
//———————————————————————//
// - Local include
#include <RAT/VertexGen_ES.hh>

// - RAT includes
#include <RAT/ESgen.hh>
#include <RAT/GLG4PosGen.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/StringUtil.hh>
#include <RAT/Log.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>

// - Geant4 includes
#include <G4ParticleDefinition.hh>
#include <G4ParticleTable.hh>
#include <G4ThreeVector.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4Event.hh>
#include <Randomize.hh>
#include <globals.hh>

// - CLHEP includes
#include <CLHEP/Vector/LorentzVector.h>

#include <sstream>
#include <CLHEP/Units/PhysicalConstants.h>
#include <cmath>
using namespace CLHEP;

namespace RAT {

VertexGen_ES::VertexGen_ES(const char *arg_dbname)
: GLG4VertexGen(arg_dbname), fNuDir(0.,0.,0.), fEnableMessenger(true),fDBName("SOLAR"),fRandomDir(false)
{
#ifdef RATDEBUG
	detail << "[VertexGen_ES]:: in constructor with arg [" << arg_dbname << "]." << newline;
#endif
	// The involved particles are always these
	fElectron = G4ParticleTable::GetParticleTable()->FindParticle("e-");
	fNue = G4ParticleTable::GetParticleTable()->FindParticle("nu_e");
	fNumu = G4ParticleTable::GetParticleTable()->FindParticle("nu_mu");
	fElectronMass = fElectron->GetPDGMass();
	fESgen = new ESgen();
#ifdef RATDEBUG
	debug << "[VertexGen_ES]:: Leaving constructor." << newline;
#endif

}

VertexGen_ES::~VertexGen_ES() {
	if (fESgen) delete fESgen;
}


void VertexGen_ES::
GeneratePrimaryVertex(G4Event* argEvent,
		G4ThreeVector& dx,
		G4double dt,
		G4ThreeVector& nudir)
{
	fNuDir = nudir;
	this->GeneratePrimaryVertex(argEvent,dx,dt);
}

void VertexGen_ES::
GeneratePrimaryVertex(G4Event* argEvent,
		G4ThreeVector& dx,
		G4double dt)
{

#ifdef RATDEBUG
debug << "[VertexGen_ES]::GeneratePrimaryVertex : Creating event with nu from  [ " << fNuDir.x() << " "<< fNuDir.y() << " " << fNuDir.z() << " ] and time [ " << dt <<" ]." << newline;
debug << "[VertexGen_ES]::GeneratePrimaryVertex : Generating ES vertex with parameters: " << newline;
debug << "                  Position : " << dx.x() << " " << dx.y() << " " << dx.z() << newline;
debug << "                  Time     : " << dt << newline;
#endif
//
// Build the primary vertex with the position and time
//
G4PrimaryVertex* vertex = new G4PrimaryVertex(dx, dt);

//
// Build the incoming neutrino direction.
// If this generator was not called from a higher level generator it is possible that the neutrino direction is not set
//

if (fNuDir.mag2() == 0.0 || fRandomDir) {
#ifdef RATDEBUG
	warn << "[VertexGen_ES]::GeneratePrimaryVertex : Found zero magnitude in direction. Generating stopped neutrino (?)."<< newline;
	warn << "[VertexGen_ES]::GeneratePrimaryVertex : Generating a random direction at each call of the generator."<< newline;
#endif
	fRandomDir = true;
	// Pick isotropic direction
	double theta = acos(2.0 * G4UniformRand() - 1.0);
	double phi = G4UniformRand() * twopi;
	fNuDir.setRThetaPhi(1.0, theta, phi);
}

// Generate elastic-scattering interaction using ESgen.
// NB: Updated to keep track of the neutrino as well as the electron
// For the moment the ESgen does not full mom_nu so it is empty
//

CLHEP::HepLorentzVector mom_electron, mom_nu;
#ifdef RATDEBUG
detail << "[VertexGen_ES]::GeneratePrimaryVertex : Calling internal event generator with direction [ "
		<< fNuDir.x() << " "<< fNuDir.y() << " " << fNuDir.z() << " ]." << newline;
#endif

fESgen->GenerateEvent( fNuDir, mom_nu, mom_electron );

// -- Create particle at vertex
// FIXME: Should we keep track of the outgoing neutrino as well?
//       If so ESgen needs to be updated to pass the new neutrino direction.


G4PrimaryParticle* electron_particle =
		new G4PrimaryParticle(fElectron,             	// particle code
				mom_electron.px(),     				// x component of momentum
				mom_electron.py(),     				// y component of momentum
				mom_electron.pz());    				// z component of momentum
electron_particle->SetMass(fElectronMass);
vertex->SetPrimary( electron_particle );

// Add the incoming neutrino as the primary
G4PrimaryParticle *neutrinoparent;
if(this->GetNuFlavor()=="nue"){
	neutrinoparent = new G4PrimaryParticle(fNue,
				mom_nu.px(),
				mom_nu.py(),
				mom_nu.pz());
}else if(this->GetNuFlavor()=="numu"){
	neutrinoparent = new G4PrimaryParticle(fNumu,
				mom_nu.px(),
				mom_nu.py(),
				mom_nu.pz());
}else{
	warn << "[VertexGen_ES] : not found parent neutrino type " << this->GetNuFlavor()
	     << newline;
	return;
}

// We DON'T Add this one to the vertex as we don't want to propagate it but alongside
// so that the information is present for extraction in Gsim
ParentPrimaryVertexInformation *vertinfo = new ParentPrimaryVertexInformation();
vertinfo->SetVertexParentParticle(neutrinoparent);

vertex->SetUserInformation(vertinfo);
argEvent->AddPrimaryVertex(vertex);




}


void VertexGen_ES::SetState(G4String newValues)
{
	// If the messenger is disallowed this option should be disallowed as well
	if (!fEnableMessenger) {
		warn << "[VertexGen_ES]::SetState : The generator is in locked state (most likely part of the solar generator)." << newline;
		warn << "[VertexGen_ES]::SetState : All options passed through messengers will be ignored." << newline;
		return;
	}

	newValues = util_strip_default(newValues); // from GLG4StringUtil
	if (newValues.length() == 0) {
		// print help and current state
		warn << "[VertexGen_ES]::SetState : Current state of this VertexGen_ES:\n"
				<< " \"" << GetState() << "\"\n" << newline;
		warn << "[VertexGen_ES]::SetState : Format of argument to VertexGen_ES::SetState: \n"
				" \"nu_dir_x nu_dir_y nu_dir_z [db_name]:db_flux:nu_flavor\"\n"
				" where fNuDir is the initial direction of the incoming neutrino.\n"
				" Does not have to be normalized.  Set to \"0. 0. 0.\" for isotropic\n"
				" neutrino direction."
				<< newline;
		return;
	}
#ifdef RATDEBUG
	detail << "[VertexGen_ES]::SetState : Setting state to ["<<newValues<<"]" << newline;
#endif
	std::istringstream is(newValues.c_str());
	double x, y, z;
	std::string rest;
	is >> x >> y >> z >> rest;
	if (is.fail()) {
		warn << "[VertexGen_ES]::SetState : Failed to extract state from input string." << newline;
		warn << "[VertexGen_ES]::SetState : Input string : [" << is.str() << "]." << newline;
		warn << "[VertexGen_ES]::SetState : Format of argument to VertexGen_ES::SetState: \n"
				" \"nu_dir_x nu_dir_y nu_dir_z [db_name]:db_flux:nu_flavor\"\n"
				" where fNuDir is the initial direction of the incoming neutrino.\n"
				" Does not have to be normalized.  Set to \"0. 0. 0.\" for isotropic\n"
				" neutrino direction."
				<< newline;
		Log::Die("[VertexGen_ES]::SetState : No DB spectrum entry nor flavor specified.",1);
		return;
	}
	// We take care of normalising the input direction here
	if ( x == 0. && y == 0. && z == 0. )
		fNuDir.set(0., 0., 0.);
	else
		fNuDir = G4ThreeVector(x, y, z).unit();

	// Now check that everything else is in "rest"
	if (rest.length() == 0) {
		warn << "[VertexGen_ES]::SetState : No DB spectrum entry nor flavor specified." << newline;
		warn << "[VertexGen_ES]::SetState : Format of argument to VertexGen_ES::SetState: \n"
				" \"nu_dir_x nu_dir_y nu_dir_z [db_name]:db_flux:nu_flavor\"\n"
				" where fNuDir is the initial direction of the incoming neutrino.\n"
				" Does not have to be normalized.  Set to \"0. 0. 0.\" for isotropic\n"
				" neutrino direction."
				<< newline;
		Log::Die("[VertexGen_ES]::SetState : No DB spectrum entry nor flavor specified.",1);
	}
	// Get separate the entries by the ':'
	std::string parseparams = ":";
	rest = strip(rest,parseparams);
	std::vector<std::string> params = util_split(rest,parseparams);
	switch (params.size()) {

	case 3:
		// First entry is the database name
		this->SetDBName(params[0]);
		this->SetFlux(params[1]);
		this->SetNuFlavor(params[2]);
		break;
	case 2:
		this->SetFlux(params[0]);
		this->SetNuFlavor(params[1]);
		break;
	default:
		warn << "[VertexGen_ES]::SetState : Detected only "<< params.size() << " neutrino state terms (2 or 3 expected)."<< newline;
		warn << "[VertexGen_ES]::SetState : Input string ["<< rest << "]."<< newline;
		warn << "[VertexGen_ES]::SetState : Format of argument to VertexGen_ES::SetState: \n"
				" \"nu_dir_x nu_dir_y nu_dir_z [db_name]:db_flux:nu_flavor\"\n"
				" where fNuDir is the initial direction of the incoming neutrino.\n"
				" Does not have to be normalized.  Set to \"0. 0. 0.\" for isotropic\n"
				" neutrino direction."
				<< newline;
		Log::Die("[VertexGen_ES]::SetState : Wrong input state string in VertexGen_ES.");
		break;
	}



}


G4String VertexGen_ES::GetState()
{
	std::ostringstream os;

	os << fNuDir.x() << "\t" << fNuDir.y() << "\t" << fNuDir.z() << std::ends;

	G4String rv(os.str());
	return rv;
}

void VertexGen_ES::SetFlux(const G4String flux) {
	if (fFlux == flux) return;
	fFlux = flux;
	fESgen->SetNuType(fFlux);
}

void VertexGen_ES::SetNuFlavor(const G4String flavor) {
	if (fNuFlavor == flavor) return;
	fNuFlavor = flavor;
	fESgen->SetNuFlavor(fNuFlavor);
}

void VertexGen_ES::SetDBName(const G4String name) {
	if (fDBName == name) return;
	fDBName = name;
	fESgen->SetDBName(fDBName);
}

void VertexGen_ES::SetEnableMessenger( bool state ) {
	fEnableMessenger = state;
}

} // namespace RAT