00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifdef GNUPRAGMA
00018 #pragma implementation
00019 #endif
00020
00021 #include "CLHEP/Vector/ThreeVector.h"
00022 #include "CLHEP/Units/PhysicalConstants.h"
00023
00024 #include <cmath>
00025
00026 namespace CLHEP {
00027
00028
00029
00030
00031
00032
00033
00034
00035 void Hep3Vector::setSpherical (
00036 double r1,
00037 double theta1,
00038 double phi1) {
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 dz = r1 * std::cos(theta1);
00050 double rho1 ( r1*std::sin(theta1));
00051 dy = rho1 * std::sin (phi1);
00052 dx = rho1 * std::cos (phi1);
00053 return;
00054 }
00055
00056 void Hep3Vector::setCylindrical (
00057 double rho1,
00058 double phi1,
00059 double z1) {
00060
00061
00062
00063
00064
00065 dz = z1;
00066 dy = rho1 * std::sin (phi1);
00067 dx = rho1 * std::cos (phi1);
00068 return;
00069 }
00070
00071 void Hep3Vector::setRhoPhiTheta (
00072 double rho1,
00073 double phi1,
00074 double theta1) {
00075 if (rho1 == 0) {
00076 std::cerr << "Hep3Vector::setRhoPhiTheta() - "
00077 << "Attempt set vector components rho, phi, theta with zero rho -- "
00078 << "zero vector is returned, ignoring theta and phi" << std::endl;
00079 dx = 0; dy = 0; dz = 0;
00080 return;
00081 }
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092 dz = rho1 / std::tan (theta1);
00093 dy = rho1 * std::sin (phi1);
00094 dx = rho1 * std::cos (phi1);
00095 return;
00096 }
00097
00098 void Hep3Vector::setRhoPhiEta (
00099 double rho1,
00100 double phi1,
00101 double eta1 ) {
00102 if (rho1 == 0) {
00103 std::cerr << "Hep3Vector::setRhoPhiEta() - "
00104 << "Attempt set vector components rho, phi, eta with zero rho -- "
00105 << "zero vector is returned, ignoring eta and phi" << std::endl;
00106 dx = 0; dy = 0; dz = 0;
00107 return;
00108 }
00109 double theta1 (2 * std::atan ( std::exp (-eta1) ));
00110 dz = rho1 / std::tan (theta1);
00111 dy = rho1 * std::sin (phi1);
00112 dx = rho1 * std::cos (phi1);
00113 return;
00114 }
00115
00116
00117
00118
00119
00120
00121
00122 int Hep3Vector::compare (const Hep3Vector & v) const {
00123 if ( dz > v.dz ) {
00124 return 1;
00125 } else if ( dz < v.dz ) {
00126 return -1;
00127 } else if ( dy > v.dy ) {
00128 return 1;
00129 } else if ( dy < v.dy ) {
00130 return -1;
00131 } else if ( dx > v.dx ) {
00132 return 1;
00133 } else if ( dx < v.dx ) {
00134 return -1;
00135 } else {
00136 return 0;
00137 }
00138 }
00139
00140
00141 bool Hep3Vector::operator > (const Hep3Vector & v) const {
00142 return (compare(v) > 0);
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
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 double Hep3Vector::howParallel (const Hep3Vector & v) const {
00169
00170 double v1v2 = std::fabs(dot(v));
00171 if ( v1v2 == 0 ) {
00172
00173 return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1;
00174 }
00175 Hep3Vector v1Xv2 ( cross(v) );
00176 double abscross = v1Xv2.mag();
00177 if ( abscross >= v1v2 ) {
00178 return 1;
00179 } else {
00180 return abscross/v1v2;
00181 }
00182 }
00183
00184 bool Hep3Vector::isParallel (const Hep3Vector & v,
00185 double epsilon) const {
00186
00187
00188
00189 static const double TOOBIG = std::pow(2.0,507);
00190 static const double SCALE = std::pow(2.0,-507);
00191 double v1v2 = std::fabs(dot(v));
00192 if ( v1v2 == 0 ) {
00193 return ( (mag2() == 0) && (v.mag2() == 0) );
00194 }
00195 if ( v1v2 >= TOOBIG ) {
00196 Hep3Vector sv1 ( *this * SCALE );
00197 Hep3Vector sv2 ( v * SCALE );
00198 Hep3Vector sv1Xsv2 = sv1.cross(sv2);
00199 double x2 = sv1Xsv2.mag2();
00200 double limit = v1v2*SCALE*SCALE;
00201 limit = epsilon*epsilon*limit*limit;
00202 return ( x2 <= limit );
00203 }
00204
00205
00206
00207 Hep3Vector v1Xv2 ( cross(v) );
00208 if ( (std::fabs (v1Xv2.dx) > TOOBIG) ||
00209 (std::fabs (v1Xv2.dy) > TOOBIG) ||
00210 (std::fabs (v1Xv2.dz) > TOOBIG) ) {
00211 return false;
00212 }
00213
00214 return ( (v1Xv2.mag2()) <= ((epsilon * v1v2) * (epsilon * v1v2)) );
00215
00216 }
00217
00218
00219 double Hep3Vector::howOrthogonal (const Hep3Vector & v) const {
00220
00221
00222 double v1v2 = std::fabs(dot(v));
00223
00224 if ( v1v2 == 0 ) {
00225 return 0;
00226 }
00227 Hep3Vector v1Xv2 ( cross(v) );
00228 double abscross = v1Xv2.mag();
00229 if ( v1v2 >= abscross ) {
00230 return 1;
00231 } else {
00232 return v1v2/abscross;
00233 }
00234
00235 }
00236
00237 bool Hep3Vector::isOrthogonal (const Hep3Vector & v,
00238 double epsilon) const {
00239
00240
00241
00242 static const double TOOBIG = std::pow(2.0,507);
00243 static const double SCALE = std::pow(2.0,-507);
00244 double v1v2 = std::fabs(dot(v));
00245
00246 if ( v1v2 >= TOOBIG ) {
00247 Hep3Vector sv1 ( *this * SCALE );
00248 Hep3Vector sv2 ( v * SCALE );
00249 Hep3Vector sv1Xsv2 = sv1.cross(sv2);
00250 double x2 = sv1Xsv2.mag2();
00251 double limit = epsilon*epsilon*x2;
00252 double y2 = v1v2*SCALE*SCALE;
00253 return ( y2*y2 <= limit );
00254 }
00255
00256
00257
00258 Hep3Vector eps_v1Xv2 ( cross(epsilon*v) );
00259 if ( (std::fabs (eps_v1Xv2.dx) > TOOBIG) ||
00260 (std::fabs (eps_v1Xv2.dy) > TOOBIG) ||
00261 (std::fabs (eps_v1Xv2.dz) > TOOBIG) ) {
00262 return true;
00263 }
00264
00265
00266
00267 return ( v1v2*v1v2 <= eps_v1Xv2.mag2() );
00268
00269 }
00270
00271 double Hep3Vector::setTolerance (double tol) {
00272
00273 double oldTolerance (tolerance);
00274 tolerance = tol;
00275 return oldTolerance;
00276 }
00277
00278
00279
00280
00281
00282
00283 double Hep3Vector::negativeInfinity() const {
00284
00285 struct Dib {
00286 union {
00287 double d;
00288 unsigned char i[8];
00289 } u;
00290 };
00291 Dib negOne;
00292 Dib posTwo;
00293 negOne.u.d = -1.0;
00294 posTwo.u.d = 2.0;
00295 Dib value;
00296 int k;
00297 for (k=0; k<8; k++) {
00298 value.u.i[k] = negOne.u.i[k] | posTwo.u.i[k];
00299 }
00300 return value.u.d;
00301 }
00302
00303 }