Geant4  10.00.p01
ThreeVector.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:$
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/ThreeVector.h"
18 #include "CLHEP/Units/PhysicalConstants.h"
19 
20 #include <cmath>
21 #include <iostream>
22 
23 namespace CLHEP {
24 
25 void Hep3Vector::setMag(double ma) {
26  double factor = mag();
27  if (factor == 0) {
28  std::cerr << "Hep3Vector::setMag() - "
29  << "zero vector can't be stretched" << std::endl;
30  }else{
31  factor = ma/factor;
32  setX(x()*factor);
33  setY(y()*factor);
34  setZ(z()*factor);
35  }
36 }
37 
38 double Hep3Vector::operator () (int i) const {
39  switch(i) {
40  case X:
41  return x();
42  case Y:
43  return y();
44  case Z:
45  return z();
46  default:
47  std::cerr << "Hep3Vector::operator () - "
48  << "Hep3Vector subscripting: bad index (" << i << ")"
49  << std::endl;
50  }
51  return 0.;
52 }
53 
54 double & Hep3Vector::operator () (int i) {
55  static double dummy;
56  switch(i) {
57  case X:
58  return dx;
59  case Y:
60  return dy;
61  case Z:
62  return dz;
63  default:
64  std::cerr
65  << "Hep3Vector::operator () - "
66  << "Hep3Vector subscripting: bad index (" << i << ")"
67  << std::endl;
68  return dummy;
69  }
70 }
71 
72 Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) {
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 
92 double Hep3Vector::pseudoRapidity() const {
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 
124 Hep3Vector & Hep3Vector::rotateX (double phi1) {
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 
134 Hep3Vector & Hep3Vector::rotateY (double phi1) {
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 
144 Hep3Vector & Hep3Vector::rotateZ (double phi1) {
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  std::cerr << "Hep3Vector::setEta() - "
223  << "Attempt to set eta of zero vector -- vector is unchanged"
224  << std::endl;
225  return;
226  }
227  std::cerr << "Hep3Vector::setEta() - "
228  << "Attempt to set eta of vector along Z axis -- will use phi = 0"
229  << std::endl;
230  r1 = std::fabs(dz);
231  } else {
232  r1 = getR();
233  phi1 = getPhi();
234  }
235  double tanHalfTheta = std::exp ( -eta1 );
236  double cosTheta1 =
237  (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
238  dz = r1 * cosTheta1;
239  double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
240  dy = rho1 * std::sin (phi1);
241  dx = rho1 * std::cos (phi1);
242  return;
243 }
244 
245 void Hep3Vector::setCylTheta (double theta1) {
246 
247  // In cylindrical coords, set theta while keeping rho and phi fixed
248 
249  if ( (dx == 0) && (dy == 0) ) {
250  if (dz == 0) {
251  std::cerr << "Hep3Vector::setCylTheta() - "
252  << "Attempt to set cylTheta of zero vector -- vector is unchanged"
253  << std::endl;
254  return;
255  }
256  if (theta1 == 0) {
257  dz = std::fabs(dz);
258  return;
259  }
260  if (theta1 == CLHEP::pi) {
261  dz = -std::fabs(dz);
262  return;
263  }
264  std::cerr << "Hep3Vector::setCylTheta() - "
265  << "Attempt set cylindrical theta of vector along Z axis "
266  << "to a non-trivial value, while keeping rho fixed -- "
267  << "will return zero vector" << std::endl;
268  dz = 0;
269  return;
270  }
271  if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
272  std::cerr << "Hep3Vector::setCylTheta() - "
273  << "Setting Cyl theta of a vector based on a value not in [0, PI]"
274  << std::endl;
275  // No special return needed if warning is ignored.
276  }
277  double phi1 (getPhi());
278  double rho1 = getRho();
279  if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
280  std::cerr << "Hep3Vector::setCylTheta() - "
281  << "Attempt to set cylindrical theta to 0 or PI "
282  << "while keeping rho fixed -- infinite Z will be computed"
283  << std::endl;
284  dz = (theta1==0) ? 1.0E72 : -1.0E72;
285  return;
286  }
287  dz = rho1 / std::tan (theta1);
288  dy = rho1 * std::sin (phi1);
289  dx = rho1 * std::cos (phi1);
290 
291 } /* setCylTheta */
292 
293 void Hep3Vector::setCylEta (double eta1) {
294 
295  // In cylindrical coords, set eta while keeping rho and phi fixed
296 
297  double theta1 = 2 * std::atan ( std::exp (-eta1) );
298 
299  //-| The remaining code is similar to setCylTheta, The reason for
300  //-| using a copy is so as to be able to change the messages in the
301  //-| ZMthrows to say eta rather than theta. Besides, we assumedly
302  //-| need not check for theta of 0 or PI.
303 
304  if ( (dx == 0) && (dy == 0) ) {
305  if (dz == 0) {
306  std::cerr << "Hep3Vector::setCylEta() - "
307  << "Attempt to set cylEta of zero vector -- vector is unchanged"
308  << std::endl;
309  return;
310  }
311  if (theta1 == 0) {
312  dz = std::fabs(dz);
313  return;
314  }
315  if (theta1 == CLHEP::pi) {
316  dz = -std::fabs(dz);
317  return;
318  }
319  std::cerr << "Hep3Vector::setCylEta() - "
320  << "Attempt set cylindrical eta of vector along Z axis "
321  << "to a non-trivial value, while keeping rho fixed -- "
322  << "will return zero vector" << std::endl;
323  dz = 0;
324  return;
325  }
326  double phi1 (getPhi());
327  double rho1 = getRho();
328  dz = rho1 / std::tan (theta1);
329  dy = rho1 * std::sin (phi1);
330  dx = rho1 * std::cos (phi1);
331 
332 } /* setCylEta */
333 
334 
335 Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
336 // if (c == 0) {
337 // std::cerr << "Hep3Vector::operator/ () - "
338 // << "Attempt to divide vector by 0 -- "
339 // << "will produce infinities and/or NANs" << std::endl;
340 // }
341  double oneOverC = 1.0/c;
342  return Hep3Vector ( v1.x() * oneOverC,
343  v1.y() * oneOverC,
344  v1.z() * oneOverC );
345 } /* v / c */
346 
347 Hep3Vector & Hep3Vector::operator/= (double c) {
348 // if (c == 0) {
349 // std::cerr << "Hep3Vector::operator/ () - "
350 // << "Attempt to do vector /= 0 -- "
351 // << "division by zero would produce infinite or NAN components"
352 // << std::endl;
353 // }
354  double oneOverC = 1.0/c;
355  dx *= oneOverC;
356  dy *= oneOverC;
357  dz *= oneOverC;
358  return *this;
359 }
360 
361 double Hep3Vector::tolerance = Hep3Vector::ToleranceTicks * 2.22045e-16;
362 
363 } // namespace CLHEP
HepLorentzVector operator/(const HepLorentzVector &w, double c)
const Hep3Vector HepYHat(0.0, 1.0, 0.0)
G4double z
Definition: TRTMaterials.hh:39
const G4double pi
G4double a
Definition: TRTMaterials.hh:39
const Hep3Vector HepXHat(1.0, 0.0, 0.0)
std::istream & operator>>(std::istream &is, HepAxisAngle &aa)
Definition: AxisAngle.cc:95
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:85
const Hep3Vector HepZHat(0.0, 0.0, 1.0)
void ZMinput3doubles(std::istream &is, const char *type, double &x, double &y, double &z)
Definition: ZMinput.cc:42