CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

SpaceVector.cc
Go to the documentation of this file.
1// -*- C++ -*-
2// ---------------------------------------------------------------------------
3//
4// This file is a part of the CLHEP - a Class Library for High Energy Physics.
5//
6// SpaceVector
7//
8// This is the implementation of those methods of the Hep3Vector class which
9// originated from the ZOOM SpaceVector class. Several groups of these methods
10// have been separated off into the following code units:
11//
12// SpaceVectorR.cc All methods involving rotation
13// SpaceVectorD.cc All methods involving angle decomposition
14// SpaceVectorP.cc Intrinsic properties and methods involving second vector
15//
16
17#ifdef GNUPRAGMA
18#pragma implementation
19#endif
20
21#include "CLHEP/Vector/defs.h"
22#include "CLHEP/Vector/ThreeVector.h"
23#include "CLHEP/Vector/ZMxpv.h"
24#include "CLHEP/Units/PhysicalConstants.h"
25
26#include <cmath>
27
28namespace CLHEP {
29
30//-*****************************
31// - 1 -
32// set (multiple components)
33// in various coordinate systems
34//
35//-*****************************
36
38 double r1,
39 double theta1,
40 double phi1) {
41 if ( r1 < 0 ) {
42 ZMthrowC (ZMxpvNegativeR(
43 "Spherical coordinates set with negative R"));
44 // No special return needed if warning is ignored.
45 }
46 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
47 ZMthrowC (ZMxpvUnusualTheta(
48 "Spherical coordinates set with theta not in [0, PI]"));
49 // No special return needed if warning is ignored.
50 }
51 dz = r1 * std::cos(theta1);
52 double rho1 ( r1*std::sin(theta1));
53 dy = rho1 * std::sin (phi1);
54 dx = rho1 * std::cos (phi1);
55 return;
56} /* setSpherical (r, theta1, phi) */
57
59 double rho1,
60 double phi1,
61 double z1) {
62 if ( rho1 < 0 ) {
63 ZMthrowC (ZMxpvNegativeR(
64 "Cylindrical coordinates supplied with negative Rho"));
65 // No special return needed if warning is ignored.
66 }
67 dz = z1;
68 dy = rho1 * std::sin (phi1);
69 dx = rho1 * std::cos (phi1);
70 return;
71} /* setCylindrical (r, phi, z) */
72
74 double rho1,
75 double phi1,
76 double theta1) {
77 if (rho1 == 0) {
78 ZMthrowC (ZMxpvZeroVector(
79 "Attempt set vector components rho, phi, theta with zero rho -- "
80 "zero vector is returned, ignoring theta and phi"));
81 dx = 0; dy = 0; dz = 0;
82 return;
83 }
84 if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
85 ZMthrowA (ZMxpvInfiniteVector(
86 "Attempt set cylindrical vector vector with finite rho and "
87 "theta along the Z axis: infinite Z would be computed"));
88 }
89 if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
90 ZMthrowC (ZMxpvUnusualTheta(
91 "Rho, phi, theta set with theta not in [0, PI]"));
92 // No special return needed if warning is ignored.
93 }
94 dz = rho1 / std::tan (theta1);
95 dy = rho1 * std::sin (phi1);
96 dx = rho1 * std::cos (phi1);
97 return;
98} /* setCyl (rho, phi, theta) */
99
101 double rho1,
102 double phi1,
103 double eta1 ) {
104 if (rho1 == 0) {
105 ZMthrowC (ZMxpvZeroVector(
106 "Attempt set vector components rho, phi, eta with zero rho -- "
107 "zero vector is returned, ignoring eta and phi"));
108 dx = 0; dy = 0; dz = 0;
109 return;
110 }
111 double theta1 (2 * std::atan ( std::exp (-eta1) ));
112 dz = rho1 / std::tan (theta1);
113 dy = rho1 * std::sin (phi1);
114 dx = rho1 * std::cos (phi1);
115 return;
116} /* setCyl (rho, phi, eta) */
117
118
119//************
120// - 3 -
121// Comparisons
122//
123//************
124
125int Hep3Vector::compare (const Hep3Vector & v) const {
126 if ( dz > v.dz ) {
127 return 1;
128 } else if ( dz < v.dz ) {
129 return -1;
130 } else if ( dy > v.dy ) {
131 return 1;
132 } else if ( dy < v.dy ) {
133 return -1;
134 } else if ( dx > v.dx ) {
135 return 1;
136 } else if ( dx < v.dx ) {
137 return -1;
138 } else {
139 return 0;
140 }
141} /* Compare */
142
143
144bool Hep3Vector::operator > (const Hep3Vector & v) const {
145 return (compare(v) > 0);
146}
147bool Hep3Vector::operator < (const Hep3Vector & v) const {
148 return (compare(v) < 0);
149}
150bool Hep3Vector::operator>= (const Hep3Vector & v) const {
151 return (compare(v) >= 0);
152}
153bool Hep3Vector::operator<= (const Hep3Vector & v) const {
154 return (compare(v) <= 0);
155}
156
157
158
159//-********
160// Nearness
161//-********
162
163// These methods all assume you can safely take mag2() of each vector.
164// Absolutely safe but slower and much uglier alternatives were
165// provided as build-time options in ZOOM SpaceVectors.
166// Also, much smaller codes were provided tht assume you can square
167// mag2() of each vector; but those return bad answers without warning
168// when components exceed 10**90.
169//
170// IsNear, HowNear, and DeltaR are found in ThreeVector.cc
171
172double Hep3Vector::howParallel (const Hep3Vector & v) const {
173 // | V1 x V2 | / | V1 dot V2 |
174 double v1v2 = std::fabs(dot(v));
175 if ( v1v2 == 0 ) {
176 // Zero is parallel to no other vector except for zero.
177 return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1;
178 }
179 Hep3Vector v1Xv2 ( cross(v) );
180 double abscross = v1Xv2.mag();
181 if ( abscross >= v1v2 ) {
182 return 1;
183 } else {
184 return abscross/v1v2;
185 }
186} /* howParallel() */
187
189 double epsilon) const {
190 // | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 | **2
191 // V1 is *this, V2 is v
192
193 static const double TOOBIG = std::pow(2.0,507);
194 static const double SCALE = std::pow(2.0,-507);
195 double v1v2 = std::fabs(dot(v));
196 if ( v1v2 == 0 ) {
197 return ( (mag2() == 0) && (v.mag2() == 0) );
198 }
199 if ( v1v2 >= TOOBIG ) {
200 Hep3Vector sv1 ( *this * SCALE );
201 Hep3Vector sv2 ( v * SCALE );
202 Hep3Vector sv1Xsv2 = sv1.cross(sv2);
203 double x2 = sv1Xsv2.mag2();
204 double limit = v1v2*SCALE*SCALE;
205 limit = epsilon*epsilon*limit*limit;
206 return ( x2 <= limit );
207 }
208
209 // At this point we know v1v2 can be squared.
210
211 Hep3Vector v1Xv2 ( cross(v) );
212 if ( (std::fabs (v1Xv2.dx) > TOOBIG) ||
213 (std::fabs (v1Xv2.dy) > TOOBIG) ||
214 (std::fabs (v1Xv2.dz) > TOOBIG) ) {
215 return false;
216 }
217
218 return ( (v1Xv2.mag2()) <= ((epsilon * v1v2) * (epsilon * v1v2)) );
219
220} /* isParallel() */
221
222
223double Hep3Vector::howOrthogonal (const Hep3Vector & v) const {
224 // | V1 dot V2 | / | V1 x V2 |
225
226 double v1v2 = std::fabs(dot(v));
227 //-| Safe because both v1 and v2 can be squared
228 if ( v1v2 == 0 ) {
229 return 0; // Even if one or both are 0, they are considered orthogonal
230 }
231 Hep3Vector v1Xv2 ( cross(v) );
232 double abscross = v1Xv2.mag();
233 if ( v1v2 >= abscross ) {
234 return 1;
235 } else {
236 return v1v2/abscross;
237 }
238
239} /* howOrthogonal() */
240
242 double epsilon) const {
243// | V1 x V2 | **2 <= epsilon **2 | V1 dot V2 | **2
244// V1 is *this, V2 is v
245
246 static const double TOOBIG = std::pow(2.0,507);
247 static const double SCALE = std::pow(2.0,-507);
248 double v1v2 = std::fabs(dot(v));
249 //-| Safe because both v1 and v2 can be squared
250 if ( v1v2 >= TOOBIG ) {
251 Hep3Vector sv1 ( *this * SCALE );
252 Hep3Vector sv2 ( v * SCALE );
253 Hep3Vector sv1Xsv2 = sv1.cross(sv2);
254 double x2 = sv1Xsv2.mag2();
255 double limit = epsilon*epsilon*x2;
256 double y2 = v1v2*SCALE*SCALE;
257 return ( y2*y2 <= limit );
258 }
259
260 // At this point we know v1v2 can be squared.
261
262 Hep3Vector eps_v1Xv2 ( cross(epsilon*v) );
263 if ( (std::fabs (eps_v1Xv2.dx) > TOOBIG) ||
264 (std::fabs (eps_v1Xv2.dy) > TOOBIG) ||
265 (std::fabs (eps_v1Xv2.dz) > TOOBIG) ) {
266 return true;
267 }
268
269 // At this point we know all the math we need can be done.
270
271 return ( v1v2*v1v2 <= eps_v1Xv2.mag2() );
272
273} /* isOrthogonal() */
274
275double Hep3Vector::setTolerance (double tol) {
276// Set the tolerance for Hep3Vectors to be considered near one another
277 double oldTolerance (tolerance);
278 tolerance = tol;
279 return oldTolerance;
280}
281
282
283//-***********************
284// Helper Methods:
285// negativeInfinity()
286//-***********************
287
289 // A byte-order-independent way to return -Infinity
290 struct Dib {
291 union {
292 double d;
293 unsigned char i[8];
294 } u;
295 };
296 Dib negOne;
297 Dib posTwo;
298 negOne.u.d = -1.0;
299 posTwo.u.d = 2.0;
300 Dib value;
301 int k;
302 for (k=0; k<8; k++) {
303 value.u.i[k] = negOne.u.i[k] | posTwo.u.i[k];
304 }
305 return value.u.d;
306}
307
308} // namespace CLHEP
309
310
#define ZMthrowC(A)
#define ZMthrowA(A)
void setRhoPhiEta(double rho, double phi, double eta)
bool operator>=(const Hep3Vector &v) const
double mag2() const
void setSpherical(double r, double theta, double phi)
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
bool operator>(const Hep3Vector &v) const
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
double negativeInfinity() const
static double setTolerance(double tol)
double mag() const
bool operator<=(const Hep3Vector &v) const
void setCylindrical(double r, double phi, double z)
int compare(const Hep3Vector &v) const
double howParallel(const Hep3Vector &v) const
double howOrthogonal(const Hep3Vector &v) const
void setRhoPhiTheta(double rho, double phi, double theta)
bool operator<(const Hep3Vector &v) const
bool isParallel(const Hep3Vector &v, double epsilon=tolerance) const