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