Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TRotation.cxx
Go to the documentation of this file.
1 // @(#)root/physics:$Id$
2 // Author: Peter Malzacher 19/06/99
3 
4 /** \class TRotation
5  \ingroup Physics
6 
7 The TRotation class describes a rotation of objects of the TVector3 class.
8 It is a 3*3 matrix of Double_t:
9 
10 ~~~
11 | xx xy xz |
12 | yx yy yz |
13 | zx zy zz |
14 ~~~
15 
16 It describes a so called active rotation, i.e. rotation of objects inside
17 a static system of coordinates. In case you want to rotate the frame and
18 want to know the coordinates of objects in the rotated system, you should
19 apply the inverse rotation to the objects. If you want to transform coordinates
20 from the rotated frame to the original frame you have to apply the direct
21 transformation.
22 
23 A rotation around a specified axis means counterclockwise rotation around
24 the positive direction of the axis.
25 
26 
27 ### Declaration, Access, Comparisons
28 
29 ~~~
30  TRotation r; // r initialized as identity
31  TRotation m(r); // m = r
32 ~~~
33 
34 There is no direct way to set the matrix elements - to ensure that
35 a TRotation object always describes a real rotation. But you can get the
36 values by the member functions XX()..ZZ() or the (,)
37 operator:
38 
39 ~~~
40  Double_t xx = r.XX(); // the same as xx=r(0,0)
41  xx = r(0,0);
42 
43  if (r==m) {...} // test for equality
44  if (r!=m) {..} // test for inequality
45  if (r.IsIdentity()) {...} // test for identity
46 ~~~
47 
48 ### Rotation around axes
49 The following matrices describe counterclockwise rotations around coordinate
50 axes
51 
52 ~~~
53  | 1 0 0 |
54 Rx(a) = | 0 cos(a) -sin(a) |
55  | 0 sin(a) cos(a) |
56 
57  | cos(a) 0 sin(a) |
58 Ry(a) = | 0 1 0 |
59  | -sin(a) 0 cos(a) |
60 
61  | cos(a) -sin(a) 0 |
62 Rz(a) = | sin(a) cos(a) 0 |
63  | 0 0 1 |
64 ~~~
65 
66 and are implemented as member functions RotateX(), RotateY() and RotateZ():
67 
68 ~~~
69  r.RotateX(TMath::Pi()); // rotation around the x-axis
70 ~~~
71 
72 ### Rotation around arbitrary axis
73 The member function Rotate() allows to rotate around an arbitrary vector
74 (not necessary a unit one) and returns the result.
75 
76 ~~~
77  r.Rotate(TMath::Pi()/3,TVector3(3,4,5));
78 ~~~
79 
80 It is possible to find a unit vector and an angle, which describe the
81 same rotation as the current one:
82 
83 ~~~
84  Double_t angle;
85  TVector3 axis;
86  r.GetAngleAxis(angle,axis);
87 ~~~
88 
89 ### Rotation of local axes
90 Member function RotateAxes() adds a rotation of local axes to
91 the current rotation and returns the result:
92 
93 ~~~
94  TVector3 newX(0,1,0);
95  TVector3 newY(0,0,1);
96  TVector3 newZ(1,0,0);
97  a.RotateAxes(newX,newY,newZ);
98 ~~~
99 
100 Member functions ThetaX(), ThetaY(), ThetaZ(),
101 PhiX(), PhiY(),PhiZ() return azimuth and polar
102 angles of the rotated axes:
103 
104 ~~~
105  Double_t tx,ty,tz,px,py,pz;
106  tx= a.ThetaX();
107  ...
108  pz= a.PhiZ();
109 ~~~
110 
111 ### Setting The Rotations
112 The member function SetToIdentity() will set the rotation object
113 to the identity (no rotation).
114 
115 With a minor caveat, the Euler angles of the rotation may be set using
116 SetXEulerAngles() or individually set with SetXPhi(),
117 SetXTheta(), and SetXPsi(). These routines set the Euler
118 angles using the X-convention which is defined by a rotation about the Z-axis,
119 about the new X-axis, and about the new Z-axis. This is the convention used
120 in Landau and Lifshitz, Goldstein and other common physics texts. The
121 Y-convention Euler angles can be set with SetYEulerAngles(),
122 SetYPhi(), SetYTheta(), and SetYPsi(). The caveat
123 is that Euler angles usually define the rotation of the new coordinate system
124 with respect to the original system, however, the TRotation class specifies
125 the rotation of the object in the original system (an active rotation). To
126 recover the usual Euler rotations (ie. rotate the system not the object), you
127 must take the inverse of the rotation.
128 
129 The member functions SetXAxis(), SetYAxis(), and
130 SetZAxis() will create a rotation which rotates the requested axis
131 of the object to be parallel to a vector. If used with one argument, the
132 rotation about that axis is arbitrary. If used with two arguments, the
133 second variable defines the XY, YZ, or ZX
134 respectively.
135 
136 
137 ### Inverse rotation
138 
139 ~~~
140  TRotation a,b;
141  ...
142  b = a.Inverse(); // b is inverse of a, a is unchanged
143  b = a.Invert(); // invert a and set b = a
144 ~~~
145 
146 ### Compound Rotations
147 The operator * has been implemented in a way that follows the
148 mathematical notation of a product of the two matrices which describe the
149 two consecutive rotations. Therefore the second rotation should be placed
150 first:
151 
152 ~~~
153  r = r2 * r1;
154 ~~~
155 
156 ### Rotation of TVector3
157 The TRotation class provides an operator * which allows to express
158 a rotation of a TVector3 analog to the mathematical notation
159 
160 ~~~
161  | x' | | xx xy xz | | x |
162  | y' | = | yx yy yz | | y |
163  | z' | | zx zy zz | | z |
164 ~~~
165 
166 e.g.:
167 
168 ~~~
169  TVector3 v(1,1,1);
170  v = r * v;
171 ~~~
172 
173 You can also use the Transform() member function or the operator *= of the
174 TVector3 class:
175 
176 ~~~
177  TVector3 v;
178  TRotation r;
179  v.Transform(r);
180  v *= r; //Attention v = r * v
181 ~~~
182 */
183 
184 #include "TRotation.h"
185 #include "TMath.h"
186 #include "TQuaternion.h"
187 
188 ClassImp(TRotation);
189 
190 #define TOLERANCE (1.0E-6)
191 
192 ////////////////////////////////////////////////////////////////////////////////
193 /// Constructor.
194 
195 TRotation::TRotation()
196 : fxx(1.0), fxy(0.0), fxz(0.0), fyx(0.0), fyy(1.0), fyz(0.0),
197  fzx(0.0), fzy(0.0), fzz(1.0) {}
198 
199 ////////////////////////////////////////////////////////////////////////////////
200 /// Constructor.
201 
202 TRotation::TRotation(const TRotation & m) : TObject(m),
203  fxx(m.fxx), fxy(m.fxy), fxz(m.fxz), fyx(m.fyx), fyy(m.fyy), fyz(m.fyz),
204  fzx(m.fzx), fzy(m.fzy), fzz(m.fzz) {}
205 
206 ////////////////////////////////////////////////////////////////////////////////
207 /// Constructor.
208 
209 TRotation::TRotation(Double_t mxx, Double_t mxy, Double_t mxz,
210  Double_t myx, Double_t myy, Double_t myz,
211  Double_t mzx, Double_t mzy, Double_t mzz)
212 : fxx(mxx), fxy(mxy), fxz(mxz), fyx(myx), fyy(myy), fyz(myz),
213  fzx(mzx), fzy(mzy), fzz(mzz) {}
214 
215 ////////////////////////////////////////////////////////////////////////////////
216 /// Dereferencing operator const.
217 
218 Double_t TRotation::operator() (int i, int j) const {
219  if (i == 0) {
220  if (j == 0) { return fxx; }
221  if (j == 1) { return fxy; }
222  if (j == 2) { return fxz; }
223  } else if (i == 1) {
224  if (j == 0) { return fyx; }
225  if (j == 1) { return fyy; }
226  if (j == 2) { return fyz; }
227  } else if (i == 2) {
228  if (j == 0) { return fzx; }
229  if (j == 1) { return fzy; }
230  if (j == 2) { return fzz; }
231  }
232 
233  Warning("operator()(i,j)", "bad indices (%d , %d)",i,j);
234 
235  return 0.0;
236 }
237 
238 ////////////////////////////////////////////////////////////////////////////////
239 /// Multiplication operator.
240 
241 TRotation TRotation::operator* (const TRotation & b) const {
242  return TRotation(fxx*b.fxx + fxy*b.fyx + fxz*b.fzx,
243  fxx*b.fxy + fxy*b.fyy + fxz*b.fzy,
244  fxx*b.fxz + fxy*b.fyz + fxz*b.fzz,
245  fyx*b.fxx + fyy*b.fyx + fyz*b.fzx,
246  fyx*b.fxy + fyy*b.fyy + fyz*b.fzy,
247  fyx*b.fxz + fyy*b.fyz + fyz*b.fzz,
248  fzx*b.fxx + fzy*b.fyx + fzz*b.fzx,
249  fzx*b.fxy + fzy*b.fyy + fzz*b.fzy,
250  fzx*b.fxz + fzy*b.fyz + fzz*b.fzz);
251 }
252 
253 ////////////////////////////////////////////////////////////////////////////////
254 /// Constructor for a rotation based on a Quaternion
255 /// if magnitude of quaternion is null, creates identity rotation
256 /// if quaternion is non-unit, creates rotation corresponding to the normalized (unit) quaternion
257 
258 TRotation::TRotation(const TQuaternion & Q) {
259 
260  double two_r2 = 2 * Q.fRealPart * Q.fRealPart;
261  double two_x2 = 2 * Q.fVectorPart.X() * Q.fVectorPart.X();
262  double two_y2 = 2 * Q.fVectorPart.Y() * Q.fVectorPart.Y();
263  double two_z2 = 2 * Q.fVectorPart.Z() * Q.fVectorPart.Z();
264  double two_xy = 2 * Q.fVectorPart.X() * Q.fVectorPart.Y();
265  double two_xz = 2 * Q.fVectorPart.X() * Q.fVectorPart.Z();
266  double two_xr = 2 * Q.fVectorPart.X() * Q.fRealPart;
267  double two_yz = 2 * Q.fVectorPart.Y() * Q.fVectorPart.Z();
268  double two_yr = 2 * Q.fVectorPart.Y() * Q.fRealPart;
269  double two_zr = 2 * Q.fVectorPart.Z() * Q.fRealPart;
270 
271  // protect agains zero quaternion
272  double mag2 = Q.QMag2();
273  if (mag2 > 0) {
274 
275  // diago + identity
276  fxx = two_r2 + two_x2;
277  fyy = two_r2 + two_y2;
278  fzz = two_r2 + two_z2;
279 
280  // line 0 column 1 and conjugate
281  fxy = two_xy - two_zr;
282  fyx = two_xy + two_zr;
283 
284  // line 0 column 2 and conjugate
285  fxz = two_xz + two_yr;
286  fzx = two_xz - two_yr;
287 
288  // line 1 column 2 and conjugate
289  fyz = two_yz - two_xr;
290  fzy = two_yz + two_xr;
291 
292  // protect agains non-unit quaternion
293  if (TMath::Abs(mag2-1) > 1e-10) {
294  fxx /= mag2;
295  fyy /= mag2;
296  fzz /= mag2;
297  fxy /= mag2;
298  fyx /= mag2;
299  fxz /= mag2;
300  fzx /= mag2;
301  fyz /= mag2;
302  fzy /= mag2;
303  }
304 
305  // diago : remove identity
306  fxx -= 1;
307  fyy -= 1;
308  fzz -= 1;
309 
310 
311  } else {
312  // Identity
313 
314  fxx = fyy = fzz = 1;
315  fxy = fyx = fxz = fzx = fyz = fzy = 0;
316 
317  }
318 
319 }
320 
321 ////////////////////////////////////////////////////////////////////////////////
322 /// Rotate along an axis.
323 
324 TRotation & TRotation::Rotate(Double_t a, const TVector3& axis) {
325  if (a != 0.0) {
326  Double_t ll = axis.Mag();
327  if (ll == 0.0) {
328  Warning("Rotate(angle,axis)"," zero axis");
329  } else {
330  Double_t sa = TMath::Sin(a), ca = TMath::Cos(a);
331  Double_t dx = axis.X()/ll, dy = axis.Y()/ll, dz = axis.Z()/ll;
332  TRotation m(
333  ca+(1-ca)*dx*dx, (1-ca)*dx*dy-sa*dz, (1-ca)*dx*dz+sa*dy,
334  (1-ca)*dy*dx+sa*dz, ca+(1-ca)*dy*dy, (1-ca)*dy*dz-sa*dx,
335  (1-ca)*dz*dx-sa*dy, (1-ca)*dz*dy+sa*dx, ca+(1-ca)*dz*dz );
336  Transform(m);
337  }
338  }
339  return *this;
340 }
341 
342 ////////////////////////////////////////////////////////////////////////////////
343 /// Rotate around x.
344 
345 TRotation & TRotation::RotateX(Double_t a) {
346  Double_t c = TMath::Cos(a);
347  Double_t s = TMath::Sin(a);
348  Double_t x = fyx, y = fyy, z = fyz;
349  fyx = c*x - s*fzx;
350  fyy = c*y - s*fzy;
351  fyz = c*z - s*fzz;
352  fzx = s*x + c*fzx;
353  fzy = s*y + c*fzy;
354  fzz = s*z + c*fzz;
355  return *this;
356 }
357 
358 ////////////////////////////////////////////////////////////////////////////////
359 /// Rotate around y.
360 
361 TRotation & TRotation::RotateY(Double_t a){
362  Double_t c = TMath::Cos(a);
363  Double_t s = TMath::Sin(a);
364  Double_t x = fzx, y = fzy, z = fzz;
365  fzx = c*x - s*fxx;
366  fzy = c*y - s*fxy;
367  fzz = c*z - s*fxz;
368  fxx = s*x + c*fxx;
369  fxy = s*y + c*fxy;
370  fxz = s*z + c*fxz;
371  return *this;
372 }
373 
374 ////////////////////////////////////////////////////////////////////////////////
375 /// Rotate around z.
376 
377 TRotation & TRotation::RotateZ(Double_t a) {
378  Double_t c = TMath::Cos(a);
379  Double_t s = TMath::Sin(a);
380  Double_t x = fxx, y = fxy, z = fxz;
381  fxx = c*x - s*fyx;
382  fxy = c*y - s*fyy;
383  fxz = c*z - s*fyz;
384  fyx = s*x + c*fyx;
385  fyy = s*y + c*fyy;
386  fyz = s*z + c*fyz;
387  return *this;
388 }
389 
390 ////////////////////////////////////////////////////////////////////////////////
391 /// Rotate axes.
392 
393 TRotation & TRotation::RotateAxes(const TVector3 &newX,
394  const TVector3 &newY,
395  const TVector3 &newZ) {
396  Double_t del = 0.001;
397  TVector3 w = newX.Cross(newY);
398 
399  if (TMath::Abs(newZ.X()-w.X()) > del ||
400  TMath::Abs(newZ.Y()-w.Y()) > del ||
401  TMath::Abs(newZ.Z()-w.Z()) > del ||
402  TMath::Abs(newX.Mag2()-1.) > del ||
403  TMath::Abs(newY.Mag2()-1.) > del ||
404  TMath::Abs(newZ.Mag2()-1.) > del ||
405  TMath::Abs(newX.Dot(newY)) > del ||
406  TMath::Abs(newY.Dot(newZ)) > del ||
407  TMath::Abs(newZ.Dot(newX)) > del) {
408  Warning("RotateAxes","bad axis vectors");
409  return *this;
410  } else {
411  return Transform(TRotation(newX.X(), newY.X(), newZ.X(),
412  newX.Y(), newY.Y(), newZ.Y(),
413  newX.Z(), newY.Z(), newZ.Z()));
414  }
415 }
416 
417 ////////////////////////////////////////////////////////////////////////////////
418 /// Return Phi.
419 
420 Double_t TRotation::PhiX() const {
421  return (fyx == 0.0 && fxx == 0.0) ? 0.0 : TMath::ATan2(fyx,fxx);
422 }
423 
424 ////////////////////////////////////////////////////////////////////////////////
425 /// Return Phi.
426 
427 Double_t TRotation::PhiY() const {
428  return (fyy == 0.0 && fxy == 0.0) ? 0.0 : TMath::ATan2(fyy,fxy);
429 }
430 
431 ////////////////////////////////////////////////////////////////////////////////
432 /// Return Phi.
433 
434 Double_t TRotation::PhiZ() const {
435  return (fyz == 0.0 && fxz == 0.0) ? 0.0 : TMath::ATan2(fyz,fxz);
436 }
437 
438 ////////////////////////////////////////////////////////////////////////////////
439 /// Return Theta.
440 
441 Double_t TRotation::ThetaX() const {
442  return TMath::ACos(fzx);
443 }
444 
445 ////////////////////////////////////////////////////////////////////////////////
446 /// Return Theta.
447 
448 Double_t TRotation::ThetaY() const {
449  return TMath::ACos(fzy);
450 }
451 
452 ////////////////////////////////////////////////////////////////////////////////
453 /// Return Theta.
454 
455 Double_t TRotation::ThetaZ() const {
456  return TMath::ACos(fzz);
457 }
458 
459 ////////////////////////////////////////////////////////////////////////////////
460 /// Rotation defined by an angle and a vector.
461 
462 void TRotation::AngleAxis(Double_t &angle, TVector3 &axis) const {
463  Double_t cosa = 0.5*(fxx+fyy+fzz-1);
464  Double_t cosa1 = 1-cosa;
465  if (cosa1 <= 0) {
466  angle = 0;
467  axis = TVector3(0,0,1);
468  } else {
469  Double_t x=0, y=0, z=0;
470  if (fxx > cosa) x = TMath::Sqrt((fxx-cosa)/cosa1);
471  if (fyy > cosa) y = TMath::Sqrt((fyy-cosa)/cosa1);
472  if (fzz > cosa) z = TMath::Sqrt((fzz-cosa)/cosa1);
473  if (fzy < fyz) x = -x;
474  if (fxz < fzx) y = -y;
475  if (fyx < fxy) z = -z;
476  angle = TMath::ACos(cosa);
477  axis = TVector3(x,y,z);
478  }
479 }
480 
481 ////////////////////////////////////////////////////////////////////////////////
482 /// Rotate using the x-convention (Landau and Lifshitz, Goldstein, &c) by
483 /// doing the explicit rotations. This is slightly less efficient than
484 /// directly applying the rotation, but makes the code much clearer. My
485 /// presumption is that this code is not going to be a speed bottle neck.
486 
487 TRotation & TRotation::SetXEulerAngles(Double_t phi,
488  Double_t theta,
489  Double_t psi) {
490  SetToIdentity();
491  RotateZ(phi);
492  RotateX(theta);
493  RotateZ(psi);
494 
495  return *this;
496 }
497 
498 ////////////////////////////////////////////////////////////////////////////////
499 /// Rotate using the y-convention.
500 
501 TRotation & TRotation::SetYEulerAngles(Double_t phi,
502  Double_t theta,
503  Double_t psi) {
504  SetToIdentity();
505  RotateZ(phi);
506  RotateY(theta);
507  RotateZ(psi);
508  return *this;
509 }
510 
511 ////////////////////////////////////////////////////////////////////////////////
512 /// Rotate using the x-convention.
513 
514 TRotation & TRotation::RotateXEulerAngles(Double_t phi,
515  Double_t theta,
516  Double_t psi) {
517  TRotation euler;
518  euler.SetXEulerAngles(phi,theta,psi);
519  return Transform(euler);
520 }
521 
522 ////////////////////////////////////////////////////////////////////////////////
523 /// Rotate using the y-convention.
524 
525 TRotation & TRotation::RotateYEulerAngles(Double_t phi,
526  Double_t theta,
527  Double_t psi) {
528  TRotation euler;
529  euler.SetYEulerAngles(phi,theta,psi);
530  return Transform(euler);
531 }
532 
533 ////////////////////////////////////////////////////////////////////////////////
534 /// Set XPhi.
535 
536 void TRotation::SetXPhi(Double_t phi) {
537  SetXEulerAngles(phi,GetXTheta(),GetXPsi());
538 }
539 
540 ////////////////////////////////////////////////////////////////////////////////
541 /// Set XTheta.
542 
543 void TRotation::SetXTheta(Double_t theta) {
544  SetXEulerAngles(GetXPhi(),theta,GetXPsi());
545 }
546 
547 ////////////////////////////////////////////////////////////////////////////////
548 /// Set XPsi.
549 
550 void TRotation::SetXPsi(Double_t psi) {
551  SetXEulerAngles(GetXPhi(),GetXTheta(),psi);
552 }
553 
554 ////////////////////////////////////////////////////////////////////////////////
555 /// Set YPhi.
556 
557 void TRotation::SetYPhi(Double_t phi) {
558  SetYEulerAngles(phi,GetYTheta(),GetYPsi());
559 }
560 
561 ////////////////////////////////////////////////////////////////////////////////
562 /// Set YTheta.
563 
564 void TRotation::SetYTheta(Double_t theta) {
565  SetYEulerAngles(GetYPhi(),theta,GetYPsi());
566 }
567 
568 ////////////////////////////////////////////////////////////////////////////////
569 /// Set YPsi.
570 
571 void TRotation::SetYPsi(Double_t psi) {
572  SetYEulerAngles(GetYPhi(),GetYTheta(),psi);
573 }
574 
575 ////////////////////////////////////////////////////////////////////////////////
576 /// Return phi angle.
577 
578 Double_t TRotation::GetXPhi(void) const {
579  Double_t finalPhi;
580 
581  Double_t s2 = 1.0 - fzz*fzz;
582  if (s2 < 0) {
583  Warning("GetPhi()"," |fzz| > 1 ");
584  s2 = 0;
585  }
586  const Double_t sinTheta = TMath::Sqrt(s2);
587 
588  if (sinTheta != 0) {
589  const Double_t cscTheta = 1/sinTheta;
590  Double_t cosAbsPhi = fzy * cscTheta;
591  if ( TMath::Abs(cosAbsPhi) > 1 ) { // NaN-proofing
592  Warning("GetPhi()","finds | cos phi | > 1");
593  cosAbsPhi = 1;
594  }
595  const Double_t absPhi = TMath::ACos(cosAbsPhi);
596  if (fzx > 0) {
597  finalPhi = absPhi;
598  } else if (fzx < 0) {
599  finalPhi = -absPhi;
600  } else if (fzy > 0) {
601  finalPhi = 0.0;
602  } else {
603  finalPhi = TMath::Pi();
604  }
605  } else { // sinTheta == 0 so |Fzz| = 1
606  const Double_t absPhi = .5 * TMath::ACos (fxx);
607  if (fxy > 0) {
608  finalPhi = -absPhi;
609  } else if (fxy < 0) {
610  finalPhi = absPhi;
611  } else if (fxx>0) {
612  finalPhi = 0.0;
613  } else {
614  finalPhi = fzz * TMath::PiOver2();
615  }
616  }
617  return finalPhi;
618 }
619 
620 ////////////////////////////////////////////////////////////////////////////////
621 /// Return YPhi.
622 
623 Double_t TRotation::GetYPhi(void) const {
624  return GetXPhi() + TMath::Pi()/2.0;
625 }
626 
627 ////////////////////////////////////////////////////////////////////////////////
628 /// Return XTheta.
629 
630 Double_t TRotation::GetXTheta(void) const {
631  return ThetaZ();
632 }
633 
634 ////////////////////////////////////////////////////////////////////////////////
635 /// Return YTheta.
636 
637 Double_t TRotation::GetYTheta(void) const {
638  return ThetaZ();
639 }
640 
641 ////////////////////////////////////////////////////////////////////////////////
642 /// Get psi angle.
643 
644 Double_t TRotation::GetXPsi(void) const {
645  double finalPsi = 0.0;
646 
647  Double_t s2 = 1.0 - fzz*fzz;
648  if (s2 < 0) {
649  Warning("GetPsi()"," |fzz| > 1 ");
650  s2 = 0;
651  }
652  const Double_t sinTheta = TMath::Sqrt(s2);
653 
654  if (sinTheta != 0) {
655  const Double_t cscTheta = 1/sinTheta;
656  Double_t cosAbsPsi = - fyz * cscTheta;
657  if ( TMath::Abs(cosAbsPsi) > 1 ) { // NaN-proofing
658  Warning("GetPsi()","| cos psi | > 1 ");
659  cosAbsPsi = 1;
660  }
661  const Double_t absPsi = TMath::ACos(cosAbsPsi);
662  if (fxz > 0) {
663  finalPsi = absPsi;
664  } else if (fxz < 0) {
665  finalPsi = -absPsi;
666  } else {
667  finalPsi = (fyz < 0) ? 0 : TMath::Pi();
668  }
669  } else { // sinTheta == 0 so |Fzz| = 1
670  Double_t absPsi = fxx;
671  if ( TMath::Abs(fxx) > 1 ) { // NaN-proofing
672  Warning("GetPsi()","| fxx | > 1 ");
673  absPsi = 1;
674  }
675  absPsi = .5 * TMath::ACos (absPsi);
676  if (fyx > 0) {
677  finalPsi = absPsi;
678  } else if (fyx < 0) {
679  finalPsi = -absPsi;
680  } else {
681  finalPsi = (fxx > 0) ? 0 : TMath::PiOver2();
682  }
683  }
684  return finalPsi;
685 }
686 
687 ////////////////////////////////////////////////////////////////////////////////
688 /// Return YPsi.
689 
690 Double_t TRotation::GetYPsi(void) const {
691  return GetXPsi() - TMath::Pi()/2;
692 }
693 
694 ////////////////////////////////////////////////////////////////////////////////
695 /// Set X axis.
696 
697 TRotation & TRotation::SetXAxis(const TVector3& axis,
698  const TVector3& xyPlane) {
699  TVector3 xAxis(xyPlane);
700  TVector3 yAxis;
701  TVector3 zAxis(axis);
702  MakeBasis(xAxis,yAxis,zAxis);
703  fxx = zAxis.X(); fyx = zAxis.Y(); fzx = zAxis.Z();
704  fxy = xAxis.X(); fyy = xAxis.Y(); fzy = xAxis.Z();
705  fxz = yAxis.X(); fyz = yAxis.Y(); fzz = yAxis.Z();
706  return *this;
707 }
708 
709 ////////////////////////////////////////////////////////////////////////////////
710 /// Set X axis.
711 
712 TRotation & TRotation::SetXAxis(const TVector3& axis) {
713  TVector3 xyPlane(0.0,1.0,0.0);
714  return SetXAxis(axis,xyPlane);
715 }
716 
717 ////////////////////////////////////////////////////////////////////////////////
718 /// Set Y axis.
719 
720 TRotation & TRotation::SetYAxis(const TVector3& axis,
721  const TVector3& yzPlane) {
722  TVector3 xAxis(yzPlane);
723  TVector3 yAxis;
724  TVector3 zAxis(axis);
725  MakeBasis(xAxis,yAxis,zAxis);
726  fxx = yAxis.X(); fyx = yAxis.Y(); fzx = yAxis.Z();
727  fxy = zAxis.X(); fyy = zAxis.Y(); fzy = zAxis.Z();
728  fxz = xAxis.X(); fyz = xAxis.Y(); fzz = xAxis.Z();
729  return *this;
730 }
731 
732 ////////////////////////////////////////////////////////////////////////////////
733 /// Set Y axis.
734 
735 TRotation & TRotation::SetYAxis(const TVector3& axis) {
736  TVector3 yzPlane(0.0,0.0,1.0);
737  return SetYAxis(axis,yzPlane);
738 }
739 
740 ////////////////////////////////////////////////////////////////////////////////
741 /// Set Z axis.
742 
743 TRotation & TRotation::SetZAxis(const TVector3& axis,
744  const TVector3& zxPlane) {
745  TVector3 xAxis(zxPlane);
746  TVector3 yAxis;
747  TVector3 zAxis(axis);
748  MakeBasis(xAxis,yAxis,zAxis);
749  fxx = xAxis.X(); fyx = xAxis.Y(); fzx = xAxis.Z();
750  fxy = yAxis.X(); fyy = yAxis.Y(); fzy = yAxis.Z();
751  fxz = zAxis.X(); fyz = zAxis.Y(); fzz = zAxis.Z();
752  return *this;
753 }
754 
755 ////////////////////////////////////////////////////////////////////////////////
756 /// Set Z axis.
757 
758 TRotation & TRotation::SetZAxis(const TVector3& axis) {
759  TVector3 zxPlane(1.0,0.0,0.0);
760  return SetZAxis(axis,zxPlane);
761 }
762 
763 ////////////////////////////////////////////////////////////////////////////////
764 /// Make the Z axis into a unit variable.
765 
766 void TRotation::MakeBasis(TVector3& xAxis,
767  TVector3& yAxis,
768  TVector3& zAxis) const {
769  Double_t zmag = zAxis.Mag();
770  if (zmag<TOLERANCE) {
771  Warning("MakeBasis(X,Y,Z)","non-zero Z Axis is required");
772  }
773  zAxis *= (1.0/zmag);
774 
775  Double_t xmag = xAxis.Mag();
776  if (xmag<TOLERANCE*zmag) {
777  xAxis = zAxis.Orthogonal();
778  xmag = 1.0;
779  }
780 
781  // Find the Y axis
782  yAxis = zAxis.Cross(xAxis)*(1.0/xmag);
783  Double_t ymag = yAxis.Mag();
784  if (ymag<TOLERANCE*zmag) {
785  yAxis = zAxis.Orthogonal();
786  } else {
787  yAxis *= (1.0/ymag);
788  }
789 
790  xAxis = yAxis.Cross(zAxis);
791 }