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