// --- Copied from AuNP.
// --- Conditions:
//              - there is always a nanoparticle in the geometry.
//              - the nanoparticle is assigned to a region called "NoDNA"
//              - NoDNA region is assigned LivermoreEMPhysics
//              - No chemistry is active


// -------------------------------------------------------------------------------------------------------|
// -------------------------------------------------------------------------------------------------------|
//                              This version of PhysicsList.cc should NOT crash                           |
// -------------------------------------------------------------------------------------------------------|
// -------------------------------------------------------------------------------------------------------|


/// \file PhysicsList.cc
/// \brief Implementation of the PhysicsList class

#include "PhysicsList.hh"

#include "DetectorConstruction.hh"

#include "G4PhysicsConstructorRegistry.hh"
#include "G4RunManager.hh"
#include "G4SystemOfUnits.hh"

#include "G4ComptonScattering.hh"
#include "G4DNAAttachment.hh"
#include "G4DNABornExcitationModel.hh"
#include "G4DNABornIonisationModel.hh"
#include "G4DNAChampionElasticModel.hh"
#include "G4DNAChemistryManager.hh"
#include "G4DNADingfelderChargeDecreaseModel.hh"
#include "G4DNADiracRMatrixExcitationModel.hh"
#include "G4DNAELSEPAElasticModel.hh"
#include "G4DNAElastic.hh"
#include "G4DNAExcitation.hh"
#include "G4DNAGenericIonsManager.hh"
#include "G4DNAIonisation.hh"
#include "G4DNAMeltonAttachmentModel.hh"
#include "G4DNAPlasmonExcitation.hh"
#include "G4DNAQuinnPlasmonExcitationModel.hh"
#include "G4DNARelativisticIonisationModel.hh"
#include "G4DNASancheExcitationModel.hh"
#include "G4DNAVibExcitation.hh"
#include "G4DummyModel.hh"
#include "G4EmConfigurator.hh"
#include "G4Gamma.hh"
#include "G4GammaConversion.hh"
#include "G4GoudsmitSaundersonMscModel.hh"
#include "G4IonFluctuations.hh"
#include "G4LivermoreBremsstrahlungModel.hh"
#include "G4LivermoreComptonModel.hh"
#include "G4LivermoreGammaConversionModel.hh"
#include "G4LivermoreIonisationModel.hh"
#include "G4LivermorePhotoElectricModel.hh"
#include "G4LossTableManager.hh"
#include "G4PenelopeBremsstrahlungModel.hh"
#include "G4PenelopeIonisationModel.hh"
#include "G4PhotoElectricEffect.hh"
#include "G4RayleighScattering.hh"
#include "G4SeltzerBergerModel.hh"
#include "G4StepLimiter.hh"
#include "G4UAtomicDeexcitation.hh"
#include "G4UniversalFluctuation.hh"
#include "G4UrbanMscModel.hh"
#include "G4VEmModel.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eIonisation.hh"
#include "G4eMultipleScattering.hh"

#include "G4hIonisation.hh"
#include "G4hMultipleScattering.hh"
#include "G4DNAIonElasticModel.hh"
#include "G4DNAMillerGreenExcitationModel.hh"
#include "G4DNABornExcitationModel.hh"
#include "G4DNARuddIonisationExtendedModel.hh"
#include "G4DNABornIonisationModel.hh"
#include "G4DNADingfelderChargeDecreaseModel.hh"
#include "G4BraggModel.hh"
#include "G4BetheBlochModel.hh"
#include "G4DNAChargeDecrease.hh"

#include "G4Neutron.hh"
#include "G4HadronElasticPhysicsHP.hh"
#include "G4HadronPhysicsQGSP_BIC_HP.hh"
#include "G4StoppingPhysics.hh"

#include "G4ProcessTable.hh"
#include "G4EmDNAChemistry_option3.hh"
#include "G4GenericIon.hh"


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

PhysicsList::PhysicsList() {
  fpDetector = dynamic_cast<const DetectorConstruction *>(
    G4RunManager::GetRunManager()->GetUserDetectorConstruction());
  defaultCutValue = 0.1 * CLHEP::nanometer;
  SetVerboseLevel(1);
  fEmDNAChemistryList = std::make_unique<G4EmDNAChemistry_option3>();
  G4EmParameters *param = G4EmParameters::Instance();
  param->SetMinEnergy(100 * eV);
  param->SetMaxEnergy(1 * GeV);

  ConstructHadronic(); // --- For neutron transport

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructParticle() {
    G4Gamma::GammaDefinition();

    G4Neutron::NeutronDefinition();

    G4Electron::ElectronDefinition();
    G4Positron::PositronDefinition();

    G4Proton::ProtonDefinition();
    G4GenericIon::GenericIonDefinition();

    auto *genericIonsManager = G4DNAGenericIonsManager::Instance();
    genericIonsManager->GetIon("alpha++");
    genericIonsManager->GetIon("alpha+");
    genericIonsManager->GetIon("helium");
    genericIonsManager->GetIon("hydrogen");

    if (fEmDNAChemistryList != nullptr) {
    fEmDNAChemistryList->ConstructParticle();
  }

    G4VModularPhysicsList::ConstructParticle();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructProcess() {

  ConstructEM();
  ConstructGeneral();

  if (fEmDNAChemistryList != nullptr) {
    fEmDNAChemistryList->ConstructProcess();
  }
  // Contruct processes of the chemistry list
  G4VModularPhysicsList::ConstructProcess();

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructEM() {
  


    ConstructElectronPhysics();

    ConstructProtonPhysics();

    ConstructGammaPhysics();

    ConstructAtomicPhysics();

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructHadronic() {
    // ------------------------------------------------------------------------
    // Hadronic physics
    // ------------------------------------------------------------------------
    RegisterPhysics( new G4HadronElasticPhysicsHP(verboseLevel) );

    RegisterPhysics( new G4HadronPhysicsQGSP_BIC_HP(verboseLevel) );

    RegisterPhysics( new G4StoppingPhysics(verboseLevel) );
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructGammaPhysics() {


  G4ProcessManager *pm = G4Gamma::Gamma()->GetProcessManager();
  G4PhotoElectricEffect *thePhotoElectricEffect = new G4PhotoElectricEffect();
  thePhotoElectricEffect->SetEmModel(new G4LivermorePhotoElectricModel(), 1);
  pm->AddDiscreteProcess(thePhotoElectricEffect);

  G4ComptonScattering *theComptonScattering = new G4ComptonScattering();
  theComptonScattering->SetEmModel(new G4LivermoreComptonModel(), 1);
  pm->AddDiscreteProcess(theComptonScattering);

  G4GammaConversion *theGammaConversion = new G4GammaConversion();
  theGammaConversion->SetEmModel(new G4LivermoreGammaConversionModel(), 1);
  pm->AddDiscreteProcess(theGammaConversion);

  G4RayleighScattering *theRayleigh = new G4RayleighScattering();
  pm->AddDiscreteProcess(theRayleigh);

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructElectronPhysics() {
  // --- ELECTRON PHYSICS --- 
  
  G4EmConfigurator *em_config = G4LossTableManager::Instance()->EmConfigurator();
  
  G4VEmModel *mod;

  G4ProcessManager *pm = G4Electron::Electron()->GetProcessManager();
  
  // === G4DNAElastic ===
  // --- --- Global registration --- ---
  const auto theDNAElasticProcess = new G4DNAElastic("e-_G4DNAElastic");
  theDNAElasticProcess->SetEmModel(new G4DummyModel(), 1); 
  pm->AddDiscreteProcess(theDNAElasticProcess);

  // --- --- Model activation in DNA region --- --- 
  mod = new G4DNAChampionElasticModel(); 
  em_config->SetExtraEmModel("e-", "e-_G4DNAElastic", mod, "DNA");

  // === G4DNAExcitation ===
  // --- --- Global registration --- ---
  const auto theDNAExcitationProcess = new G4DNAExcitation("e-_G4DNAExcitation");
  theDNAExcitationProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAExcitationProcess);
  
  // --- --- Model activation in DNA region --- --- 
  mod = new G4DNABornExcitationModel(); 
  em_config->SetExtraEmModel("e-", "e-_G4DNAExcitation", mod, "DNA");

  // === G4DNAIonisation ===
  // --- --- Global registration --- ---
  const auto theDNAIonizationProcess = new G4DNAIonisation("e-_G4DNAIonisation");
  theDNAIonizationProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAIonizationProcess);

  // --- --- Model activation in DNA region --- ---
  mod = new G4DNABornIonisationModel();
  em_config->SetExtraEmModel("e-", "e-_G4DNAIonisation", mod, "DNA");


  // === G4DNAAttachment ===
  // --- --- Global registration --- ---
  const auto theDNAAttachmentProcess = new G4DNAAttachment("e-_G4DNAAttachment");
  theDNAAttachmentProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAAttachmentProcess);

  // --- --- Model activation in DNA region --- ---
  mod = new G4DNAMeltonAttachmentModel();
  em_config->SetExtraEmModel("e-", "e-_G4DNAAttachment", mod, "DNA");

  // --- G4DNAVibExcitation ---
  // --- --- Global registration --- ---
  const auto theDNAVibExcProcess = new G4DNAVibExcitation("e-_G4DNAVibExcitation");
  theDNAVibExcProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAVibExcProcess);

  // --- --- Model activation in DNA region --- ---
  mod = new G4DNASancheExcitationModel();
  em_config->SetExtraEmModel("e-", "e-_G4DNAVibExcitation", mod, "DNA");

  // --- G4DNABremsstrahlung ---
  // --- --- Global registration --- ---
  const auto theDNABremProcess = new G4eBremsstrahlung("e-_G4DNABremsstrahlung");
  theDNABremProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNABremProcess);

  // --- --- Model activation in DNA region --- ---
  mod = new G4SeltzerBergerModel();
  em_config->SetExtraEmModel("e-", "e-_G4DNABremsstrahlung", mod, "DNA");

  // === === === === === BEGIN LIVERMORE PHYSICS MODELS === === === === ===

  // === G4eMultipleScattering ===
  // --- --- Global dummy model registration --- ---
  const auto msc = new G4eMultipleScattering();
  msc ->SetEmModel(new G4UrbanMscModel(), 1); 
  pm->AddDiscreteProcess(msc);

  // --- --- Model inactivated in DNA region --- ---      
  mod = new G4UrbanMscModel();
  // mod = new G4GoudsmitSaundersonMscModel(); // --- alternate model
  mod->SetActivationLowEnergyLimit(1 * GeV);
  em_config->SetExtraEmModel("e-", "msc", mod, "DNA");

  // === G4eIonisation ===
  // --- --- Global dummy model registration --- ---
  const auto eIoni = new G4eIonisation("e-_G4LivermoreIoni");
  eIoni->SetEmModel(new G4LivermoreIonisationModel(), 1);
  pm->AddDiscreteProcess(eIoni);

  // --- --- Model inactivated in DNA region --- ---
  mod = new G4LivermoreIonisationModel();
  mod->SetActivationLowEnergyLimit(1 * GeV);
  em_config->SetExtraEmModel("e-", "e-_G4LivermoreIoni", mod, "DNA");

  // === G4eBremsstrahlung ===
  // --- --- Global dummy model registration --- ---
  const auto eBrem = new G4eBremsstrahlung("e-_G4LivermoreBrem");
  eBrem->SetEmModel(new G4LivermoreBremsstrahlungModel(), 1);
  pm->AddDiscreteProcess(eBrem);

  // --- --- Model inactivated in DNA region --- ---
  mod = new G4LivermoreBremsstrahlungModel();
  mod->SetActivationLowEnergyLimit(1 * GeV);
  em_config->SetExtraEmModel("e-", "e-_G4LivermoreBrem", mod, "DNA");

  // --- G4StepLimiter ---
  const auto steplimit = new G4StepLimiter();
  pm->AddDiscreteProcess(steplimit);
    


  //// Capture of low-energy e-
  // G4ElectronCapture* ecap = new G4ElectronCapture("NoDNA",10.*eV);
  // pm->AddDiscreteProcess(ecap); 
}
    
    
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructProtonPhysics() {

  // --- PROTON PHYSICS ---

  G4EmConfigurator *em_config = G4LossTableManager::Instance()->EmConfigurator();
  
  G4VEmModel *mod;

  G4ProcessManager *pm = G4Proton::Proton()->GetProcessManager();
  // === G4DNAElastic ===
  // --- --- Global registration --- ---
  const auto theDNAElasticProcess = new G4DNAElastic("proton_G4DNAElastic");
  theDNAElasticProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAElasticProcess);

  // --- --- Model activation in DNA region --- ---
  mod = new G4DNAIonElasticModel();
  em_config->SetExtraEmModel("proton", "proton_G4DNAElastic", mod, "DNA");

  // === G4DNAExcitation === 
  // --- --- Global registration --- ---
  const auto theDNAExcitationProcess = new G4DNAExcitation("proton_G4DNAExcitation");
  theDNAExcitationProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAExcitationProcess);

  // --- --- Model activation in DNA region --- --- 
  mod = new G4DNAMillerGreenExcitationModel();
  em_config->SetExtraEmModel("proton","proton_G4DNAExcitation", mod, "DNA");
  mod = new G4DNABornExcitationModel();
  em_config->SetExtraEmModel("proton","proton_G4DNAExcitation", mod, "DNA");

  // === G4DNAIonisation ===
  // --- --- Global registration --- ---
  const auto theDNAIonizationProcess = new G4DNAIonisation("proton_G4DNAIonisation");
  theDNAIonizationProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAIonizationProcess);


  // --- --- Model activation in DNA region --- ---
  mod = new G4DNARuddIonisationExtendedModel();
  em_config->SetExtraEmModel("proton","proton_G4DNAIonisation", mod, "DNA");
  mod = new G4DNABornIonisationModel();
  em_config->SetExtraEmModel("proton","proton_G4DNAIonisation", mod, "DNA");

  // === G4DNAChargeDecrease ===
  // --- --- Global registration --- ---
  const auto theDNAChargeDecreaseProcess = new G4DNAChargeDecrease("proton_G4DNAChargeDecrease");
  theDNAChargeDecreaseProcess->SetEmModel(new G4DummyModel(), 1);
  pm->AddDiscreteProcess(theDNAChargeDecreaseProcess);

  // --- --- Model activation in DNA region --- ---
  mod = new G4DNADingfelderChargeDecreaseModel();
  em_config->SetExtraEmModel("proton", "proton_G4DNAChargeDecrease", mod, "DNA");

  // === === === === === Begin Livermore Physics Models === === === === ===

  // === G4hMultipleScattering ===
  // --- --- Global dummy model registration --- ---
  const auto msc = new G4hMultipleScattering();
  msc->SetEmModel(new G4UrbanMscModel(), 1);
  pm->AddDiscreteProcess(msc);

  // --- --- Model inactivation in DNA region --- ---
  mod = new G4UrbanMscModel();
  mod->SetActivationLowEnergyLimit(1 * GeV);
  em_config->SetExtraEmModel("proton", "msc", mod, "DNA");

  // === G4hIonisation ===
  
  // --- --- Global dummy model registration --- ---
  const auto hIoni = new G4hIonisation("proton_G4LivermoreIoni");
  hIoni->SetEmModel(new G4BraggModel(), 1);
  hIoni->SetEmModel(new G4BetheBlochModel(), 2);
  pm->AddDiscreteProcess(hIoni);

    // --- --- Model inactivation in DNA region --- ---
  mod = new G4BraggModel();
  mod->SetActivationLowEnergyLimit(1 * GeV);
  em_config->SetExtraEmModel("proton", "proton_G4LivermoreIoni", mod, "DNA");
  mod = new G4BetheBlochModel();
  mod->SetActivationLowEnergyLimit(1 * GeV);
  em_config->SetExtraEmModel("proton", "proton_G4LivermoreIoni", mod, "DNA");

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void PhysicsList::ConstructAtomicPhysics() {

  G4bool atomicDeexcitationActive = true;

  G4VAtomDeexcitation *de = new G4UAtomicDeexcitation();
  G4LossTableManager::Instance()->SetAtomDeexcitation(de);
  de->SetFluo(atomicDeexcitationActive);
  de->SetPIXE(atomicDeexcitationActive);
  de->SetAuger(atomicDeexcitationActive);
  de->SetAugerCascade(atomicDeexcitationActive);

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....


void PhysicsList::SetCuts() {
  if (verboseLevel > 0) {
    G4cout << "PhysicsList::SetCuts:";
    G4cout << "CutLength : " << G4BestUnit(defaultCutValue, "Length") << G4endl;
  }

  SetCutValue(fcutForGamma, "gamma");
  SetCutValue(fcutForElectron, "e-");
  SetCutValue(fcutForPositron, "e+");

  if (verboseLevel > 0) {
    DumpCutValuesTable();
  }
}