CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: ThreeVector.cc,v 1.3 2003/08/13 20:00:14 garren Exp $ 00003 // --------------------------------------------------------------------------- 00004 // 00005 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00006 // 00007 // This is the implementation of the Hep3Vector class. 00008 // 00009 // See also ThreeVectorR.cc for implementation of Hep3Vector methods which 00010 // would couple in all the HepRotation methods. 00011 // 00012 00013 #ifdef GNUPRAGMA 00014 #pragma implementation 00015 #endif 00016 00017 #include "CLHEP/Vector/defs.h" 00018 #include "CLHEP/Vector/ThreeVector.h" 00019 #include "CLHEP/Vector/ZMxpv.h" 00020 #include "CLHEP/Units/PhysicalConstants.h" 00021 00022 #include <cmath> 00023 #include <iostream> 00024 00025 namespace CLHEP { 00026 00027 void Hep3Vector::setMag(double ma) { 00028 double factor = mag(); 00029 if (factor == 0) { 00030 ZMthrowA ( ZMxpvZeroVector ( 00031 "Hep3Vector::setMag : zero vector can't be stretched")); 00032 }else{ 00033 factor = ma/factor; 00034 setX(x()*factor); 00035 setY(y()*factor); 00036 setZ(z()*factor); 00037 } 00038 } 00039 00040 double Hep3Vector::operator () (int i) const { 00041 switch(i) { 00042 case X: 00043 return x(); 00044 case Y: 00045 return y(); 00046 case Z: 00047 return z(); 00048 default: 00049 std::cerr << "Hep3Vector subscripting: bad index (" << i << ")" 00050 << std::endl; 00051 } 00052 return 0.; 00053 } 00054 00055 double & Hep3Vector::operator () (int i) { 00056 static double dummy; 00057 switch(i) { 00058 case X: 00059 return dx; 00060 case Y: 00061 return dy; 00062 case Z: 00063 return dz; 00064 default: 00065 std::cerr 00066 << "Hep3Vector subscripting: bad index (" << i << ")" 00067 << std::endl; 00068 return dummy; 00069 } 00070 } 00071 00072 Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) { 00073 // NewUzVector must be normalized ! 00074 00075 double u1 = NewUzVector.x(); 00076 double u2 = NewUzVector.y(); 00077 double u3 = NewUzVector.z(); 00078 double up = u1*u1 + u2*u2; 00079 00080 if (up>0) { 00081 up = sqrt(up); 00082 double px = dx, py = dy, pz = dz; 00083 dx = (u1*u3*px - u2*py)/up + u1*pz; 00084 dy = (u2*u3*px + u1*py)/up + u2*pz; 00085 dz = -up*px + u3*pz; 00086 } 00087 else if (u3 < 0.) { dx = -dx; dz = -dz; } // phi=0 teta=pi 00088 else {}; 00089 return *this; 00090 } 00091 00092 double Hep3Vector::pseudoRapidity() const { 00093 double m = mag(); 00094 if ( m== 0 ) return 0.0; 00095 if ( m== z() ) return 1.0E72; 00096 if ( m== -z() ) return -1.0E72; 00097 return 0.5*log( (m+z())/(m-z()) ); 00098 } 00099 00100 std::ostream & operator<< (std::ostream & os, const Hep3Vector & v) { 00101 return os << "(" << v.x() << "," << v.y() << "," << v.z() << ")"; 00102 } 00103 00104 void ZMinput3doubles ( std::istream & is, const char * type, 00105 double & x, double & y, double & z ); 00106 00107 std::istream & operator>>(std::istream & is, Hep3Vector & v) { 00108 double x, y, z; 00109 ZMinput3doubles ( is, "Hep3Vector", x, y, z ); 00110 v.set(x, y, z); 00111 return is; 00112 } // operator>>() 00113 00114 const Hep3Vector HepXHat(1.0, 0.0, 0.0); 00115 const Hep3Vector HepYHat(0.0, 1.0, 0.0); 00116 const Hep3Vector HepZHat(0.0, 0.0, 1.0); 00117 00118 //------------------- 00119 // 00120 // New methods introduced when ZOOM PhysicsVectors was merged in: 00121 // 00122 //------------------- 00123 00124 Hep3Vector & Hep3Vector::rotateX (double phi) { 00125 double sinphi = sin(phi); 00126 double cosphi = cos(phi); 00127 double ty; 00128 ty = dy * cosphi - dz * sinphi; 00129 dz = dz * cosphi + dy * sinphi; 00130 dy = ty; 00131 return *this; 00132 } /* rotateX */ 00133 00134 Hep3Vector & Hep3Vector::rotateY (double phi) { 00135 double sinphi = sin(phi); 00136 double cosphi = cos(phi); 00137 double tz; 00138 tz = dz * cosphi - dx * sinphi; 00139 dx = dx * cosphi + dz * sinphi; 00140 dz = tz; 00141 return *this; 00142 } /* rotateY */ 00143 00144 Hep3Vector & Hep3Vector::rotateZ (double phi) { 00145 double sinphi = sin(phi); 00146 double cosphi = cos(phi); 00147 double tx; 00148 tx = dx * cosphi - dy * sinphi; 00149 dy = dy * cosphi + dx * sinphi; 00150 dx = tx; 00151 return *this; 00152 } /* rotateZ */ 00153 00154 bool Hep3Vector::isNear(const Hep3Vector & v, double epsilon) const { 00155 double limit = dot(v)*epsilon*epsilon; 00156 return ( (*this - v).mag2() <= limit ); 00157 } /* isNear() */ 00158 00159 double Hep3Vector::howNear(const Hep3Vector & v ) const { 00160 // | V1 - V2 | **2 / V1 dot V2, up to 1 00161 double d = (*this - v).mag2(); 00162 double vdv = dot(v); 00163 if ( (vdv > 0) && (d < vdv) ) { 00164 return sqrt (d/vdv); 00165 } else if ( (vdv == 0) && (d == 0) ) { 00166 return 0; 00167 } else { 00168 return 1; 00169 } 00170 } /* howNear */ 00171 00172 double Hep3Vector::deltaPhi (const Hep3Vector & v2) const { 00173 double dphi = v2.getPhi() - getPhi(); 00174 if ( dphi > CLHEP::pi ) { 00175 dphi -= CLHEP::twopi; 00176 } else if ( dphi <= -CLHEP::pi ) { 00177 dphi += CLHEP::twopi; 00178 } 00179 return dphi; 00180 } /* deltaPhi */ 00181 00182 double Hep3Vector::deltaR ( const Hep3Vector & v ) const { 00183 double a = eta() - v.eta(); 00184 double b = deltaPhi(v); 00185 return sqrt ( a*a + b*b ); 00186 } /* deltaR */ 00187 00188 double Hep3Vector::cosTheta(const Hep3Vector & q) const { 00189 double arg; 00190 double ptot2 = mag2()*q.mag2(); 00191 if(ptot2 <= 0) { 00192 arg = 0.0; 00193 }else{ 00194 arg = dot(q)/sqrt(ptot2); 00195 if(arg > 1.0) arg = 1.0; 00196 if(arg < -1.0) arg = -1.0; 00197 } 00198 return arg; 00199 } 00200 00201 double Hep3Vector::cos2Theta(const Hep3Vector & q) const { 00202 double arg; 00203 double ptot2 = mag2(); 00204 double qtot2 = q.mag2(); 00205 if ( ptot2 == 0 || qtot2 == 0 ) { 00206 arg = 1.0; 00207 }else{ 00208 double pdq = dot(q); 00209 arg = (pdq/ptot2) * (pdq/qtot2); 00210 // More naive methods overflow on vectors which can be squared 00211 // but can't be raised to the 4th power. 00212 if(arg > 1.0) arg = 1.0; 00213 } 00214 return arg; 00215 } 00216 00217 void Hep3Vector::setEta (double eta) { 00218 double phi = 0; 00219 double r; 00220 if ( (dx == 0) && (dy == 0) ) { 00221 if (dz == 0) { 00222 ZMthrowC (ZMxpvZeroVector( 00223 "Attempt to set eta of zero vector -- vector is unchanged")); 00224 return; 00225 } 00226 ZMthrowC (ZMxpvZeroVector( 00227 "Attempt to set eta of vector along Z axis -- will use phi = 0")); 00228 r = fabs(dz); 00229 } else { 00230 r = getR(); 00231 phi = getPhi(); 00232 } 00233 double tanHalfTheta = exp ( -eta ); 00234 double cosTheta = 00235 (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta); 00236 dz = r * cosTheta; 00237 double rho = r*sqrt(1 - cosTheta*cosTheta); 00238 dy = rho * sin (phi); 00239 dx = rho * cos (phi); 00240 return; 00241 } 00242 00243 void Hep3Vector::setCylTheta (double theta) { 00244 00245 // In cylindrical coords, set theta while keeping rho and phi fixed 00246 00247 if ( (dx == 0) && (dy == 0) ) { 00248 if (dz == 0) { 00249 ZMthrowC (ZMxpvZeroVector( 00250 "Attempt to set cylTheta of zero vector -- vector is unchanged")); 00251 return; 00252 } 00253 if (theta == 0) { 00254 dz = fabs(dz); 00255 return; 00256 } 00257 if (theta == CLHEP::pi) { 00258 dz = -fabs(dz); 00259 return; 00260 } 00261 ZMthrowC (ZMxpvZeroVector( 00262 "Attempt set cylindrical theta of vector along Z axis " 00263 "to a non-trivial value, while keeping rho fixed -- " 00264 "will return zero vector")); 00265 dz = 0; 00266 return; 00267 } 00268 if ( (theta < 0) || (theta > CLHEP::pi) ) { 00269 ZMthrowC (ZMxpvUnusualTheta( 00270 "Setting Cyl theta of a vector based on a value not in [0, PI]")); 00271 // No special return needed if warning is ignored. 00272 } 00273 double phi (getPhi()); 00274 double rho = getRho(); 00275 if ( (theta == 0) || (theta == CLHEP::pi) ) { 00276 ZMthrowC (ZMxpvInfiniteVector( 00277 "Attempt to set cylindrical theta to 0 or PI " 00278 "while keeping rho fixed -- infinite Z will be computed")); 00279 dz = (theta==0) ? 1.0E72 : -1.0E72; 00280 return; 00281 } 00282 dz = rho / tan (theta); 00283 dy = rho * sin (phi); 00284 dx = rho * cos (phi); 00285 00286 } /* setCylTheta */ 00287 00288 void Hep3Vector::setCylEta (double eta) { 00289 00290 // In cylindrical coords, set eta while keeping rho and phi fixed 00291 00292 double theta = 2 * atan ( exp (-eta) ); 00293 00294 //-| The remaining code is similar to setCylTheta, The reason for 00295 //-| using a copy is so as to be able to change the messages in the 00296 //-| ZMthrows to say eta rather than theta. Besides, we assumedly 00297 //-| need not check for theta of 0 or PI. 00298 00299 if ( (dx == 0) && (dy == 0) ) { 00300 if (dz == 0) { 00301 ZMthrowC (ZMxpvZeroVector( 00302 "Attempt to set cylEta of zero vector -- vector is unchanged")); 00303 return; 00304 } 00305 if (theta == 0) { 00306 dz = fabs(dz); 00307 return; 00308 } 00309 if (theta == CLHEP::pi) { 00310 dz = -fabs(dz); 00311 return; 00312 } 00313 ZMthrowC (ZMxpvZeroVector( 00314 "Attempt set cylindrical eta of vector along Z axis " 00315 "to a non-trivial value, while keeping rho fixed -- " 00316 "will return zero vector")); 00317 dz = 0; 00318 return; 00319 } 00320 double phi (getPhi()); 00321 double rho = getRho(); 00322 dz = rho / tan (theta); 00323 dy = rho * sin (phi); 00324 dx = rho * cos (phi); 00325 00326 } /* setCylEta */ 00327 00328 00329 Hep3Vector operator/ ( const Hep3Vector & v1, double c ) { 00330 if (c == 0) { 00331 ZMthrowA ( ZMxpvInfiniteVector ( 00332 "Attempt to divide vector by 0 -- " 00333 "will produce infinities and/or NANs")); 00334 } 00335 double oneOverC = 1.0/c; 00336 return Hep3Vector ( v1.x() * oneOverC, 00337 v1.y() * oneOverC, 00338 v1.z() * oneOverC ); 00339 } /* v / c */ 00340 00341 Hep3Vector & Hep3Vector::operator/= (double c) { 00342 if (c == 0) { 00343 ZMthrowA (ZMxpvInfiniteVector( 00344 "Attempt to do vector /= 0 -- " 00345 "division by zero would produce infinite or NAN components")); 00346 } 00347 double oneOverC = 1.0/c; 00348 dx *= oneOverC; 00349 dy *= oneOverC; 00350 dz *= oneOverC; 00351 return *this; 00352 } 00353 00354 double Hep3Vector::tolerance = Hep3Vector::ToleranceTicks * 2.22045e-16; 00355 00356 } // namespace CLHEP 00357 00358 00359