#include <RAT/PhysicsList.hh>
#include <RAT/PhysicsListMessenger.hh>

#include <RAT/EmPhysicsListV3.hh>
#include <RAT/EmPhysicsListV4.hh>
#include <RAT/EmPhysicsListV5.hh>

#include <RAT/HadronicPhysicsListV4.hh>
#include <RAT/HadronicPhysicsList_borex.hh>

#include <RAT/Log.hh>

#include <G4LossTableManager.hh>
#include <G4UnitsTable.hh>
#include <G4ProcessManager.hh>

// ConstructParticle
#include <G4BosonConstructor.hh>
#include <G4LeptonConstructor.hh>
#include <G4BaryonConstructor.hh>
#include <G4MesonConstructor.hh>
#include <G4IonConstructor.hh> // Deuteron/Alpha
#include <G4ShortLivedConstructor.hh> //enuclear processes

// ConstructProcess
#include <G4EmProcessOptions.hh>

// AddParameterisation
#include <G4FastSimulationManagerProcess.hh>

// SelectEMPhysicsList
#include <G4EmStandardPhysics.hh>
#include <G4EmLivermorePhysics.hh>
#include <G4EmPenelopePhysics.hh>

// ConstructENuclear
#include <G4PhotoNuclearProcess.hh>
#include <G4TheoFSGenerator.hh>
#include <G4GammaParticipants.hh>
#include <G4GammaConversionToMuons.hh>
#include <G4QGSModel.hh>
#include <G4QGSMFragmentation.hh>
#include <G4ExcitedStringDecay.hh>
#include <G4GeneratorPrecompoundInterface.hh>
#include <G4ElectronNuclearProcess.hh>
#include <G4PositronNuclearProcess.hh>
#include <G4ElectroVDNuclearModel.hh>
#include <G4MuonNuclearProcess.hh>
#include <G4CascadeInterface.hh>
#include <G4MuonVDNuclearModel.hh>

// AddDecay
#include <G4Decay.hh>
#include <G4RadioactiveDecay.hh>

// AddAttenuation
#include <G4OpAbsorption.hh>
#include <G4OpBoundaryProcess.hh>
#include <G4OpMieHG.hh>
#include <RAT/Cerenkov.hh>
#include <RAT/SNOMANCerenkov.hh>
#include <RAT/OpRayleigh.hh>
#include <RAT/SNOMANOpRayleigh.hh>

#include <RAT/GLG4Scint.hh>

// SetCuts
#include <G4Region.hh>
#include <G4ProductionCuts.hh>
#include <G4Gamma.hh>
#include <G4Electron.hh>
#include <G4Positron.hh>
#include <G4Proton.hh>

#include <RAT/TrackKiller.hh>
// General
#include <CLHEP/Units/SystemOfUnits.h>
using namespace CLHEP;
#include <G4Version.hh>

bool RAT::PhysicsList::fsOmitHadronic = false;
bool RAT::PhysicsList::fsOmitMuonic = false;

RAT::PhysicsList::PhysicsList() : G4VModularPhysicsList()
{
  g4vuplInstanceID = subInstanceManager.CreateSubInstance();

  G4LossTableManager::Instance();

  double currentDefaultCut   = 40.0*um;
  fTimeCut = -1.0; // Equivalent to off
  fGammaCut         = currentDefaultCut;
  fElectronCut      = currentDefaultCut;
  fPositronCut      = currentDefaultCut;
  fProtonCut        = currentDefaultCut;
  fMuonCut          = currentDefaultCut;

  fPhyListMessenger = new PhysicsListMessenger(this);

  SetVerboseLevel(0);

  // Default EM physics
  fEMName = G4String("V5");
  fEMPhysicsList = new EmPhysicsListV5(0);
  // Default Hadronic physics
  fHadronicName = G4String("V4");
  fHadronicPhysicsList = new HadronicPhysicsListV4();
  // Default Cerenkov process
  fCerenkovName = G4String("Cerenkov");
  // Default OpRayleigh process
  fOpRayleighName = G4String("OpRayleigh");

  //Default ALL processes turned on
  fOmitAll = false;
  fOmitOpticalAttenuation = false;
  fOmitOpticalAbsorption = false;
  fOmitOpticalRayleigh = false;
  fOmitOpticalMie = true;
  fsOmitHadronic = false;
  fsOmitMuonic = false;
  fOmitHadronicListOnly = false;
  fOmitCerenkov = false;
  fOmitOpticalBoundaryEffects = false;
}

RAT::PhysicsList::~PhysicsList()
{
  delete fEMPhysicsList;
  delete fHadronicPhysicsList;
  delete fPhyListMessenger;
}

////////////////////////////////////////////////////////////////////////
/// ConstructParticles - Specific To RAT
///
/// Currently construct all particles, doesn't seem to be any penalty, so
/// should construct all possible - PGJ
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::ConstructParticle()
{
  G4BosonConstructor pBosonConstructor;
  pBosonConstructor.ConstructParticle();

  G4LeptonConstructor pLeptonConstructor;
  pLeptonConstructor.ConstructParticle();

  G4BaryonConstructor pBaryonConstructor;
  pBaryonConstructor.ConstructParticle();

  G4MesonConstructor pMesonConstructor;
  pMesonConstructor.ConstructParticle();

  G4ShortLivedConstructor pShortLivedConstructor;
  pShortLivedConstructor.ConstructParticle();

  G4IonConstructor pIonConstructor;
  pIonConstructor.ConstructParticle();
}

////////////////////////////////////////////////////////////////////////
/// ConstructProcess
///
/// Construct the physical processes.
/// RAT requires specialist Attenuation, parameterisation, decays, em and
/// hadronic processes.
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::ConstructProcess()
{
  if( fOmitAll )
    {
      info << "PhysicsList::ConstructProcess: Omitting ALL physics." << newline;
      return;
    }
  AddTransportation(); // Standard G4 transportation
  AddParameterisation();
  fEMPhysicsList->ConstructProcess();

  if( fsOmitHadronic )
    info << "PhysicsList: all hadronic processes omitted" << newline;
  else if( fOmitHadronicListOnly )
  {
    info << "PhysicsList: only hadronic physics list processes omitted" << newline;
    AddENuclear();
  }
  else
  {
    fHadronicPhysicsList->ConstructProcess();
    AddENuclear();
  }

  AddDecay();
  AddAttenuationAndCerenkov();
}

////////////////////////////////////////////////////////////////////////
/// AddParameterisation - Specific to RAT
///
/// Choose where the fast simulation process should act (pre, post or along step)
/// RAT chooses postStep. Also chooses to enable this only for opticalphotons.
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::AddParameterisation()
{
  G4FastSimulationManagerProcess* theFastSimulationManagerProcess = new G4FastSimulationManagerProcess();
  theParticleIterator->reset();
  while( (*theParticleIterator)() )
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      // Parameterisation only exists for optical photons PHIL
      if( particle->GetParticleName() != "opticalphoton" )
        continue;

      pmanager->AddProcess(theFastSimulationManagerProcess, -1, -1, 1);
    }
}

////////////////////////////////////////////////////////////////////////
/// SelectEMPhysicsList
///
/// Allows a standard em physics list to be used instead of the specialised
/// RAT EM physics list
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::SelectEMPhysicsList(const G4String& name)
{
  if( verboseLevel > 0 )
    RAT::info << "PhysicsList::SelectEMPhysicsList: Selecting " << name << " EM List\n";

  if( name == fEMName )
    return;
  fEMName = name;

  if( fEMName == "V1" )
    Log::Die( "PhysicsList::SelectEMPhysicsList, V1 list is no longer compatible." );
  else if( fEMName == "V2" )
    Log::Die( "PhysicsList::SelectEMPhysicsList, V2 list is no longer compatible." );
  else if( fEMName == "V3" )
    {
      delete fEMPhysicsList;
      fEMPhysicsList = new EmPhysicsListV3(0);
    }
  else if( fEMName == "V4" )
    {
      delete fEMPhysicsList;
      fEMPhysicsList = new EmPhysicsListV4(0);
    }
  else if( fEMName == "V5" )
    {
      delete fEMPhysicsList;
      fEMPhysicsList = new EmPhysicsListV5(0);
    }
  else if( fEMName == "standard" )
    {
      delete fEMPhysicsList;
      fEMPhysicsList = new G4EmStandardPhysics();
    }
  else if( fEMName == "penelope" )
    {
      delete fEMPhysicsList;
      fEMPhysicsList = new G4EmPenelopePhysics();
    }
  else if( fEMName == "livermore" )
    {
      delete fEMPhysicsList;
      fEMPhysicsList = new G4EmLivermorePhysics();
    }
  else
    RAT::warn << "PhysicsList::SelectEMPhysicsList: " << fEMName << " is not defined.\n";
}

////////////////////////////////////////////////////////////////////////
/// SelectHadronicPhysicsList
///
/// Allows the hadronic physics to be selected
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::SelectHadronicPhysicsList(const G4String& name)
{
  if( verboseLevel > 0 )
    RAT::info << "PhysicsList::SelectHadronicPhysicsList: Selecting " << name << " Hadronic List\n";

  if( name == fHadronicName )
    return;
  fHadronicName = name;

  if( fHadronicName == "V1" )
    Log::Die( "PhysicsList::SelectHadronicPhysicsList, V1 list is no longer compatible." );
  else if( fHadronicName == "V2" )
    Log::Die( "PhysicsList::SelectHadronicPhysicsList, V2 list is no longer compatible." );
  else if( fHadronicName == "V3" )
    Log::Die( "PhysicsList::SelectHadronicPhysicsList, V3 list is no longer compatible." );
  else if( fHadronicName == "V4" )
    {
      delete fHadronicPhysicsList;
      fHadronicPhysicsList = new HadronicPhysicsListV4();
    }
  else if( fHadronicName == "borex" )
    {
      RAT::warn << "PhysicsList::SelectHadronicPhysicsList: " << fHadronicName << " selected (V4 is default) \n";
      delete fHadronicPhysicsList;
      fHadronicPhysicsList = new HadronicPhysicsList_borex();
    }
  else
    RAT::warn << "PhysicsList::SelectHadronicPhysicsList: " << fHadronicName << " is not defined options V4 or borex available \n";
}

////////////////////////////////////////////////////////////////////////
/// AddENuclear - Specific to RAT
///
/// This is inelastic nuclear scattering of gammas, e, mu. Includes option
/// to omit the muonic processes.
/// Do we need e/mu nuclear interactions?
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::AddENuclear()
{
  //a flat variant of the G4's electronuclear builder because of the ratdb muon flags
  G4String particleName = "";

  theParticleIterator->reset();
  while ((*theParticleIterator)())
    {
      G4ParticleDefinition* particle = theParticleIterator->value();
      G4ProcessManager* pmanager = particle->GetProcessManager();
      particleName = particle->GetParticleName();
      if (particleName == "gamma")
        {
          G4PhotoNuclearProcess* g_nuc = new G4PhotoNuclearProcess();
          g_nuc->RegisterMe( new G4CascadeInterface() );
          pmanager->AddDiscreteProcess( g_nuc );
#if G4VERSION_NUMBER < 950 // Only geant4.9.4 or lower include this
          //ftf?
          G4TheoFSGenerator* theModel = new G4TheoFSGenerator(); //above 3.0 GeV
          G4QGSModel<G4GammaParticipants>* theStringModel = new G4QGSModel<G4GammaParticipants>;
          G4QGSMFragmentation* theFragmentation = new G4QGSMFragmentation();
          G4ExcitedStringDecay* theStringDecay = new G4ExcitedStringDecay(theFragmentation);
          theStringModel->SetFragmentationModel(theStringDecay);
          G4GeneratorPrecompoundInterface* theCascade = new G4GeneratorPrecompoundInterface;
          theModel->SetTransport(theCascade);
          theModel->SetHighEnergyGenerator(theStringModel);
          theModel->SetMinEnergy(3.0*GeV);
          theModel->SetMaxEnergy(100.0*GeV);
          g_nuc->RegisterMe(theModel);
#endif
          pmanager->AddDiscreteProcess(g_nuc);
          pmanager->AddDiscreteProcess(new G4GammaConversionToMuons);
        }
      else if (particleName == "e-")
        {
          G4ElectronNuclearProcess* e_nuc = new G4ElectronNuclearProcess();
          e_nuc->RegisterMe(new G4ElectroVDNuclearModel);
          pmanager->AddDiscreteProcess(e_nuc);
        }
      else if (particleName == "e+")
        {
          G4PositronNuclearProcess* ep_nuc = new G4PositronNuclearProcess();
          ep_nuc->RegisterMe(new G4ElectroVDNuclearModel);
          pmanager->AddDiscreteProcess(ep_nuc);
        }
      else if (particleName == "mu-" && !fsOmitMuonic )
        {
          G4MuonNuclearProcess* mu_nuc = new G4MuonNuclearProcess();
          mu_nuc->RegisterMe( new G4MuonVDNuclearModel );
          pmanager->AddDiscreteProcess( mu_nuc );
        }
      else if (particleName == "mu+" && !fsOmitMuonic )
        {
          G4MuonNuclearProcess* mu_nuc = new G4MuonNuclearProcess();
          mu_nuc->RegisterMe( new G4MuonVDNuclearModel );
          pmanager->AddDiscreteProcess( mu_nuc );
        }
    }
}

////////////////////////////////////////////////////////////////////////
/// AddDecay - specific to RAT
///
/// Allows particles to decay e.g. in flight
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::AddDecay()
{
  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
  G4Decay* decayProcess = new G4Decay();
  theParticleIterator->reset();
  G4ParticleDefinition* particle=0;

  while( (*theParticleIterator)() )
    {
      particle = theParticleIterator->value();
      if( decayProcess->IsApplicable(*particle) )
        {
          ph->RegisterProcess(decayProcess, particle);
        }
    }

  G4RadioactiveDecay* theRadioactiveDecay = new G4RadioactiveDecay();
  G4GenericIon* ion = G4GenericIon::GenericIon();
  theRadioactiveDecay->SetAnalogueMonteCarlo(true); //default is true
  theRadioactiveDecay->SetFBeta(true); //fast beta mode, default is false

  G4ProcessManager* pmanager = ion->GetProcessManager();
  pmanager->AddProcess(theRadioactiveDecay);
  pmanager->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
  pmanager->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
}

////////////////////////////////////////////////////////////////////////
/// AddAttenuation - specific to RAT
///
/// Add the specialist attenuation and cerenkov processes
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::AddAttenuationAndCerenkov()
{
  G4double alphaMass = G4Alpha::Alpha()->GetPDGMass();
  G4double neutronMass = G4Neutron::Neutron()->GetPDGMass();

  GLG4Scint* theDefaultScintProcess = new GLG4Scint();
  GLG4Scint* theNeutronScintProcess = new GLG4Scint("neutron", 0.9*neutronMass); //The 0.9*neutronMass is the lower mass limit
  GLG4Scint* theAlphaScintProcess = new GLG4Scint("alpha", 0.9*alphaMass);

  theDefaultScintProcess->SetVerboseLevel( verboseLevel );
  theAlphaScintProcess->SetVerboseLevel( verboseLevel );
  theNeutronScintProcess->SetVerboseLevel( verboseLevel );

  // Boundary process deals with optical photons crossing media boundaries, does it reflect refract etc...


  G4ProcessManager* pmanager = G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
  // NOTE: GLG4Scint is special, it is not a true G4 physics process so do not add it here - PGJ
  if( !fOmitOpticalAbsorption && !fOmitOpticalAttenuation )
    pmanager->AddDiscreteProcess( new G4OpAbsorption() );
  else
    info << "PhysicsList::AddAttenuationAndCerenkov: Omitting optical absorption" << newline;
  //  OpRayleigh implements Rayleigh scattering.
  //  G4OpRayleigh is not used for the following two reasons:
  //    1) It doesn't even try to work for anything other than water.
  //    2) It doesn't actually work for water, either.
  if( !fOmitOpticalRayleigh && !fOmitOpticalAttenuation )
    {
      if( fOpRayleighName == "SNOMAN" )
        {
          info << "PhysicsList::AddAttenuationAndCerenkov: Using the SNOMAN Rayleigh model." << newline;
          pmanager->AddDiscreteProcess( new SNOMANOpRayleigh() );
        }
      else
        pmanager->AddDiscreteProcess( new OpRayleigh() );
    }
  else
    info << "PhysicsList::AddAttenuationAndCerenkov: Omitting optical Rayleigh scattering" << newline;
  if( !fOmitOpticalMie && !fOmitOpticalAttenuation )
    pmanager->AddDiscreteProcess( new G4OpMieHG() );
  else
    info << "PhysicsList::AddAttenuationAndCerenkov Omitting optical Mie scattering" << newline;
  if( !fOmitOpticalBoundaryEffects )
    pmanager->AddDiscreteProcess( new G4OpBoundaryProcess() );
  else
    info << "PhysicsList::AddAttenuationAndCerenkov: Omitting crossing media boundary processes (reflection/refraction)" << newline;
  if( !fOmitCerenkov )
    {
      Cerenkov* theCerenkovProcess;
      if( fCerenkovName == "SNOMAN" )
        {
          info << "PhysicsList::AddAttenuationAndCerenkov: Using the SNOMAN Cerenkov model." << newline;
          theCerenkovProcess = new SNOMANCerenkov("Cerenkov");
        }
      else
        theCerenkovProcess = new Cerenkov("Cerenkov");

      theCerenkovProcess->BuildThePhysicsTable();
      theCerenkovProcess->SetTrackSecondariesFirst(true);
      theCerenkovProcess->SetMaxNumPhotonsPerStep(4);

      // Add the Cerenkov process to all applicable particles (charged)
      theParticleIterator->reset();
      while((*theParticleIterator)())
      {
        G4ParticleDefinition* particle = theParticleIterator->value();
        G4ProcessManager* pmanager = particle->GetProcessManager();

        if( theCerenkovProcess->IsApplicable(*particle) )
        {
          pmanager->AddProcess(theCerenkovProcess);
          pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
        }
      }
    }
  else
    info << "PhysicsList: Cerenkov processes omitted" << newline;
}
////////////////////////////////////////////////////////////////////////
/// SetCuts
///
/// Set the production cuts, with special values for the scint region
/// PHIL Does the Scint region need to be special, what about bg sim?
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::SetCuts()
{
  G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(100.*eV, 100.*GeV);

  SetCutValue( fGammaCut, "gamma");
  SetCutValue( fElectronCut, "e-");
  SetCutValue( fPositronCut, "e+");
  SetCutValue( fProtonCut, "proton");
  SetCutValue( fMuonCut, "mu-" );
  SetCutValue( fMuonCut, "mu+" );

  if( verboseLevel > 0 )
    DumpCutValuesTable();

  // Now the special track killer if needed
  if( fTimeCut > 0.0 )
    {
      theParticleIterator->reset();
      TrackKiller* theTrackKiller = new TrackKiller();
      while( (*theParticleIterator)() )
        {
          G4ParticleDefinition* particle = theParticleIterator->value();
          G4ProcessManager* pmanager = particle->GetProcessManager();
          pmanager->AddDiscreteProcess( theTrackKiller );
        }
    }
}

////////////////////////////////////////////////////////////////////////
/// SetCut - specific to RAT
///
/// Allows a cut be selected from the messenger
////////////////////////////////////////////////////////////////////////

void RAT::PhysicsList::SetCut(const G4String& cutName, const double value)
{
  if( cutName == "Time" )
    fTimeCut = value;
  else if( cutName == "Gamma" )
    fGammaCut = value;
  else if( cutName == "Electron" )
    {
      fElectronCut = value;
      fPositronCut = value;
    }
  else if( cutName == "Muon" )
    fMuonCut = value;
  else if( cutName == "All" )
    {
      fGammaCut = value;
      fElectronCut = value;
      fPositronCut = value;
      fProtonCut = value;
      fMuonCut = value;
    }
}

////////////////////////////////////////////////////////////////////////
/// Omit Muonic/Hadronic Processes  - specific to RAT
///
/// The omits Muonic/Hadronic Processes
////////////////////////////////////////////////////////////////////////
void RAT::PhysicsList::OmitProcesses( G4String processName, bool omit )
{
  if( processName == "all" )
    fOmitAll = omit;
  else if( processName == "muonic")
    fsOmitMuonic = omit;
  else if(processName == "hadronic")
    fsOmitHadronic = omit;
  else if (processName == "hadronicListOnly")
    fOmitHadronicListOnly = omit;
  else if (processName == "cerenkov")
    fOmitCerenkov = omit;
  else if( processName == "attenuation" )
    fOmitOpticalAttenuation = omit;
  else if( processName == "absorption" )
    fOmitOpticalAbsorption = omit;
  else if( processName == "rayleigh" )
    fOmitOpticalRayleigh = omit;
  else if( processName == "mie" )
    fOmitOpticalMie = omit;
  else if( processName == "boundary" )
    fOmitOpticalBoundaryEffects = omit;
  // This will automatically keep the processes on unless specified otherwise
  else
    RAT::Log::Die( "PhysicsList::OmitProcesses: ERROR > Invalid Command " + processName );
}