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 
28 namespace 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 
125 int 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 
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 bool 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 
172 double 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 
223 double 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 
275 double 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 
bool operator>(const Hep3Vector &v) const
Definition: SpaceVector.cc:144
bool operator>=(const Hep3Vector &v) const
Definition: SpaceVector.cc:150
double negativeInfinity() const
Definition: SpaceVector.cc:288
void setCylindrical(double r, double phi, double z)
Definition: SpaceVector.cc:58
int compare(const Hep3Vector &v) const
Definition: SpaceVector.cc:125
double dot(const Hep3Vector &) const
bool operator<(const Hep3Vector &v) const
Definition: SpaceVector.cc:147
#define ZMthrowC(A)
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:241
bool operator<=(const Hep3Vector &v) const
Definition: SpaceVector.cc:153
void setRhoPhiEta(double rho, double phi, double eta)
Definition: SpaceVector.cc:100
double mag2() const
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:223
#define ZMthrowA(A)
void setRhoPhiTheta(double rho, double phi, double theta)
Definition: SpaceVector.cc:73
void setSpherical(double r, double theta, double phi)
Definition: SpaceVector.cc:37
bool isParallel(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:188
Hep3Vector cross(const Hep3Vector &) const
static double setTolerance(double tol)
Definition: SpaceVector.cc:275
double mag() const
double howParallel(const Hep3Vector &v) const
Definition: SpaceVector.cc:172