CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

ThreeVector.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: ThreeVector.cc,v 1.3 2003/08/13 20:00:14 garren Exp $
3 // ---------------------------------------------------------------------------
4 //
5 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
6 //
7 // This is the implementation of the Hep3Vector class.
8 //
9 // See also ThreeVectorR.cc for implementation of Hep3Vector methods which
10 // would couple in all the HepRotation methods.
11 //
12 
13 #ifdef GNUPRAGMA
14 #pragma implementation
15 #endif
16 
17 #include "CLHEP/Vector/defs.h"
18 #include "CLHEP/Vector/ThreeVector.h"
19 #include "CLHEP/Vector/ZMxpv.h"
20 #include "CLHEP/Units/PhysicalConstants.h"
21 
22 #include <cmath>
23 #include <iostream>
24 
25 namespace CLHEP {
26 
27 void Hep3Vector::setMag(double ma) {
28  double factor = mag();
29  if (factor == 0) {
30  ZMthrowA ( ZMxpvZeroVector (
31  "Hep3Vector::setMag : zero vector can't be stretched"));
32  }else{
33  factor = ma/factor;
34  setX(x()*factor);
35  setY(y()*factor);
36  setZ(z()*factor);
37  }
38 }
39 
40 double Hep3Vector::operator () (int i) const {
41  switch(i) {
42  case X:
43  return x();
44  case Y:
45  return y();
46  case Z:
47  return z();
48  default:
49  std::cerr << "Hep3Vector subscripting: bad index (" << i << ")"
50  << std::endl;
51  }
52  return 0.;
53 }
54 
55 double & Hep3Vector::operator () (int i) {
56  static double dummy;
57  switch(i) {
58  case X:
59  return dx;
60  case Y:
61  return dy;
62  case Z:
63  return dz;
64  default:
65  std::cerr
66  << "Hep3Vector subscripting: bad index (" << i << ")"
67  << std::endl;
68  return dummy;
69  }
70 }
71 
73  // NewUzVector must be normalized !
74 
75  double u1 = NewUzVector.x();
76  double u2 = NewUzVector.y();
77  double u3 = NewUzVector.z();
78  double up = u1*u1 + u2*u2;
79 
80  if (up>0) {
81  up = std::sqrt(up);
82  double px = dx, py = dy, pz = dz;
83  dx = (u1*u3*px - u2*py)/up + u1*pz;
84  dy = (u2*u3*px + u1*py)/up + u2*pz;
85  dz = -up*px + u3*pz;
86  }
87  else if (u3 < 0.) { dx = -dx; dz = -dz; } // phi=0 teta=pi
88  else {};
89  return *this;
90 }
91 
93  double m1 = mag();
94  if ( m1== 0 ) return 0.0;
95  if ( m1== z() ) return 1.0E72;
96  if ( m1== -z() ) return -1.0E72;
97  return 0.5*std::log( (m1+z())/(m1-z()) );
98 }
99 
100 std::ostream & operator<< (std::ostream & os, const Hep3Vector & v) {
101  return os << "(" << v.x() << "," << v.y() << "," << v.z() << ")";
102 }
103 
104 void ZMinput3doubles ( std::istream & is, const char * type,
105  double & x, double & y, double & z );
106 
107 std::istream & operator>>(std::istream & is, Hep3Vector & v) {
108  double x, y, z;
109  ZMinput3doubles ( is, "Hep3Vector", x, y, z );
110  v.set(x, y, z);
111  return is;
112 } // operator>>()
113 
114 const Hep3Vector HepXHat(1.0, 0.0, 0.0);
115 const Hep3Vector HepYHat(0.0, 1.0, 0.0);
116 const Hep3Vector HepZHat(0.0, 0.0, 1.0);
117 
118 //-------------------
119 //
120 // New methods introduced when ZOOM PhysicsVectors was merged in:
121 //
122 //-------------------
123 
125  double sinphi = std::sin(phi1);
126  double cosphi = std::cos(phi1);
127  double ty;
128  ty = dy * cosphi - dz * sinphi;
129  dz = dz * cosphi + dy * sinphi;
130  dy = ty;
131  return *this;
132 } /* rotateX */
133 
135  double sinphi = std::sin(phi1);
136  double cosphi = std::cos(phi1);
137  double tz;
138  tz = dz * cosphi - dx * sinphi;
139  dx = dx * cosphi + dz * sinphi;
140  dz = tz;
141  return *this;
142 } /* rotateY */
143 
145  double sinphi = std::sin(phi1);
146  double cosphi = std::cos(phi1);
147  double tx;
148  tx = dx * cosphi - dy * sinphi;
149  dy = dy * cosphi + dx * sinphi;
150  dx = tx;
151  return *this;
152 } /* rotateZ */
153 
154 bool Hep3Vector::isNear(const Hep3Vector & v, double epsilon) const {
155  double limit = dot(v)*epsilon*epsilon;
156  return ( (*this - v).mag2() <= limit );
157 } /* isNear() */
158 
159 double Hep3Vector::howNear(const Hep3Vector & v ) const {
160  // | V1 - V2 | **2 / V1 dot V2, up to 1
161  double d = (*this - v).mag2();
162  double vdv = dot(v);
163  if ( (vdv > 0) && (d < vdv) ) {
164  return std::sqrt (d/vdv);
165  } else if ( (vdv == 0) && (d == 0) ) {
166  return 0;
167  } else {
168  return 1;
169  }
170 } /* howNear */
171 
172 double Hep3Vector::deltaPhi (const Hep3Vector & v2) const {
173  double dphi = v2.getPhi() - getPhi();
174  if ( dphi > CLHEP::pi ) {
175  dphi -= CLHEP::twopi;
176  } else if ( dphi <= -CLHEP::pi ) {
177  dphi += CLHEP::twopi;
178  }
179  return dphi;
180 } /* deltaPhi */
181 
182 double Hep3Vector::deltaR ( const Hep3Vector & v ) const {
183  double a = eta() - v.eta();
184  double b = deltaPhi(v);
185  return std::sqrt ( a*a + b*b );
186 } /* deltaR */
187 
188 double Hep3Vector::cosTheta(const Hep3Vector & q) const {
189  double arg;
190  double ptot2 = mag2()*q.mag2();
191  if(ptot2 <= 0) {
192  arg = 0.0;
193  }else{
194  arg = dot(q)/std::sqrt(ptot2);
195  if(arg > 1.0) arg = 1.0;
196  if(arg < -1.0) arg = -1.0;
197  }
198  return arg;
199 }
200 
201 double Hep3Vector::cos2Theta(const Hep3Vector & q) const {
202  double arg;
203  double ptot2 = mag2();
204  double qtot2 = q.mag2();
205  if ( ptot2 == 0 || qtot2 == 0 ) {
206  arg = 1.0;
207  }else{
208  double pdq = dot(q);
209  arg = (pdq/ptot2) * (pdq/qtot2);
210  // More naive methods overflow on vectors which can be squared
211  // but can't be raised to the 4th power.
212  if(arg > 1.0) arg = 1.0;
213  }
214  return arg;
215 }
216 
217 void Hep3Vector::setEta (double eta1) {
218  double phi1 = 0;
219  double r1;
220  if ( (dx == 0) && (dy == 0) ) {
221  if (dz == 0) {
222  ZMthrowC (ZMxpvZeroVector(
223  "Attempt to set eta of zero vector -- vector is unchanged"));
224  return;
225  }
226  ZMthrowC (ZMxpvZeroVector(
227  "Attempt to set eta of vector along Z axis -- will use phi = 0"));
228  r1 = std::fabs(dz);
229  } else {
230  r1 = getR();
231  phi1 = getPhi();
232  }
233  double tanHalfTheta = std::exp ( -eta1 );
234  double cosTheta1 =
235  (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
236  dz = r1 * cosTheta1;
237  double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
238  dy = rho1 * std::sin (phi1);
239  dx = rho1 * std::cos (phi1);
240  return;
241 }
242 
243 void Hep3Vector::setCylTheta (double theta1) {
244 
245  // In cylindrical coords, set theta while keeping rho and phi fixed
246 
247  if ( (dx == 0) && (dy == 0) ) {
248  if (dz == 0) {
249  ZMthrowC (ZMxpvZeroVector(
250  "Attempt to set cylTheta of zero vector -- vector is unchanged"));
251  return;
252  }
253  if (theta1 == 0) {
254  dz = std::fabs(dz);
255  return;
256  }
257  if (theta1 == CLHEP::pi) {
258  dz = -std::fabs(dz);
259  return;
260  }
261  ZMthrowC (ZMxpvZeroVector(
262  "Attempt set cylindrical theta of vector along Z axis "
263  "to a non-trivial value, while keeping rho fixed -- "
264  "will return zero vector"));
265  dz = 0;
266  return;
267  }
268  if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
269  ZMthrowC (ZMxpvUnusualTheta(
270  "Setting Cyl theta of a vector based on a value not in [0, PI]"));
271  // No special return needed if warning is ignored.
272  }
273  double phi1 (getPhi());
274  double rho1 = getRho();
275  if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
276  ZMthrowC (ZMxpvInfiniteVector(
277  "Attempt to set cylindrical theta to 0 or PI "
278  "while keeping rho fixed -- infinite Z will be computed"));
279  dz = (theta1==0) ? 1.0E72 : -1.0E72;
280  return;
281  }
282  dz = rho1 / std::tan (theta1);
283  dy = rho1 * std::sin (phi1);
284  dx = rho1 * std::cos (phi1);
285 
286 } /* setCylTheta */
287 
288 void Hep3Vector::setCylEta (double eta1) {
289 
290  // In cylindrical coords, set eta while keeping rho and phi fixed
291 
292  double theta1 = 2 * std::atan ( std::exp (-eta1) );
293 
294  //-| The remaining code is similar to setCylTheta, The reason for
295  //-| using a copy is so as to be able to change the messages in the
296  //-| ZMthrows to say eta rather than theta. Besides, we assumedly
297  //-| need not check for theta of 0 or PI.
298 
299  if ( (dx == 0) && (dy == 0) ) {
300  if (dz == 0) {
301  ZMthrowC (ZMxpvZeroVector(
302  "Attempt to set cylEta of zero vector -- vector is unchanged"));
303  return;
304  }
305  if (theta1 == 0) {
306  dz = std::fabs(dz);
307  return;
308  }
309  if (theta1 == CLHEP::pi) {
310  dz = -std::fabs(dz);
311  return;
312  }
313  ZMthrowC (ZMxpvZeroVector(
314  "Attempt set cylindrical eta of vector along Z axis "
315  "to a non-trivial value, while keeping rho fixed -- "
316  "will return zero vector"));
317  dz = 0;
318  return;
319  }
320  double phi1 (getPhi());
321  double rho1 = getRho();
322  dz = rho1 / std::tan (theta1);
323  dy = rho1 * std::sin (phi1);
324  dx = rho1 * std::cos (phi1);
325 
326 } /* setCylEta */
327 
328 
329 Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
330  if (c == 0) {
331  ZMthrowA ( ZMxpvInfiniteVector (
332  "Attempt to divide vector by 0 -- "
333  "will produce infinities and/or NANs"));
334  }
335  double oneOverC = 1.0/c;
336  return Hep3Vector ( v1.x() * oneOverC,
337  v1.y() * oneOverC,
338  v1.z() * oneOverC );
339 } /* v / c */
340 
342  if (c == 0) {
343  ZMthrowA (ZMxpvInfiniteVector(
344  "Attempt to do vector /= 0 -- "
345  "division by zero would produce infinite or NAN components"));
346  }
347  double oneOverC = 1.0/c;
348  dx *= oneOverC;
349  dy *= oneOverC;
350  dz *= oneOverC;
351  return *this;
352 }
353 
355 
356 } // namespace CLHEP
357 
358 
359 
CLHEP::Hep3Vector::cos2Theta
double cos2Theta() const
CLHEP::Hep3Vector::setEta
void setEta(double p)
Definition: ThreeVector.cc:217
CLHEP::Hep3Vector::deltaR
double deltaR(const Hep3Vector &v) const
Definition: ThreeVector.cc:182
CLHEP::Hep3Vector::setX
void setX(double)
CLHEP::Hep3Vector::X
@ X
Definition: Geometry/CLHEP/Vector/ThreeVector.h:47
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::getPhi
double getPhi() const
a
@ a
Definition: testCategories.cc:125
CLHEP::Hep3Vector::operator/=
Hep3Vector & operator/=(double)
Definition: ThreeVector.cc:341
CLHEP::HepYHat
const Hep3Vector HepYHat
Definition: Geometry/CLHEP/Vector/ThreeVector.h:425
CLHEP::Hep3Vector::eta
double eta() const
b
@ b
Definition: testCategories.cc:125
CLHEP::Hep3Vector::tolerance
static double tolerance
Definition: Geometry/CLHEP/Vector/ThreeVector.h:400
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
CLHEP::Hep3Vector::cosTheta
double cosTheta() const
CLHEP::Hep3Vector::isNear
bool isNear(const Hep3Vector &, double epsilon=tolerance) const
Definition: ThreeVector.cc:154
CLHEP::Hep3Vector::rotateY
Hep3Vector & rotateY(double)
Definition: ThreeVector.cc:134
CLHEP::Hep3Vector::setCylTheta
void setCylTheta(double)
Definition: ThreeVector.cc:243
CLHEP::Hep3Vector::operator()
double operator()(int) const
Definition: ThreeVector.cc:40
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::Hep3Vector::getRho
double getRho() const
CLHEP::Hep3Vector::rotateX
Hep3Vector & rotateX(double)
Definition: ThreeVector.cc:124
CLHEP::Hep3Vector::dy
double dy
Definition: Geometry/CLHEP/Vector/ThreeVector.h:396
CLHEP::Hep3Vector::x
double x() const
CLHEP::operator/
HepLorentzVector operator/(const HepLorentzVector &, double a)
Definition: LorentzVector.cc:162
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::dz
double dz
Definition: Geometry/CLHEP/Vector/ThreeVector.h:397
CLHEP
Definition: ClhepVersion.h:13
type
Signatures of Hep3Vector::rotate For equivalent ZOOM axis There is no harm in leaving this axis CLHEP has implemented a first forming an identity then rotating that by axis and I leave the CLHEP code alone people are of course free to use the ZOOM originated method with signature which I believe will be faster Return types for rotateZ CLHEP and PhysicsVectors each have these three and they are identical except that the ZOOM version returns a reference to while in CLHEP they return void Having methods that alter an object return a reference to that object is convenient for certain chained and costs nothing I don t wish to potentially break ZOOM user code for no good so I have made these CLHEP method conform to this convention There are a couple of other CLHEP rotate and which use the void return type
Definition: minorMergeIssues.doc:113
CLHEP::Hep3Vector::deltaPhi
double deltaPhi(const Hep3Vector &v2) const
Definition: ThreeVector.cc:172
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::ToleranceTicks
@ ToleranceTicks
Definition: Geometry/CLHEP/Vector/ThreeVector.h:299
CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
CLHEP::Hep3Vector::pseudoRapidity
double pseudoRapidity() const
Definition: ThreeVector.cc:92
Hep3Vector
Issues Concerning the PhysicsVectors CLHEP Vector Merge The merge of ZOOM PhysicsVdectors and the CLHEP Vector package is completed The purpose of this document is to list the major issues that affected the merge of these and where relevant describe the resolutions More detailed documents describe more minor issues General Approach As agreed at the June CLHEP the approach is to combine the features of each ZOOM class with the corresponding CLHEP class expanding the interface to create a single lingua franca of what a Hep3Vector(for example) means. We are not forming SpaceVector as an class derived from Hep3Vector and enhancing it in that way. Another rule imposed by the agreement is to avoid using the Exceptions package(even though that will later go into CLHEP for other uses). A desirable goal is to avoid cluttering the interface and enlarging the code linked in when ordinary CLHEP Vector functionallity is used. To this end
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::howNear
double howNear(const Hep3Vector &v) const
Definition: ThreeVector.cc:159
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::HepXHat
const Hep3Vector HepXHat
Definition: Matrix/CLHEP/Vector/ThreeVector.h:425
CLHEP::Hep3Vector::setCylEta
void setCylEta(double p)
Definition: ThreeVector.cc:288
CLHEP::Hep3Vector
Definition: Geometry/CLHEP/Vector/ThreeVector.h:41
CLHEP::Hep3Vector::Z
@ Z
Definition: Geometry/CLHEP/Vector/ThreeVector.h:47
ZMthrowC
#define ZMthrowC(A)
Definition: Geometry/CLHEP/Vector/ZMxpv.h:132
CLHEP::operator>>
std::istream & operator>>(std::istream &is, HepAxisAngle &aa)
Definition: AxisAngle.cc:96
CLHEP::Hep3Vector::dx
double dx
Definition: Geometry/CLHEP/Vector/ThreeVector.h:395
CLHEP::HepZHat
const Hep3Vector HepZHat
Definition: Geometry/CLHEP/Vector/ThreeVector.h:425
CLHEP::Hep3Vector::rotateZ
Hep3Vector & rotateZ(double)
Definition: ThreeVector.cc:144
CLHEP::Hep3Vector::Y
@ Y
Definition: Geometry/CLHEP/Vector/ThreeVector.h:47
i
long i
Definition: JamesRandomSeeding.txt:27
CLHEP::Hep3Vector::y
double y() const
CLHEP::ZMinput3doubles
void ZMinput3doubles(std::istream &is, const char *type, double &x, double &y, double &z)
Definition: ZMinput.cc:37
CLHEP::Hep3Vector::setMag
void setMag(double)
Definition: ThreeVector.cc:27
CLHEP::Hep3Vector::getR
double getR() const
CLHEP::operator<<
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt: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