| Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*-
2 // ------------------------------------------- 2 // ---------------------------------------------------------------------------
3 3
4 #include <cmath> 4 #include <cmath>
5 #include <float.h> 5 #include <float.h>
6 #include <iostream> 6 #include <iostream>
7 #include "CLHEP/Geometry/BasicVector3D.h" 7 #include "CLHEP/Geometry/BasicVector3D.h"
8 8
9 namespace HepGeom { 9 namespace HepGeom {
10 //------------------------------------------ 10 //--------------------------------------------------------------------------
11 template<> 11 template<>
12 float BasicVector3D<float>::pseudoRapidity() 12 float BasicVector3D<float>::pseudoRapidity() const {
13 float ma = mag(), dz = z(); 13 float ma = mag(), dz = z();
14 if (ma == 0) return 0; 14 if (ma == 0) return 0;
15 if (ma == dz) return FLT_MAX; 15 if (ma == dz) return FLT_MAX;
16 if (ma == -dz) return -FLT_MAX; 16 if (ma == -dz) return -FLT_MAX;
17 return 0.5*std::log((ma+dz)/(ma-dz)); 17 return 0.5*std::log((ma+dz)/(ma-dz));
18 } 18 }
19 19
20 //------------------------------------------ 20 //--------------------------------------------------------------------------
21 template<> 21 template<>
22 void BasicVector3D<float>::setEta(float a) { 22 void BasicVector3D<float>::setEta(float a) {
23 double ma = mag(); 23 double ma = mag();
24 if (ma == 0) return; 24 if (ma == 0) return;
25 double tanHalfTheta = std::exp(-a); 25 double tanHalfTheta = std::exp(-a);
26 double tanHalfTheta2 = tanHalfTheta * tanH 26 double tanHalfTheta2 = tanHalfTheta * tanHalfTheta;
27 double cosTheta1 = (1 - tanHalfTheta2 27 double cosTheta1 = (1 - tanHalfTheta2) / (1 + tanHalfTheta2);
28 double rh = ma * std::sqrt(1 - 28 double rh = ma * std::sqrt(1 - cosTheta1*cosTheta1);
29 double ph = phi(); 29 double ph = phi();
30 set(rh*std::cos(ph), rh*std::sin(ph), ma*c 30 set(rh*std::cos(ph), rh*std::sin(ph), ma*cosTheta1);
31 } 31 }
32 32
33 //------------------------------------------ 33 //--------------------------------------------------------------------------
34 template<> 34 template<>
35 float BasicVector3D<float>::angle(const Basi 35 float BasicVector3D<float>::angle(const BasicVector3D<float> & v) const {
36 double cosa = 0; 36 double cosa = 0;
37 double ptot = mag()*v.mag(); 37 double ptot = mag()*v.mag();
38 if(ptot > 0) { 38 if(ptot > 0) {
39 cosa = dot(v)/ptot; 39 cosa = dot(v)/ptot;
40 if(cosa > 1) cosa = 1; 40 if(cosa > 1) cosa = 1;
41 if(cosa < -1) cosa = -1; 41 if(cosa < -1) cosa = -1;
42 } 42 }
43 return std::acos(cosa); 43 return std::acos(cosa);
44 } 44 }
45 45
46 //------------------------------------------ 46 //--------------------------------------------------------------------------
47 template<> 47 template<>
48 BasicVector3D<float> & BasicVector3D<float>: 48 BasicVector3D<float> & BasicVector3D<float>::rotateX(float a) {
49 double sina = std::sin(a), cosa = std::cos 49 double sina = std::sin(a), cosa = std::cos(a), dy = y(), dz = z();
50 setY(dy*cosa-dz*sina); 50 setY(dy*cosa-dz*sina);
51 setZ(dz*cosa+dy*sina); 51 setZ(dz*cosa+dy*sina);
52 return *this; 52 return *this;
53 } 53 }
54 54
55 //------------------------------------------ 55 //--------------------------------------------------------------------------
56 template<> 56 template<>
57 BasicVector3D<float> & BasicVector3D<float>: 57 BasicVector3D<float> & BasicVector3D<float>::rotateY(float a) {
58 double sina = std::sin(a), cosa = std::cos 58 double sina = std::sin(a), cosa = std::cos(a), dz = z(), dx = x();
59 setZ(dz*cosa-dx*sina); 59 setZ(dz*cosa-dx*sina);
60 setX(dx*cosa+dz*sina); 60 setX(dx*cosa+dz*sina);
61 return *this; 61 return *this;
62 } 62 }
63 63
64 //------------------------------------------ 64 //--------------------------------------------------------------------------
65 template<> 65 template<>
66 BasicVector3D<float> & BasicVector3D<float>: 66 BasicVector3D<float> & BasicVector3D<float>::rotateZ(float a) {
67 double sina = std::sin(a), cosa = std::cos 67 double sina = std::sin(a), cosa = std::cos(a), dx = x(), dy = y();
68 setX(dx*cosa-dy*sina); 68 setX(dx*cosa-dy*sina);
69 setY(dy*cosa+dx*sina); 69 setY(dy*cosa+dx*sina);
70 return *this; 70 return *this;
71 } 71 }
72 72
73 //------------------------------------------ 73 //--------------------------------------------------------------------------
74 template<> 74 template<>
75 BasicVector3D<float> & 75 BasicVector3D<float> &
76 BasicVector3D<float>::rotate(float a, const 76 BasicVector3D<float>::rotate(float a, const BasicVector3D<float> & v) {
77 if (a == 0) return *this; 77 if (a == 0) return *this;
78 double cx = v.x(), cy = v.y(), cz = v.z(); 78 double cx = v.x(), cy = v.y(), cz = v.z();
79 double ll = std::sqrt(cx*cx + cy*cy + cz*c 79 double ll = std::sqrt(cx*cx + cy*cy + cz*cz);
80 if (ll == 0) { 80 if (ll == 0) {
81 std::cerr << "BasicVector<float>::rotate 81 std::cerr << "BasicVector<float>::rotate() : zero axis" << std::endl;
82 return *this; 82 return *this;
83 } 83 }
84 double cosa = std::cos(a), sina = std::sin 84 double cosa = std::cos(a), sina = std::sin(a);
85 cx /= ll; cy /= ll; cz /= ll; 85 cx /= ll; cy /= ll; cz /= ll;
86 86
87 double xx = cosa + (1-cosa)*cx*cx; 87 double xx = cosa + (1-cosa)*cx*cx;
88 double xy = (1-cosa)*cx*cy - sina*c 88 double xy = (1-cosa)*cx*cy - sina*cz;
89 double xz = (1-cosa)*cx*cz + sina*c 89 double xz = (1-cosa)*cx*cz + sina*cy;
90 90
91 double yx = (1-cosa)*cy*cx + sina*c 91 double yx = (1-cosa)*cy*cx + sina*cz;
92 double yy = cosa + (1-cosa)*cy*cy; 92 double yy = cosa + (1-cosa)*cy*cy;
93 double yz = (1-cosa)*cy*cz - sina*c 93 double yz = (1-cosa)*cy*cz - sina*cx;
94 94
95 double zx = (1-cosa)*cz*cx - sina*c 95 double zx = (1-cosa)*cz*cx - sina*cy;
96 double zy = (1-cosa)*cz*cy + sina*c 96 double zy = (1-cosa)*cz*cy + sina*cx;
97 double zz = cosa + (1-cosa)*cz*cz; 97 double zz = cosa + (1-cosa)*cz*cz;
98 98
99 cx = x(); cy = y(); cz = z(); 99 cx = x(); cy = y(); cz = z();
100 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, 100 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz);
101 return *this; 101 return *this;
102 } 102 }
103 103
104 //------------------------------------------ 104 //--------------------------------------------------------------------------
105 std::ostream & 105 std::ostream &
106 operator<<(std::ostream & os, const BasicVec 106 operator<<(std::ostream & os, const BasicVector3D<float> & a)
107 { 107 {
108 return os << "(" << a.x() << "," << a.y() 108 return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")";
109 } 109 }
110 110
111 //------------------------------------------ 111 //--------------------------------------------------------------------------
112 std::istream & 112 std::istream &
113 operator>> (std::istream & is, BasicVector3D 113 operator>> (std::istream & is, BasicVector3D<float> & a)
114 { 114 {
115 // Required format is ( a, b, c ) that is, 115 // Required format is ( a, b, c ) that is, three numbers, preceded by
116 // (, followed by ), and separated by comm 116 // (, followed by ), and separated by commas. The three numbers are
117 // taken as x, y, z. 117 // taken as x, y, z.
118 118
119 float x, y, z; 119 float x, y, z;
120 char c; 120 char c;
121 121
122 is >> std::ws >> c; 122 is >> std::ws >> c;
123 // ws is defined to invoke eatwhite(istrea 123 // ws is defined to invoke eatwhite(istream & )
124 // see (Stroustrup gray book) page 333 and 124 // see (Stroustrup gray book) page 333 and 345.
125 if (is.fail() || c != '(' ) { 125 if (is.fail() || c != '(' ) {
126 std::cerr 126 std::cerr
127 << "Could not find required opening parenthe 127 << "Could not find required opening parenthesis "
128 << "in input of a BasicVector3D<float>" 128 << "in input of a BasicVector3D<float>"
129 << std::endl; 129 << std::endl;
130 return is; 130 return is;
131 } 131 }
132 132
133 is >> x >> std::ws >> c; 133 is >> x >> std::ws >> c;
134 if (is.fail() || c != ',' ) { 134 if (is.fail() || c != ',' ) {
135 std::cerr 135 std::cerr
136 << "Could not find x value and required trai 136 << "Could not find x value and required trailing comma "
137 << "in input of a BasicVector3D<float>" 137 << "in input of a BasicVector3D<float>"
138 << std::endl; 138 << std::endl;
139 return is; 139 return is;
140 } 140 }
141 141
142 is >> y >> std::ws >> c; 142 is >> y >> std::ws >> c;
143 if (is.fail() || c != ',' ) { 143 if (is.fail() || c != ',' ) {
144 std::cerr 144 std::cerr
145 << "Could not find y value and required trai 145 << "Could not find y value and required trailing comma "
146 << "in input of a BasicVector3D<float>" 146 << "in input of a BasicVector3D<float>"
147 << std::endl; 147 << std::endl;
148 return is; 148 return is;
149 } 149 }
150 150
151 is >> z >> std::ws >> c; 151 is >> z >> std::ws >> c;
152 if (is.fail() || c != ')' ) { 152 if (is.fail() || c != ')' ) {
153 std::cerr 153 std::cerr
154 << "Could not find z value and required clos 154 << "Could not find z value and required close parenthesis "
155 << "in input of a BasicVector3D<float>" 155 << "in input of a BasicVector3D<float>"
156 << std::endl; 156 << std::endl;
157 return is; 157 return is;
158 } 158 }
159 159
160 a.setX(x); 160 a.setX(x);
161 a.setY(y); 161 a.setY(y);
162 a.setZ(z); 162 a.setZ(z);
163 return is; 163 return is;
164 } 164 }
165 165
166 //------------------------------------------ 166 //--------------------------------------------------------------------------
167 template<> 167 template<>
168 double BasicVector3D<double>::pseudoRapidity 168 double BasicVector3D<double>::pseudoRapidity() const {
169 double ma = mag(), dz = z(); 169 double ma = mag(), dz = z();
170 if (ma == 0) return 0; 170 if (ma == 0) return 0;
171 if (ma == dz) return DBL_MAX; 171 if (ma == dz) return DBL_MAX;
172 if (ma == -dz) return -DBL_MAX; 172 if (ma == -dz) return -DBL_MAX;
173 return 0.5*std::log((ma+dz)/(ma-dz)); 173 return 0.5*std::log((ma+dz)/(ma-dz));
174 } 174 }
175 175
176 //------------------------------------------ 176 //--------------------------------------------------------------------------
177 template<> 177 template<>
178 void BasicVector3D<double>::setEta(double a) 178 void BasicVector3D<double>::setEta(double a) {
179 double ma = mag(); 179 double ma = mag();
180 if (ma == 0) return; 180 if (ma == 0) return;
181 double tanHalfTheta = std::exp(-a); 181 double tanHalfTheta = std::exp(-a);
182 double tanHalfTheta2 = tanHalfTheta * tanH 182 double tanHalfTheta2 = tanHalfTheta * tanHalfTheta;
183 double cosTheta1 = (1 - tanHalfTheta2 183 double cosTheta1 = (1 - tanHalfTheta2) / (1 + tanHalfTheta2);
184 double rh = ma * std::sqrt(1 - 184 double rh = ma * std::sqrt(1 - cosTheta1*cosTheta1);
185 double ph = phi(); 185 double ph = phi();
186 set(rh*std::cos(ph), rh*std::sin(ph), ma*c 186 set(rh*std::cos(ph), rh*std::sin(ph), ma*cosTheta1);
187 } 187 }
188 188
189 //------------------------------------------ 189 //--------------------------------------------------------------------------
190 template<> 190 template<>
191 double BasicVector3D<double>::angle(const Ba 191 double BasicVector3D<double>::angle(const BasicVector3D<double> & v) const {
192 double cosa = 0; 192 double cosa = 0;
193 double ptot = mag()*v.mag(); 193 double ptot = mag()*v.mag();
194 if(ptot > 0) { 194 if(ptot > 0) {
195 cosa = dot(v)/ptot; 195 cosa = dot(v)/ptot;
196 if(cosa > 1) cosa = 1; 196 if(cosa > 1) cosa = 1;
197 if(cosa < -1) cosa = -1; 197 if(cosa < -1) cosa = -1;
198 } 198 }
199 return std::acos(cosa); 199 return std::acos(cosa);
200 } 200 }
201 201
202 //------------------------------------------ 202 //--------------------------------------------------------------------------
203 template<> 203 template<>
204 BasicVector3D<double> & BasicVector3D<double 204 BasicVector3D<double> & BasicVector3D<double>::rotateX(double a) {
205 double sina = std::sin(a), cosa = std::cos 205 double sina = std::sin(a), cosa = std::cos(a), dy = y(), dz = z();
206 setY(dy*cosa-dz*sina); 206 setY(dy*cosa-dz*sina);
207 setZ(dz*cosa+dy*sina); 207 setZ(dz*cosa+dy*sina);
208 return *this; 208 return *this;
209 } 209 }
210 210
211 //------------------------------------------ 211 //--------------------------------------------------------------------------
212 template<> 212 template<>
213 BasicVector3D<double> & BasicVector3D<double 213 BasicVector3D<double> & BasicVector3D<double>::rotateY(double a) {
214 double sina = std::sin(a), cosa = std::cos 214 double sina = std::sin(a), cosa = std::cos(a), dz = z(), dx = x();
215 setZ(dz*cosa-dx*sina); 215 setZ(dz*cosa-dx*sina);
216 setX(dx*cosa+dz*sina); 216 setX(dx*cosa+dz*sina);
217 return *this; 217 return *this;
218 } 218 }
219 219
220 //------------------------------------------ 220 //--------------------------------------------------------------------------
221 template<> 221 template<>
222 BasicVector3D<double> & BasicVector3D<double 222 BasicVector3D<double> & BasicVector3D<double>::rotateZ(double a) {
223 double sina = std::sin(a), cosa = std::cos 223 double sina = std::sin(a), cosa = std::cos(a), dx = x(), dy = y();
224 setX(dx*cosa-dy*sina); 224 setX(dx*cosa-dy*sina);
225 setY(dy*cosa+dx*sina); 225 setY(dy*cosa+dx*sina);
226 return *this; 226 return *this;
227 } 227 }
228 228
229 //------------------------------------------ 229 //--------------------------------------------------------------------------
230 template<> 230 template<>
231 BasicVector3D<double> & 231 BasicVector3D<double> &
232 BasicVector3D<double>::rotate(double a, cons 232 BasicVector3D<double>::rotate(double a, const BasicVector3D<double> & v) {
233 if (a == 0) return *this; 233 if (a == 0) return *this;
234 double cx = v.x(), cy = v.y(), cz = v.z(); 234 double cx = v.x(), cy = v.y(), cz = v.z();
235 double ll = std::sqrt(cx*cx + cy*cy + cz*c 235 double ll = std::sqrt(cx*cx + cy*cy + cz*cz);
236 if (ll == 0) { 236 if (ll == 0) {
237 std::cerr << "BasicVector<double>::rotat 237 std::cerr << "BasicVector<double>::rotate() : zero axis" << std::endl;
238 return *this; 238 return *this;
239 } 239 }
240 double cosa = std::cos(a), sina = std::sin 240 double cosa = std::cos(a), sina = std::sin(a);
241 cx /= ll; cy /= ll; cz /= ll; 241 cx /= ll; cy /= ll; cz /= ll;
242 242
243 double xx = cosa + (1-cosa)*cx*cx; 243 double xx = cosa + (1-cosa)*cx*cx;
244 double xy = (1-cosa)*cx*cy - sina*c 244 double xy = (1-cosa)*cx*cy - sina*cz;
245 double xz = (1-cosa)*cx*cz + sina*c 245 double xz = (1-cosa)*cx*cz + sina*cy;
246 246
247 double yx = (1-cosa)*cy*cx + sina*c 247 double yx = (1-cosa)*cy*cx + sina*cz;
248 double yy = cosa + (1-cosa)*cy*cy; 248 double yy = cosa + (1-cosa)*cy*cy;
249 double yz = (1-cosa)*cy*cz - sina*c 249 double yz = (1-cosa)*cy*cz - sina*cx;
250 250
251 double zx = (1-cosa)*cz*cx - sina*c 251 double zx = (1-cosa)*cz*cx - sina*cy;
252 double zy = (1-cosa)*cz*cy + sina*c 252 double zy = (1-cosa)*cz*cy + sina*cx;
253 double zz = cosa + (1-cosa)*cz*cz; 253 double zz = cosa + (1-cosa)*cz*cz;
254 254
255 cx = x(); cy = y(); cz = z(); 255 cx = x(); cy = y(); cz = z();
256 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, 256 set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz);
257 return *this; 257 return *this;
258 } 258 }
259 259
260 //------------------------------------------ 260 //--------------------------------------------------------------------------
261 std::ostream & 261 std::ostream &
262 operator<< (std::ostream & os, const BasicVe 262 operator<< (std::ostream & os, const BasicVector3D<double> & a)
263 { 263 {
264 return os << "(" << a.x() << "," << a.y() 264 return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")";
265 } 265 }
266 266
267 //------------------------------------------ 267 //--------------------------------------------------------------------------
268 std::istream & 268 std::istream &
269 operator>> (std::istream & is, BasicVector3D 269 operator>> (std::istream & is, BasicVector3D<double> & a)
270 { 270 {
271 // Required format is ( a, b, c ) that is, 271 // Required format is ( a, b, c ) that is, three numbers, preceded by
272 // (, followed by ), and separated by comm 272 // (, followed by ), and separated by commas. The three numbers are
273 // taken as x, y, z. 273 // taken as x, y, z.
274 274
275 double x, y, z; 275 double x, y, z;
276 char c; 276 char c;
277 277
278 is >> std::ws >> c; 278 is >> std::ws >> c;
279 // ws is defined to invoke eatwhite(istrea 279 // ws is defined to invoke eatwhite(istream & )
280 // see (Stroustrup gray book) page 333 and 280 // see (Stroustrup gray book) page 333 and 345.
281 if (is.fail() || c != '(' ) { 281 if (is.fail() || c != '(' ) {
282 std::cerr 282 std::cerr
283 << "Could not find required opening parenthe 283 << "Could not find required opening parenthesis "
284 << "in input of a BasicVector3D<double>" 284 << "in input of a BasicVector3D<double>"
285 << std::endl; 285 << std::endl;
286 return is; 286 return is;
287 } 287 }
288 288
289 is >> x >> std::ws >> c; 289 is >> x >> std::ws >> c;
290 if (is.fail() || c != ',' ) { 290 if (is.fail() || c != ',' ) {
291 std::cerr 291 std::cerr
292 << "Could not find x value and required trai 292 << "Could not find x value and required trailing comma "
293 << "in input of a BasicVector3D<double>" 293 << "in input of a BasicVector3D<double>"
294 << std::endl; 294 << std::endl;
295 return is; 295 return is;
296 } 296 }
297 297
298 is >> y >> std::ws >> c; 298 is >> y >> std::ws >> c;
299 if (is.fail() || c != ',' ) { 299 if (is.fail() || c != ',' ) {
300 std::cerr 300 std::cerr
301 << "Could not find y value and required trai 301 << "Could not find y value and required trailing comma "
302 << "in input of a BasicVector3D<double>" 302 << "in input of a BasicVector3D<double>"
303 << std::endl; 303 << std::endl;
304 return is; 304 return is;
305 } 305 }
306 306
307 is >> z >> std::ws >> c; 307 is >> z >> std::ws >> c;
308 if (is.fail() || c != ')' ) { 308 if (is.fail() || c != ')' ) {
309 std::cerr 309 std::cerr
310 << "Could not find z value and required clos 310 << "Could not find z value and required close parenthesis "
311 << "in input of a BasicVector3D<double>" 311 << "in input of a BasicVector3D<double>"
312 << std::endl; 312 << std::endl;
313 return is; 313 return is;
314 } 314 }
315 315
316 a.setX(x); 316 a.setX(x);
317 a.setY(y); 317 a.setY(y);
318 a.setZ(z); 318 a.setZ(z);
319 return is; 319 return is;
320 } 320 }
321 } /* namespace HepGeom */ 321 } /* namespace HepGeom */
322 322