Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- 2 // ------------------------------------------- 2 // --------------------------------------------------------------------------- 3 // 3 // 4 // This file is a part of the CLHEP - a Class 4 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 5 // 5 // 6 // SpaceVector 6 // SpaceVector 7 // 7 // 8 // This is the implementation of those methods 8 // This is the implementation of those methods of the Hep3Vector class which 9 // originated from the ZOOM SpaceVector class. 9 // originated from the ZOOM SpaceVector class. Several groups of these methods 10 // have been separated off into the following 10 // have been separated off into the following code units: 11 // 11 // 12 // SpaceVectorR.cc All methods involving rota 12 // SpaceVectorR.cc All methods involving rotation 13 // SpaceVectorD.cc All methods involving angl 13 // SpaceVectorD.cc All methods involving angle decomposition 14 // SpaceVectorP.cc Intrinsic properties and m 14 // SpaceVectorP.cc Intrinsic properties and methods involving second vector 15 // 15 // 16 16 >> 17 #ifdef GNUPRAGMA >> 18 #pragma implementation >> 19 #endif >> 20 17 #include "CLHEP/Vector/ThreeVector.h" 21 #include "CLHEP/Vector/ThreeVector.h" 18 #include "CLHEP/Units/PhysicalConstants.h" 22 #include "CLHEP/Units/PhysicalConstants.h" 19 23 20 #include <cmath> 24 #include <cmath> 21 25 22 namespace CLHEP { 26 namespace CLHEP { 23 27 24 //-***************************** 28 //-***************************** 25 // - 1 - 29 // - 1 - 26 // set (multiple components) 30 // set (multiple components) 27 // in various coordinate systems 31 // in various coordinate systems 28 // 32 // 29 //-***************************** 33 //-***************************** 30 34 31 void Hep3Vector::setSpherical ( 35 void Hep3Vector::setSpherical ( 32 double r1, 36 double r1, 33 double theta1, 37 double theta1, 34 double phi1) { 38 double phi1) { 35 // if ( r1 < 0 ) { 39 // if ( r1 < 0 ) { 36 // std::cerr << "Hep3Vector::setSpherical() 40 // std::cerr << "Hep3Vector::setSpherical() - " 37 // << "Spherical coordinates set with neg 41 // << "Spherical coordinates set with negative R" << std::endl; 38 // // No special return needed if warning i 42 // // No special return needed if warning is ignored. 39 // } 43 // } 40 // if ( (theta1 < 0) || (theta1 > CLHEP::pi) 44 // if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) { 41 // std::cerr << "Hep3Vector::setSpherical() 45 // std::cerr << "Hep3Vector::setSpherical() - " 42 // << "Spherical coordinates set with the 46 // << "Spherical coordinates set with theta not in [0, PI]" << std::endl; 43 // // No special return needed if warning is 47 // // No special return needed if warning is ignored. 44 // } 48 // } >> 49 dz = r1 * std::cos(theta1); 45 double rho1 ( r1*std::sin(theta1)); 50 double rho1 ( r1*std::sin(theta1)); 46 setZ(r1 * std::cos(theta1)); << 51 dy = rho1 * std::sin (phi1); 47 setY(rho1 * std::sin (phi1)); << 52 dx = rho1 * std::cos (phi1); 48 setX(rho1 * std::cos (phi1)); << 49 return; 53 return; 50 } /* setSpherical (r, theta1, phi1) */ 54 } /* setSpherical (r, theta1, phi1) */ 51 55 52 void Hep3Vector::setCylindrical ( 56 void Hep3Vector::setCylindrical ( 53 double rho1, 57 double rho1, 54 double phi1, 58 double phi1, 55 double z1) { 59 double z1) { 56 // if ( rho1 < 0 ) { 60 // if ( rho1 < 0 ) { 57 // std::cerr << "Hep3Vector::setCylindrical 61 // std::cerr << "Hep3Vector::setCylindrical() - " 58 // << "Cylindrical coordinates supplied w 62 // << "Cylindrical coordinates supplied with negative Rho" << std::endl; 59 // // No special return needed if warning i 63 // // No special return needed if warning is ignored. 60 // } 64 // } 61 setZ(z1); << 65 dz = z1; 62 setY(rho1 * std::sin (phi1)); << 66 dy = rho1 * std::sin (phi1); 63 setX(rho1 * std::cos (phi1)); << 67 dx = rho1 * std::cos (phi1); 64 return; 68 return; 65 } /* setCylindrical (r, phi, z) */ 69 } /* setCylindrical (r, phi, z) */ 66 70 67 void Hep3Vector::setRhoPhiTheta ( 71 void Hep3Vector::setRhoPhiTheta ( 68 double rho1, 72 double rho1, 69 double phi1, 73 double phi1, 70 double theta1) { 74 double theta1) { 71 if (rho1 == 0) { 75 if (rho1 == 0) { 72 std::cerr << "Hep3Vector::setRhoPhiTheta() 76 std::cerr << "Hep3Vector::setRhoPhiTheta() - " 73 << "Attempt set vector components rho, p 77 << "Attempt set vector components rho, phi, theta with zero rho -- " 74 << "zero vector is returned, ignoring th 78 << "zero vector is returned, ignoring theta and phi" << std::endl; 75 set(0.0, 0.0, 0.0); << 79 dx = 0; dy = 0; dz = 0; 76 return; 80 return; 77 } 81 } 78 // if ( (theta1 == 0) || (theta1 == CLHEP::pi 82 // if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) { 79 // std::cerr << "Hep3Vector::setRhoPhiTheta 83 // std::cerr << "Hep3Vector::setRhoPhiTheta() - " 80 // << "Attempt set cylindrical vector vec 84 // << "Attempt set cylindrical vector vector with finite rho and " 81 // << "theta along the Z axis: infinite 85 // << "theta along the Z axis: infinite Z would be computed" << std::endl; 82 // } 86 // } 83 // if ( (theta1 < 0) || (theta1 > CLHEP::pi) 87 // if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) { 84 // std::cerr << "Hep3Vector::setRhoPhiTheta 88 // std::cerr << "Hep3Vector::setRhoPhiTheta() - " 85 // << "Rho, phi, theta set with theta not 89 // << "Rho, phi, theta set with theta not in [0, PI]" << std::endl; 86 // // No special return needed if warning is 90 // // No special return needed if warning is ignored. 87 // } 91 // } 88 setZ(rho1 / std::tan (theta1)); << 92 dz = rho1 / std::tan (theta1); 89 setY(rho1 * std::sin (phi1)); << 93 dy = rho1 * std::sin (phi1); 90 setX(rho1 * std::cos (phi1)); << 94 dx = rho1 * std::cos (phi1); 91 return; 95 return; 92 } /* setCyl (rho, phi, theta) */ 96 } /* setCyl (rho, phi, theta) */ 93 97 94 void Hep3Vector::setRhoPhiEta ( 98 void Hep3Vector::setRhoPhiEta ( 95 double rho1, 99 double rho1, 96 double phi1, 100 double phi1, 97 double eta1 ) { 101 double eta1 ) { 98 if (rho1 == 0) { 102 if (rho1 == 0) { 99 std::cerr << "Hep3Vector::setRhoPhiEta() - 103 std::cerr << "Hep3Vector::setRhoPhiEta() - " 100 << "Attempt set vector components rho, p 104 << "Attempt set vector components rho, phi, eta with zero rho -- " 101 << "zero vector is returned, ignoring et 105 << "zero vector is returned, ignoring eta and phi" << std::endl; 102 set(0.0, 0.0, 0.0); << 106 dx = 0; dy = 0; dz = 0; 103 return; 107 return; 104 } 108 } 105 double theta1 (2 * std::atan ( std::exp (-et 109 double theta1 (2 * std::atan ( std::exp (-eta1) )); 106 setZ(rho1 / std::tan (theta1)); << 110 dz = rho1 / std::tan (theta1); 107 setY(rho1 * std::sin (phi1)); << 111 dy = rho1 * std::sin (phi1); 108 setX(rho1 * std::cos (phi1)); << 112 dx = rho1 * std::cos (phi1); 109 return; 113 return; 110 } /* setCyl (rho, phi, eta) */ 114 } /* setCyl (rho, phi, eta) */ 111 115 112 //************ 116 //************ 113 // - 3 - 117 // - 3 - 114 // Comparisons 118 // Comparisons 115 // 119 // 116 //************ 120 //************ 117 121 118 int Hep3Vector::compare (const Hep3Vector & v) 122 int Hep3Vector::compare (const Hep3Vector & v) const { 119 if ( z() > v.z() ) { << 123 if ( dz > v.dz ) { 120 return 1; 124 return 1; 121 } else if ( z() < v.z() ) { << 125 } else if ( dz < v.dz ) { 122 return -1; 126 return -1; 123 } else if ( y() > v.y() ) { << 127 } else if ( dy > v.dy ) { 124 return 1; 128 return 1; 125 } else if ( y() < v.y() ) { << 129 } else if ( dy < v.dy ) { 126 return -1; 130 return -1; 127 } else if ( x() > v.x() ) { << 131 } else if ( dx > v.dx ) { 128 return 1; 132 return 1; 129 } else if ( x() < v.x() ) { << 133 } else if ( dx < v.dx ) { 130 return -1; 134 return -1; 131 } else { 135 } else { 132 return 0; 136 return 0; 133 } 137 } 134 } /* Compare */ 138 } /* Compare */ 135 139 136 140 137 bool Hep3Vector::operator > (const Hep3Vector 141 bool Hep3Vector::operator > (const Hep3Vector & v) const { 138 return (compare(v) > 0); 142 return (compare(v) > 0); 139 } 143 } 140 bool Hep3Vector::operator < (const Hep3Vector 144 bool Hep3Vector::operator < (const Hep3Vector & v) const { 141 return (compare(v) < 0); 145 return (compare(v) < 0); 142 } 146 } 143 bool Hep3Vector::operator>= (const Hep3Vector 147 bool Hep3Vector::operator>= (const Hep3Vector & v) const { 144 return (compare(v) >= 0); 148 return (compare(v) >= 0); 145 } 149 } 146 bool Hep3Vector::operator<= (const Hep3Vector 150 bool Hep3Vector::operator<= (const Hep3Vector & v) const { 147 return (compare(v) <= 0); 151 return (compare(v) <= 0); 148 } 152 } 149 153 150 154 151 //-******** 155 //-******** 152 // Nearness 156 // Nearness 153 //-******** 157 //-******** 154 158 155 // These methods all assume you can safely tak 159 // These methods all assume you can safely take mag2() of each vector. 156 // Absolutely safe but slower and much uglier 160 // Absolutely safe but slower and much uglier alternatives were 157 // provided as build-time options in ZOOM Spac 161 // provided as build-time options in ZOOM SpaceVectors. 158 // Also, much smaller codes were provided tht 162 // Also, much smaller codes were provided tht assume you can square 159 // mag2() of each vector; but those return bad 163 // mag2() of each vector; but those return bad answers without warning 160 // when components exceed 10**90. 164 // when components exceed 10**90. 161 // 165 // 162 // IsNear, HowNear, and DeltaR are found in Th 166 // IsNear, HowNear, and DeltaR are found in ThreeVector.cc 163 167 164 double Hep3Vector::howParallel (const Hep3Vect 168 double Hep3Vector::howParallel (const Hep3Vector & v) const { 165 // | V1 x V2 | / | V1 dot V2 | 169 // | V1 x V2 | / | V1 dot V2 | 166 double v1v2 = std::fabs(dot(v)); 170 double v1v2 = std::fabs(dot(v)); 167 if ( v1v2 == 0 ) { 171 if ( v1v2 == 0 ) { 168 // Zero is parallel to no other vector exc 172 // Zero is parallel to no other vector except for zero. 169 return ( (mag2() == 0) && (v.mag2() == 0) 173 return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1; 170 } 174 } 171 Hep3Vector v1Xv2 ( cross(v) ); 175 Hep3Vector v1Xv2 ( cross(v) ); 172 double abscross = v1Xv2.mag(); 176 double abscross = v1Xv2.mag(); 173 if ( abscross >= v1v2 ) { 177 if ( abscross >= v1v2 ) { 174 return 1; 178 return 1; 175 } else { 179 } else { 176 return abscross/v1v2; 180 return abscross/v1v2; 177 } 181 } 178 } /* howParallel() */ 182 } /* howParallel() */ 179 183 180 bool Hep3Vector::isParallel (const Hep3Vector 184 bool Hep3Vector::isParallel (const Hep3Vector & v, 181 double epsilon) 185 double epsilon) const { 182 // | V1 x V2 | **2 <= epsilon **2 | V1 dot 186 // | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 | **2 183 // V1 is *this, V2 is v 187 // V1 is *this, V2 is v 184 188 185 static const double TOOBIG = std::pow(2.0,50 189 static const double TOOBIG = std::pow(2.0,507); 186 static const double SCALE = std::pow(2.0,-5 190 static const double SCALE = std::pow(2.0,-507); 187 double v1v2 = std::fabs(dot(v)); 191 double v1v2 = std::fabs(dot(v)); 188 if ( v1v2 == 0 ) { 192 if ( v1v2 == 0 ) { 189 return ( (mag2() == 0) && (v.mag2() == 0) 193 return ( (mag2() == 0) && (v.mag2() == 0) ); 190 } 194 } 191 if ( v1v2 >= TOOBIG ) { 195 if ( v1v2 >= TOOBIG ) { 192 Hep3Vector sv1 ( *this * SCALE ); 196 Hep3Vector sv1 ( *this * SCALE ); 193 Hep3Vector sv2 ( v * SCALE ); 197 Hep3Vector sv2 ( v * SCALE ); 194 Hep3Vector sv1Xsv2 = sv1.cross(sv2); 198 Hep3Vector sv1Xsv2 = sv1.cross(sv2); 195 double x2 = sv1Xsv2.mag2(); 199 double x2 = sv1Xsv2.mag2(); 196 double limit = v1v2*SCALE*SCALE; 200 double limit = v1v2*SCALE*SCALE; 197 limit = epsilon*epsilon*limit*limit; 201 limit = epsilon*epsilon*limit*limit; 198 return ( x2 <= limit ); 202 return ( x2 <= limit ); 199 } 203 } 200 204 201 // At this point we know v1v2 can be squared 205 // At this point we know v1v2 can be squared. 202 206 203 Hep3Vector v1Xv2 ( cross(v) ); 207 Hep3Vector v1Xv2 ( cross(v) ); 204 if ( (std::fabs (v1Xv2.x()) > TOOBIG) || << 208 if ( (std::fabs (v1Xv2.dx) > TOOBIG) || 205 (std::fabs (v1Xv2.y()) > TOOBIG) || << 209 (std::fabs (v1Xv2.dy) > TOOBIG) || 206 (std::fabs (v1Xv2.z()) > TOOBIG) ) { << 210 (std::fabs (v1Xv2.dz) > TOOBIG) ) { 207 return false; 211 return false; 208 } 212 } 209 213 210 return ( (v1Xv2.mag2()) <= ((epsilon * v1v2) 214 return ( (v1Xv2.mag2()) <= ((epsilon * v1v2) * (epsilon * v1v2)) ); 211 215 212 } /* isParallel() */ 216 } /* isParallel() */ 213 217 214 218 215 double Hep3Vector::howOrthogonal (const Hep3Ve 219 double Hep3Vector::howOrthogonal (const Hep3Vector & v) const { 216 // | V1 dot V2 | / | V1 x V2 | 220 // | V1 dot V2 | / | V1 x V2 | 217 221 218 double v1v2 = std::fabs(dot(v)); 222 double v1v2 = std::fabs(dot(v)); 219 //-| Safe because both v1 and v2 can be squa 223 //-| Safe because both v1 and v2 can be squared 220 if ( v1v2 == 0 ) { 224 if ( v1v2 == 0 ) { 221 return 0; // Even if one or both are 0, th 225 return 0; // Even if one or both are 0, they are considered orthogonal 222 } 226 } 223 Hep3Vector v1Xv2 ( cross(v) ); 227 Hep3Vector v1Xv2 ( cross(v) ); 224 double abscross = v1Xv2.mag(); 228 double abscross = v1Xv2.mag(); 225 if ( v1v2 >= abscross ) { 229 if ( v1v2 >= abscross ) { 226 return 1; 230 return 1; 227 } else { 231 } else { 228 return v1v2/abscross; 232 return v1v2/abscross; 229 } 233 } 230 234 231 } /* howOrthogonal() */ 235 } /* howOrthogonal() */ 232 236 233 bool Hep3Vector::isOrthogonal (const Hep3Vecto 237 bool Hep3Vector::isOrthogonal (const Hep3Vector & v, 234 double epsilon) const { 238 double epsilon) const { 235 // | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 239 // | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 | **2 236 // V1 is *this, V2 is v 240 // V1 is *this, V2 is v 237 241 238 static const double TOOBIG = std::pow(2.0,50 242 static const double TOOBIG = std::pow(2.0,507); 239 static const double SCALE = std::pow(2.0,-50 243 static const double SCALE = std::pow(2.0,-507); 240 double v1v2 = std::fabs(dot(v)); 244 double v1v2 = std::fabs(dot(v)); 241 //-| Safe because both v1 and v2 can b 245 //-| Safe because both v1 and v2 can be squared 242 if ( v1v2 >= TOOBIG ) { 246 if ( v1v2 >= TOOBIG ) { 243 Hep3Vector sv1 ( *this * SCALE ); 247 Hep3Vector sv1 ( *this * SCALE ); 244 Hep3Vector sv2 ( v * SCALE ); 248 Hep3Vector sv2 ( v * SCALE ); 245 Hep3Vector sv1Xsv2 = sv1.cross(sv2); 249 Hep3Vector sv1Xsv2 = sv1.cross(sv2); 246 double x2 = sv1Xsv2.mag2(); 250 double x2 = sv1Xsv2.mag2(); 247 double limit = epsilon*epsilon*x2; 251 double limit = epsilon*epsilon*x2; 248 double y2 = v1v2*SCALE*SCALE; 252 double y2 = v1v2*SCALE*SCALE; 249 return ( y2*y2 <= limit ); 253 return ( y2*y2 <= limit ); 250 } 254 } 251 255 252 // At this point we know v1v2 can be squared 256 // At this point we know v1v2 can be squared. 253 257 254 Hep3Vector eps_v1Xv2 ( cross(epsilon*v) ); 258 Hep3Vector eps_v1Xv2 ( cross(epsilon*v) ); 255 if ( (std::fabs (eps_v1Xv2.x()) > TOOBIG) | << 259 if ( (std::fabs (eps_v1Xv2.dx) > TOOBIG) || 256 (std::fabs (eps_v1Xv2.y()) > TOOBIG) | << 260 (std::fabs (eps_v1Xv2.dy) > TOOBIG) || 257 (std::fabs (eps_v1Xv2.z()) > TOOBIG) ) << 261 (std::fabs (eps_v1Xv2.dz) > TOOBIG) ) { 258 return true; 262 return true; 259 } 263 } 260 264 261 // At this point we know all the math we nee 265 // At this point we know all the math we need can be done. 262 266 263 return ( v1v2*v1v2 <= eps_v1Xv2.mag2() ); 267 return ( v1v2*v1v2 <= eps_v1Xv2.mag2() ); 264 268 265 } /* isOrthogonal() */ 269 } /* isOrthogonal() */ 266 270 267 double Hep3Vector::setTolerance (double tol) { 271 double Hep3Vector::setTolerance (double tol) { 268 // Set the tolerance for Hep3Vectors to be con 272 // Set the tolerance for Hep3Vectors to be considered near one another 269 double oldTolerance (tolerance); 273 double oldTolerance (tolerance); 270 tolerance = tol; 274 tolerance = tol; 271 return oldTolerance; 275 return oldTolerance; 272 } 276 } 273 277 274 //-*********************** 278 //-*********************** 275 // Helper Methods: 279 // Helper Methods: 276 // negativeInfinity() 280 // negativeInfinity() 277 //-*********************** 281 //-*********************** 278 282 279 double Hep3Vector::negativeInfinity() const { 283 double Hep3Vector::negativeInfinity() const { 280 // A byte-order-independent way to return -I 284 // A byte-order-independent way to return -Infinity 281 struct Dib { 285 struct Dib { 282 union { 286 union { 283 double d; 287 double d; 284 unsigned char i[8]; 288 unsigned char i[8]; 285 } u; 289 } u; 286 }; 290 }; 287 Dib negOne; 291 Dib negOne; 288 Dib posTwo; 292 Dib posTwo; 289 negOne.u.d = -1.0; 293 negOne.u.d = -1.0; 290 posTwo.u.d = 2.0; 294 posTwo.u.d = 2.0; 291 Dib value; 295 Dib value; 292 int k; 296 int k; 293 for (k=0; k<8; k++) { 297 for (k=0; k<8; k++) { 294 value.u.i[k] = negOne.u.i[k] | posTwo.u.i[ 298 value.u.i[k] = negOne.u.i[k] | posTwo.u.i[k]; 295 } 299 } 296 return value.u.d; 300 return value.u.d; 297 } 301 } 298 302 299 } // namespace CLHEP 303 } // namespace CLHEP 300 304