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