Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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