00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifdef GNUPRAGMA
00014 #pragma implementation
00015 #endif
00016
00017 #include "CLHEP/Vector/ThreeVector.h"
00018 #include "CLHEP/Units/PhysicalConstants.h"
00019
00020 #include <cmath>
00021 #include <iostream>
00022
00023 namespace CLHEP {
00024
00025 void Hep3Vector::setMag(double ma) {
00026 double factor = mag();
00027 if (factor == 0) {
00028 std::cerr << "Hep3Vector::setMag() - "
00029 << "zero vector can't be stretched" << std::endl;
00030 }else{
00031 factor = ma/factor;
00032 setX(x()*factor);
00033 setY(y()*factor);
00034 setZ(z()*factor);
00035 }
00036 }
00037
00038 double Hep3Vector::operator () (int i) const {
00039 switch(i) {
00040 case X:
00041 return x();
00042 case Y:
00043 return y();
00044 case Z:
00045 return z();
00046 default:
00047 std::cerr << "Hep3Vector::operator () - "
00048 << "Hep3Vector subscripting: bad index (" << i << ")"
00049 << std::endl;
00050 }
00051 return 0.;
00052 }
00053
00054 double & Hep3Vector::operator () (int i) {
00055 static double dummy;
00056 switch(i) {
00057 case X:
00058 return dx;
00059 case Y:
00060 return dy;
00061 case Z:
00062 return dz;
00063 default:
00064 std::cerr
00065 << "Hep3Vector::operator () - "
00066 << "Hep3Vector subscripting: bad index (" << i << ")"
00067 << std::endl;
00068 return dummy;
00069 }
00070 }
00071
00072 Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) {
00073
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 = std::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; }
00088 else {};
00089 return *this;
00090 }
00091
00092 double Hep3Vector::pseudoRapidity() const {
00093 double m1 = mag();
00094 if ( m1== 0 ) return 0.0;
00095 if ( m1== z() ) return 1.0E72;
00096 if ( m1== -z() ) return -1.0E72;
00097 return 0.5*std::log( (m1+z())/(m1-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 }
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
00121
00122
00123
00124 Hep3Vector & Hep3Vector::rotateX (double phi1) {
00125 double sinphi = std::sin(phi1);
00126 double cosphi = std::cos(phi1);
00127 double ty;
00128 ty = dy * cosphi - dz * sinphi;
00129 dz = dz * cosphi + dy * sinphi;
00130 dy = ty;
00131 return *this;
00132 }
00133
00134 Hep3Vector & Hep3Vector::rotateY (double phi1) {
00135 double sinphi = std::sin(phi1);
00136 double cosphi = std::cos(phi1);
00137 double tz;
00138 tz = dz * cosphi - dx * sinphi;
00139 dx = dx * cosphi + dz * sinphi;
00140 dz = tz;
00141 return *this;
00142 }
00143
00144 Hep3Vector & Hep3Vector::rotateZ (double phi1) {
00145 double sinphi = std::sin(phi1);
00146 double cosphi = std::cos(phi1);
00147 double tx;
00148 tx = dx * cosphi - dy * sinphi;
00149 dy = dy * cosphi + dx * sinphi;
00150 dx = tx;
00151 return *this;
00152 }
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 }
00158
00159 double Hep3Vector::howNear(const Hep3Vector & v ) const {
00160
00161 double d = (*this - v).mag2();
00162 double vdv = dot(v);
00163 if ( (vdv > 0) && (d < vdv) ) {
00164 return std::sqrt (d/vdv);
00165 } else if ( (vdv == 0) && (d == 0) ) {
00166 return 0;
00167 } else {
00168 return 1;
00169 }
00170 }
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 }
00181
00182 double Hep3Vector::deltaR ( const Hep3Vector & v ) const {
00183 double a = eta() - v.eta();
00184 double b = deltaPhi(v);
00185 return std::sqrt ( a*a + b*b );
00186 }
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)/std::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
00211
00212 if(arg > 1.0) arg = 1.0;
00213 }
00214 return arg;
00215 }
00216
00217 void Hep3Vector::setEta (double eta1) {
00218 double phi1 = 0;
00219 double r1;
00220 if ( (dx == 0) && (dy == 0) ) {
00221 if (dz == 0) {
00222 std::cerr << "Hep3Vector::setEta() - "
00223 << "Attempt to set eta of zero vector -- vector is unchanged"
00224 << std::endl;
00225 return;
00226 }
00227 std::cerr << "Hep3Vector::setEta() - "
00228 << "Attempt to set eta of vector along Z axis -- will use phi = 0"
00229 << std::endl;
00230 r1 = std::fabs(dz);
00231 } else {
00232 r1 = getR();
00233 phi1 = getPhi();
00234 }
00235 double tanHalfTheta = std::exp ( -eta1 );
00236 double cosTheta1 =
00237 (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
00238 dz = r1 * cosTheta1;
00239 double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
00240 dy = rho1 * std::sin (phi1);
00241 dx = rho1 * std::cos (phi1);
00242 return;
00243 }
00244
00245 void Hep3Vector::setCylTheta (double theta1) {
00246
00247
00248
00249 if ( (dx == 0) && (dy == 0) ) {
00250 if (dz == 0) {
00251 std::cerr << "Hep3Vector::setCylTheta() - "
00252 << "Attempt to set cylTheta of zero vector -- vector is unchanged"
00253 << std::endl;
00254 return;
00255 }
00256 if (theta1 == 0) {
00257 dz = std::fabs(dz);
00258 return;
00259 }
00260 if (theta1 == CLHEP::pi) {
00261 dz = -std::fabs(dz);
00262 return;
00263 }
00264 std::cerr << "Hep3Vector::setCylTheta() - "
00265 << "Attempt set cylindrical theta of vector along Z axis "
00266 << "to a non-trivial value, while keeping rho fixed -- "
00267 << "will return zero vector" << std::endl;
00268 dz = 0;
00269 return;
00270 }
00271 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
00272 std::cerr << "Hep3Vector::setCylTheta() - "
00273 << "Setting Cyl theta of a vector based on a value not in [0, PI]"
00274 << std::endl;
00275
00276 }
00277 double phi1 (getPhi());
00278 double rho1 = getRho();
00279 if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
00280 std::cerr << "Hep3Vector::setCylTheta() - "
00281 << "Attempt to set cylindrical theta to 0 or PI "
00282 << "while keeping rho fixed -- infinite Z will be computed"
00283 << std::endl;
00284 dz = (theta1==0) ? 1.0E72 : -1.0E72;
00285 return;
00286 }
00287 dz = rho1 / std::tan (theta1);
00288 dy = rho1 * std::sin (phi1);
00289 dx = rho1 * std::cos (phi1);
00290
00291 }
00292
00293 void Hep3Vector::setCylEta (double eta1) {
00294
00295
00296
00297 double theta1 = 2 * std::atan ( std::exp (-eta1) );
00298
00299
00300
00301
00302
00303
00304 if ( (dx == 0) && (dy == 0) ) {
00305 if (dz == 0) {
00306 std::cerr << "Hep3Vector::setCylEta() - "
00307 << "Attempt to set cylEta of zero vector -- vector is unchanged"
00308 << std::endl;
00309 return;
00310 }
00311 if (theta1 == 0) {
00312 dz = std::fabs(dz);
00313 return;
00314 }
00315 if (theta1 == CLHEP::pi) {
00316 dz = -std::fabs(dz);
00317 return;
00318 }
00319 std::cerr << "Hep3Vector::setCylEta() - "
00320 << "Attempt set cylindrical eta of vector along Z axis "
00321 << "to a non-trivial value, while keeping rho fixed -- "
00322 << "will return zero vector" << std::endl;
00323 dz = 0;
00324 return;
00325 }
00326 double phi1 (getPhi());
00327 double rho1 = getRho();
00328 dz = rho1 / std::tan (theta1);
00329 dy = rho1 * std::sin (phi1);
00330 dx = rho1 * std::cos (phi1);
00331
00332 }
00333
00334
00335 Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
00336
00337
00338
00339
00340
00341 double oneOverC = 1.0/c;
00342 return Hep3Vector ( v1.x() * oneOverC,
00343 v1.y() * oneOverC,
00344 v1.z() * oneOverC );
00345 }
00346
00347 Hep3Vector & Hep3Vector::operator/= (double c) {
00348
00349
00350
00351
00352
00353
00354 double oneOverC = 1.0/c;
00355 dx *= oneOverC;
00356 dy *= oneOverC;
00357 dz *= oneOverC;
00358 return *this;
00359 }
00360
00361 double Hep3Vector::tolerance = Hep3Vector::ToleranceTicks * 2.22045e-16;
00362
00363 }