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