ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SpaceVector.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SpaceVector.cc
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 
23 
24 #include <cmath>
25 
26 namespace CLHEP {
27 
28 //-*****************************
29 // - 1 -
30 // set (multiple components)
31 // in various coordinate systems
32 //
33 //-*****************************
34 
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 
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 
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 
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 
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 
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 
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