fix for CLHEP-141 supplied by Evgueni Tcherniaev

ChangeLog

+
+2017-03-21 Lynn Garren <garren@fnal.gov>
+
+  * Vector/src/RotationA.cc:
+    fix for CLHEP-141 supplied by Evgueni Tcherniaev
+  * Vector/test/testRotationAxis.cc:
+    test for CLHEP-141 supplied by Evgueni Tcherniaev
+
 ==============================
 30.11.16 Release CLHEP-2.3.4.3
 ==============================

Vector/ChangeLog

+
+2017-03-21 Lynn Garren <garren@fnal.gov>
+
+  * src/RotationA.cc:
+    fix for CLHEP-141 supplied by Evgueni Tcherniaev
+  * test/testRotationAxis.cc:
+    test for CLHEP-141 supplied by Evgueni Tcherniaev
+
 ==============================
 30.11.16 Release CLHEP-2.3.4.3
 ==============================

Vector/src/RotationA.cc

@@ -64,8 +64,6 @@ HepRotation::HepRotation ( const HepAxisAngle & ax )
   set ( ax.axis(), ax.delta() );
 }
 
-
-
 double    HepRotation::delta() const {
 
   double cosdelta = (rxx + ryy + rzz - 1.0) / 2.0;
@@ -81,20 +79,42 @@ double    HepRotation::delta() const {
 
 Hep3Vector HepRotation::axis () const {
 
-  // Determine 2*std::sin(delta) times the u components (I call this uX, uY, Uz)
-  // Normalization is not needed; it will be done when returning the 3-Vector
-
-  double  Uz = ryx - rxy;
-  double  Uy = rxz - rzx;
-  double  Ux = rzy - ryz;
-
-  if ( (Uz==0) && (Uy==0) && (Ux==0) ) {
-    if        ( rzz>0 ) {
-      return Hep3Vector(0,0,1);
-    } else if ( ryy>0 ) {
-      return Hep3Vector(0,1,0);
+  const double eps = 1e-15;
+
+  double Ux = rzy - ryz;
+  double Uy = rxz - rzx;
+  double Uz = ryx - rxy;
+  if (std::abs(Ux) < eps && std::abs(Uy) < eps && std::abs(Uz) < eps) {
+
+    double cosdelta = (rxx + ryy + rzz - 1.0) / 2.0;
+    if (cosdelta > 0.0) return Hep3Vector(0,0,1); // angle = 0, any axis is good
+
+    double xx = (rxx + 1)/2;
+    double yy = (ryy + 1)/2;
+    double zz = (rzz + 1)/2;
+    double xy = (rxy + ryx)/4;
+    double xz = (rxz + rzx)/4;
+    double yz = (ryz + rzy)/4;
+    double x, y, z;
+
+    if (xx > ryy && xx > rzz) {
+      x = std::sqrt(xx);
+      if (rzy - ryz < 0) x = -x;
+      y = xy/x;
+      z = xz/x;
+      return  Hep3Vector( x, y, z ).unit();
+    } else if (yy > zz) {
+      y = std::sqrt(yy);
+      if (rxz - rzx < 0) y = -y;
+      x = xy/y;
+      z = yz/y;
+      return  Hep3Vector( x, y, z ).unit();
     } else {
-      return Hep3Vector(1,0,0);
+      z = std::sqrt(zz);
+      if (ryx - rxy < 0) z = -z;
+      x = xz/z;
+      y = yz/z;
+      return  Hep3Vector( x, y, z ).unit();
     }
   } else {
     return  Hep3Vector( Ux, Uy, Uz ).unit();