Created attachment 15 [details] Patch Hello! I submit a (3-to-5-fold) bug report (+ possible solution) for a problem reported at the last Geant4 workshop in Tokyo. Performing some validations of the radioactive decay module (i.e. comparing it with the alternative simulation code MGGPOD and theoretically expected values) revealed a significant problem concerning the electron capture decay channel and the subsequent de-excitation of the generated atom: the ratios of the G4-simulated lines are not correct. The problem reveals itself, if one tries to simulate the decay of 67-Ga, which decays via electron capture to the meta-stable 67-Zn. Looking e.g. at the simulated line complex around 300 keV should reveal three line complexes, one around 300.3 keV (M-shell EC), 301.3 keV (L-shell EC), and 309.8 keV (K-shell EC). (If in your simulation code you do not stop the simulation BEFORE the decay of the meta-stable state, you should add 93.3 keV to your observed values). According to the Geant4 radioactive decay data files ("z31.a67"), the expected ratios should be: K:L ~8.8, L:M ~5.6. This K:(L+M) ratio although corresponds roughly with measured values (e.g. on-board TGRS satellite). The simulation with geant4.9.1-p01 however results in e.g. a K:L value around 4:3, with several lines missing or at the wrong position! The cause for the problem can be found in the de-excitation handling in the class G4NuclearDecayChannel.cc - a closer look will reveal at least 3 connected bugs resulting in the wrong line ratios. Problem 1: No local energy deposits handled Comparing the de-excitation handling in G4NuclearDecayChannel::DecayIt() with the same functionality in e.g. G4LowEnergyPhotoEffect reveals that the decay class does not take care of local deposit, e.g. from the unhandled Auger-Electrons (they are unhandled because the Auger data is only loaded for materials defined during detector construction, and thus not for all possible elements which might decay). In the attached patch the de-ecitation handling from the G4LowEnergyPhotoEffect is adapted for the decay class – including loacl energy deposits. Problem 2: Mixing of different "shellIndex" and "shellID" indices In the "case" loop in G4NuclearDecayChannel::DecayIt() a random sub-shell is chosen, ranging from 1 to 9. This is then the input for the function call: atomDeex->GenerateParticles(aZ,eshell); However, GenerateParticles takes not the shell index but the shell ID defined by the underlying EADL data, i.e. K=1; L=3,5,6; M=8,10,11,13,14 This shell ID can be retrieved via: const G4AtomicShell* shell = transitionManager->Shell(aZ, eshell); G4int shellId = shell->ShellId(); In addition the startindex has to be 0 not one: case KshellEC: { eShell = 0; // <-- } break; Code validity question 1: "case" loop in G4NuclearDecayChannel::DecayIt() The "case" loop in G4NuclearDecayChannel::DecayIt() selects a random sub-shell of the atom. This takes into account that the electrons in the sub-shells have different binding energies, but it does not take into account that the shells have different electron occupancies! I did not fix this! Code validity question 2: Again "case" loop in G4NuclearDecayChannel::DecayIt() The "case" loop does not protect against non-existing shells, e.g. for an atom with not completely filled M-shell the random generator should only sample those sub-shells which really exist!! I did not fix this! Problem 3: G4VDecayChannel interface does not allow for propagating local energy deposits The next BIG problem is the propagation of the local energy deposit through the "DecayIt" functions back to the G4RadioactiveDecay class: Currently now such functionality is provided!! This constitutes a severe design issue in the G4VDecayChannel interface which is used by many other classes in Geant4! In order to minimize the changes I decided to change ALL DecayIt functions, i.e. that of the base class and of all derived classes from: virtual G4DecayProducts* DecayIt(G4double parentMass = -1.0) = 0; to: virtual G4DecayProducts* DecayIt(G4double parentMass = -1.0, G4double* localDeposit = 0) = 0; I admit this change does not represent the best programming standards, so if you have a better solution please use it instead The next step is to propagate the local deposits in the RadioactiveDecay class through the "DoDecay"-function and through all calls of "DecayIt". After all those changes, decaying Ga-67 leads exactly to the observed line ratios for the line complex around 300 keV!! Problem fixed :) However I have several additional enhancement requests: Enhancement request 1: It is currently not possible to simulate Auger-Electrons, since the Auger data is only loaded for those elements, which exist during detector construction. In order to generate Auger electrons, the atomic de-excitation class needs a facility, which automatically loads the data if the de-excitation of a non-existing elements is requested! Enhancement request 2: For interface symmetry reasons it would be nice to also allow the user to activate Auger-Electron and Fluorescence via functions such as void ActivateAuger(G4bool); in the class G4LowEnergyPhotoEffect The modified source file are attached as a patch. Simply switch to your $G4INSTALL/source directory, and do: cat <path to>/Decay.patch | patch -p1 Please let me know when you have implemented the changes, I volunteer to test them on my data. Thanks, Andreas Zoglauer
Created attachment 16 [details] Original, wrong line ratios Attached is an image of the original, wrong line ratios generated with geant4.9.1.p01
Created attachment 17 [details] Correct line ratios after patch The attached image shows the correct line ratios after applying the patch
Hi again, Not sure if anybody is listening here, but... The same problem appears to be in G4PhotonEvaporation for decays where an immediate internal conversion from emitted gamma ray to shell electron appears. However, it is significantly more complicated to propagate the local energy deposit out of this class. Thus I have not yet succeeded, since this involved changing > 30 classes, and gave me some nice seg faults... ;-( Thus, this requires a non-negligible design iteration, which I do not want to do by myself without enough time and insight into Geant4... I hope that somebody of you can do it --- at least please let me know how you intend to proceed, since I have to rely on getting the correct decay lines out of Geant4... Thanks for taking care of this problem. Andreas
I am only in cc in this problem report, but I wish to say at least that somebody is listening... I fully share your evaluation that the RadioactiveDecay requires a significant design iteration: I have been saying the same for years. Also PhotonEvaporation requires a significant design investment: I designed it myself in October-November 1998, but the original design has undergone substantial modifications in the next years without any change management nor any design review. The process of improving the software quality and physical behaviour of this subdomain requires substantial investment of manpower and adequate architectural expertise supported by a sound software process: neither is easy to find. Best wishes, Maria Grazia Pia
Hi Andreas, Thanks for reporting the problem, the detailed analysis and fix. I will look into this and evaluate your fix. As Maria Grazia has pointed out, the radioactive decay and the associated codes were developed many years ago and are in need of a design review and upgrade. But at the moment there is no budget for such activities. All of us should try to convience ESA, who funded the original development, to support this. You mentioned the MGGPOD code, as the author of the original decay code in MGGPOD, I am quite interested in any comparison you performed between these two codes. I never managed to find the time to conduct a proper comparison of them. I have missed your initial report on this problem. I rely on the email alert for any problem being reported here. For some reason, I didn't receive the email in February. Fan
Hi Fan, Thanks for the answer. > You mentioned the MGGPOD code, as the author of the original decay code in > MGGPOD, I am quite interested in any comparison you performed between these two > codes. I never managed to find the time to conduct a proper comparison of them. I did some preliminary tests for the Geant4 Space Users Workshop in Tokyo. Here is a link to my talk: http://www.astro.isas.jaxa.jp/conference/g4space5/slide/14/ActivationSimulation.pdf We have basically implemented a first version of activation simulation with Geant4. And I was basically trying to understand the differences between Geant4 simulations, MGGPOD simulations and measurements. I am not yet done with it --- there are still many more differences which I do not yet understand. I will return to this analysis early April, but have to do some more urgent work in the mean time. Ciao & thanks for your efforts, Andreas
Hello Fan & Maria, I have been able to continue debugging the radioactive decay mechanisms in Geant4. The changes to Geant4.9.2-beta1 helped a lot, nevertheless there are several problems remaining, especially concerning meta-stable isotopes: Some of them can be demonstrated by looking at the following decay chain: 73-Ga → 73m-Ge (66.7 keV) → 73m-Ge (13.3 keV) → 73-Ge (1) If there is an internal conversion (electron is ejected instead of gamma-ray), then the energy conservation is sometimes broken: If there is a discrete gamma transition, then the function "void G4DiscreteGammaTransition::SelectGamma()", selects the transition and, if an internal conversion electron is emitted, it randomly selects a shell and SUBTRACTS the binding energy from the ejected electron. The shell number is stored and later retrieved by G4DecayProducts *G4NuclearDecayChannel::DecayIt (G4double theParentMass) to produce the Auger-Electrons. The problem now is, that during G4FragmentVector* G4VGammaDeexcitation::DoChain() MULTIPLE conversions can happen, however, only the LAST shell index is stored, and in the DecayIt function and only the LAST emptied shell is filled, resulting in missing energy! (2) Geant4 correctly recognizes that the 73-Ga decay ends in a metastable isotope at 66.7 keV. However, the second meta-stable state at 13.3 keV is not reached, although the data tables contain the information along with the correct half life. (3) There is only limited user intervention possible to postpone meta-stable decays to the next G4Track, since they occur with the same step. --> (1), (2), and (3) can be fixed by only allowing one decay or one nuclear transition per Geant4-step. (4) Not all meta-stable states are listed in the RadioactiveDecay data sets, but only in the photon-evaporation data sets. (5) The G4ParticleDefinition::PDGLifeTime() is not present for all metastable states → (4) & (5) I think are related. I added some quick fixes in G4RadioactiveDecay which circumvent the problem, but don't really solve it. Enhancement requests: (1) The function deexcitation->SetMaxHalfLife(...) should be accessible by the user from the interface of the G4RadioactiveDecay-class. In the attached patch for Geant4.9.2-b01 all the above problems a fixed *privisorily* My changes start with // ← acz and end with // →. Apply the patch like a normal "Kernel" patch: Simply switch to your $G4INSTALL/source directory, and do: cat <path to>/Decay.patch | patch -p1 Please review them carefully, since I am not familiar enough with the code, to determine if they have any negative side effects. Some patches are really only quick fixes. Since you have written the original code, you probably know best how to fix them. However, running some activation simulations with a Germanium detector utilizing the above fixes, gives results rather close to an alternate activation simulation tool, MGGPOD :) Ciao Andreas
Created attachment 30 [details] Patch for geant4.9.2-b01 Patch for geant4.9.2-b01 to fix the problems mentioned in last comment.
In the original post: Problems 1 and 2 have been fixed as of Geant4 9.3. Problem 3 has been deemed a design improvement and will be handled as a user request, as will enhancement requests 1 (Auger simulation) and 2 (interface symmetry with G4LowEnergyPhotoEffect). In Comment #7: Items 1 (energy non-conservation) and 2 (meta-stable state not reached) will be re-assigned in a separate bug report. Items 3 (postponement of meta-stable decays), 4 (inclusion of meta-stable states in database) and 5 (missing meta-stable state lifetimes) will be treated as design improvements and handled as user requests.
