| Summary: | Rounding error in BetaDecayIt() method of G4NuclearDecayChannel | ||
|---|---|---|---|
| Product: | Geant4 | Reporter: | Ben Morgan <Ben.Morgan> |
| Component: | processes/hadronic/models | Assignee: | dennis.herbert.wright |
| Status: | RESOLVED FIXED | ||
| Severity: | major | ||
| Priority: | P3 | ||
| Version: | 7.1 | ||
| Hardware: | PC | ||
| OS: | Linux | ||
Thanks, Ben. Your fix will go in the next patch or release. Just to note, the rounding error issue has not been fixed in release 8.1. I'm sorry, I made your change to a similar, but commented out, part of the code and thus it had no effect. The change is now made in the right part of the code and will be included in the next patch. Thanks for pointing it out. |
Hi, The calculation of the momentum of the recoiling nucleus on Line 647 of G4NuclearDecayChannel (in the BetaDecayIt() method): daughtermomentum[1] = std::sqrt(daughterenergy[1]*daughterenergy[1] + 2.0*daughterenergy[1] * daughtermass[1]); can result in a 'nan' and subsequent lock of an application because the quantity whose square root is taken can go negative due to rounding errors from previous calculations. A suggested fix would be to check that the quantity in the square root is positive before taking the square root, and take it as zero otherwise: G4double recoilmomentumsquared = daughterenergy[1]*daughterenergy[1] + 2.0*daughterenergy[1] * daughtermass[1]; if( recoilmomentumsquared < 0.0) { recoilmomentumsquared = 0.0; } daughtermomentum[1] = std::sqrt(recoilmomentumsquared); I've only seen the quantity 'recoilmomentumsquared' go to slightly negative values (typically of order 10^-8), so I don't think taking it as zero will violate conservation of energy too much :). Ben Morgan.