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