Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- >> 2 // $Id:$ 2 // ------------------------------------------- 3 // --------------------------------------------------------------------------- 3 // 4 // 4 // This file is a part of the CLHEP - a Class 5 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 5 // 6 // 6 // This is the implementation of the Hep3Vecto 7 // This is the implementation of the Hep3Vector class. 7 // 8 // 8 // See also ThreeVectorR.cc for implementation 9 // See also ThreeVectorR.cc for implementation of Hep3Vector methods which 9 // would couple in all the HepRotation methods 10 // would couple in all the HepRotation methods. 10 // 11 // 11 12 >> 13 #ifdef GNUPRAGMA >> 14 #pragma implementation >> 15 #endif >> 16 12 #include "CLHEP/Vector/ThreeVector.h" 17 #include "CLHEP/Vector/ThreeVector.h" 13 #include "CLHEP/Units/PhysicalConstants.h" 18 #include "CLHEP/Units/PhysicalConstants.h" 14 19 15 #include <cmath> 20 #include <cmath> 16 #include <iostream> 21 #include <iostream> 17 22 18 namespace CLHEP { 23 namespace CLHEP { 19 24 20 void Hep3Vector::setMag(double ma) { 25 void Hep3Vector::setMag(double ma) { 21 double factor = mag(); 26 double factor = mag(); 22 if (factor == 0) { 27 if (factor == 0) { 23 std::cerr << "Hep3Vector::setMag() - " 28 std::cerr << "Hep3Vector::setMag() - " 24 << "zero vector can't be stretch 29 << "zero vector can't be stretched" << std::endl; 25 }else{ 30 }else{ 26 factor = ma/factor; 31 factor = ma/factor; 27 setX(x()*factor); 32 setX(x()*factor); 28 setY(y()*factor); 33 setY(y()*factor); 29 setZ(z()*factor); 34 setZ(z()*factor); 30 } 35 } 31 } 36 } 32 37 33 Hep3Vector & Hep3Vector::rotateUz(const Hep3Ve 38 Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) { 34 // NewUzVector must be normalized ! 39 // NewUzVector must be normalized ! 35 40 36 double u1 = NewUzVector.x(); 41 double u1 = NewUzVector.x(); 37 double u2 = NewUzVector.y(); 42 double u2 = NewUzVector.y(); 38 double u3 = NewUzVector.z(); 43 double u3 = NewUzVector.z(); 39 double up = u1*u1 + u2*u2; 44 double up = u1*u1 + u2*u2; 40 45 41 if (up > 0) { << 46 if (up>0) { 42 up = std::sqrt(up); << 47 up = std::sqrt(up); 43 double px = (u1 * u3 * x() - u2 * y()) / u << 48 double px = dx, py = dy, pz = dz; 44 double py = (u2 * u3 * x() + u1 * y()) / u << 49 dx = (u1*u3*px - u2*py)/up + u1*pz; 45 double pz = -up * x() + u3 * z(); << 50 dy = (u2*u3*px + u1*py)/up + u2*pz; 46 set(px, py, pz); << 51 dz = -up*px + u3*pz; 47 } else if (u3 < 0.) { << 52 } 48 setX(-x()); << 53 else if (u3 < 0.) { dx = -dx; dz = -dz; } // phi=0 teta=pi 49 setZ(-z()); << 54 else {}; 50 } // phi=0 teta=pi << 51 << 52 return *this; 55 return *this; 53 } 56 } 54 57 55 double Hep3Vector::pseudoRapidity() const { 58 double Hep3Vector::pseudoRapidity() const { 56 double m1 = mag(); 59 double m1 = mag(); 57 if ( m1== 0 ) return 0.0; 60 if ( m1== 0 ) return 0.0; 58 if ( m1== z() ) return 1.0E72; 61 if ( m1== z() ) return 1.0E72; 59 if ( m1== -z() ) return -1.0E72; 62 if ( m1== -z() ) return -1.0E72; 60 return 0.5*std::log( (m1+z())/(m1-z()) ); 63 return 0.5*std::log( (m1+z())/(m1-z()) ); 61 } 64 } 62 65 63 std::ostream & operator<< (std::ostream & os, 66 std::ostream & operator<< (std::ostream & os, const Hep3Vector & v) { 64 return os << "(" << v.x() << "," << v.y() << 67 return os << "(" << v.x() << "," << v.y() << "," << v.z() << ")"; 65 } 68 } 66 69 67 void ZMinput3doubles ( std::istream & is, cons 70 void ZMinput3doubles ( std::istream & is, const char * type, 68 double & x, double & y, 71 double & x, double & y, double & z ); 69 72 70 std::istream & operator>>(std::istream & is, H 73 std::istream & operator>>(std::istream & is, Hep3Vector & v) { 71 double x, y, z; 74 double x, y, z; 72 ZMinput3doubles ( is, "Hep3Vector", x, y, z 75 ZMinput3doubles ( is, "Hep3Vector", x, y, z ); 73 v.set(x, y, z); 76 v.set(x, y, z); 74 return is; 77 return is; 75 } // operator>>() 78 } // operator>>() 76 79 77 const Hep3Vector HepXHat(1.0, 0.0, 0.0); 80 const Hep3Vector HepXHat(1.0, 0.0, 0.0); 78 const Hep3Vector HepYHat(0.0, 1.0, 0.0); 81 const Hep3Vector HepYHat(0.0, 1.0, 0.0); 79 const Hep3Vector HepZHat(0.0, 0.0, 1.0); 82 const Hep3Vector HepZHat(0.0, 0.0, 1.0); 80 83 81 //------------------- 84 //------------------- 82 // 85 // 83 // New methods introduced when ZOOM PhysicsVec 86 // New methods introduced when ZOOM PhysicsVectors was merged in: 84 // 87 // 85 //------------------- 88 //------------------- 86 89 87 Hep3Vector & Hep3Vector::rotateX (double phi1) 90 Hep3Vector & Hep3Vector::rotateX (double phi1) { 88 double sinphi = std::sin(phi1); 91 double sinphi = std::sin(phi1); 89 double cosphi = std::cos(phi1); 92 double cosphi = std::cos(phi1); 90 double ty = y() * cosphi - z() * sinphi; << 93 double ty; 91 double tz = z() * cosphi + y() * sinphi; << 94 ty = dy * cosphi - dz * sinphi; 92 setY(ty); << 95 dz = dz * cosphi + dy * sinphi; 93 setZ(tz); << 96 dy = ty; 94 return *this; 97 return *this; 95 } /* rotateX */ 98 } /* rotateX */ 96 99 97 Hep3Vector & Hep3Vector::rotateY (double phi1) 100 Hep3Vector & Hep3Vector::rotateY (double phi1) { 98 double sinphi = std::sin(phi1); 101 double sinphi = std::sin(phi1); 99 double cosphi = std::cos(phi1); 102 double cosphi = std::cos(phi1); 100 double tx = x() * cosphi + z() * sinphi; << 103 double tz; 101 double tz = z() * cosphi - x() * sinphi; << 104 tz = dz * cosphi - dx * sinphi; 102 setX(tx); << 105 dx = dx * cosphi + dz * sinphi; 103 setZ(tz); << 106 dz = tz; 104 return *this; 107 return *this; 105 } /* rotateY */ 108 } /* rotateY */ 106 109 107 Hep3Vector & Hep3Vector::rotateZ (double phi1) 110 Hep3Vector & Hep3Vector::rotateZ (double phi1) { 108 double sinphi = std::sin(phi1); 111 double sinphi = std::sin(phi1); 109 double cosphi = std::cos(phi1); 112 double cosphi = std::cos(phi1); 110 double tx = x() * cosphi - y() * sinphi; << 113 double tx; 111 double ty = y() * cosphi + x() * sinphi; << 114 tx = dx * cosphi - dy * sinphi; 112 setX(tx); << 115 dy = dy * cosphi + dx * sinphi; 113 setY(ty); << 116 dx = tx; 114 return *this; 117 return *this; 115 } /* rotateZ */ 118 } /* rotateZ */ 116 119 117 bool Hep3Vector::isNear(const Hep3Vector & v, 120 bool Hep3Vector::isNear(const Hep3Vector & v, double epsilon) const { 118 double limit = dot(v)*epsilon*epsilon; 121 double limit = dot(v)*epsilon*epsilon; 119 return ( (*this - v).mag2() <= limit ); 122 return ( (*this - v).mag2() <= limit ); 120 } /* isNear() */ 123 } /* isNear() */ 121 124 122 double Hep3Vector::howNear(const Hep3Vector & 125 double Hep3Vector::howNear(const Hep3Vector & v ) const { 123 // | V1 - V2 | **2 / V1 dot V2, up to 1 126 // | V1 - V2 | **2 / V1 dot V2, up to 1 124 double d = (*this - v).mag2(); 127 double d = (*this - v).mag2(); 125 double vdv = dot(v); 128 double vdv = dot(v); 126 if ( (vdv > 0) && (d < vdv) ) { 129 if ( (vdv > 0) && (d < vdv) ) { 127 return std::sqrt (d/vdv); 130 return std::sqrt (d/vdv); 128 } else if ( (vdv == 0) && (d == 0) ) { 131 } else if ( (vdv == 0) && (d == 0) ) { 129 return 0; 132 return 0; 130 } else { 133 } else { 131 return 1; 134 return 1; 132 } 135 } 133 } /* howNear */ 136 } /* howNear */ 134 137 135 double Hep3Vector::deltaPhi (const Hep3Vector 138 double Hep3Vector::deltaPhi (const Hep3Vector & v2) const { 136 double dphi = v2.getPhi() - getPhi(); 139 double dphi = v2.getPhi() - getPhi(); 137 if ( dphi > CLHEP::pi ) { 140 if ( dphi > CLHEP::pi ) { 138 dphi -= CLHEP::twopi; 141 dphi -= CLHEP::twopi; 139 } else if ( dphi <= -CLHEP::pi ) { 142 } else if ( dphi <= -CLHEP::pi ) { 140 dphi += CLHEP::twopi; 143 dphi += CLHEP::twopi; 141 } 144 } 142 return dphi; 145 return dphi; 143 } /* deltaPhi */ 146 } /* deltaPhi */ 144 147 145 double Hep3Vector::deltaR ( const Hep3Vector & 148 double Hep3Vector::deltaR ( const Hep3Vector & v ) const { 146 double a = eta() - v.eta(); 149 double a = eta() - v.eta(); 147 double b = deltaPhi(v); 150 double b = deltaPhi(v); 148 return std::sqrt ( a*a + b*b ); 151 return std::sqrt ( a*a + b*b ); 149 } /* deltaR */ 152 } /* deltaR */ 150 153 151 double Hep3Vector::cosTheta(const Hep3Vector & 154 double Hep3Vector::cosTheta(const Hep3Vector & q) const { 152 double arg; 155 double arg; 153 double ptot2 = mag2()*q.mag2(); 156 double ptot2 = mag2()*q.mag2(); 154 if(ptot2 <= 0) { 157 if(ptot2 <= 0) { 155 arg = 0.0; 158 arg = 0.0; 156 }else{ 159 }else{ 157 arg = dot(q)/std::sqrt(ptot2); 160 arg = dot(q)/std::sqrt(ptot2); 158 if(arg > 1.0) arg = 1.0; 161 if(arg > 1.0) arg = 1.0; 159 if(arg < -1.0) arg = -1.0; 162 if(arg < -1.0) arg = -1.0; 160 } 163 } 161 return arg; 164 return arg; 162 } 165 } 163 166 164 double Hep3Vector::cos2Theta(const Hep3Vector 167 double Hep3Vector::cos2Theta(const Hep3Vector & q) const { 165 double arg; 168 double arg; 166 double ptot2 = mag2(); 169 double ptot2 = mag2(); 167 double qtot2 = q.mag2(); 170 double qtot2 = q.mag2(); 168 if ( ptot2 == 0 || qtot2 == 0 ) { 171 if ( ptot2 == 0 || qtot2 == 0 ) { 169 arg = 1.0; 172 arg = 1.0; 170 }else{ 173 }else{ 171 double pdq = dot(q); 174 double pdq = dot(q); 172 arg = (pdq/ptot2) * (pdq/qtot2); 175 arg = (pdq/ptot2) * (pdq/qtot2); 173 // More naive methods overflow on vect 176 // More naive methods overflow on vectors which can be squared 174 // but can't be raised to the 4th powe 177 // but can't be raised to the 4th power. 175 if(arg > 1.0) arg = 1.0; 178 if(arg > 1.0) arg = 1.0; 176 } 179 } 177 return arg; 180 return arg; 178 } 181 } 179 182 180 void Hep3Vector::setEta (double eta1) { 183 void Hep3Vector::setEta (double eta1) { 181 double phi1 = 0; 184 double phi1 = 0; 182 double r1; 185 double r1; 183 if ( (x() == 0) && (y() == 0) ) { << 186 if ( (dx == 0) && (dy == 0) ) { 184 if (z() == 0) { << 187 if (dz == 0) { 185 std::cerr << "Hep3Vector::setEta() - " 188 std::cerr << "Hep3Vector::setEta() - " 186 << "Attempt to set eta of zero 189 << "Attempt to set eta of zero vector -- vector is unchanged" 187 << std::endl; 190 << std::endl; 188 return; 191 return; 189 } 192 } 190 std::cerr << "Hep3Vector::setEta() - " 193 std::cerr << "Hep3Vector::setEta() - " 191 << "Attempt to set eta of vector a 194 << "Attempt to set eta of vector along Z axis -- will use phi = 0" 192 << std::endl; 195 << std::endl; 193 r1 = std::fabs(z()); << 196 r1 = std::fabs(dz); 194 } else { 197 } else { 195 r1 = getR(); 198 r1 = getR(); 196 phi1 = getPhi(); 199 phi1 = getPhi(); 197 } 200 } 198 double tanHalfTheta = std::exp ( -eta1 ); 201 double tanHalfTheta = std::exp ( -eta1 ); 199 double cosTheta1 = 202 double cosTheta1 = 200 (1 - tanHalfTheta*tanHalfTheta) / (1 + 203 (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta); >> 204 dz = r1 * cosTheta1; 201 double rho1 = r1*std::sqrt(1 - cosTheta1*cos 205 double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1); 202 setZ(r1 * cosTheta1); << 206 dy = rho1 * std::sin (phi1); 203 setY(rho1 * std::sin (phi1)); << 207 dx = rho1 * std::cos (phi1); 204 setX(rho1 * std::cos (phi1)); << 205 return; 208 return; 206 } 209 } 207 210 208 void Hep3Vector::setCylTheta (double theta1) { 211 void Hep3Vector::setCylTheta (double theta1) { 209 212 210 // In cylindrical coords, set theta while ke 213 // In cylindrical coords, set theta while keeping rho and phi fixed 211 214 212 if ( (x() == 0) && (y() == 0) ) { << 215 if ( (dx == 0) && (dy == 0) ) { 213 if (z() == 0) { << 216 if (dz == 0) { 214 std::cerr << "Hep3Vector::setCylTheta() 217 std::cerr << "Hep3Vector::setCylTheta() - " 215 << "Attempt to set cylTheta of 218 << "Attempt to set cylTheta of zero vector -- vector is unchanged" 216 << std::endl; 219 << std::endl; 217 return; 220 return; 218 } 221 } 219 if (theta1 == 0) { 222 if (theta1 == 0) { 220 setZ(std::fabs(z())); << 223 dz = std::fabs(dz); 221 return; 224 return; 222 } 225 } 223 if (theta1 == CLHEP::pi) { 226 if (theta1 == CLHEP::pi) { 224 setZ(-std::fabs(z())); << 227 dz = -std::fabs(dz); 225 return; 228 return; 226 } 229 } 227 std::cerr << "Hep3Vector::setCylTheta() - 230 std::cerr << "Hep3Vector::setCylTheta() - " 228 << "Attempt set cylindrical theta of vec 231 << "Attempt set cylindrical theta of vector along Z axis " 229 << "to a non-trivial value, while keepin 232 << "to a non-trivial value, while keeping rho fixed -- " 230 << "will return zero vector" << std::end 233 << "will return zero vector" << std::endl; 231 setZ(0.0); << 234 dz = 0; 232 return; 235 return; 233 } 236 } 234 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) 237 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) { 235 std::cerr << "Hep3Vector::setCylTheta() - 238 std::cerr << "Hep3Vector::setCylTheta() - " 236 << "Setting Cyl theta of a vector based 239 << "Setting Cyl theta of a vector based on a value not in [0, PI]" 237 << std::endl; 240 << std::endl; 238 // No special return needed if warning 241 // No special return needed if warning is ignored. 239 } 242 } 240 double phi1 (getPhi()); 243 double phi1 (getPhi()); 241 double rho1 = getRho(); 244 double rho1 = getRho(); 242 if ( (theta1 == 0) || (theta1 == CLHEP::pi) 245 if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) { 243 std::cerr << "Hep3Vector::setCylTheta() - 246 std::cerr << "Hep3Vector::setCylTheta() - " 244 << "Attempt to set cylindrical theta to 247 << "Attempt to set cylindrical theta to 0 or PI " 245 << "while keeping rho fixed -- infinite 248 << "while keeping rho fixed -- infinite Z will be computed" 246 << std::endl; 249 << std::endl; 247 setZ((theta1==0) ? 1.0E72 : -1.0E72); << 250 dz = (theta1==0) ? 1.0E72 : -1.0E72; 248 return; 251 return; 249 } 252 } 250 setZ(rho1 / std::tan (theta1)); << 253 dz = rho1 / std::tan (theta1); 251 setY(rho1 * std::sin (phi1)); << 254 dy = rho1 * std::sin (phi1); 252 setX(rho1 * std::cos (phi1)); << 255 dx = rho1 * std::cos (phi1); 253 256 254 } /* setCylTheta */ 257 } /* setCylTheta */ 255 258 256 void Hep3Vector::setCylEta (double eta1) { 259 void Hep3Vector::setCylEta (double eta1) { 257 260 258 // In cylindrical coords, set eta while keep 261 // In cylindrical coords, set eta while keeping rho and phi fixed 259 262 260 double theta1 = 2 * std::atan ( std::exp (-e 263 double theta1 = 2 * std::atan ( std::exp (-eta1) ); 261 264 262 //-| The remaining code is similar to 265 //-| The remaining code is similar to setCylTheta, The reason for 263 //-| using a copy is so as to be able 266 //-| using a copy is so as to be able to change the messages in the 264 //-| ZMthrows to say eta rather than t 267 //-| ZMthrows to say eta rather than theta. Besides, we assumedly 265 //-| need not check for theta of 0 or 268 //-| need not check for theta of 0 or PI. 266 269 267 if ( (x() == 0) && (y() == 0) ) { << 270 if ( (dx == 0) && (dy == 0) ) { 268 if (z() == 0) { << 271 if (dz == 0) { 269 std::cerr << "Hep3Vector::setCylEta() - 272 std::cerr << "Hep3Vector::setCylEta() - " 270 << "Attempt to set cylEta of zero vect 273 << "Attempt to set cylEta of zero vector -- vector is unchanged" 271 << std::endl; 274 << std::endl; 272 return; 275 return; 273 } 276 } 274 if (theta1 == 0) { 277 if (theta1 == 0) { 275 setZ(std::fabs(z())); << 278 dz = std::fabs(dz); 276 return; 279 return; 277 } 280 } 278 if (theta1 == CLHEP::pi) { 281 if (theta1 == CLHEP::pi) { 279 setZ(-std::fabs(z())); << 282 dz = -std::fabs(dz); 280 return; 283 return; 281 } 284 } 282 std::cerr << "Hep3Vector::setCylEta() - " 285 std::cerr << "Hep3Vector::setCylEta() - " 283 << "Attempt set cylindrical eta of vecto 286 << "Attempt set cylindrical eta of vector along Z axis " 284 << "to a non-trivial value, while keepin 287 << "to a non-trivial value, while keeping rho fixed -- " 285 << "will return zero vector" << std::end 288 << "will return zero vector" << std::endl; 286 setZ(0.0); << 289 dz = 0; 287 return; 290 return; 288 } 291 } 289 double phi1 (getPhi()); 292 double phi1 (getPhi()); 290 double rho1 = getRho(); 293 double rho1 = getRho(); 291 setZ(rho1 / std::tan (theta1)); << 294 dz = rho1 / std::tan (theta1); 292 setY(rho1 * std::sin (phi1)); << 295 dy = rho1 * std::sin (phi1); 293 setX(rho1 * std::cos (phi1)); << 296 dx = rho1 * std::cos (phi1); 294 297 295 } /* setCylEta */ 298 } /* setCylEta */ 296 299 297 300 298 Hep3Vector operator/ ( const Hep3Vector & v1, << 301 Hep3Vector operator/ ( const Hep3Vector & v1, double c ) { 299 // if (c == 0) { 302 // if (c == 0) { 300 // std::cerr << "Hep3Vector::operator/ () - 303 // std::cerr << "Hep3Vector::operator/ () - " 301 // << "Attempt to divide vector by 0 -- " 304 // << "Attempt to divide vector by 0 -- " 302 // << "will produce infinities and/or NAN 305 // << "will produce infinities and/or NANs" << std::endl; 303 // } 306 // } 304 return v1 * (1.0/c); << 307 double oneOverC = 1.0/c; >> 308 return Hep3Vector ( v1.x() * oneOverC, >> 309 v1.y() * oneOverC, >> 310 v1.z() * oneOverC ); 305 } /* v / c */ 311 } /* v / c */ 306 312 307 Hep3Vector & Hep3Vector::operator/= (double c) 313 Hep3Vector & Hep3Vector::operator/= (double c) { 308 // if (c == 0) { 314 // if (c == 0) { 309 // std::cerr << "Hep3Vector::operator/ () - 315 // std::cerr << "Hep3Vector::operator/ () - " 310 // << "Attempt to do vector /= 0 -- " 316 // << "Attempt to do vector /= 0 -- " 311 // << "division by zero would produce inf 317 // << "division by zero would produce infinite or NAN components" 312 // << std::endl; 318 // << std::endl; 313 // } 319 // } 314 *this *= 1.0/c; << 320 double oneOverC = 1.0/c; >> 321 dx *= oneOverC; >> 322 dy *= oneOverC; >> 323 dz *= oneOverC; 315 return *this; 324 return *this; 316 } 325 } 317 326 318 double Hep3Vector::tolerance = Hep3Vector::Tol 327 double Hep3Vector::tolerance = Hep3Vector::ToleranceTicks * 2.22045e-16; 319 328 320 } // namespace CLHEP 329 } // namespace CLHEP 321 330