1 // -*- C++ -*- 2 // $Id: BasicVector3D.cc,v 1.3 2003/08/13 20:00:11 garren Exp $ 3 // --------------------------------------------------------------------------- 4 5 #include <math.h> 6 #include <float.h> 7 #include <iostream> 8 #include "CLHEP/Geometry/defs.h" 9 #include "CLHEP/Geometry/BasicVector3D.h" 10 11 namespace HepGeom { 12 //-------------------------------------------------------------------------- 13 template<> pseudoRapidity() const14 float BasicVector3D<float>::pseudoRapidity() const { 15 float ma = mag(), dz = z(); 16 if (ma == 0) return 0; 17 if (ma == dz) return FLT_MAX; 18 if (ma == -dz) return -FLT_MAX; 19 return 0.5*std::log((ma+dz)/(ma-dz)); 20 } 21 22 //-------------------------------------------------------------------------- 23 template<> setEta(float a)24 void BasicVector3D<float>::setEta(float a) { 25 double ma = mag(); 26 if (ma == 0) return; 27 double tanHalfTheta = std::exp(-a); 28 double tanHalfTheta2 = tanHalfTheta * tanHalfTheta; 29 double cosTheta1 = (1 - tanHalfTheta2) / (1 + tanHalfTheta2); 30 double rh = ma * std::sqrt(1 - cosTheta1*cosTheta1); 31 double ph = phi(); 32 set(rh*std::cos(ph), rh*std::sin(ph), ma*cosTheta1); 33 } 34 35 //-------------------------------------------------------------------------- 36 template<> angle(const BasicVector3D<float> & v) const37 float BasicVector3D<float>::angle(const BasicVector3D<float> & v) const { 38 double cosa = 0; 39 double ptot = mag()*v.mag(); 40 if(ptot > 0) { 41 cosa = dot(v)/ptot; 42 if(cosa > 1) cosa = 1; 43 if(cosa < -1) cosa = -1; 44 } 45 return std::acos(cosa); 46 } 47 48 //-------------------------------------------------------------------------- 49 template<> rotateX(float a)50 BasicVector3D<float> & BasicVector3D<float>::rotateX(float a) { 51 double sina = std::sin(a), cosa = std::cos(a), dy = y(), dz = z(); 52 setY(dy*cosa-dz*sina); 53 setZ(dz*cosa+dy*sina); 54 return *this; 55 } 56 57 //-------------------------------------------------------------------------- 58 template<> rotateY(float a)59 BasicVector3D<float> & BasicVector3D<float>::rotateY(float a) { 60 double sina = std::sin(a), cosa = std::cos(a), dz = z(), dx = x(); 61 setZ(dz*cosa-dx*sina); 62 setX(dx*cosa+dz*sina); 63 return *this; 64 } 65 66 //-------------------------------------------------------------------------- 67 template<> rotateZ(float a)68 BasicVector3D<float> & BasicVector3D<float>::rotateZ(float a) { 69 double sina = std::sin(a), cosa = std::cos(a), dx = x(), dy = y(); 70 setX(dx*cosa-dy*sina); 71 setY(dy*cosa+dx*sina); 72 return *this; 73 } 74 75 //-------------------------------------------------------------------------- 76 template<> 77 BasicVector3D<float> & rotate(float a,const BasicVector3D<float> & v)78 BasicVector3D<float>::rotate(float a, const BasicVector3D<float> & v) { 79 if (a == 0) return *this; 80 double cx = v.x(), cy = v.y(), cz = v.z(); 81 double ll = std::sqrt(cx*cx + cy*cy + cz*cz); 82 if (ll == 0) { 83 std::cerr << "BasicVector<float>::rotate() : zero axis" << std::endl; 84 return *this; 85 } 86 double cosa = std::cos(a), sina = std::sin(a); 87 cx /= ll; cy /= ll; cz /= ll; 88 89 double xx = cosa + (1-cosa)*cx*cx; 90 double xy = (1-cosa)*cx*cy - sina*cz; 91 double xz = (1-cosa)*cx*cz + sina*cy; 92 93 double yx = (1-cosa)*cy*cx + sina*cz; 94 double yy = cosa + (1-cosa)*cy*cy; 95 double yz = (1-cosa)*cy*cz - sina*cx; 96 97 double zx = (1-cosa)*cz*cx - sina*cy; 98 double zy = (1-cosa)*cz*cy + sina*cx; 99 double zz = cosa + (1-cosa)*cz*cz; 100 101 cx = x(); cy = y(); cz = z(); 102 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz); 103 return *this; 104 } 105 106 //-------------------------------------------------------------------------- 107 std::ostream & operator <<(std::ostream & os,const BasicVector3D<float> & a)108 operator<<(std::ostream & os, const BasicVector3D<float> & a) 109 { 110 return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")"; 111 } 112 113 //-------------------------------------------------------------------------- 114 std::istream & operator >>(std::istream & is,BasicVector3D<float> & a)115 operator>> (std::istream & is, BasicVector3D<float> & a) 116 { 117 // Required format is ( a, b, c ) that is, three numbers, preceded by 118 // (, followed by ), and separated by commas. The three numbers are 119 // taken as x, y, z. 120 121 float x, y, z; 122 char c; 123 124 is >> std::ws >> c; 125 // ws is defined to invoke eatwhite(istream & ) 126 // see (Stroustrup gray book) page 333 and 345. 127 if (is.fail() || c != '(' ) { 128 std::cerr 129 << "Could not find required opening parenthesis " 130 << "in input of a BasicVector3D<float>" 131 << std::endl; 132 return is; 133 } 134 135 is >> x >> std::ws >> c; 136 if (is.fail() || c != ',' ) { 137 std::cerr 138 << "Could not find x value and required trailing comma " 139 << "in input of a BasicVector3D<float>" 140 << std::endl; 141 return is; 142 } 143 144 is >> y >> std::ws >> c; 145 if (is.fail() || c != ',' ) { 146 std::cerr 147 << "Could not find y value and required trailing comma " 148 << "in input of a BasicVector3D<float>" 149 << std::endl; 150 return is; 151 } 152 153 is >> z >> std::ws >> c; 154 if (is.fail() || c != ')' ) { 155 std::cerr 156 << "Could not find z value and required close parenthesis " 157 << "in input of a BasicVector3D<float>" 158 << std::endl; 159 return is; 160 } 161 162 a.setX(x); 163 a.setY(y); 164 a.setZ(z); 165 return is; 166 } 167 168 //-------------------------------------------------------------------------- 169 template<> pseudoRapidity() const170 double BasicVector3D<double>::pseudoRapidity() const { 171 double ma = mag(), dz = z(); 172 if (ma == 0) return 0; 173 if (ma == dz) return DBL_MAX; 174 if (ma == -dz) return -DBL_MAX; 175 return 0.5*std::log((ma+dz)/(ma-dz)); 176 } 177 178 //-------------------------------------------------------------------------- 179 template<> setEta(double a)180 void BasicVector3D<double>::setEta(double a) { 181 double ma = mag(); 182 if (ma == 0) return; 183 double tanHalfTheta = std::exp(-a); 184 double tanHalfTheta2 = tanHalfTheta * tanHalfTheta; 185 double cosTheta1 = (1 - tanHalfTheta2) / (1 + tanHalfTheta2); 186 double rh = ma * std::sqrt(1 - cosTheta1*cosTheta1); 187 double ph = phi(); 188 set(rh*std::cos(ph), rh*std::sin(ph), ma*cosTheta1); 189 } 190 191 //-------------------------------------------------------------------------- 192 template<> angle(const BasicVector3D<double> & v) const193 double BasicVector3D<double>::angle(const BasicVector3D<double> & v) const { 194 double cosa = 0; 195 double ptot = mag()*v.mag(); 196 if(ptot > 0) { 197 cosa = dot(v)/ptot; 198 if(cosa > 1) cosa = 1; 199 if(cosa < -1) cosa = -1; 200 } 201 return std::acos(cosa); 202 } 203 204 //-------------------------------------------------------------------------- 205 template<> rotateX(double a)206 BasicVector3D<double> & BasicVector3D<double>::rotateX(double a) { 207 double sina = std::sin(a), cosa = std::cos(a), dy = y(), dz = z(); 208 setY(dy*cosa-dz*sina); 209 setZ(dz*cosa+dy*sina); 210 return *this; 211 } 212 213 //-------------------------------------------------------------------------- 214 template<> rotateY(double a)215 BasicVector3D<double> & BasicVector3D<double>::rotateY(double a) { 216 double sina = std::sin(a), cosa = std::cos(a), dz = z(), dx = x(); 217 setZ(dz*cosa-dx*sina); 218 setX(dx*cosa+dz*sina); 219 return *this; 220 } 221 222 //-------------------------------------------------------------------------- 223 template<> rotateZ(double a)224 BasicVector3D<double> & BasicVector3D<double>::rotateZ(double a) { 225 double sina = std::sin(a), cosa = std::cos(a), dx = x(), dy = y(); 226 setX(dx*cosa-dy*sina); 227 setY(dy*cosa+dx*sina); 228 return *this; 229 } 230 231 //-------------------------------------------------------------------------- 232 template<> 233 BasicVector3D<double> & rotate(double a,const BasicVector3D<double> & v)234 BasicVector3D<double>::rotate(double a, const BasicVector3D<double> & v) { 235 if (a == 0) return *this; 236 double cx = v.x(), cy = v.y(), cz = v.z(); 237 double ll = std::sqrt(cx*cx + cy*cy + cz*cz); 238 if (ll == 0) { 239 std::cerr << "BasicVector<double>::rotate() : zero axis" << std::endl; 240 return *this; 241 } 242 double cosa = std::cos(a), sina = std::sin(a); 243 cx /= ll; cy /= ll; cz /= ll; 244 245 double xx = cosa + (1-cosa)*cx*cx; 246 double xy = (1-cosa)*cx*cy - sina*cz; 247 double xz = (1-cosa)*cx*cz + sina*cy; 248 249 double yx = (1-cosa)*cy*cx + sina*cz; 250 double yy = cosa + (1-cosa)*cy*cy; 251 double yz = (1-cosa)*cy*cz - sina*cx; 252 253 double zx = (1-cosa)*cz*cx - sina*cy; 254 double zy = (1-cosa)*cz*cy + sina*cx; 255 double zz = cosa + (1-cosa)*cz*cz; 256 257 cx = x(); cy = y(); cz = z(); 258 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz); 259 return *this; 260 } 261 262 //-------------------------------------------------------------------------- 263 std::ostream & operator <<(std::ostream & os,const BasicVector3D<double> & a)264 operator<< (std::ostream & os, const BasicVector3D<double> & a) 265 { 266 return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")"; 267 } 268 269 //-------------------------------------------------------------------------- 270 std::istream & operator >>(std::istream & is,BasicVector3D<double> & a)271 operator>> (std::istream & is, BasicVector3D<double> & a) 272 { 273 // Required format is ( a, b, c ) that is, three numbers, preceded by 274 // (, followed by ), and separated by commas. The three numbers are 275 // taken as x, y, z. 276 277 double x, y, z; 278 char c; 279 280 is >> std::ws >> c; 281 // ws is defined to invoke eatwhite(istream & ) 282 // see (Stroustrup gray book) page 333 and 345. 283 if (is.fail() || c != '(' ) { 284 std::cerr 285 << "Could not find required opening parenthesis " 286 << "in input of a BasicVector3D<double>" 287 << std::endl; 288 return is; 289 } 290 291 is >> x >> std::ws >> c; 292 if (is.fail() || c != ',' ) { 293 std::cerr 294 << "Could not find x value and required trailing comma " 295 << "in input of a BasicVector3D<double>" 296 << std::endl; 297 return is; 298 } 299 300 is >> y >> std::ws >> c; 301 if (is.fail() || c != ',' ) { 302 std::cerr 303 << "Could not find y value and required trailing comma " 304 << "in input of a BasicVector3D<double>" 305 << std::endl; 306 return is; 307 } 308 309 is >> z >> std::ws >> c; 310 if (is.fail() || c != ')' ) { 311 std::cerr 312 << "Could not find z value and required close parenthesis " 313 << "in input of a BasicVector3D<double>" 314 << std::endl; 315 return is; 316 } 317 318 a.setX(x); 319 a.setY(y); 320 a.setZ(z); 321 return is; 322 } 323 } /* namespace HepGeom */ 324