CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

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