Geant4  10.01.p01
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 
10 namespace 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 */
z
G4double z
Definition: TRTMaterials.hh:39
a
G4double a
Definition: TRTMaterials.hh:39
HepGeom::operator>>
std::istream & operator>>(std::istream &is, BasicVector3D< float > &a)
Definition: BasicVector3D.cc:114
FLT_MAX
#define FLT_MAX
Definition: templates.hh:99
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83