Implementing Single Scattering into custom physics list

Physics Processes, Models and Cross Sections
1

Hello,

I currently have a project where I am investing electron escape yields of tritium decay in spherical dust particles of tungsten. I have made a custom physics list which is based primarily off the MicroElec example provided by Geant4, but I have modified it so that below 7 keV, all processes are modelled by MicroElec and above 7 keV standard processes are used.

Before I was using multiple scattering (msc) above 7 keV, however I wanted to try and implement single scattering using G4eCoulombScattering. However, when I run I see that the process is registered but I do not see any scattering at all. I have tested multiple energies, diameters of dust etc but I cant see any scattering, even is ionization may dominate. Therefore I suspect that my implementation is wrong, given this custom list has now become quite messy.

Here is my physics list I am using:

#include “G4SystemOfUnits.hh”

#include “physicsMicroElecNew.hh”

// Geant4-MicroElec MODELS

#include “G4MicroElecElastic.hh”
#include “G4MicroElecElasticModel_new.hh”
#include “G4MicroElecInelastic.hh”
#include “G4MicroElecInelasticModel_new.hh”
#include “G4MicroElecLOPhononModel.hh”
#include “G4MicroElecLOPhononScattering.hh”
#include “G4MicroElecSurface.hh”

//

#include “G4BetheBlochModel.hh”
#include “G4BraggIonModel.hh”
#include “G4BraggModel.hh”
#include “G4ComptonScattering.hh”
#include “G4CoulombScattering.hh”
#include “G4DummyModel.hh”
#include “G4EmConfigurator.hh”
#include “G4EmStandardPhysics_option4.hh”
#include “G4IonFluctuations.hh”
#include “G4LivermoreComptonModel.hh”
#include “G4LivermorePhotoElectricModel.hh”
#include “G4LivermoreRayleighModel.hh”
#include “G4LossTableManager.hh”
#include “G4MicroElecCapture.hh”
#include “G4MollerBhabhaModel.hh”
#include “G4PenelopeBremsstrahlungModel.hh”
#include “G4PenelopeComptonModel.hh”
#include “G4PenelopeGammaConversionModel.hh”
#include “G4PenelopeIonisationModel.hh”
#include “G4PenelopePhotoElectricModel.hh”
#include “G4PenelopeRayleighModel.hh”
#include “G4PhotoElectricEffect.hh”
#include “G4RayleighScattering.hh”
#include “G4SeltzerBergerModel.hh”
#include “G4UAtomicDeexcitation.hh”
#include “G4UniversalFluctuation.hh”
#include “G4UrbanMscModel.hh”
#include “G4VEmModel.hh”
#include “G4eBremsstrahlung.hh”
#include “G4eCoulombScatteringModel.hh”
#include “G4eDPWACoulombScatteringModel.hh”
#include “G4eIonisation.hh”
#include “G4eMultipleScattering.hh”
#include “G4hIonisation.hh”
#include “G4hMultipleScattering.hh”
#include “G4ionIonisation.hh”

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

MicroElecPhysics::MicroElecPhysics() : G4VUserPhysicsList()
{
defaultCutValue = 1 * micrometer;
cutForGamma = defaultCutValue;
cutForElectron = defaultCutValue;
cutForPositron = defaultCutValue;
cutForProton = defaultCutValue;

// G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(5 * eV, 10 * MeV);

SetVerboseLevel(1);
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

MicroElecPhysics::~MicroElecPhysics() {}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructParticle()
{
ConstructBosons();
ConstructLeptons();
ConstructBarions();
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructBosons()
{
// gamma
G4Gamma::GammaDefinition();
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructLeptons()
{
// leptons
G4Electron::ElectronDefinition();
G4Positron::PositronDefinition();
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructBarions()
{
// baryons
G4Proton::ProtonDefinition();
G4GenericIon::GenericIonDefinition();
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructProcess()
{
AddTransportation();
ConstructEM();
ConstructGeneral();
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructEM()
{
G4EmParameters* param = G4EmParameters::Instance();
// param->SetDefaults();
param->SetBuildCSDARange(true);
// param->SetMscStepLimitType(fUseSafetyPlus);
// param->SetMscStepLimitType(fUseDistanceToBoundary);
param->SetMscStepLimitType(fUseSafety);
param->RegionsMicroElec();
// physicList ISS
param->SetDefaults();
param->SetMinEnergy(0.1 * eV);
param->SetMaxEnergy(10 * TeV);
param->SetLowestElectronEnergy(
0 * eV); //<— Energie de cut dans le vide!!! A fixer 0eV pour ne pas fausser les SEY
param->SetNumberOfBinsPerDecade(20);
param->ActivateAngularGeneratorForIonisation(true);
param->SetAugerCascade(true);
// param->SetUseMottCorrection(true); //*/

auto particleIterator = GetParticleIterator();
particleIterator->reset();

while ((particleIterator)()) {
G4ParticleDefinition particle = particleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();

// *********************************
// 1) Processes for the World region
// *********************************

if (particleName == "e-") {
  // STANDARD msc is active in the world
  G4eMultipleScattering* msc = new G4eMultipleScattering();
  msc->AddEmModel(1, new G4UrbanMscModel());
  pmanager->AddProcess(msc, -1, 1, -1);

  // STANDARD ionisation is active in the world
  G4eIonisation* eion = new G4eIonisation();
  pmanager->AddProcess(eion, -1, 2, 2);

  G4eBremsstrahlung* brem = new G4eBremsstrahlung();
  pmanager->AddProcess(brem, -1, 3, 3);

  G4CoulombScattering* cs = new G4CoulombScattering();
  cs->SetEmModel(new G4eCoulombScatteringModel());
  pmanager->AddDiscreteProcess(cs);

  // MicroElec elastic is not active in the world
  // G4MicroElecElasticCorrected* theMicroElecElasticProcess = new
  // G4MicroElecElasticCorrected("e-_G4MicroElecElastic");
  // theMicroElecElasticProcess->SetEmModel(new G4DummyModel(),1);
  // pmanager->AddDiscreteProcess(theMicroElecElasticProcess);

  G4MicroElecElastic* theMicroElecElasticProcess =
    new G4MicroElecElastic("e-_G4MicroElecElastic");
  theMicroElecElasticProcess->SetEmModel(new G4DummyModel(), 1);
  // G4MicroElecElasticModel_new* mod = new G4MicroElecElasticModel_new();
  // theMicroElecElasticProcess->AddEmModel(0,mod);
  pmanager->AddDiscreteProcess(theMicroElecElasticProcess);

  // MicroElec ionisation is not active in the world
  /*G4MicroElecInelastic* microelecioni = new G4MicroElecInelastic("e-_G4MicroElecInelastic");
  microelecioni->SetEmModel(new G4DummyModel(),1);
  pmanager->AddDiscreteProcess(microelecioni);*/

  G4MicroElecInelastic* microelecioni = new G4MicroElecInelastic("e-_G4Dielectrics");
  microelecioni->SetEmModel(new G4DummyModel(), 1);
  pmanager->AddDiscreteProcess(microelecioni);

  // Phonons for SiO2

  G4MicroElecLOPhononScattering* opticalPhonon =
    new G4MicroElecLOPhononScattering("e-_G4LOPhononScattering");
  opticalPhonon->SetEmModel(new G4DummyModel(), 1);
  pmanager->AddDiscreteProcess(opticalPhonon);

  /*G4LOPhononScattering* LO60 = new G4LOPhononScattering("e-_G4LO60");
  LO60->SetEmModel(new G4DummyModel(), 1);
  pmanager->AddDiscreteProcess(LO60);//*/

  G4MicroElecSurface* MicroElecSurf = new G4MicroElecSurface("e-_G4MicroElecSurface");
  MicroElecSurf->SetProcessManager(pmanager);
  pmanager->AddDiscreteProcess(MicroElecSurf);  //*/

  G4MicroElecCapture* ecap =
    new G4MicroElecCapture("Target", 0.9 * eV);  //<--- Piges pour Al2O3
  pmanager->AddDiscreteProcess(ecap);  //*/
}
else if (particleName == "proton") {
  // STANDARD msc is active in the world
  /*G4hMultipleScattering* msc = new G4hMultipleScattering();
  msc->AddEmModel(1, new G4UrbanMscModel());
  pmanager->AddProcess(msc, -1, 1, -1);*/

  // STANDARD ionisation is active in the world
  G4hIonisation* hion = new G4hIonisation();
  pmanager->AddProcess(hion, -1, 2, 2);

  //// Dielectric ionisation is not active in the world
  // G4MicroElecInelastic* dielectricioni = new G4MicroElecInelastic("p_G4Dielectrics");
  // dielectricioni->SetEmModel(new G4DummyModel(), 1);
  //// dielectricioni->SetEmModel(new G4DummyModel(),2);
  // pmanager->AddDiscreteProcess(dielectricioni);

  G4MicroElecCapture* pcap = new G4MicroElecCapture("Target", 100.1 * eV);
  pmanager->AddDiscreteProcess(pcap);  //*/
}
else if (particleName == "alpha") {
  // STANDARD ionisation is active in the world
  G4ionIonisation* hion = new G4ionIonisation();
  pmanager->AddProcess(hion, -1, 2, 2);
  //// Dielectric ionisation is not active in the world
  // G4MicroElecInelastic* dielectricioni = new G4MicroElecInelastic("alpha_G4Dielectrics");
  // dielectricioni->SetEmModel(new G4DummyModel(), 1);
  // dielectricioni->SetEmModel(new G4DummyModel(), 2);
  // pmanager->AddDiscreteProcess(dielectricioni);
}
else if (particleName == "GenericIon") {
  // STANDARD msc is active in the world
  /*G4hMultipleScattering* msc = new G4hMultipleScattering();
  msc->AddEmModel(1, new G4UrbanMscModel());
  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, -1);*/

  /*G4CoulombScattering* cs = new G4CoulombScattering();
  cs->AddEmModel(0, new G4IonCoulombScatteringModel());
  cs->SetBuildTableFlag(false);
  pmanager->AddDiscreteProcess(cs);*/

  // STANDARD ionisation is active in the world
  G4ionIonisation* hion = new G4ionIonisation();
  pmanager->AddProcess(hion, -1, 2, 2);

  //// Dielectric ionisation is not actived in the world
  // G4MicroElecInelastic* dielectricioni = new G4MicroElecInelastic("ion_G4Dielectrics");
  // dielectricioni->SetEmModel(new G4DummyModel(), 1);
  // dielectricioni->SetEmModel(new G4DummyModel(), 2);
  // pmanager->AddDiscreteProcess(dielectricioni);
}

}

// **************************************
// 2) Define processes for Target region
// **************************************

// STANDARD EM processes should be inactivated when corresponding MicroElec processes are used
// - STANDARD EM e- processes are inactivated below 100 MeV
// - STANDARD EM proton & ion processes are inactivated below standEnergyLimit
//
G4EmConfigurator* em_config = G4LossTableManager::Instance()->EmConfigurator();

G4VEmModel* mod;
// *** e-

// —> STANDARD EM processes are inactivated below 100 MeV

G4DummyModel* dummyMsc = new G4DummyModel();
em_config->SetExtraEmModel(“e-”, “msc”, dummyMsc, “Target”, 0 * eV, 10 * TeV);

//// Then set up single scattering
////G4CoulombScattering* ss = new G4CoulombScattering();
// G4VEmModel* coulombModel = nullptr;
// if (param->UseMottCorrection()) {
// coulombModel = new G4eDPWACoulombScatteringModel();
// }
// else {
// coulombModel = new G4eCoulombScatteringModel();
// }
// coulombModel->SetActivationLowEnergyLimit(7 * keV);
////ss->SetEmModel(coulombModel);
// em_config->SetExtraEmModel(“e-”, “CoulombScat”, coulombModel, “Target”, 7 * keV, 100 * MeV);

// G4UrbanMscModel* msc = new G4UrbanMscModel();
// msc->SetActivationLowEnergyLimit(7 * keV);
// em_config->SetExtraEmModel(“e-”, “msc”, msc, “Target”, 7 *keV, 100 * MeV);

// mod = new G4MollerBhabhaModel();
// mod->SetActivationLowEnergyLimit(10 * keV);
// em_config->SetExtraEmModel(“e-”, “eIoni”, mod, “Target”, 10 * keV, 10 * TeV,
// new G4UniversalFluctuation());

mod = new G4PenelopeIonisationModel();
mod->SetActivationLowEnergyLimit(7 * keV);
mod->SetHighEnergyLimit(1 * MeV);
em_config->SetExtraEmModel(“e-”, “eIoni”, mod, “Target”, 7 * keV, 1 * MeV,
new G4UniversalFluctuation());

mod = new G4eCoulombScatteringModel();
mod->SetActivationLowEnergyLimit(7 * keV);
em_config->SetExtraEmModel(“e-”, “CoulombScat”, mod, “Target”, 7 * keV, 1 * TeV);

// brem
mod = new G4SeltzerBergerModel();
mod->SetActivationLowEnergyLimit(7 * keV);
mod->SetHighEnergyLimit(1 * GeV);
em_config->SetExtraEmModel(“e-”, “eBrem”, mod, “Target”, 7 * keV, 1 * GeV);

// —> MicroElec processes activated

mod = new G4MicroElecElasticModel_new();
em_config->SetExtraEmModel(“e-”, “e-_G4MicroElecElastic”, mod, “Target”, 0.1 * eV, 7 * keV);

mod = new G4MicroElecInelasticModel_new();
em_config->SetExtraEmModel(“e-”, “e-_G4Dielectrics”, mod, “Target”, 0.1 * eV, 7 * keV);

// G4double hw = 0.15*eV;

// Old phonon

/mod = new LOPhononModel(0.153eV,false);
em_config->SetExtraEmModel(“e-”, “e-_G4LOPhononScattering”, mod, “Target”, 0.153*eV, 10 * MeV);

mod = new LOPhononModel(0.063eV,false);
em_config->SetExtraEmModel(“e-”, “e-_G4LO60”, mod, “Target”, 0.06eV, 10 * MeV);//*/

// Phonons LO pour sio2 et al2o3

mod = new G4MicroElecLOPhononModel();
// mod->SetActivationLowEnergyLimit(5 * eV);
em_config->SetExtraEmModel(“e-”, “e-_G4LOPhononScattering”, mod, “Target”, 0.1 * eV,
7 * keV); //*/

// *** proton ----------------------------------------------------------

// —> STANDARD EM processes inactivated below standEnergyLimit

// STANDARD msc is still active
// Inactivate following STANDARD processes

// il faut desactiver Bragg puisque notre modle descend en-dessous de 50 keV
/mod = new G4BraggModel();
mod->SetActivationHighEnergyLimit(50keV);
em_config->SetExtraEmModel(“proton”,“hIoni”,mod,“Target”,0.0,2MeV, new G4IonFluctuations());/

mod = new G4BetheBlochModel();
mod->SetActivationLowEnergyLimit(10 * MeV);
em_config->SetExtraEmModel(“proton”, “hIoni”, mod, “Target”, 2 * MeV, 10 * TeV,
new G4IonFluctuations());

// —> Dielectric processes activated

mod = new G4MicroElecInelasticModel_new();
mod->SetActivationLowEnergyLimit(100 * eV);
em_config->SetExtraEmModel(“proton”, “p_G4Dielectrics”, mod, “Target”, 100 * eV, 7 * keV);
// em_config->SetExtraEmModel(“proton”,“p_G4Dielectrics”,new G4DummyModel,“Target”,10MeV,10TeV);

//*/

// *** alpha ----------------------------------------------------------
mod = new G4BetheBlochModel();
mod->SetActivationLowEnergyLimit(10 * MeV);
em_config->SetExtraEmModel(“alpha”, “ionIoni”, mod, “Target”, 10 * MeV, 10 * TeV,
new G4IonFluctuations());

/mod = new G4MicroElecInelasticModel_new();
//mod->SetActivationLowEnergyLimit(100eV);
em_config->SetExtraEmModel(“alpha”,“alpha_G4Dielectrics”,mod,“Target”,0.0,10MeV);///

// *** ion ----------------------------------------------------------

// —> STANDARD EM processes inactivated below standEnergyLimit

// STANDARD msc is still active
// Inactivate following STANDARD processes

/mod = new G4BraggIonModel();
mod->SetActivationHighEnergyLimit(50keV);
em_config->SetExtraEmModel(“GenericIon”,“ionIoni”,mod,“Target”,0.0,2MeV, new
G4IonFluctuations());/

mod = new G4BetheBlochModel();
mod->SetActivationLowEnergyLimit(10 * MeV);
em_config->SetExtraEmModel(“GenericIon”, “ionIoni”, mod, “Target”, 10 * MeV, 10 * TeV,
new G4IonFluctuations());

// —> Dielectric processes activated
mod = new G4MicroElecInelasticModel_new();
mod->SetActivationLowEnergyLimit(100 * eV);
em_config->SetExtraEmModel(“GenericIon”, “ion_G4Dielectrics”, mod, “Target”, 0.0, 7 * keV);
// em_config->SetExtraEmModel(“GenericIon”,“ion_G4Dielectrics”,new
// G4DummyModel,“Target”,10GeV,10TeV);

// *** gamma ----------------------------------------------------------

G4VEmModel* comptonModel = new G4LivermoreComptonModel();
comptonModel->SetLowEnergyLimit(7 * keV);
em_config->SetExtraEmModel(“gamma”, “compt”, comptonModel, “Target”, 7 * keV, 10 * GeV);

G4VEmModel* photoModel = new G4LivermorePhotoElectricModel();
photoModel->SetLowEnergyLimit(7 * keV);
em_config->SetExtraEmModel(“gamma”, “phot”, photoModel, “Target”, 7 * keV, 10 * GeV);

G4VEmModel* convModel = new G4PenelopeGammaConversionModel();
convModel->SetLowEnergyLimit(7 * keV);
em_config->SetExtraEmModel(“gamma”, “conv”, convModel, “Target”, 7 * keV, 10 * GeV);

G4VEmModel* rayleighModel = new G4LivermoreRayleighModel();
rayleighModel->SetLowEnergyLimit(7 * keV);
rayleighModel->SetHighEnergyLimit(10 * GeV);
em_config->SetExtraEmModel(“gamma”, “Rayleigh”, rayleighModel, “Target”, 7 * keV, 10 * GeV);

// Deexcitation
//
G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
G4LossTableManager::Instance()->SetAtomDeexcitation(de);
de->SetFluo(true);
de->SetAuger(true);
de->SetPIXE(true);
de->InitialiseForNewRun();

// G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(5eV, 100.0GeV);
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

void MicroElecPhysics::ConstructGeneral() {}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

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

// set cut values for gamma at first and for e- second and next for e+,
// because some processes for e+/e- need cut values for gamma
SetCutValue(cutForGamma, “gamma”);
SetCutValue(cutForElectron, “e-”);
SetCutValue(cutForPositron, “e+”);
SetCutValue(cutForProton, “proton”);

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

1 Like