SpaceVector.cc 8.12 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// -*- C++ -*-
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
//
// SpaceVector
//
// This is the implementation of those methods of the Hep3Vector class which
// originated from the ZOOM SpaceVector class.  Several groups of these methods
// have been separated off into the following code units:
//
// SpaceVectorR.cc	All methods involving rotation
// SpaceVectorD.cc	All methods involving angle decomposition
// SpaceVectorP.cc	Intrinsic properties and methods involving second vector
//

17
#include "CLHEP/Vector/defs.h"
18
19
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/ZMxpv.h"
Lynn Garren's avatar
Lynn Garren committed
20
#include "CLHEP/Units/PhysicalConstants.h"
21
22
23
24
25
26
27
28
29
30
31
32
33

#include <cmath>

namespace CLHEP  {

//-*****************************
//           - 1 -
// set (multiple components)
// in various coordinate systems
//
//-*****************************

void Hep3Vector::setSpherical (
34
35
36
37
		double r1,
                double theta1,
                double phi1) {
  if ( r1 < 0 ) {
38
39
40
41
    ZMthrowC (ZMxpvNegativeR(
      "Spherical coordinates set with negative   R"));
    // No special return needed if warning is ignored.
  }
42
  if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
43
44
45
46
    ZMthrowC (ZMxpvUnusualTheta(
      "Spherical coordinates set with theta not in [0, PI]"));
	// No special return needed if warning is ignored.
  }
47
  double rho1 ( r1*std::sin(theta1));
48
49
50
  setZ(r1 * std::cos(theta1));
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
51
  return;
52
} /* setSpherical (r, theta1, phi) */
53
54

void Hep3Vector::setCylindrical (
55
56
57
58
 		double rho1,
                double phi1,
                double z1) {
  if ( rho1 < 0 ) {
59
60
61
62
    ZMthrowC (ZMxpvNegativeR(
      "Cylindrical coordinates supplied with negative Rho"));
    // No special return needed if warning is ignored.
  }
63
64
65
  setZ(z1);
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
66
67
68
69
  return;
} /* setCylindrical (r, phi, z) */

void Hep3Vector::setRhoPhiTheta (
70
71
72
73
 		double rho1,
                double phi1,
                double theta1) {
  if (rho1 == 0) {
74
75
76
    ZMthrowC (ZMxpvZeroVector(
      "Attempt set vector components rho, phi, theta with zero rho -- "
      "zero vector is returned, ignoring theta and phi"));
77
    set(0.0, 0.0, 0.0);
78
79
    return;
  }
80
  if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
81
82
83
84
    ZMthrowA (ZMxpvInfiniteVector(
      "Attempt set cylindrical vector vector with finite rho and "
      "theta along the Z axis:  infinite Z would be computed"));
  }
85
  if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
86
87
88
89
    ZMthrowC (ZMxpvUnusualTheta(
      "Rho, phi, theta set with theta not in [0, PI]"));
	// No special return needed if warning is ignored.
  }
90
91
92
  setZ(rho1 / std::tan (theta1));
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
93
94
95
96
  return;
} /* setCyl (rho, phi, theta) */

void Hep3Vector::setRhoPhiEta (
97
98
99
100
 		double rho1,
                double phi1,
                double eta1 ) {
  if (rho1 == 0) {
101
102
103
    ZMthrowC (ZMxpvZeroVector(
      "Attempt set vector components rho, phi, eta with zero rho -- "
      "zero vector is returned, ignoring eta and phi"));
104
    set(0.0, 0.0, 0.0);
105
106
    return;
  }
107
  double theta1 (2 * std::atan ( std::exp (-eta1) ));
108
109
110
  setZ(rho1 / std::tan (theta1));
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
111
112
113
114
115
116
117
118
119
120
121
  return;
} /* setCyl (rho, phi, eta) */


//************
//    - 3 - 
// Comparisons
//
//************

int Hep3Vector::compare (const Hep3Vector & v) const {
122
  if       ( z() > v.z() ) {
123
    return 1;
124
  } else if ( z() < v.z() ) {
125
    return -1;
126
  } else if ( y() > v.y() ) {
127
    return 1;
128
  } else if ( y() < v.y() ) {
129
    return -1;
130
  } else if ( x() > v.x() ) {
131
    return 1;
132
  } else if ( x() < v.x() ) {
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
    return -1;
  } else {
    return 0;
  }
} /* Compare */


bool Hep3Vector::operator > (const Hep3Vector & v) const {
	return (compare(v)  > 0);
}
bool Hep3Vector::operator < (const Hep3Vector & v) const {
	return (compare(v)  < 0);
}
bool Hep3Vector::operator>= (const Hep3Vector & v) const {
	return (compare(v) >= 0);
}
bool Hep3Vector::operator<= (const Hep3Vector & v) const {
	return (compare(v) <= 0);
}



//-********
// Nearness
//-********

// These methods all assume you can safely take mag2() of each vector.
// Absolutely safe but slower and much uglier alternatives were 
// provided as build-time options in ZOOM SpaceVectors.
// Also, much smaller codes were provided tht assume you can square
// mag2() of each vector; but those return bad answers without warning
// when components exceed 10**90. 
//
// IsNear, HowNear, and DeltaR are found in ThreeVector.cc

double Hep3Vector::howParallel (const Hep3Vector & v) const {
  // | V1 x V2 | / | V1 dot V2 |
170
  double v1v2 = std::fabs(dot(v));
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
  if ( v1v2 == 0 ) {
    // Zero is parallel to no other vector except for zero.
    return ( (mag2() == 0) && (v.mag2() == 0) ) ? 0 : 1;
  }
  Hep3Vector v1Xv2 ( cross(v) );
  double abscross = v1Xv2.mag();
  if ( abscross >= v1v2 ) {
    return 1;
  } else {
    return abscross/v1v2;
  }
} /* howParallel() */

bool Hep3Vector::isParallel (const Hep3Vector & v,
                              double epsilon) const {
  // | V1 x V2 | **2  <= epsilon **2 | V1 dot V2 | **2
  // V1 is *this, V2 is v

189
190
191
  static const double TOOBIG = std::pow(2.0,507);
  static const double SCALE  = std::pow(2.0,-507);
  double v1v2 = std::fabs(dot(v));
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
  if ( v1v2 == 0 ) {
    return ( (mag2() == 0) && (v.mag2() == 0) );
  }
  if ( v1v2 >= TOOBIG ) {
    Hep3Vector sv1 ( *this * SCALE );
    Hep3Vector sv2 ( v * SCALE );
    Hep3Vector sv1Xsv2 = sv1.cross(sv2);
    double x2 = sv1Xsv2.mag2();
    double limit = v1v2*SCALE*SCALE;
    limit = epsilon*epsilon*limit*limit;
    return ( x2 <= limit );
  }

  // At this point we know v1v2 can be squared.

  Hep3Vector v1Xv2 ( cross(v) );
208
209
210
  if (  (std::fabs (v1Xv2.x()) > TOOBIG) ||
        (std::fabs (v1Xv2.y()) > TOOBIG) ||
        (std::fabs (v1Xv2.z()) > TOOBIG) ) {
211
212
213
214
215
216
217
218
219
220
221
    return false;
  }
                    
  return ( (v1Xv2.mag2()) <= ((epsilon * v1v2) * (epsilon * v1v2)) );

} /* isParallel() */


double Hep3Vector::howOrthogonal (const Hep3Vector & v) const {
  // | V1 dot V2 | / | V1 x V2 | 

222
  double v1v2 = std::fabs(dot(v));
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
	//-| Safe because both v1 and v2 can be squared
  if ( v1v2 == 0 ) {
    return 0;	// Even if one or both are 0, they are considered orthogonal
  }
  Hep3Vector v1Xv2 ( cross(v) );
  double abscross = v1Xv2.mag();
  if ( v1v2 >= abscross ) {
    return 1;
  } else {
    return v1v2/abscross;
  }

} /* howOrthogonal() */

bool Hep3Vector::isOrthogonal (const Hep3Vector & v,
			     double epsilon) const {
// | V1 x V2 | **2  <= epsilon **2 | V1 dot V2 | **2
// V1 is *this, V2 is v

242
243
244
  static const double TOOBIG = std::pow(2.0,507);
  static const double SCALE = std::pow(2.0,-507);
  double v1v2 = std::fabs(dot(v));
245
246
247
248
249
250
251
252
253
254
255
256
257
258
        //-| Safe because both v1 and v2 can be squared
  if ( v1v2 >= TOOBIG ) {
    Hep3Vector sv1 ( *this * SCALE );
    Hep3Vector sv2 ( v * SCALE );
    Hep3Vector sv1Xsv2 = sv1.cross(sv2);
    double x2 = sv1Xsv2.mag2();
    double limit = epsilon*epsilon*x2;
    double y2 = v1v2*SCALE*SCALE;
    return ( y2*y2 <= limit );
  }

  // At this point we know v1v2 can be squared.

  Hep3Vector eps_v1Xv2 ( cross(epsilon*v) );
259
260
261
  if (  (std::fabs (eps_v1Xv2.x()) > TOOBIG) ||
        (std::fabs (eps_v1Xv2.y()) > TOOBIG) ||
        (std::fabs (eps_v1Xv2.z()) > TOOBIG) ) {
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
    return true;
  }

  // At this point we know all the math we need can be done.

  return ( v1v2*v1v2 <= eps_v1Xv2.mag2() );

} /* isOrthogonal() */

double Hep3Vector::setTolerance (double tol) {
// Set the tolerance for Hep3Vectors to be considered near one another
  double oldTolerance (tolerance);
  tolerance = tol;
  return oldTolerance;
}


//-***********************
// Helper Methods:
//	negativeInfinity()
//-***********************

double Hep3Vector::negativeInfinity() const {
  // A byte-order-independent way to return -Infinity
  struct Dib {
    union {
      double d;
      unsigned char i[8];
    } u;
  };
  Dib negOne;
  Dib posTwo;
  negOne.u.d = -1.0;
  posTwo.u.d =  2.0;
  Dib value;
  int k;
  for (k=0; k<8; k++) {
    value.u.i[k] = negOne.u.i[k] | posTwo.u.i[k];
  }
  return value.u.d;
}

}  // namespace CLHEP