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 
CLHEP::Hep3Vector::howOrthogonal
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:223
CLHEP::Hep3Vector::compare
int compare(const Hep3Vector &v) const
Definition: SpaceVector.cc:125
CLHEP::Hep3Vector::operator<=
bool operator<=(const Hep3Vector &v) const
Definition: SpaceVector.cc:153
CLHEP::Hep3Vector::howParallel
double howParallel(const Hep3Vector &v) const
Definition: SpaceVector.cc:172
CLHEP::Hep3Vector::tolerance
static double tolerance
Definition: Geometry/CLHEP/Vector/ThreeVector.h:400
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::Hep3Vector::dy
double dy
Definition: Geometry/CLHEP/Vector/ThreeVector.h:396
CLHEP::Hep3Vector::setRhoPhiTheta
void setRhoPhiTheta(double rho, double phi, double theta)
Definition: SpaceVector.cc:73
CLHEP::Hep3Vector::isParallel
bool isParallel(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:188
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::Hep3Vector::dz
double dz
Definition: Geometry/CLHEP/Vector/ThreeVector.h:397
CLHEP
Definition: ClhepVersion.h:13
CLHEP::Hep3Vector::mag
double mag() const
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::Hep3Vector::operator>=
bool operator>=(const Hep3Vector &v) const
Definition: SpaceVector.cc:150
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector
Definition: Geometry/CLHEP/Vector/ThreeVector.h:41
CLHEP::Hep3Vector::operator<
bool operator<(const Hep3Vector &v) const
Definition: SpaceVector.cc:147
ZMthrowC
#define ZMthrowC(A)
Definition: Geometry/CLHEP/Vector/ZMxpv.h:132
CLHEP::Hep3Vector::dx
double dx
Definition: Geometry/CLHEP/Vector/ThreeVector.h:395
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::Hep3Vector::setCylindrical
void setCylindrical(double r, double phi, double z)
Definition: SpaceVector.cc:58
CLHEP::Hep3Vector::operator>
bool operator>(const Hep3Vector &v) const
Definition: SpaceVector.cc:144
CLHEP::Hep3Vector::setSpherical
void setSpherical(double r, double theta, double phi)
Definition: SpaceVector.cc:37
ZMthrowA
it has advantages For I leave the ZMthrows but substitute I replaced ZMthrow with ZMthrowA in this package ZMthrowA
Definition: keyMergeIssues.doc:69
CLHEP::Hep3Vector::negativeInfinity
double negativeInfinity() const
Definition: SpaceVector.cc:288
k
long k
Definition: JamesRandomSeeding.txt:29
CLHEP::Hep3Vector::setRhoPhiEta
void setRhoPhiEta(double rho, double phi, double eta)
Definition: SpaceVector.cc:100
CLHEP::Hep3Vector::setTolerance
static double setTolerance(double tol)
Definition: SpaceVector.cc:275
CLHEP::Hep3Vector::isOrthogonal
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:241