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"
9 #include "CLHEP/Geometry/BasicVector3D.h"
10 
11 namespace HepGeom {
12   //--------------------------------------------------------------------------
13   template<>
pseudoRapidity() const14   float BasicVector3D<float>::pseudoRapidity() const {
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<>
setEta(float a)24   void BasicVector3D<float>::setEta(float a) {
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<>
angle(const BasicVector3D<float> & v) const37   float BasicVector3D<float>::angle(const BasicVector3D<float> & v) const {
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<>
rotateX(float a)50   BasicVector3D<float> & BasicVector3D<float>::rotateX(float a) {
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<>
rotateY(float a)59   BasicVector3D<float> & BasicVector3D<float>::rotateY(float a) {
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<>
rotateZ(float a)68   BasicVector3D<float> & BasicVector3D<float>::rotateZ(float a) {
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<>
77   BasicVector3D<float> &
rotate(float a,const BasicVector3D<float> & v)78   BasicVector3D<float>::rotate(float a, const BasicVector3D<float> & v) {
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 &
operator <<(std::ostream & os,const BasicVector3D<float> & a)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 &
operator >>(std::istream & is,BasicVector3D<float> & a)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<>
pseudoRapidity() const170   double BasicVector3D<double>::pseudoRapidity() const {
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<>
setEta(double a)180   void BasicVector3D<double>::setEta(double a) {
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<>
angle(const BasicVector3D<double> & v) const193   double BasicVector3D<double>::angle(const BasicVector3D<double> & v) const {
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<>
rotateX(double a)206   BasicVector3D<double> & BasicVector3D<double>::rotateX(double a) {
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<>
rotateY(double a)215   BasicVector3D<double> & BasicVector3D<double>::rotateY(double a) {
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<>
rotateZ(double a)224   BasicVector3D<double> & BasicVector3D<double>::rotateZ(double a) {
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<>
233   BasicVector3D<double> &
rotate(double a,const BasicVector3D<double> & v)234   BasicVector3D<double>::rotate(double a, const BasicVector3D<double> & v) {
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 &
operator <<(std::ostream & os,const BasicVector3D<double> & a)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 &
operator >>(std::istream & is,BasicVector3D<double> & a)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 */
324