00001 #ifndef GeometryVector_newBasic3DVector_h
00002 #define GeometryVector_newBasic3DVector_h
00003
00004 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
00005 #include "DataFormats/GeometryVector/interface/Theta.h"
00006 #include "DataFormats/GeometryVector/interface/Phi.h"
00007 #include "DataFormats/GeometryVector/interface/PreciseFloatType.h"
00008 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
00009 #include "DataFormats/Math/interface/SSEVec.h"
00010 #include <iosfwd>
00011 #include <cmath>
00012
00013 namespace detailsBasic3DVector {
00014 inline float __attribute__((always_inline)) __attribute__ ((pure))
00015 eta(float x, float y, float z) { float t(z/std::sqrt(x*x+y*y)); return ::asinhf(t);}
00016 inline double __attribute__((always_inline)) __attribute__ ((pure))
00017 eta(double x, double y, double z) { double t(z/std::sqrt(x*x+y*y)); return ::asinh(t);}
00018 inline long double __attribute__((always_inline)) __attribute__ ((pure))
00019 eta(long double x, long double y, long double z) { long double t(z/std::sqrt(x*x+y*y)); return ::asinhl(t);}
00020 }
00021
00022
00023 template < typename T>
00024 class Basic3DVector {
00025 public:
00026
00027 typedef T ScalarType;
00028 typedef Geom::Cylindrical2Cartesian<T> Cylindrical;
00029 typedef Geom::Spherical2Cartesian<T> Spherical;
00030 typedef Spherical Polar;
00031
00036 Basic3DVector() {}
00037
00039 Basic3DVector( const Basic3DVector & p) :
00040 v(p.v) {}
00041
00043 template <class U>
00044 Basic3DVector( const Basic3DVector<U> & p) :
00045 v(p.v) {}
00046
00047
00049 Basic3DVector( const Basic2DVector<T> & p) :
00050 v(p.x(),p.y(),0) {}
00051
00052
00061 template <class OtherPoint>
00062 explicit Basic3DVector( const OtherPoint& p) :
00063 v(p.x(),p.y(),p.z()) {}
00064
00065
00066
00067 template<class U>
00068 Basic3DVector(mathSSE::Vec4<U> const& iv) : v(iv){}
00069
00071 Basic3DVector( const T& x, const T& y, const T& z) :
00072 v(x,y,z){}
00073
00078 template <typename U>
00079 Basic3DVector( const Geom::Theta<U>& theta,
00080 const Geom::Phi<U>& phi, const T& r) {
00081 Polar p( theta.value(), phi.value(), r);
00082 v.o.theX = p.x(); v.o.theY = p.y(); v.o.theZ = p.z();
00083 }
00084
00086 T x() const { return v.o.theX;}
00087
00089 T y() const { return v.o.theY;}
00090
00092 T z() const { return v.o.theZ;}
00093
00094 Basic2DVector<T> xy() const { return v.xy();}
00095
00096
00097 bool operator==(const Basic3DVector& rh) const {
00098 return v==rh.v;
00099 }
00100
00102 T mag2() const { return ::dot(v,v);}
00103
00105 T mag() const { return std::sqrt( mag2());}
00106
00108 T perp2() const { return x()*x() + y()*y();}
00109
00111 T perp() const { return std::sqrt( perp2());}
00112
00114 T transverse() const { return perp();}
00115
00120 T barePhi() const {return std::atan2(y(),x());}
00121 Geom::Phi<T> phi() const {return Geom::Phi<T>(barePhi());}
00122
00127 T bareTheta() const {return std::atan2(perp(),z());}
00128 Geom::Theta<T> theta() const {return Geom::Theta<T>(std::atan2(perp(),z()));}
00129
00134
00135 T eta() const { return detailsBasic3DVector::eta(x(),y(),z());}
00136
00140 Basic3DVector unit() const {
00141 T my_mag = mag2();
00142 if (my_mag==0) return *this;
00143 my_mag = T(1)/std::sqrt(my_mag);
00144 return *this * my_mag;
00145 }
00146
00149 template <class U>
00150 Basic3DVector& operator+= ( const Basic3DVector<U>& p) {
00151 v = v + p.v;
00152 return *this;
00153 }
00154
00157 template <class U>
00158 Basic3DVector& operator-= ( const Basic3DVector<U>& p) {
00159 v = v - p.v;
00160 return *this;
00161 }
00162
00164 Basic3DVector operator-() const { return Basic3DVector(-v);}
00165
00167 Basic3DVector& operator*= ( T t) {
00168 v = t*v;
00169 return *this;
00170 }
00171
00173 Basic3DVector& operator/= ( T t) {
00174 t = T(1)/t;
00175 v = t*v;
00176 return *this;
00177 }
00178
00180 T dot( const Basic3DVector& rh) const {
00181 return ::dot(v,rh.v);
00182 }
00183
00189 template <class U>
00190 typename PreciseFloatType<T,U>::Type dot( const Basic3DVector<U>& lh) const {
00191 return Basic3DVector<typename PreciseFloatType<T,U>::Type>(*this)
00192 .dot(Basic3DVector<typename PreciseFloatType<T,U>::Type>(lh));
00193 }
00194
00196 Basic3DVector cross( const Basic3DVector& lh) const {
00197 return ::cross(v,lh.v);
00198 }
00199
00200
00206 template <class U>
00207 Basic3DVector<typename PreciseFloatType<T,U>::Type>
00208 cross( const Basic3DVector<U>& lh) const {
00209 return Basic3DVector<typename PreciseFloatType<T,U>::Type>(*this)
00210 .cross(Basic3DVector<typename PreciseFloatType<T,U>::Type>(lh));
00211 }
00212
00213 public:
00214 mathSSE::Vec4<T> v;
00215 } __attribute__ ((aligned (16)));
00216
00217
00218 namespace geometryDetails {
00219 std::ostream & print3D(std::ostream& s, double x, double y, double z);
00220 }
00221
00223 template <class T>
00224 inline std::ostream & operator<<( std::ostream& s, const Basic3DVector<T>& v) {
00225 return geometryDetails::print3D(s, v.x(),v.y(), v.z());
00226 }
00227
00228
00230 template <class T>
00231 inline Basic3DVector<T>
00232 operator+( const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
00233 return a.v+b.v;
00234 }
00235 template <class T>
00236 inline Basic3DVector<T>
00237 operator-( const Basic3DVector<T>& a, const Basic3DVector<T>& b) {
00238 return a.v-b.v;
00239 }
00240
00241 template <class T, class U>
00242 inline Basic3DVector<typename PreciseFloatType<T,U>::Type>
00243 operator+( const Basic3DVector<T>& a, const Basic3DVector<U>& b) {
00244 typedef Basic3DVector<typename PreciseFloatType<T,U>::Type> RT;
00245 return RT(a).v+RT(b).v;
00246 }
00247
00248 template <class T, class U>
00249 inline Basic3DVector<typename PreciseFloatType<T,U>::Type>
00250 operator-( const Basic3DVector<T>& a, const Basic3DVector<U>& b) {
00251 typedef Basic3DVector<typename PreciseFloatType<T,U>::Type> RT;
00252 return RT(a).v-RT(b).v;
00253 }
00254
00256 template <class T>
00257 inline T operator*( const Basic3DVector<T>& v1, const Basic3DVector<T>& v2) {
00258 return v1.dot(v2);
00259 }
00260
00262 template <class T, class U>
00263 inline typename PreciseFloatType<T,U>::Type operator*( const Basic3DVector<T>& v1,
00264 const Basic3DVector<U>& v2) {
00265 return v1.dot(v2);
00266 }
00267
00271 template <class T>
00272 inline Basic3DVector<T> operator*( const Basic3DVector<T>& v, T t) {
00273 return v.v*t;
00274 }
00275
00277 template <class T>
00278 inline Basic3DVector<T> operator*(T t, const Basic3DVector<T>& v) {
00279 return v.v*t;
00280 }
00281
00282 template <class T, typename S>
00283 inline Basic3DVector<T> operator*(S t, const Basic3DVector<T>& v) {
00284 return static_cast<T>(t)*v;
00285 }
00286
00287 template <class T, typename S>
00288 inline Basic3DVector<T> operator*(const Basic3DVector<T>& v, S t) {
00289 return static_cast<T>(t)*v;
00290 }
00291
00292
00296 template <class T, typename S>
00297 inline Basic3DVector<T> operator/( const Basic3DVector<T>& v, S s) {
00298 T t = T(1)/s;
00299 return v*t;
00300 }
00301
00302
00303 typedef Basic3DVector<float> Basic3DVectorF;
00304 typedef Basic3DVector<double> Basic3DVectorD;
00305
00306
00307
00308 #include "Basic3DVectorLD.h"
00309
00310 #endif // GeometryVector_Basic3DVector_h
00311
00312