Problem 2123 - Excited nucleus sometimes fails to decay
Summary: Excited nucleus sometimes fails to decay
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: processes/hadronic/models/de_excitation/photon_evaporation (show other problems)
Version: 10.5
Hardware: All All
: P4 normal
Assignee: Vladimir.Ivantchenko
URL:
Depends on:
Blocks:
 
Reported: 2019-01-28 01:18 CET by Clifford Dugal
Modified: 2019-04-08 16:31 CEST (History)
1 user (show)

See Also:


Attachments
New Photon evaporation file for At217 (3.35 KB, text/plain)
2019-02-01 17:34 CET, Laurent Desorgher
Details

Note You need to log in before you can comment on or make changes to this problem.
Description Clifford Dugal 2019-01-28 01:18:50 CET
I've found that sometimes At-217 at 310.3 keV excitation is not able to decay, depending on the randomly selected final energy of the de-excitation calculation.

It appears that the problem is in G4PhotonEvaporation::GenerateGamma:

---
    // final discrete level
    if(fLevelManager) {
      if(efinal < fLevelEnergyMax) {
	//G4cout << "Efinal= " << efinal << "  idx= " << fIndex << G4endl;
	fIndex = fLevelManager->NearestLevelIndex(efinal, fIndex);
	efinal = fLevelManager->LevelEnergy(fIndex);
	// protection - take level below
	if(efinal >= eexc && 0 < fIndex) {
	  --fIndex;
	  efinal = fLevelManager->LevelEnergy(fIndex);
	} 
	nucleus->SetFloatingLevelNumber(fLevelManager->FloatingLevel(fIndex));

	// not allowed to have final energy above max energy
	// if G4LevelManager exist
      } else {
	efinal = fLevelEnergyMax;
	fIndex = fLevelManager->NearestLevelIndex(efinal, fIndex);
      }
    }
---

For the case that fails to decay, the "protection" if statement above sees the following:

  efinal = 0.31030000000000002
  eexc   = 0.31030000001192093

In this case, eexc is numerically just slightly larger than efinal, so the "protection" code forcing a de-excitation is not run, no de-excitation occurs, and At-217 at 310.3 keV excitation appears not decay further.
Comment 1 Vladimir.Ivantchenko 2019-01-30 17:16:07 CET
Hello,

thank you for pointing to the problem. I had a look and my impression of today is different from yours. I see, that At-217 level 310.3 has no decay channels: 

PhotonEvaporation5.3/z85.a217 

If you have some extra information about this level, please, let us know. I will try to understand, what we should do for such levels. 

VI
Comment 2 Clifford Dugal 2019-01-30 21:15:38 CET
I'm new to Geant4, and I only looked into why At-217 at 310.3 keV sometimes doesn't decay. It's defined to be a daughter of Fr-221, so occasionally I observed my longer decay chain not reach a stable nuclide (I'm using something close to the rdecay01 example). I didn't look at a "good" decay for reference, so perhaps the photon evaporation file is really what needs to be adjusted for full correctness.

Most of the time (>94.7%), Geant4 has At-217 at 310.3 keV decay with an initial gamma of 38.23, 92.18, 210.05, or 310.30 keV. The energy is selected in the G4PhotonEvaporation::GenerateGamma function, despite there not being a decay channel in PhotonEvaporation5.3/z85.a217.

It appears to me that the "protection" code in G4PhotonEvaporation::GenerateGamma was intended to force a decay to the nearest lower energy level in the case that the selected final level matched (or even exceeded) the initial level. It seems like numerical precision issues made the if statement not work as intended, so the forced decay did not occur.

While I don't know exactly how this should be handled, I think that decays for excited nuclei should never get stuck. I would much rather have any reasonable de-excitation (even if it's not perfectly correct) and continue with the decay chain than have the remainder of the decay chain omitted.
Comment 3 Laurent Desorgher 2019-02-01 17:34:56 CET
Created attachment 539 [details]
New Photon evaporation file  for At217

Please find attached to this comment  a new version of the photoevaporation file for At217 with the 310.3 level decay corrected.

Please let us now if this solve your issue


Best regards


Laurent Desorgher
Comment 4 Clifford Dugal 2019-02-04 18:08:23 CET
The attached file doesn't quite work (at least with Geant4 10.4 at work). It is missing the gamma de-excitation channels for levels 12, 14, and 15, and that causes a G4Exception:

===

### G4LevelReader: broken transition 0 from level 12 to 13 for isotope Z= 85 A= 217 - use ground level

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : had014
      issued by : G4LevelReader::LevelManager(..)
 G4LevelReader: wrong data file for Z= 85 A= 217 level #13 has index 0
Check G4LEVELGAMMADATA

===

Reverting back to what was in the original z85.a217 file for those broken levels fixes it.

Once the file is fixed, At-217 at 310.3 keV excitation does appear to consistently decay. The specific instance that revealed the issue to me is now addressed. Thanks!


Making sure that this isn't an issue elsewhere, I tried Rn-217 at 295 keV excitation. It decays over 90.7% of the time, but otherwise has the same issue of not decaying. It was rather easy for me to find another broken case, so I'm thinking that this is probably a wide-spread issue.

Even if you manage to fix all of the photon evaporation data files now, I think it's bad practise to keep the partially-silently-failing code in Geant4. I think Geant4 should consistently do one of two things:

    1. Always decay excited atoms, even if gamma de-excitation channels are missing.

    2. Never decay excited atoms with missing gamma de-excitation channels. Instead, throw an exception.

Currently, Geant4 is mostly taking the approach of 1. (except when the code silently fails), and I prefer that approach for my current work. I suppose Geant4 could also have switchable behaviour between the above two options.
Comment 5 Vladimir.Ivantchenko 2019-04-08 16:31:00 CEST
Hello,

For the case, when there is no data for decay probabilities, the fix is introduced: decay will be performed to be closest level.

The fix is validated and will be available with the next public version and in the patch for 10.5.

VI