eic-opticks/qudarap/qutil.h

0001 #pragma once
0002
0003
0004 #if defined(__CUDACC__) || defined(__CUDABE__)
0005    #define QUTIL_METHOD __device__
0006 #else
0007    #define QUTIL_METHOD
0008 #endif
0009
0010
0011 struct qutil
0012 {
0013     QUTIL_METHOD static void rotateUz(float3& d, const float3& u );
0014 };
0015
0016
0017 /**
0018 qutil::rotateUz
0019 -----------------
0020
0021 This rotates the reference frame of a vector such that the original Z-axis will lie in the
0022 direction of *u*. Many rotations would accomplish this; the one selected
0023 uses *u* as its third column and is given by the below matrix.
0024
0025 The below CUDA implementation follows the CLHEP implementation used by Geant4::
0026
0027      // geant4.10.00.p01/source/externals/clhep/src/ThreeVector.cc
0028      72 Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) {
0029      73   // NewUzVector must be normalized !
0030      74
0031      75   double u1 = NewUzVector.x();
0032      76   double u2 = NewUzVector.y();
0033      77   double u3 = NewUzVector.z();
0034      78   double up = u1*u1 + u2*u2;
0035      79
0036      80   if (up>0) {
0037      81       up = std::sqrt(up);
0038      82       double px = dx,  py = dy,  pz = dz;
0039      83       dx = (u1*u3*px - u2*py)/up + u1*pz;
0040      84       dy = (u2*u3*px + u1*py)/up + u2*pz;
0041      85       dz =    -up*px +             u3*pz;
0042      86     }
0043      87   else if (u3 < 0.) { dx = -dx; dz = -dz; }      // phi=0  teta=pi
0044      88   else {};
0045      89   return *this;
0046      90 }
0047
0048 This implements rotation of (px,py,pz) vector into (dx,dy,dz)
0049 using the below rotation matrix, the columns of which must be
0050 orthogonal unit vectors.::
0051
0052                 |  u.x * u.z / up   -u.y / up    u.x  |
0053         d  =    |  u.y * u.z / up   +u.x / up    u.y  |      p
0054                 |   -up               0.         u.z  |
0055
0056 Taking dot products between and within columns shows that to
0057 be the case for normalized u. See oxrap/rotateUz.h for the algebra.
0058
0059 Special cases:
0060
0061 u = [0,0,1] (up=0.)
0062    does nothing, effectively identity matrix
0063
0064 u = [0,0,-1] (up=0., u.z<0. )
0065    flip x, and z which is a rotation of pi/2 about y
0066
0067                |   -1    0     0   |
0068       d =      |    0    1     0   |   p
0069                |    0    0    -1   |
0070
0071 **/
0072
0073 inline QUTIL_METHOD void qutil::rotateUz(float3& d, const float3& u )
0074 {
0075     float up = u.x*u.x + u.y*u.y ;
0076     if (up>0.f)
0077     {
0078         up = sqrt(up);
0079         float px = d.x ;
0080         float py = d.y ;
0081         float pz = d.z ;
0082         d.x = (u.x*u.z*px - u.y*py)/up + u.x*pz;
0083         d.y = (u.y*u.z*px + u.x*py)/up + u.y*pz;
0084         d.z =    -up*px +                u.z*pz;
0085     }
0086     else if (u.z < 0.f )
0087     {
0088         d.x = -d.x;
0089         d.z = -d.z;
0090     }
0091 }
0092
0093
0094