Geant4  10.01.p03
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 Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) {
39  // NewUzVector must be normalized !
40 
41  double u1 = NewUzVector.x();
42  double u2 = NewUzVector.y();
43  double u3 = NewUzVector.z();
44  double up = u1*u1 + u2*u2;
45 
46  if (up>0) {
47  up = std::sqrt(up);
48  double px = dx, py = dy, pz = dz;
49  dx = (u1*u3*px - u2*py)/up + u1*pz;
50  dy = (u2*u3*px + u1*py)/up + u2*pz;
51  dz = -up*px + u3*pz;
52  }
53  else if (u3 < 0.) { dx = -dx; dz = -dz; } // phi=0 teta=pi
54  else {};
55  return *this;
56 }
57 
58 double Hep3Vector::pseudoRapidity() const {
59  double m1 = mag();
60  if ( m1== 0 ) return 0.0;
61  if ( m1== z() ) return 1.0E72;
62  if ( m1== -z() ) return -1.0E72;
63  return 0.5*std::log( (m1+z())/(m1-z()) );
64 }
65 
66 std::ostream & operator<< (std::ostream & os, const Hep3Vector & v) {
67  return os << "(" << v.x() << "," << v.y() << "," << v.z() << ")";
68 }
69 
70 void ZMinput3doubles ( std::istream & is, const char * type,
71  double & x, double & y, double & z );
72 
73 std::istream & operator>>(std::istream & is, Hep3Vector & v) {
74  double x, y, z;
75  ZMinput3doubles ( is, "Hep3Vector", x, y, z );
76  v.set(x, y, z);
77  return is;
78 } // operator>>()
79 
80 const Hep3Vector HepXHat(1.0, 0.0, 0.0);
81 const Hep3Vector HepYHat(0.0, 1.0, 0.0);
82 const Hep3Vector HepZHat(0.0, 0.0, 1.0);
83 
84 //-------------------
85 //
86 // New methods introduced when ZOOM PhysicsVectors was merged in:
87 //
88 //-------------------
89 
90 Hep3Vector & Hep3Vector::rotateX (double phi1) {
91  double sinphi = std::sin(phi1);
92  double cosphi = std::cos(phi1);
93  double ty;
94  ty = dy * cosphi - dz * sinphi;
95  dz = dz * cosphi + dy * sinphi;
96  dy = ty;
97  return *this;
98 } /* rotateX */
99 
100 Hep3Vector & Hep3Vector::rotateY (double phi1) {
101  double sinphi = std::sin(phi1);
102  double cosphi = std::cos(phi1);
103  double tz;
104  tz = dz * cosphi - dx * sinphi;
105  dx = dx * cosphi + dz * sinphi;
106  dz = tz;
107  return *this;
108 } /* rotateY */
109 
110 Hep3Vector & Hep3Vector::rotateZ (double phi1) {
111  double sinphi = std::sin(phi1);
112  double cosphi = std::cos(phi1);
113  double tx;
114  tx = dx * cosphi - dy * sinphi;
115  dy = dy * cosphi + dx * sinphi;
116  dx = tx;
117  return *this;
118 } /* rotateZ */
119 
120 bool Hep3Vector::isNear(const Hep3Vector & v, double epsilon) const {
121  double limit = dot(v)*epsilon*epsilon;
122  return ( (*this - v).mag2() <= limit );
123 } /* isNear() */
124 
125 double Hep3Vector::howNear(const Hep3Vector & v ) const {
126  // | V1 - V2 | **2 / V1 dot V2, up to 1
127  double d = (*this - v).mag2();
128  double vdv = dot(v);
129  if ( (vdv > 0) && (d < vdv) ) {
130  return std::sqrt (d/vdv);
131  } else if ( (vdv == 0) && (d == 0) ) {
132  return 0;
133  } else {
134  return 1;
135  }
136 } /* howNear */
137 
138 double Hep3Vector::deltaPhi (const Hep3Vector & v2) const {
139  double dphi = v2.getPhi() - getPhi();
140  if ( dphi > CLHEP::pi ) {
141  dphi -= CLHEP::twopi;
142  } else if ( dphi <= -CLHEP::pi ) {
143  dphi += CLHEP::twopi;
144  }
145  return dphi;
146 } /* deltaPhi */
147 
148 double Hep3Vector::deltaR ( const Hep3Vector & v ) const {
149  double a = eta() - v.eta();
150  double b = deltaPhi(v);
151  return std::sqrt ( a*a + b*b );
152 } /* deltaR */
153 
154 double Hep3Vector::cosTheta(const Hep3Vector & q) const {
155  double arg;
156  double ptot2 = mag2()*q.mag2();
157  if(ptot2 <= 0) {
158  arg = 0.0;
159  }else{
160  arg = dot(q)/std::sqrt(ptot2);
161  if(arg > 1.0) arg = 1.0;
162  if(arg < -1.0) arg = -1.0;
163  }
164  return arg;
165 }
166 
167 double Hep3Vector::cos2Theta(const Hep3Vector & q) const {
168  double arg;
169  double ptot2 = mag2();
170  double qtot2 = q.mag2();
171  if ( ptot2 == 0 || qtot2 == 0 ) {
172  arg = 1.0;
173  }else{
174  double pdq = dot(q);
175  arg = (pdq/ptot2) * (pdq/qtot2);
176  // More naive methods overflow on vectors which can be squared
177  // but can't be raised to the 4th power.
178  if(arg > 1.0) arg = 1.0;
179  }
180  return arg;
181 }
182 
183 void Hep3Vector::setEta (double eta1) {
184  double phi1 = 0;
185  double r1;
186  if ( (dx == 0) && (dy == 0) ) {
187  if (dz == 0) {
188  std::cerr << "Hep3Vector::setEta() - "
189  << "Attempt to set eta of zero vector -- vector is unchanged"
190  << std::endl;
191  return;
192  }
193  std::cerr << "Hep3Vector::setEta() - "
194  << "Attempt to set eta of vector along Z axis -- will use phi = 0"
195  << std::endl;
196  r1 = std::fabs(dz);
197  } else {
198  r1 = getR();
199  phi1 = getPhi();
200  }
201  double tanHalfTheta = std::exp ( -eta1 );
202  double cosTheta1 =
203  (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
204  dz = r1 * cosTheta1;
205  double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
206  dy = rho1 * std::sin (phi1);
207  dx = rho1 * std::cos (phi1);
208  return;
209 }
210 
211 void Hep3Vector::setCylTheta (double theta1) {
212 
213  // In cylindrical coords, set theta while keeping rho and phi fixed
214 
215  if ( (dx == 0) && (dy == 0) ) {
216  if (dz == 0) {
217  std::cerr << "Hep3Vector::setCylTheta() - "
218  << "Attempt to set cylTheta of zero vector -- vector is unchanged"
219  << std::endl;
220  return;
221  }
222  if (theta1 == 0) {
223  dz = std::fabs(dz);
224  return;
225  }
226  if (theta1 == CLHEP::pi) {
227  dz = -std::fabs(dz);
228  return;
229  }
230  std::cerr << "Hep3Vector::setCylTheta() - "
231  << "Attempt set cylindrical theta of vector along Z axis "
232  << "to a non-trivial value, while keeping rho fixed -- "
233  << "will return zero vector" << std::endl;
234  dz = 0;
235  return;
236  }
237  if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
238  std::cerr << "Hep3Vector::setCylTheta() - "
239  << "Setting Cyl theta of a vector based on a value not in [0, PI]"
240  << std::endl;
241  // No special return needed if warning is ignored.
242  }
243  double phi1 (getPhi());
244  double rho1 = getRho();
245  if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
246  std::cerr << "Hep3Vector::setCylTheta() - "
247  << "Attempt to set cylindrical theta to 0 or PI "
248  << "while keeping rho fixed -- infinite Z will be computed"
249  << std::endl;
250  dz = (theta1==0) ? 1.0E72 : -1.0E72;
251  return;
252  }
253  dz = rho1 / std::tan (theta1);
254  dy = rho1 * std::sin (phi1);
255  dx = rho1 * std::cos (phi1);
256 
257 } /* setCylTheta */
258 
259 void Hep3Vector::setCylEta (double eta1) {
260 
261  // In cylindrical coords, set eta while keeping rho and phi fixed
262 
263  double theta1 = 2 * std::atan ( std::exp (-eta1) );
264 
265  //-| The remaining code is similar to setCylTheta, The reason for
266  //-| using a copy is so as to be able to change the messages in the
267  //-| ZMthrows to say eta rather than theta. Besides, we assumedly
268  //-| need not check for theta of 0 or PI.
269 
270  if ( (dx == 0) && (dy == 0) ) {
271  if (dz == 0) {
272  std::cerr << "Hep3Vector::setCylEta() - "
273  << "Attempt to set cylEta of zero vector -- vector is unchanged"
274  << std::endl;
275  return;
276  }
277  if (theta1 == 0) {
278  dz = std::fabs(dz);
279  return;
280  }
281  if (theta1 == CLHEP::pi) {
282  dz = -std::fabs(dz);
283  return;
284  }
285  std::cerr << "Hep3Vector::setCylEta() - "
286  << "Attempt set cylindrical eta of vector along Z axis "
287  << "to a non-trivial value, while keeping rho fixed -- "
288  << "will return zero vector" << std::endl;
289  dz = 0;
290  return;
291  }
292  double phi1 (getPhi());
293  double rho1 = getRho();
294  dz = rho1 / std::tan (theta1);
295  dy = rho1 * std::sin (phi1);
296  dx = rho1 * std::cos (phi1);
297 
298 } /* setCylEta */
299 
300 
301 Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
302 // if (c == 0) {
303 // std::cerr << "Hep3Vector::operator/ () - "
304 // << "Attempt to divide vector by 0 -- "
305 // << "will produce infinities and/or NANs" << std::endl;
306 // }
307  double oneOverC = 1.0/c;
308  return Hep3Vector ( v1.x() * oneOverC,
309  v1.y() * oneOverC,
310  v1.z() * oneOverC );
311 } /* v / c */
312 
313 Hep3Vector & Hep3Vector::operator/= (double c) {
314 // if (c == 0) {
315 // std::cerr << "Hep3Vector::operator/ () - "
316 // << "Attempt to do vector /= 0 -- "
317 // << "division by zero would produce infinite or NAN components"
318 // << std::endl;
319 // }
320  double oneOverC = 1.0/c;
321  dx *= oneOverC;
322  dy *= oneOverC;
323  dz *= oneOverC;
324  return *this;
325 }
326 
327 double Hep3Vector::tolerance = Hep3Vector::ToleranceTicks * 2.22045e-16;
328 
329 } // 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
static const G4double tolerance
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
static const G4double factor
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