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