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