CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // --------------------------------------------------------------------------- 00003 // 00004 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00005 // 00006 // SpaceVector 00007 // 00008 // This is the implementation of those methods of the Hep3Vector class which 00009 // originated from the ZOOM SpaceVector class. Several groups of these methods 00010 // have been separated off into the following code units: 00011 // 00012 // SpaceVectorR.cc All methods involving rotation 00013 // SpaceVectorD.cc All methods involving angle decomposition 00014 // SpaceVectorP.cc Intrinsic properties and methods involving second vector 00015 // 00016 00017 #ifdef GNUPRAGMA 00018 #pragma implementation 00019 #endif 00020 00021 #include "CLHEP/Vector/defs.h" 00022 #include "CLHEP/Vector/ThreeVector.h" 00023 #include "CLHEP/Vector/ZMxpv.h" 00024 #include "CLHEP/Units/PhysicalConstants.h" 00025 00026 #include <cmath> 00027 00028 namespace CLHEP { 00029 00030 //-***************************** 00031 // - 1 - 00032 // set (multiple components) 00033 // in various coordinate systems 00034 // 00035 //-***************************** 00036 00037 void Hep3Vector::setSpherical ( 00038 double r1, 00039 double theta1, 00040 double phi1) { 00041 if ( r1 < 0 ) { 00042 ZMthrowC (ZMxpvNegativeR( 00043 "Spherical coordinates set with negative R")); 00044 // No special return needed if warning is ignored. 00045 } 00046 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) { 00047 ZMthrowC (ZMxpvUnusualTheta( 00048 "Spherical coordinates set with theta not in [0, PI]")); 00049 // No special return needed if warning is ignored. 00050 } 00051 dz = r1 * std::cos(theta1); 00052 double rho1 ( r1*std::sin(theta1)); 00053 dy = rho1 * std::sin (phi1); 00054 dx = rho1 * std::cos (phi1); 00055 return; 00056 } /* setSpherical (r, theta1, phi) */ 00057 00058 void Hep3Vector::setCylindrical ( 00059 double rho1, 00060 double phi1, 00061 double z1) { 00062 if ( rho1 < 0 ) { 00063 ZMthrowC (ZMxpvNegativeR( 00064 "Cylindrical coordinates supplied with negative Rho")); 00065 // No special return needed if warning is ignored. 00066 } 00067 dz = z1; 00068 dy = rho1 * std::sin (phi1); 00069 dx = rho1 * std::cos (phi1); 00070 return; 00071 } /* setCylindrical (r, phi, z) */ 00072 00073 void Hep3Vector::setRhoPhiTheta ( 00074 double rho1, 00075 double phi1, 00076 double theta1) { 00077 if (rho1 == 0) { 00078 ZMthrowC (ZMxpvZeroVector( 00079 "Attempt set vector components rho, phi, theta with zero rho -- " 00080 "zero vector is returned, ignoring theta and phi")); 00081 dx = 0; dy = 0; dz = 0; 00082 return; 00083 } 00084 if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) { 00085 ZMthrowA (ZMxpvInfiniteVector( 00086 "Attempt set cylindrical vector vector with finite rho and " 00087 "theta along the Z axis: infinite Z would be computed")); 00088 } 00089 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) { 00090 ZMthrowC (ZMxpvUnusualTheta( 00091 "Rho, phi, theta set with theta not in [0, PI]")); 00092 // No special return needed if warning is ignored. 00093 } 00094 dz = rho1 / std::tan (theta1); 00095 dy = rho1 * std::sin (phi1); 00096 dx = rho1 * std::cos (phi1); 00097 return; 00098 } /* setCyl (rho, phi, theta) */ 00099 00100 void Hep3Vector::setRhoPhiEta ( 00101 double rho1, 00102 double phi1, 00103 double eta1 ) { 00104 if (rho1 == 0) { 00105 ZMthrowC (ZMxpvZeroVector( 00106 "Attempt set vector components rho, phi, eta with zero rho -- " 00107 "zero vector is returned, ignoring eta and phi")); 00108 dx = 0; dy = 0; dz = 0; 00109 return; 00110 } 00111 double theta1 (2 * std::atan ( std::exp (-eta1) )); 00112 dz = rho1 / std::tan (theta1); 00113 dy = rho1 * std::sin (phi1); 00114 dx = rho1 * std::cos (phi1); 00115 return; 00116 } /* setCyl (rho, phi, eta) */ 00117 00118 00119 //************ 00120 // - 3 - 00121 // Comparisons 00122 // 00123 //************ 00124 00125 int Hep3Vector::compare (const Hep3Vector & v) const { 00126 if ( dz > v.dz ) { 00127 return 1; 00128 } else if ( dz < v.dz ) { 00129 return -1; 00130 } else if ( dy > v.dy ) { 00131 return 1; 00132 } else if ( dy < v.dy ) { 00133 return -1; 00134 } else if ( dx > v.dx ) { 00135 return 1; 00136 } else if ( dx < v.dx ) { 00137 return -1; 00138 } else { 00139 return 0; 00140 } 00141 } /* Compare */ 00142 00143 00144 bool Hep3Vector::operator > (const Hep3Vector & v) const { 00145 return (compare(v) > 0); 00146 } 00147 bool Hep3Vector::operator < (const Hep3Vector & v) const { 00148 return (compare(v) < 0); 00149 } 00150 bool Hep3Vector::operator>= (const Hep3Vector & v) const { 00151 return (compare(v) >= 0); 00152 } 00153 bool Hep3Vector::operator<= (const Hep3Vector & v) const { 00154 return (compare(v) <= 0); 00155 } 00156 00157 00158 00159 //-******** 00160 // Nearness 00161 //-******** 00162 00163 // These methods all assume you can safely take mag2() of each vector. 00164 // Absolutely safe but slower and much uglier alternatives were 00165 // provided as build-time options in ZOOM SpaceVectors. 00166 // Also, much smaller codes were provided tht assume you can square 00167 // mag2() of each vector; but those return bad answers without warning 00168 // when components exceed 10**90. 00169 // 00170 // IsNear, HowNear, and DeltaR are found in ThreeVector.cc 00171 00172 double Hep3Vector::howParallel (const Hep3Vector & v) const { 00173 // | V1 x V2 | / | V1 dot V2 | 00174 double v1v2 = std::fabs(dot(v)); 00175 if ( v1v2 == 0 ) { 00176 // Zero is parallel to no other vector except for zero. 00177 return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1; 00178 } 00179 Hep3Vector v1Xv2 ( cross(v) ); 00180 double abscross = v1Xv2.mag(); 00181 if ( abscross >= v1v2 ) { 00182 return 1; 00183 } else { 00184 return abscross/v1v2; 00185 } 00186 } /* howParallel() */ 00187 00188 bool Hep3Vector::isParallel (const Hep3Vector & v, 00189 double epsilon) const { 00190 // | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 | **2 00191 // V1 is *this, V2 is v 00192 00193 static const double TOOBIG = std::pow(2.0,507); 00194 static const double SCALE = std::pow(2.0,-507); 00195 double v1v2 = std::fabs(dot(v)); 00196 if ( v1v2 == 0 ) { 00197 return ( (mag2() == 0) && (v.mag2() == 0) ); 00198 } 00199 if ( v1v2 >= TOOBIG ) { 00200 Hep3Vector sv1 ( *this * SCALE ); 00201 Hep3Vector sv2 ( v * SCALE ); 00202 Hep3Vector sv1Xsv2 = sv1.cross(sv2); 00203 double x2 = sv1Xsv2.mag2(); 00204 double limit = v1v2*SCALE*SCALE; 00205 limit = epsilon*epsilon*limit*limit; 00206 return ( x2 <= limit ); 00207 } 00208 00209 // At this point we know v1v2 can be squared. 00210 00211 Hep3Vector v1Xv2 ( cross(v) ); 00212 if ( (std::fabs (v1Xv2.dx) > TOOBIG) || 00213 (std::fabs (v1Xv2.dy) > TOOBIG) || 00214 (std::fabs (v1Xv2.dz) > TOOBIG) ) { 00215 return false; 00216 } 00217 00218 return ( (v1Xv2.mag2()) <= ((epsilon * v1v2) * (epsilon * v1v2)) ); 00219 00220 } /* isParallel() */ 00221 00222 00223 double Hep3Vector::howOrthogonal (const Hep3Vector & v) const { 00224 // | V1 dot V2 | / | V1 x V2 | 00225 00226 double v1v2 = std::fabs(dot(v)); 00227 //-| Safe because both v1 and v2 can be squared 00228 if ( v1v2 == 0 ) { 00229 return 0; // Even if one or both are 0, they are considered orthogonal 00230 } 00231 Hep3Vector v1Xv2 ( cross(v) ); 00232 double abscross = v1Xv2.mag(); 00233 if ( v1v2 >= abscross ) { 00234 return 1; 00235 } else { 00236 return v1v2/abscross; 00237 } 00238 00239 } /* howOrthogonal() */ 00240 00241 bool Hep3Vector::isOrthogonal (const Hep3Vector & v, 00242 double epsilon) const { 00243 // | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 | **2 00244 // V1 is *this, V2 is v 00245 00246 static const double TOOBIG = std::pow(2.0,507); 00247 static const double SCALE = std::pow(2.0,-507); 00248 double v1v2 = std::fabs(dot(v)); 00249 //-| Safe because both v1 and v2 can be squared 00250 if ( v1v2 >= TOOBIG ) { 00251 Hep3Vector sv1 ( *this * SCALE ); 00252 Hep3Vector sv2 ( v * SCALE ); 00253 Hep3Vector sv1Xsv2 = sv1.cross(sv2); 00254 double x2 = sv1Xsv2.mag2(); 00255 double limit = epsilon*epsilon*x2; 00256 double y2 = v1v2*SCALE*SCALE; 00257 return ( y2*y2 <= limit ); 00258 } 00259 00260 // At this point we know v1v2 can be squared. 00261 00262 Hep3Vector eps_v1Xv2 ( cross(epsilon*v) ); 00263 if ( (std::fabs (eps_v1Xv2.dx) > TOOBIG) || 00264 (std::fabs (eps_v1Xv2.dy) > TOOBIG) || 00265 (std::fabs (eps_v1Xv2.dz) > TOOBIG) ) { 00266 return true; 00267 } 00268 00269 // At this point we know all the math we need can be done. 00270 00271 return ( v1v2*v1v2 <= eps_v1Xv2.mag2() ); 00272 00273 } /* isOrthogonal() */ 00274 00275 double Hep3Vector::setTolerance (double tol) { 00276 // Set the tolerance for Hep3Vectors to be considered near one another 00277 double oldTolerance (tolerance); 00278 tolerance = tol; 00279 return oldTolerance; 00280 } 00281 00282 00283 //-*********************** 00284 // Helper Methods: 00285 // negativeInfinity() 00286 //-*********************** 00287 00288 double Hep3Vector::negativeInfinity() const { 00289 // A byte-order-independent way to return -Infinity 00290 struct Dib { 00291 union { 00292 double d; 00293 unsigned char i[8]; 00294 } u; 00295 }; 00296 Dib negOne; 00297 Dib posTwo; 00298 negOne.u.d = -1.0; 00299 posTwo.u.d = 2.0; 00300 Dib value; 00301 int k; 00302 for (k=0; k<8; k++) { 00303 value.u.i[k] = negOne.u.i[k] | posTwo.u.i[k]; 00304 } 00305 return value.u.d; 00306 } 00307 00308 } // namespace CLHEP 00309 00310