CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: BasicVector3D.cc,v 1.3 2003/08/13 20:00:11 garren Exp $ 00003 // --------------------------------------------------------------------------- 00004 00005 #include <math.h> 00006 #include <float.h> 00007 #include <iostream> 00008 #include "CLHEP/Geometry/defs.h" 00009 #include "CLHEP/Geometry/BasicVector3D.h" 00010 00011 namespace HepGeom { 00012 //-------------------------------------------------------------------------- 00013 template<> 00014 float BasicVector3D<float>::pseudoRapidity() const { 00015 float ma = mag(), dz = z(); 00016 if (ma == 0) return 0; 00017 if (ma == dz) return FLT_MAX; 00018 if (ma == -dz) return -FLT_MAX; 00019 return 0.5*log((ma+dz)/(ma-dz)); 00020 } 00021 00022 //-------------------------------------------------------------------------- 00023 template<> 00024 void BasicVector3D<float>::setEta(float a) { 00025 double ma = mag(); 00026 if (ma == 0) return; 00027 double tanHalfTheta = exp(-a); 00028 double tanHalfTheta2 = tanHalfTheta * tanHalfTheta; 00029 double cosTheta = (1 - tanHalfTheta2) / (1 + tanHalfTheta2); 00030 double rh = ma * sqrt(1 - cosTheta*cosTheta); 00031 double ph = phi(); 00032 set(rh*cos(ph), rh*sin(ph), ma*cosTheta); 00033 } 00034 00035 //-------------------------------------------------------------------------- 00036 template<> 00037 float BasicVector3D<float>::angle(const BasicVector3D<float> & v) const { 00038 double cosa = 0; 00039 double ptot = mag()*v.mag(); 00040 if(ptot > 0) { 00041 cosa = dot(v)/ptot; 00042 if(cosa > 1) cosa = 1; 00043 if(cosa < -1) cosa = -1; 00044 } 00045 return acos(cosa); 00046 } 00047 00048 //-------------------------------------------------------------------------- 00049 template<> 00050 BasicVector3D<float> & BasicVector3D<float>::rotateX(float a) { 00051 double sina = sin(a), cosa = cos(a), dy = y(), dz = z(); 00052 setY(dy*cosa-dz*sina); 00053 setZ(dz*cosa+dy*sina); 00054 return *this; 00055 } 00056 00057 //-------------------------------------------------------------------------- 00058 template<> 00059 BasicVector3D<float> & BasicVector3D<float>::rotateY(float a) { 00060 double sina = sin(a), cosa = cos(a), dz = z(), dx = x(); 00061 setZ(dz*cosa-dx*sina); 00062 setX(dx*cosa+dz*sina); 00063 return *this; 00064 } 00065 00066 //-------------------------------------------------------------------------- 00067 template<> 00068 BasicVector3D<float> & BasicVector3D<float>::rotateZ(float a) { 00069 double sina = sin(a), cosa = cos(a), dx = x(), dy = y(); 00070 setX(dx*cosa-dy*sina); 00071 setY(dy*cosa+dx*sina); 00072 return *this; 00073 } 00074 00075 //-------------------------------------------------------------------------- 00076 template<> 00077 BasicVector3D<float> & 00078 BasicVector3D<float>::rotate(float a, const BasicVector3D<float> & v) { 00079 if (a == 0) return *this; 00080 double cx = v.x(), cy = v.y(), cz = v.z(); 00081 double ll = sqrt(cx*cx + cy*cy + cz*cz); 00082 if (ll == 0) { 00083 std::cerr << "BasicVector<float>::rotate() : zero axis" << std::endl; 00084 return *this; 00085 } 00086 double cosa = cos(a), sina = sin(a); 00087 cx /= ll; cy /= ll; cz /= ll; 00088 00089 double xx = cosa + (1-cosa)*cx*cx; 00090 double xy = (1-cosa)*cx*cy - sina*cz; 00091 double xz = (1-cosa)*cx*cz + sina*cy; 00092 00093 double yx = (1-cosa)*cy*cx + sina*cz; 00094 double yy = cosa + (1-cosa)*cy*cy; 00095 double yz = (1-cosa)*cy*cz - sina*cx; 00096 00097 double zx = (1-cosa)*cz*cx - sina*cy; 00098 double zy = (1-cosa)*cz*cy + sina*cx; 00099 double zz = cosa + (1-cosa)*cz*cz; 00100 00101 cx = x(); cy = y(); cz = z(); 00102 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz); 00103 return *this; 00104 } 00105 00106 //-------------------------------------------------------------------------- 00107 std::ostream & 00108 operator<<(std::ostream & os, const BasicVector3D<float> & a) 00109 { 00110 return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")"; 00111 } 00112 00113 //-------------------------------------------------------------------------- 00114 std::istream & 00115 operator>> (std::istream & is, BasicVector3D<float> & a) 00116 { 00117 // Required format is ( a, b, c ) that is, three numbers, preceded by 00118 // (, followed by ), and separated by commas. The three numbers are 00119 // taken as x, y, z. 00120 00121 float x, y, z; 00122 char c; 00123 00124 is >> std::ws >> c; 00125 // ws is defined to invoke eatwhite(istream & ) 00126 // see (Stroustrup gray book) page 333 and 345. 00127 if (is.fail() || c != '(' ) { 00128 std::cerr 00129 << "Could not find required opening parenthesis " 00130 << "in input of a BasicVector3D<float>" 00131 << std::endl; 00132 return is; 00133 } 00134 00135 is >> x >> std::ws >> c; 00136 if (is.fail() || c != ',' ) { 00137 std::cerr 00138 << "Could not find x value and required trailing comma " 00139 << "in input of a BasicVector3D<float>" 00140 << std::endl; 00141 return is; 00142 } 00143 00144 is >> y >> std::ws >> c; 00145 if (is.fail() || c != ',' ) { 00146 std::cerr 00147 << "Could not find y value and required trailing comma " 00148 << "in input of a BasicVector3D<float>" 00149 << std::endl; 00150 return is; 00151 } 00152 00153 is >> z >> std::ws >> c; 00154 if (is.fail() || c != ')' ) { 00155 std::cerr 00156 << "Could not find z value and required close parenthesis " 00157 << "in input of a BasicVector3D<float>" 00158 << std::endl; 00159 return is; 00160 } 00161 00162 a.setX(x); 00163 a.setY(y); 00164 a.setZ(z); 00165 return is; 00166 } 00167 00168 //-------------------------------------------------------------------------- 00169 template<> 00170 double BasicVector3D<double>::pseudoRapidity() const { 00171 double ma = mag(), dz = z(); 00172 if (ma == 0) return 0; 00173 if (ma == dz) return DBL_MAX; 00174 if (ma == -dz) return -DBL_MAX; 00175 return 0.5*log((ma+dz)/(ma-dz)); 00176 } 00177 00178 //-------------------------------------------------------------------------- 00179 template<> 00180 void BasicVector3D<double>::setEta(double a) { 00181 double ma = mag(); 00182 if (ma == 0) return; 00183 double tanHalfTheta = exp(-a); 00184 double tanHalfTheta2 = tanHalfTheta * tanHalfTheta; 00185 double cosTheta = (1 - tanHalfTheta2) / (1 + tanHalfTheta2); 00186 double rh = ma * sqrt(1 - cosTheta*cosTheta); 00187 double ph = phi(); 00188 set(rh*cos(ph), rh*sin(ph), ma*cosTheta); 00189 } 00190 00191 //-------------------------------------------------------------------------- 00192 template<> 00193 double BasicVector3D<double>::angle(const BasicVector3D<double> & v) const { 00194 double cosa = 0; 00195 double ptot = mag()*v.mag(); 00196 if(ptot > 0) { 00197 cosa = dot(v)/ptot; 00198 if(cosa > 1) cosa = 1; 00199 if(cosa < -1) cosa = -1; 00200 } 00201 return acos(cosa); 00202 } 00203 00204 //-------------------------------------------------------------------------- 00205 template<> 00206 BasicVector3D<double> & BasicVector3D<double>::rotateX(double a) { 00207 double sina = sin(a), cosa = cos(a), dy = y(), dz = z(); 00208 setY(dy*cosa-dz*sina); 00209 setZ(dz*cosa+dy*sina); 00210 return *this; 00211 } 00212 00213 //-------------------------------------------------------------------------- 00214 template<> 00215 BasicVector3D<double> & BasicVector3D<double>::rotateY(double a) { 00216 double sina = sin(a), cosa = cos(a), dz = z(), dx = x(); 00217 setZ(dz*cosa-dx*sina); 00218 setX(dx*cosa+dz*sina); 00219 return *this; 00220 } 00221 00222 //-------------------------------------------------------------------------- 00223 template<> 00224 BasicVector3D<double> & BasicVector3D<double>::rotateZ(double a) { 00225 double sina = sin(a), cosa = cos(a), dx = x(), dy = y(); 00226 setX(dx*cosa-dy*sina); 00227 setY(dy*cosa+dx*sina); 00228 return *this; 00229 } 00230 00231 //-------------------------------------------------------------------------- 00232 template<> 00233 BasicVector3D<double> & 00234 BasicVector3D<double>::rotate(double a, const BasicVector3D<double> & v) { 00235 if (a == 0) return *this; 00236 double cx = v.x(), cy = v.y(), cz = v.z(); 00237 double ll = sqrt(cx*cx + cy*cy + cz*cz); 00238 if (ll == 0) { 00239 std::cerr << "BasicVector<double>::rotate() : zero axis" << std::endl; 00240 return *this; 00241 } 00242 double cosa = cos(a), sina = sin(a); 00243 cx /= ll; cy /= ll; cz /= ll; 00244 00245 double xx = cosa + (1-cosa)*cx*cx; 00246 double xy = (1-cosa)*cx*cy - sina*cz; 00247 double xz = (1-cosa)*cx*cz + sina*cy; 00248 00249 double yx = (1-cosa)*cy*cx + sina*cz; 00250 double yy = cosa + (1-cosa)*cy*cy; 00251 double yz = (1-cosa)*cy*cz - sina*cx; 00252 00253 double zx = (1-cosa)*cz*cx - sina*cy; 00254 double zy = (1-cosa)*cz*cy + sina*cx; 00255 double zz = cosa + (1-cosa)*cz*cz; 00256 00257 cx = x(); cy = y(); cz = z(); 00258 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz); 00259 return *this; 00260 } 00261 00262 //-------------------------------------------------------------------------- 00263 std::ostream & 00264 operator<< (std::ostream & os, const BasicVector3D<double> & a) 00265 { 00266 return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")"; 00267 } 00268 00269 //-------------------------------------------------------------------------- 00270 std::istream & 00271 operator>> (std::istream & is, BasicVector3D<double> & a) 00272 { 00273 // Required format is ( a, b, c ) that is, three numbers, preceded by 00274 // (, followed by ), and separated by commas. The three numbers are 00275 // taken as x, y, z. 00276 00277 double x, y, z; 00278 char c; 00279 00280 is >> std::ws >> c; 00281 // ws is defined to invoke eatwhite(istream & ) 00282 // see (Stroustrup gray book) page 333 and 345. 00283 if (is.fail() || c != '(' ) { 00284 std::cerr 00285 << "Could not find required opening parenthesis " 00286 << "in input of a BasicVector3D<double>" 00287 << std::endl; 00288 return is; 00289 } 00290 00291 is >> x >> std::ws >> c; 00292 if (is.fail() || c != ',' ) { 00293 std::cerr 00294 << "Could not find x value and required trailing comma " 00295 << "in input of a BasicVector3D<double>" 00296 << std::endl; 00297 return is; 00298 } 00299 00300 is >> y >> std::ws >> c; 00301 if (is.fail() || c != ',' ) { 00302 std::cerr 00303 << "Could not find y value and required trailing comma " 00304 << "in input of a BasicVector3D<double>" 00305 << std::endl; 00306 return is; 00307 } 00308 00309 is >> z >> std::ws >> c; 00310 if (is.fail() || c != ')' ) { 00311 std::cerr 00312 << "Could not find z value and required close parenthesis " 00313 << "in input of a BasicVector3D<double>" 00314 << std::endl; 00315 return is; 00316 } 00317 00318 a.setX(x); 00319 a.setY(y); 00320 a.setZ(z); 00321 return is; 00322 } 00323 } /* namespace HepGeom */