// VertexGen_PrimakoffRosen.cc
// Contact person:Jeanne Wilson <j.r.wilson@qmul.ac.uk>
// See VertexGen_PrimakoffRosen.hh for more details
//———————————————————————//
#include <CLHEP/Units/PhysicalConstants.h>
using namespace CLHEP;

#include <RAT/VertexGen_PrimakoffRosen.hh>
#include <RAT/Log.hh>
#include <Randomize.hh>
#include <G4Event.hh>
#include <G4ParticleTable.hh>
#include <G4IonTable.hh>
#include <G4PrimaryParticle.hh>
#include <G4PrimaryVertex.hh>
#include <G4ThreeVector.hh>
#include <G4ParticleMomentum.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>
#include <G4IonTable.hh>
#include <G4Ions.hh>

#define G4std  std
namespace RAT {

VertexGen_PrimakoffRosen::VertexGen_PrimakoffRosen(const char *arg_dbname): GLG4VertexGen(arg_dbname)
{
	// constructor
	fELimUniLow  = 0;
	fELimUniHigh  = 0;
	fELimTempLow  = 0;
	fELimTempHigh  = 0;
	// we are working with electrons
	fPDef = G4ParticleTable::GetParticleTable()->FindParticle("e-");
}

///-------------------------------------------------------------------------
VertexGen_PrimakoffRosen::~VertexGen_PrimakoffRosen()
{
	// destructor
}

///-------------------------------------------------------------------------
void VertexGen_PrimakoffRosen::GeneratePrimaryVertex(G4Event *event, G4ThreeVector &dx, G4double dt)
{
	G4PrimaryVertex* vertex= new G4PrimaryVertex( dx, dt );

	G4double E1, E2, btheta;
	// which mode of decay (set particle definition inside this loop
	// to allow flexibility for adding positron modes etc.
	if (fMode==0){
	   // choose electron energies (total) and separation angle for
	   // neutrinoless double beta decay
	   ZeroNu2nBasic(E1,E2,btheta);
	}
        else if (fMode==2){
	   // chooseelectron energies (total) and separation angle for
	   // two neutrino double beta decay
	   // Option to limit the energies for this decay
	   G4double Elo = this->EMinimum();
	   G4double Ehi = this->EMaximum();
	   G4double sumE = -999;		// kinetic energy
	   while (!(sumE>Elo&&sumE<Ehi)){
		 TwoNu2nBasic(E1,E2,btheta);
		 sumE = E1+E2-2*fPDef->GetPDGMass();
	   }
	   debug << "Simulating 2nu event with " << sumE << " " << E1 << " " << E2 << "\n";
	   // after generating we must reset any temporary limits on the energy range
	   fELimTempLow = fELimUniLow;
	   fELimTempHigh = fELimUniHigh;
	}
        else {
	   // We shouldn't get to here - there is a problem
	   warn << "Problem in VertexGen_Double: no known mode for decay \n";
	}

	// Generate a random direction for the first beta
  	G4double theta, costheta, sintheta, phi;
	G4double PI = acos(-1);
  	costheta = G4UniformRand()*2. - 1.;		// Default isotropic events
  	sintheta = sqrt(1.-costheta*costheta);
  	theta = acos(costheta);
  	phi   = G4UniformRand()*2.*PI;
	G4ParticleMomentum momdir = G4ParticleMomentum(sin(phi)*sintheta, cos(phi)*sintheta, costheta);

        // Calculate the momentum for this first beta
	G4double pmom, px, py, pz;
	G4double E_Mass = fPDef->GetPDGMass();
	pmom = sqrt(std::pow((E1),2)-std::pow(E_Mass,2));
        px = pmom*momdir.x();
        py = pmom*momdir.y();
        pz = pmom*momdir.z();

	// Do I care about polarisation?

	// Add this first beta to the vertex
	G4PrimaryParticle* beta1 = new G4PrimaryParticle(fPDef,px,py,pz);
        vertex->SetPrimary(beta1);

	// Calculate the direction of the second beta using the angle wrt the first beta
	G4ParticleMomentum momdir2 = G4ParticleMomentum(sin(phi)*sin(theta+btheta), cos(phi)*sin(theta+btheta), cos(theta+btheta));

	// We've set the angle between the two vectors but they can be random in phi
	// so rotate around beta1 direction
  	G4double phi2   = G4UniformRand()*2.*PI;
	momdir2.rotate(phi2,G4ThreeVector(sin(phi)*sintheta, cos(phi)*sintheta,  costheta));
	// Calculate the momentum
	pmom = sqrt(std::pow((E2),2)-std::pow(E_Mass,2));
        px = pmom*momdir2.x();
        py = pmom*momdir2.y();
        pz = pmom*momdir2.z();

	// Do I care about polarisation?

	// Add this second beta to the vertex
	G4PrimaryParticle* beta2 = new G4PrimaryParticle(fPDef,px,py,pz);
	vertex->SetPrimary(beta2);

	// Which isotope are we generating
	int Z = 0;
	int A = 0;
	if(fIsotope == "Nd150"){
		Z = 60;
		A = 150;
	}else if(fIsotope == "Nd148"){
		Z = 60;
		A = 148;
	}else if(fIsotope == "Te130"){
		Z = 52;
		A = 130;
	}else{
		Log::Die("Unknown Isotope in VertexGen_PrimakoffRosen : " + fIsotope);
	}
    G4IonTable* ionTable = G4IonTable::GetIonTable();
    G4ParticleDefinition* particleDef = ionTable->GetIon(Z, A, 0);

    // Add the parent isotope information (stationary, no excitation)
    G4PrimaryParticle *isotopeparent = new G4PrimaryParticle(particleDef,0,0,0);
    // 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(isotopeparent);
	vertex->SetUserInformation(vertinfo);

	// And add the vertex to the event
	event->AddPrimaryVertex(vertex);
}

///-------------------------------------------------------------------------
void VertexGen_PrimakoffRosen::SetState(G4String newValues)
{
	if (newValues.length() == 0) {
	   // print help and current state
	   info	<< "Current state of this VertexGen_PrimakoffRosen:\n"
     	  	<< " \"" << GetState() << "\"\n";
	   info	<< "Format of argument to VertexGen_PrimakoffRosen::SetState: \n"
		<< " \"isotope mode (Elo) (Ehi)\"\n"
		<< " isotope = parent isotope (Nd150, Nd148) \n"
		<< " mode = decay mode (0=neutrinoless, 2=2neutrino) \n"
		<< " Elo Ehi = optional limits on the energy range of emitted particles \n";
	   return;
	}

	std::istringstream is(newValues.c_str());
	G4std::string isotope;
	// Read in the isotope
	is >> isotope;
	if (is.fail() || isotope.length()==0){
	   Log::Die("VertexGen_PrimakoffRosen: Incorrect vertex setting " + newValues);
	}
	// check that isotope is in database - not sure what happens if this fails - tidy!
	DBLinkPtr ltemp = DB::Get()->GetLink("DBG", isotope);
	try {
	  ltemp->GetD("q_value");
	  fLdbdb  = ltemp;
	  fIsotope = isotope;
	  info << "PrimakoffRosen Vertex: Setting isotope " << isotope << "\n";
	}
	catch(DBNotFoundError &e) {
	  Log::Die("VertexGen_PrimakoffRosen: can't find database file for " + isotope);
	}

	// Read in the mode
	int mode;
	is >> mode;
	// Should be 0 or 2
	if (mode==0||mode==2){
	   fMode = mode;
	}
        else {
	   fMode = 0;
	   warn	<< "VertexGen_PrimakoffRosen: Don't understand mode " << mode
		<< " defaulting to 0nuBB decay \n";
 	}

	// Has the user supplied energy limits too?
	float Elo, Ehi;
	is >> Elo >> Ehi;
	if (is.fail()){
	   // obviously not, so make the universal limits huge
	   fELimUniLow = 0;
	   fELimUniHigh = 9999999;
	}
        else {
	   // user has selected universal limits - these only make sense for 2nu mode
	   if (mode==2){
	      if (Ehi>Elo&&Ehi>0&&Elo>=0){
		 fELimUniLow = Elo;
		 fELimUniHigh = Ehi;
		 info << "PrimakoffRosen Vertex: Setting energy limits "
		      << Elo << " - " << Ehi << "\n";
	      }
              else {
		 warn << "PrimakoffRosen Vertex: Nonsensical energy limits "
		      << Elo << " - " << Ehi << " not applied \n";
		 fELimUniLow = 0;
		 fELimUniHigh = 9999999;
	      }
	    }
            else {
	      warn << "PrimakoffRosen Vertex: Energy limits can't be applied to "
		   << " neutrinoless mode \n";
	      fELimUniLow = 0;
	      fELimUniHigh = 9999999;
	    }
	}
	// temporary limits should be the same as the universal ones
	fELimTempLow = fELimUniLow;
	fELimTempHigh = fELimUniHigh;

	return;
}

///-------------------------------------------------------------------------
G4String VertexGen_PrimakoffRosen::GetState()
{
	// State setting is the isotope name and mode	- return it
	G4std::ostringstream os;
	os << fIsotope << " " << fMode << G4std::ends;
	G4String rv(os.str());
	return rv;
}

///-------------------------------------------------------------------------

void VertexGen_PrimakoffRosen::ZeroNu2nBasic(G4double& E1, G4double& E2, G4double& btheta){

	// Calculate/Select the energies of the two betas and the angle between them
	// See Tretyak&Zdesenko, "Tables of Double Beta Decay Data", Atomic Data and Nuclear
	// Data Tables, 61, 43-90 (1995) for equations.

	// we will work in terms of electron mass
	G4double E_Mass = fPDef->GetPDGMass();

	G4double t0, t1, t2, shot, beta12, cosphi;

	t0 = (fLdbdb->GetD("q_value")*MeV)/E_Mass;	// work with energy in units of beta mass

	// sample an energy for the first beta
	do {
	   t1 = G4UniformRand()*t0;
	   // shot is selected between zero and maximum of distribution
	   shot = G4UniformRand()*std::pow((t0/2.0 + 1.0),4);
	}
        while( shot > std::pow((t1+1.0),2)*std::pow((t0+1.0-t1),2) );

	E1 = (t1+1)*E_Mass;	// +1 for total energy, t1 is kinetic

	// rest of energy goes to the second beta
	t2 = t0 - t1;
	E2 = (t2+1)*E_Mass;	// +1 for total energy, t2 is kinetic

	// need angular distribution of this vector with respect to other one
	beta12 = sqrt(t1*(t1+2.0)*t2*(t2+2.0))/((t1+1)*(t2+1));
	do {
	   cosphi = 2.0*G4UniformRand()-1.0;
	   // again shot is selected between zero and maximum of distribution
	   // maximum is when costheta = -1 -> (1+beta12) term.
	   shot = G4UniformRand()*(1.0+beta12);
	}
        while (shot > (1.0-beta12*cosphi));

	// this gives us the angle between the two betas
	btheta = acos(cosphi);
	return;
}

///-------------------------------------------------------------------------

void VertexGen_PrimakoffRosen::TwoNu2nBasic(G4double& E1, G4double& E2, G4double& btheta){

	// Calculate/Select the energies of the two betas and the angle between them
	// See Tretyak&Zdesenko, "Tables of Double Beta Decay Data", Atomic Data and Nuclear
  	// Data Tables, 61, 43-90 (1995) for equations.

	// we will work in terms of electron mass
	G4double E_Mass = fPDef->GetPDGMass();

	G4double t0, ts, t1, t2, shot, qfval, sval, beta12, cosphi;

	t0 = fLdbdb->GetD("q_value")/E_Mass;	// work with energy in units of beta mass
	qfval = fLdbdb->GetD("maxPD")/E_Mass;	// this is the maximum of the distribution

	// sample an energy for the total beta energy
	do {
	   ts = G4UniformRand()*t0;
	   shot = G4UniformRand()*qfval;
	   sval = (std::pow(ts,4)+10*std::pow(ts,3)+40*std::pow(ts,2)+60*ts+30)
					*ts*std::pow((t0-ts),5);
	   if (sval > qfval){
	      warn << "VertexGen_PrimakoffRosen::TwoNu2nBasic - warning: \n"
		   << "Shot exceeded assumed max.!\n";
	   }
	}
        while(shot > sval);

	// sample an energy for the first beta
	do {
	   t1 = G4UniformRand()*ts;
	   shot = G4UniformRand()*std::pow((ts/2.0+1.0),4);
	   sval = std::pow((t1+1.),2)*std::pow((ts-t1+1.),2);
	} while(shot > sval);

	E1 = (t1+1)*E_Mass;	// +1 for total energy, t1 is kinetic

	// rest of energy goes to the second beta
	t2 = ts - t1;
	E2 = (t2+1)*E_Mass;	// +1 for total energy, t2 is kinetic

	// need angular distribution of this vector with respect to other one
	beta12=sqrt(t1*(t1+2.0)*t2*(t2+2.0))/((t1+1)*(t2+1));
	do {
	   cosphi = 2.0*G4UniformRand()-1.0;
	   shot = G4UniformRand()*(1.0+beta12);
	}
        while(shot > (1.0-beta12*cosphi));

	// this gives us the angle between the two betas
	btheta = acos(cosphi);
	return;
}

///-------------------------------------------------------------------------
bool VertexGen_PrimakoffRosen::ELimitable()
{
 	if (fMode==2){
	   return true;
	}
        else {
	   return false;
	}
}

///-------------------------------------------------------------------------
void VertexGen_PrimakoffRosen::LimitEnergies(float Elo, float Ehi)
{
	// Set the limits for the generated energy range
	// first check this makes sense
	if (Elo>fLdbdb->GetD("q_value")||((Ehi-Elo)<=0)){
	   warn	<< "VertexGen_PrimakoffRosen: temporary energy limits " << Elo << " - " << Ehi
	        << " don't make sense, not applied \n"
		<< "Q value is " << fLdbdb->GetD("q_value")
		<< "\n";
	   return;
	}
	fELimTempLow = Elo;
	fELimTempHigh = Ehi;
	debug << "Applying temporary energy limits " << fELimTempLow << " - " << fELimTempHigh << "\n";
	return;
}

///-------------------------------------------------------------------------
float VertexGen_PrimakoffRosen::EMaximum()
{
	// return the maximum possible energy (ie. Qvalue) unless a limit applies
	float Q = float(fLdbdb->GetD("q_value"));
	if (Q<fELimTempHigh&&Q<fELimUniHigh){
	   return Q;
	}
        else if (fELimTempHigh< fELimUniHigh){
	   return fELimTempHigh;
	}
        else {
	   return fELimUniHigh;
	}
}

///-------------------------------------------------------------------------
float VertexGen_PrimakoffRosen::EMinimum()
{
	// return the minimum possible energy (which is 0 unless limits have been applied)
	if (fELimTempLow<0&&fELimUniLow<0){
 	   return 0.;
	}
        else if(fELimTempLow<fELimUniLow){
	   return fELimUniLow;
	}
        else{
	   return fELimTempLow;
	}
}

///-------------------------------------------------------------------------
}