Problem 1931 - G4Torus.cc numerically unstable when rho and Rtor nearly equal
Summary: G4Torus.cc numerically unstable when rho and Rtor nearly equal
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: geometry/solids (show other problems)
Version: 10.3
Hardware: All All
: P4 normal
Assignee: Evgueni.Tcherniaev
URL:
Depends on:
Blocks:
 
Reported: 2016-12-16 17:04 CET by Helmut.Burkhardt
Modified: 2016-12-20 11:28 CET (History)
0 users

See Also:


Attachments
Improved numerical precision using differences in radius and hypot to avoid numerical cancellations between large squared numbers (48.06 KB, text/plain)
2016-12-16 17:04 CET, Helmut.Burkhardt
Details

Note You need to log in before you can comment on or make changes to this problem.
Description Helmut.Burkhardt 2016-12-16 17:04:34 CET
Created attachment 436 [details]
Improved numerical precision using differences in radius and hypot to avoid numerical cancellations between large squared numbers

The calculation of transverse distances called “pt" or “it" (often using pt2, it2 for distance squared) is done in G4Torus.cc in several places using code like

rho2 = p.x()*p.x() + p.y()*p.y();
rho = std::sqrt(rho2);
pt2 = rho2+p.z()*p.z() +fRtor * (fRtor-2*rho);
pt2 = std::max(pt2, 0.0);
pt = std::sqrt(pt2) ;

resulting in poor accuracy (and in some cases incorrect pt=0) when fRtor and rho are similar in magnitude.

A much better precision is obtained using instead :

rho = std::hypot(p.x(),p.y());
pt  = std::hypot(p.z(),rho-fRtor);

or basically using (a-b)^2 instead of $a^2+b^2-2ab$ avoiding cancellations between large squared numbers.

Improved G4Torus.cc attached.
Comment 1 Evgueni.Tcherniaev 2016-12-20 11:28:14 CET
Dear Helmut, thank you for the improvements. The code has been committed to the Geant4
repository and will be included	in the next patch release.