Geant4 Version: v11.3.2 (v11 onwards)
Operating System: Linux Alma9
Compiler/Version: g++ (GCC) 14.2.0
CMake Version: cmake version 3.30.6
Hi,
I’m hoping to get some help with changes to the production of secondaries in v11 of geant4.
The changes in BraggModel from v10 to v11 mean that the production of delta rays is suppressed for alphas below ~2.0 MeV, leading to a two orders of magnitude reduction in the number of secondary electrons produced compared to v10 (specifically comparing to geant4-10-07).
I’m using QGSP_BIC_HP and G4EMPenelopePhysics in my physics lists.
This seems to come from changes to G4BraggModel::ComputeCrossSectionPerElectron which calculates the energy cut as const G4double cutEnergy = std::max(cut, lowestKinEnergy*massRate), where massRate is the ratio of the alpha to proton mass (3.9726), scaling up the lowest kinetic energy and giving a cut floor of 397eV, rather than the defined production cut. This means below ~1.8 MeV no discrete delta-rays are produced by the alpha track.
This can be reproduced in a basic simulation of an alpha in a gaseous nitrogen target with:
void RunAction::BeginOfRunAction(const G4Run* run) {
G4EmCalculator calc;
G4Material* mat = G4Material::GetMaterial("N2");
G4cout << "MFP of delta rays produced by alphas in N2" << G4endl;
std::vector<G4double> energies = { 0.1*MeV, 0.5*MeV, 1.0*MeV, 2.0*MeV, 5.0*MeV};
for (G4double e : energies) {
G4double mfp = calc.GetMeanFreePath(e, "alpha", "ionIoni", mat->GetName());
G4cout << " E = " <<e/MeV << " MeV MFP = " << mfp/mm << " mm" << G4endl;
}
}
This returns:
E = 0.1 MeV MFP = 1.79769e+308 mm
E = 0.5 MeV MFP = 1.79769e+308 mm
E = 1 MeV MFP = 1.79769e+308 mm
E = 2 MeV MFP = 0.354217 mm
E = 5 MeV MFP = 0.116996 mm
Where low energy alphas produce delta rays with infinite mean free path.
I traced this back to the energy cut by adding this in my SteppingAction:
void SteppingAction::UserSteppingAction(const G4Step* step) {
G4Track *track = step->GetTrack();
static G4bool printed = false;
if (!printed && track->GetParticleDefinition()->GetParticleName() == "alpha") {
printed = true;
G4ProcessManager* pm = G4Alpha::Alpha()->GetProcessManager();
for (G4int i = 0; i < pm->GetProcessList()->size(); i++) {
G4VProcess* proc = (*pm->GetProcessList())[i];
if (proc->GetProcessName() == "ionIoni") {
G4VEnergyLossProcess* ioni = dynamic_cast<G4VEnergyLossProcess*>(proc);
G4BraggIonModel* bim = dynamic_cast<G4BraggIonModel*>(ioni->GetModelByIndex(0));
if (bim) {
G4Material* mat = track->GetMaterial();
G4double KE = 1.0*MeV;
G4double mass = G4Alpha::Alpha()->GetPDGMass();
G4double ratio = CLHEP::electron_mass_c2 / mass;
G4double q = G4Alpha::Alpha()->GetPDGCharge()/CLHEP::eplus;
G4double chargeSq = q * q;
G4double lowestKE = bim->LowEnergyLimit();
G4double massRate = mass / CLHEP::proton_mass_c2;
G4double tau = KE / mass;
G4double tmax = 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.0) /
(1.0 + 2.0*(tau + 1.0)*ratio + ratio*ratio);
G4double cutEnergy = std::max(1*eV, lowestKE * massRate);
G4double maxEnergy = std::min(tmax, 1*MeV);
G4double beta2 = KE*(KE + 2.0*mass) / ((KE+mass)*(KE+mass));
G4cout << "----- ComputeCrossSectionPerElectron -----" << G4endl;
G4cout << "mass = " << mass/MeV << " MeV" << G4endl;
G4cout << "ratio = " << ratio << G4endl;
G4cout << "chargeSquare = " << chargeSq << G4endl;
G4cout << "massRate = " << massRate << G4endl;
G4cout << "lowestKE = " << lowestKE/MeV << " MeV" << G4endl;
G4cout << "tau = " << tau << G4endl;
G4cout << "tmax = " << tmax/eV << " eV" << G4endl;
G4cout << "cutEnergy = " << cutEnergy/eV << " eV" << G4endl;
G4cout << "maxEnergy = " << maxEnergy/eV << " eV" << G4endl;
G4cout << "beta2 = " << beta2 << G4endl;
if (cutEnergy < maxEnergy) {
G4double cross = (maxEnergy - cutEnergy)/(cutEnergy*maxEnergy)
- beta2*std::log(maxEnergy/cutEnergy)/tmax;
cross *= CLHEP::twopi_mc2_rcl2 * chargeSq / beta2;
G4cout << " xsec (manual) = " << cross/cm2 << " cm2" << G4endl;
}
}
}
}
which returns:
mass = 3727.38 MeV
ratio = 0.000137093
chargeSquare = 4
massRate = 3.9726
lowestKE = 0.0001 MeV
tau = 0.000268285
tmax = 548.297 eV
cutEnergy = 397.26 eV
maxEnergy = 548.297 eV
beta2 = 0.000536354
xsec (manual) = 1.31785e-18 cm2
where the energy cut is 397.26 eV rather than the 20 eV I set in my PhysicsList. Calling G4double xsec = bim->ComputeCrossSectionPerAtom(G4Alpha::Alpha(), 1*MeV, 7.0, 0.0, 1*eV, 1*MeV); directly also gives 0 cross-section below 2 MeV.
This is particularly significant in gaseous detector simulations where I need to accurately count the number of secondary electrons produced. Is there something I should be calling differently to track the electrons produced by my low energy alphas?
Thanks!
Lex
