// DecayChain_Gen.cc
// Contact person: Eric Vazquez-Jauregui <ericvj@fisica.unam.mx>
// See DecayChain_Gen.hh for more details
//———————————————————————//

#include <RAT/DecayChain_Gen.hh>
#include <RAT/DecayChain.hh>

#include <RAT/GLG4PosGen.hh>
#include <RAT/GLG4TimeGen.hh>
#include <RAT/Factory.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>
#include <RAT/Log.hh>

#include <G4Event.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4IonTable.hh>
#include <G4Electron.hh>
#include <G4Positron.hh>
#include <G4Gamma.hh>
#include <G4Alpha.hh>
#include <G4ThreeVector.hh>
#include <G4UnitsTable.hh>

#include <cstring>

#undef DEBUG

namespace RAT{

DecayChain_Gen::DecayChain_Gen(const G4String & name) : GLG4Gen(name), fStateStr(""),fDecayChain(0)
{
  // As in the combo generator, use a default time generator if the
  // user does not supply one.
  fTimeGen = new GLG4TimeGen_Poisson();
  fInMiddle=false;
  fInAlphaDecay=false;
  fInGammaDecay=false;
  fScreenCorr=false;
}

DecayChain_Gen::~DecayChain_Gen()
{ }

void DecayChain_Gen::GenerateEvent(G4Event* event)
{
  if (fDecayChain == 0) {
	  warn << "RAT::DecayChain_Gen::GenerateEvent: No decay chain found; "
	         << "please use /generator/add decaychain ISOTOPE:POSITION[:TIME][:midchain/alpha/gamma]\n";
	  return;
  }

  // Generate the position of the isotope.  Note that, for now, we
  // don't change the position of the isotope as it decays.
  G4ThreeVector position;
  fPosGen->GeneratePosition(position);

  // Generate the chain of decay products for this isotope.
  fDecayChain->GenerateFullChain();
  G4int nPrimaries = fDecayChain->GetNGenerated();

#ifdef DEBUG
  debug << "RAT::DecayChain_Gen::GenerateEvent: "
	<< nPrimaries << " primaries in decay chain "
  	<< fDecayChain->GetChainName() << "\n";
#endif

  for ( G4int iPrimary = 0; iPrimary < nPrimaries; iPrimary++ ) {
	   // particle type
	   G4int pid = fDecayChain->GetEventID(iPrimary); // not the PDG code!

	   G4ParticleDefinition* particleDef = 0;
	   if ( pid >= 100000000 ) {
	      G4int A = (pid-100000000)/1000;
	      G4int Z = (pid-100000000) - A*1000;
	      G4double excitationEnergy = 0.0; // assume all in ground state
	      G4IonTable* theIonTable = G4IonTable::GetIonTable();
	      particleDef = theIonTable->GetIon(Z, A, excitationEnergy);
	   }
	   else {
	      if ( pid == DecayBeta ) {
	         particleDef = G4Electron::Electron();
	      }
              else if ( pid == DecayEC ) {
                 particleDef = G4Positron::Positron();
              }
	      else if ( pid == DecayGamma ) {
	         particleDef = G4Gamma::Gamma();
	      }
	      else if( pid == DecayAlpha ) {
	         particleDef = G4Alpha::Alpha();
	      }
	   }
	   if ( particleDef == 0 ) {
	      warn << "RAT::DecayChain_Gen::GenerateEvent: "
		          << "didn't know how to handle particle with ID " << pid << "\n";
	      continue;
	   }

	   // Get the particle's information (momentum and time)
	   DecayChain::ParticleInfo_t particleInfo = fDecayChain->GetParticleInfo(iPrimary);
	   G4double time = NextTime() + fDecayChain->GetEventTime(iPrimary);
    	// generate a vertex with a primary particle
	   G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time);
	   G4PrimaryParticle* particle = new G4PrimaryParticle(particleDef,
							                                     particleInfo.vector.px(),
							                                     particleInfo.vector.py(),
							                                     particleInfo.vector.pz());
	   particle->SetMass(particleDef->GetPDGMass()); // Apparently this is useful in IBD code.
	   vertex->SetPrimary(particle);
	   event->AddPrimaryVertex(vertex);

// Set the parent information for this vertex (but don't generate parent!)
           G4IonTable *theIonTable = G4IonTable::GetIonTable();
           G4ParticleDefinition *parentDef = theIonTable->GetIon(fDecayChain->GetParentCharge(iPrimary),fDecayChain->GetParentMass(iPrimary),0.);
           G4PrimaryParticle *parent = new G4PrimaryParticle(parentDef,0,0,0);


           ParentPrimaryVertexInformation *vertinfo = new ParentPrimaryVertexInformation();
           vertinfo->SetVertexParentParticle(parent);
           vertex->SetUserInformation(vertinfo);

#ifdef DEBUG
	debug << "RAT::DecayChain_Gen::GenerateEvent: "
	      << "Primary " << iPrimary+1 << " of " << nPrimaries
	      << ", pid=" << pid
	      << ", name=" << particleDef->GetParticleName() << "\n";
	debug << "    time=" << G4BestUnit(time,"Time")
	      << ", position=" << G4BestUnit(position,"Length")
	      << ", momentum=" << G4BestUnit(particleInfo.vector,"Energy")
	      << "\n";
#endif

    } // for each primary
}

void DecayChain_Gen::ResetTime(double offset)
{
  double eventTime = fTimeGen->GenerateEventTime();
  fNextTime = eventTime + offset;
#ifdef DEBUG
  debug << "RAT::DecayChain_Gen::ResetTime:"
	      << " eventTime=" << G4BestUnit(eventTime,"Time")
   	   << ", offset=" << G4BestUnit(offset,"Time")
	      << ", nextTime=" << G4BestUnit(fNextTime,"Time")
	      << "\n";
#endif
}

void DecayChain_Gen::SetState(G4String state)
{
#ifdef DEBUG
  debug << "RAT::DecayChain_Gen::SetState called with state='"
	      << state << "'" << "\n";
#endif

  // Break the argument to the this generator into sub-strings
  // separated by ":".
  state = util_strip_default(state);
  vector<string> parts = util_split(state, ":");
  size_t nArgs = parts.size();

#ifdef DEBUG
  debug << "RAT::DecayChain_Gen::SetState: nArgs = " << nArgs << "\n";
#endif

  try {
   if ( nArgs >= 5 ){
           // The fifth argument allows one to choose to use a screening
           // correction to the Fermi Function
       string ScreenCorr = parts[4];
       if (ScreenCorr=="screen"){
          fScreenCorr=true;
          }
       if (ScreenCorr=="none"){
          fScreenCorr=false;
          }
       }

   if ( nArgs >= 4 ){
           // The fourth argument allows one to start midchain
       string InMiddle = parts[3];
       if (InMiddle=="midchain"){
          info << "RAT::DecayChain_Gen: starting chain at Isotope\n";
          fInMiddle=true;
          }
       if (InMiddle=="alpha"){
          info << "RAT::DecayChain_Gen: alpha decay of Isotope\n";
          fInAlphaDecay=true;
          }
       else if (InMiddle=="gamma"){
          info << "RAT::DecayChain_Gen: gamma decay of Isotope\n";
          fInGammaDecay=true;
          }
        }

   if ( nArgs >= 3 ) {
	   // The third argument is an optional time generator
	   delete fTimeGen; fTimeGen = 0; // In case of exception in next line
	   fTimeGen = GlobalFactory<GLG4TimeGen>::New(parts[2]);
	}

   if ( nArgs >= 2 )	{
	   // The first argument is the isotope that starts the decay
	   // chain.
	   string isotope = parts[0];

	   // Don't bother building a new decay chain if we already have
	   // one for that isotope.
	   if ( fDecayChain == 0  || isotope != fDecayChain->GetChainName() ) {
	      delete fDecayChain;
	      fDecayChain = new DecayChain(isotope);
#ifdef DEBUG
	      fDecayChain->SetVerbose(true);
#endif
         if (fInMiddle) fDecayChain->SetInMiddleChain(true);
         if (fInAlphaDecay) fDecayChain->SetAlphaDecayStart(true);
         if (fInGammaDecay) fDecayChain->SetGammaDecayStart(true);
         if (fScreenCorr) fDecayChain->SetScreenCorr(true);

	      bool found = fDecayChain->ReadInputFile(fDecayChain->GetChainName());
         info << "found: "<<found<<", fDecayChain->GetChainName(): "<<fDecayChain->GetChainName() << "\n";
	      if (!found) {
		 warn << "RAT::DecayChain_Gen::SetState: couldn't find data for isotope "
	              << isotope << "\n";
		 delete fDecayChain; fDecayChain=0;
              }
	      info << "RAT::DecayChain_Gen::SetState: successfully created decay chain for "
		         << fDecayChain->GetChainName() << "\n";
#ifdef DEBUG
	      fDecayChain->Show();
#endif
	    }

	    // The second argument is a position generator.
	    delete fPosGen; fPosGen = 0;
	    fPosGen = GlobalFactory<GLG4PosGen>::New(parts[1]);
	}
   else	{
	  G4Exception(__FILE__, "Invalid Parameter", FatalException, ("decaychain generator syntax error: '"+
								      state+
								      "' does not have a position generator").c_str());
	}

   fStateStr = state; // Save for later call to GetState()
    } catch (FactoryUnknownID &unknown) {
      warn << "Unknown generator \"" << unknown.id << "\"\n";
   }
}

G4String DecayChain_Gen::GetState() const
{
  return fStateStr;
}

void DecayChain_Gen::SetTimeState(G4String state)
{
  if (fTimeGen)
      fTimeGen->SetState(state);
  else
      warn << "DecayChain_Gen error: Cannot set time state, no time generator selected\n";
}

G4String DecayChain_Gen::GetTimeState() const
{
  if (fTimeGen)
     return fTimeGen->GetState();
  else
     return G4String("DecayChain_Gen error: no time generator selected");
}

void DecayChain_Gen::SetPosState(G4String state)
{
  if (fPosGen)
     fPosGen->SetState(state);
  else
     warn << "DecayChain_Gen error: Cannot set position state, no position generator selected\n";
}

G4String DecayChain_Gen::GetPosState() const
{
  if (fPosGen)
     return fPosGen->GetState();
  else
     return G4String("DecayChain_Gen error: no pos generator selected");
}


}