00001 #ifndef Geom_oldTkRotation_H
00002 #define Geom_oldTkRotation_H
00003
00004 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
00005 #include "DataFormats/GeometryVector/interface/Basic3DVector.h"
00006 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00007
00008
00009
00010
00011 #include <iosfwd>
00012
00013 template <class T> class TkRotation;
00014
00015 template <class T>
00016 std::ostream & operator<<( std::ostream& s, const TkRotation<T>& r);
00017
00018 namespace geometryDetails {
00019 void TkRotationErr1();
00020 void TkRotationErr2();
00021
00022 }
00023
00024
00028 template <class T>
00029 class TkRotation {
00030 public:
00031 typedef Vector3DBase< T, GlobalTag> GlobalVector;
00032
00033 TkRotation() :
00034 R11( 1), R12( 0), R13( 0),
00035 R21( 0), R22( 1), R23( 0),
00036 R31( 0), R32( 0), R33( 1) {}
00037
00038 TkRotation( T xx, T xy, T xz, T yx, T yy, T yz, T zx, T zy, T zz) :
00039 R11(xx), R12(xy), R13(xz),
00040 R21(yx), R22(yy), R23(yz),
00041 R31(zx), R32(zy), R33(zz) {}
00042
00043 TkRotation( const T* p) :
00044 R11(p[0]), R12(p[1]), R13(p[2]),
00045 R21(p[3]), R22(p[4]), R23(p[5]),
00046 R31(p[6]), R32(p[7]), R33(p[8]) {}
00047
00048 TkRotation( const GlobalVector & aX, const GlobalVector & aY) {
00049
00050 GlobalVector uX = aX.unit();
00051 GlobalVector uY = aY.unit();
00052 GlobalVector uZ(uX.cross(uY));
00053
00054 R11 = uX.x(); R12 = uX.y(); R13 = uX.z();
00055 R21 = uY.x(); R22 = uY.y(); R23 = uY.z();
00056 R31 = uZ.x(); R32 = uZ.y(); R33 = uZ.z();
00057
00058 }
00059
00064 TkRotation( const GlobalVector & aX, const GlobalVector & aY,
00065 const GlobalVector & aZ) :
00066 R11( aX.x()), R12( aX.y()), R13( aX.z()),
00067 R21( aY.x()), R22( aY.y()), R23( aY.z()),
00068 R31( aZ.x()), R32( aZ.y()), R33( aZ.z()) {}
00069
00070
00079 TkRotation( const Basic3DVector<T>& axis, T phi) :
00080 R11( cos(phi) ), R12( sin(phi)), R13( 0),
00081 R21( -sin(phi)), R22( cos(phi)), R23( 0),
00082 R31( 0), R32( 0), R33( 1) {
00083
00084
00085 TkRotation tmpRotz ( cos(axis.phi()), sin(axis.phi()), 0.,
00086 -sin(axis.phi()), cos(axis.phi()), 0.,
00087 0., 0., 1. );
00088
00089 TkRotation tmpRoty ( cos(axis.theta()), 0.,-sin(axis.theta()),
00090 0., 1., 0.,
00091 sin(axis.theta()), 0., cos(axis.theta()) );
00092 (*this)*=tmpRoty;
00093 (*this)*=tmpRotz;
00094
00095
00096
00097 *this = (tmpRoty*tmpRotz).multiplyInverse(*this);
00098
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 template <typename U>
00134 TkRotation( const TkRotation<U>& a) :
00135 R11(a.xx()), R12(a.xy()), R13(a.xz()),
00136 R21(a.yx()), R22(a.yy()), R23(a.yz()),
00137 R31(a.zx()), R32(a.zy()), R33(a.zz()) {}
00138
00139 TkRotation transposed() const {
00140 return TkRotation( R11, R21, R31,
00141 R12, R22, R32,
00142 R13, R23, R33);
00143 }
00144
00145 Basic3DVector<T> operator*( const Basic3DVector<T>& v) const {
00146 return Basic3DVector<T>( R11*v.x() + R12*v.y() + R13*v.z(),
00147 R21*v.x() + R22*v.y() + R23*v.z(),
00148 R31*v.x() + R32*v.y() + R33*v.z());
00149 }
00150
00151 Basic3DVector<T> multiplyInverse( const Basic3DVector<T>& v) const {
00152 return Basic3DVector<T>( R11*v.x() + R21*v.y() + R31*v.z(),
00153 R12*v.x() + R22*v.y() + R32*v.z(),
00154 R13*v.x() + R23*v.y() + R33*v.z());
00155 }
00156
00157 #ifndef CMS_NO_TEMPLATE_MEMBERS
00158 template <class Scalar>
00159 Basic3DVector<Scalar> multiplyInverse( const Basic3DVector<Scalar>& v) const {
00160 return Basic3DVector<Scalar>( R11*v.x() + R21*v.y() + R31*v.z(),
00161 R12*v.x() + R22*v.y() + R32*v.z(),
00162 R13*v.x() + R23*v.y() + R33*v.z());
00163 }
00164 #endif
00165
00166 Basic3DVector<T> operator*( const Basic2DVector<T>& v) const {
00167 return Basic3DVector<T>( R11*v.x() + R12*v.y(),
00168 R21*v.x() + R22*v.y(),
00169 R31*v.x() + R32*v.y());
00170 }
00171 Basic3DVector<T> multiplyInverse( const Basic2DVector<T>& v) const {
00172 return Basic3DVector<T>( R11*v.x() + R21*v.y(),
00173 R12*v.x() + R22*v.y(),
00174 R13*v.x() + R23*v.y());
00175 }
00176
00177
00178
00179 TkRotation operator*( const TkRotation& b) const {
00180 return TkRotation(R11*b.R11 + R12*b.R21 + R13*b.R31,
00181 R11*b.R12 + R12*b.R22 + R13*b.R32,
00182 R11*b.R13 + R12*b.R23 + R13*b.R33,
00183 R21*b.R11 + R22*b.R21 + R23*b.R31,
00184 R21*b.R12 + R22*b.R22 + R23*b.R32,
00185 R21*b.R13 + R22*b.R23 + R23*b.R33,
00186 R31*b.R11 + R32*b.R21 + R33*b.R31,
00187 R31*b.R12 + R32*b.R22 + R33*b.R32,
00188 R31*b.R13 + R32*b.R23 + R33*b.R33);
00189 }
00190
00191 TkRotation multiplyInverse( const TkRotation& b) const {
00192 return TkRotation(R11*b.R11 + R21*b.R21 + R31*b.R31,
00193 R11*b.R12 + R21*b.R22 + R31*b.R32,
00194 R11*b.R13 + R21*b.R23 + R31*b.R33,
00195 R12*b.R11 + R22*b.R21 + R32*b.R31,
00196 R12*b.R12 + R22*b.R22 + R32*b.R32,
00197 R12*b.R13 + R22*b.R23 + R32*b.R33,
00198 R13*b.R11 + R23*b.R21 + R33*b.R31,
00199 R13*b.R12 + R23*b.R22 + R33*b.R32,
00200 R13*b.R13 + R23*b.R23 + R33*b.R33);
00201 }
00202
00203 TkRotation& operator*=( const TkRotation& b) {
00204 return *this = operator * (b);
00205 }
00206
00207
00208 TkRotation& transform(const TkRotation& b) {
00209 return *this = b.operator * (*this);
00210 }
00211
00212 TkRotation & rotateAxes(const Basic3DVector<T>& newX,
00213 const Basic3DVector<T>& newY,
00214 const Basic3DVector<T>& newZ) {
00215 T del = 0.001;
00216
00217 if (
00218
00219
00220
00221
00222
00223
00224
00225 fabs(newX.mag2()-1.) > del ||
00226 fabs(newY.mag2()-1.) > del ||
00227 fabs(newZ.mag2()-1.) > del ||
00228 fabs(newX.dot(newY)) > del ||
00229 fabs(newY.dot(newZ)) > del ||
00230 fabs(newZ.dot(newX)) > del) {
00231 geometryDetails::TkRotationErr2();
00232 return *this;
00233 } else {
00234 return transform(TkRotation(newX.x(), newY.x(), newZ.x(),
00235 newX.y(), newY.y(), newZ.y(),
00236 newX.z(), newY.z(), newZ.z()));
00237 }
00238 }
00239
00240 Basic3DVector<T> x() const { return Basic3DVector<T>(xx(),xy(),xz());}
00241 Basic3DVector<T> y() const { return Basic3DVector<T>(yx(),yy(),yz());}
00242 Basic3DVector<T> z() const { return Basic3DVector<T>(zx(),zy(),zz());}
00243
00244
00245 T const &xx() const { return R11;}
00246 T const &xy() const { return R12;}
00247 T const &xz() const { return R13;}
00248 T const &yx() const { return R21;}
00249 T const &yy() const { return R22;}
00250 T const &yz() const { return R23;}
00251 T const &zx() const { return R31;}
00252 T const &zy() const { return R32;}
00253 T const &zz() const { return R33;}
00254
00255 private:
00256
00257 T R11, R12, R13;
00258 T R21, R22, R23;
00259 T R31, R32, R33;
00260 };
00261
00262
00263 template<>
00264 std::ostream & operator<< <float>( std::ostream& s, const TkRotation<float>& r);
00265
00266 template<>
00267 std::ostream & operator<< <double>( std::ostream& s, const TkRotation<double>& r);
00268
00269
00270 template <class T, class U>
00271 inline Basic3DVector<U> operator*( const TkRotation<T>& r, const Basic3DVector<U>& v) {
00272 return Basic3DVector<U>( r.xx()*v.x() + r.xy()*v.y() + r.xz()*v.z(),
00273 r.yx()*v.x() + r.yy()*v.y() + r.yz()*v.z(),
00274 r.zx()*v.x() + r.zy()*v.y() + r.zz()*v.z());
00275 }
00276
00277 template <class T, class U>
00278 inline TkRotation<typename PreciseFloatType<T,U>::Type>
00279 operator*( const TkRotation<T>& a, const TkRotation<U>& b) {
00280 typedef TkRotation<typename PreciseFloatType<T,U>::Type> RT;
00281 return RT( a.xx()*b.xx() + a.xy()*b.yx() + a.xz()*b.zx(),
00282 a.xx()*b.xy() + a.xy()*b.yy() + a.xz()*b.zy(),
00283 a.xx()*b.xz() + a.xy()*b.yz() + a.xz()*b.zz(),
00284 a.yx()*b.xx() + a.yy()*b.yx() + a.yz()*b.zx(),
00285 a.yx()*b.xy() + a.yy()*b.yy() + a.yz()*b.zy(),
00286 a.yx()*b.xz() + a.yy()*b.yz() + a.yz()*b.zz(),
00287 a.zx()*b.xx() + a.zy()*b.yx() + a.zz()*b.zx(),
00288 a.zx()*b.xy() + a.zy()*b.yy() + a.zz()*b.zy(),
00289 a.zx()*b.xz() + a.zy()*b.yz() + a.zz()*b.zz());
00290 }
00291
00292 #endif
00293
00294
00295
00296