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