Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
REveTrans.cxx
Go to the documentation of this file.
1 // @(#)root/eve7:$Id$
2 // Authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007, 2018
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2019, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include <ROOT/REveTrans.hxx>
13 #include <ROOT/REveTypes.hxx>
14 
15 #include "TBuffer.h"
16 #include "TClass.h"
17 #include "TMath.h"
18 
19 #include "Riostream.h"
20 
21 #include <cctype>
22 
23 #define F00 0
24 #define F01 4
25 #define F02 8
26 #define F03 12
27 
28 #define F10 1
29 #define F11 5
30 #define F12 9
31 #define F13 13
32 
33 #define F20 2
34 #define F21 6
35 #define F22 10
36 #define F23 14
37 
38 #define F30 3
39 #define F31 7
40 #define F32 11
41 #define F33 15
42 
43 using namespace ROOT::Experimental;
44 namespace REX = ROOT::Experimental;
45 
46 /** \class REveTrans
47 \ingroup REve
48 REveTrans is a 4x4 transformation matrix for homogeneous coordinates
49 stored internally in a column-major order to allow direct usage by
50 GL. The element type is Double32_t as statically the floats would
51 be precise enough but continuous operations on the matrix must
52 retain precision of column vectors.
53 
54 Cartan angles are stored in fA[1-3] (+z, -y, +x). They are
55 recalculated on demand.
56 
57 Direct element access (first two should be used with care):
58  - operator[i] direct access to elements, i:0->15
59  - CM(i,j) element 4*j + i; i,j:0->3 { CM ~ c-matrix }
60  - operator(i,j) element 4*(j-1) + i - 1 i,j:1->4
61 
62 Column-vector access:
63 USet Get/SetBaseVec(), Get/SetPos() and Arr[XYZT]() methods.
64 
65 For all methods taking the matrix indices:
66 1->X, 2->Y, 3->Z; 4->Position (if applicable). 0 reserved for time.
67 
68 Shorthands in method-names:
69 LF ~ LocalFrame; PF ~ ParentFrame; IP ~ InPlace
70 */
71 
72 ////////////////////////////////////////////////////////////////////////////////
73 /// Default constructor.
74 
75 REveTrans::REveTrans() :
76  TObject(),
77  fA1(0), fA2(0), fA3(0), fAsOK(kFALSE),
78  fUseTrans (kTRUE),
79  fEditTrans(kFALSE),
80  fEditRotation(kTRUE),
81  fEditScale(kTRUE)
82 {
83  UnitTrans();
84 }
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// Constructor.
88 
89 REveTrans::REveTrans(const REveTrans& t) :
90  TObject(),
91  fA1(t.fA1), fA2(t.fA2), fA3(t.fA3), fAsOK(t.fAsOK),
92  fUseTrans (t.fUseTrans),
93  fEditTrans(t.fEditTrans),
94  fEditRotation(kTRUE),
95  fEditScale(kTRUE)
96 {
97  SetTrans(t, kFALSE);
98 }
99 
100 ////////////////////////////////////////////////////////////////////////////////
101 /// Constructor.
102 
103 REveTrans::REveTrans(const Double_t arr[16]) :
104  TObject(),
105  fA1(0), fA2(0), fA3(0), fAsOK(kFALSE),
106  fUseTrans (kTRUE),
107  fEditTrans(kFALSE),
108  fEditRotation(kTRUE),
109  fEditScale(kTRUE)
110 {
111  SetFromArray(arr);
112 }
113 
114 ////////////////////////////////////////////////////////////////////////////////
115 /// Constructor.
116 
117 REveTrans::REveTrans(const Float_t arr[16]) :
118  TObject(),
119  fA1(0), fA2(0), fA3(0), fAsOK(kFALSE),
120  fUseTrans (kTRUE),
121  fEditTrans(kFALSE),
122  fEditRotation(kTRUE),
123  fEditScale(kTRUE)
124 {
125  SetFromArray(arr);
126 }
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 /// Reset matrix to unity.
130 
131 void REveTrans::UnitTrans()
132 {
133  memset(fM, 0, 16*sizeof(Double_t));
134  fM[F00] = fM[F11] = fM[F22] = fM[F33] = 1;
135  fA1 = fA2 = fA3 = 0;
136  fAsOK = kTRUE;
137 }
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// Reset matrix to zero, only the perspective scaling is set to w
141 /// (1 by default).
142 
143 void REveTrans::ZeroTrans(Double_t w)
144 {
145  memset(fM, 0, 16*sizeof(Double_t));
146  fM[F33] = w;
147  fA1 = fA2 = fA3 = 0;
148  fAsOK = kFALSE;
149 }
150 
151 ////////////////////////////////////////////////////////////////////////////////
152 /// Reset rotation part of the matrix to unity.
153 
154 void REveTrans::UnitRot()
155 {
156  memset(fM, 0, 12*sizeof(Double_t));
157  fM[F00] = fM[F11] = fM[F22] = 1;
158  fA1 = fA2 = fA3 = 0;
159  fAsOK = kTRUE;
160 }
161 
162 ////////////////////////////////////////////////////////////////////////////////
163 /// Set matrix from another,
164 
165 void REveTrans::SetTrans(const REveTrans& t, Bool_t copyAngles)
166 {
167  memcpy(fM, t.fM, sizeof(fM));
168  if (copyAngles && t.fAsOK) {
169  fAsOK = kTRUE;
170  fA1 = t.fA1; fA2 = t.fA2; fA3 = t.fA3;
171  } else {
172  fAsOK = kFALSE;
173  }
174 }
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 /// Set matrix from Double_t array.
178 
179 void REveTrans::SetFromArray(const Double_t arr[16])
180 {
181  for(Int_t i=0; i<16; ++i) fM[i] = arr[i];
182  fAsOK = kFALSE;
183 }
184 
185 ////////////////////////////////////////////////////////////////////////////////
186 /// Set matrix from Float_t array.
187 
188 void REveTrans::SetFromArray(const Float_t arr[16])
189 {
190  for(Int_t i=0; i<16; ++i) fM[i] = arr[i];
191  fAsOK = kFALSE;
192 }
193 
194 ////////////////////////////////////////////////////////////////////////////////
195 /// Setup the matrix as an elementary rotation.
196 /// Optimized versions of left/right multiplication with an elementary
197 /// rotation matrix are implemented in RotatePF/RotateLF.
198 /// Expects identity matrix.
199 
200 void REveTrans::SetupRotation(Int_t i, Int_t j, Double_t f)
201 {
202  if(i == j) return;
203  REveTrans& t = *this;
204  t(i,i) = t(j,j) = TMath::Cos(f);
205  Double_t s = TMath::Sin(f);
206  t(i,j) = -s; t(j,i) = s;
207  fAsOK = kFALSE;
208 }
209 
210 ////////////////////////////////////////////////////////////////////////////////
211 /// A function for creating a rotation matrix that rotates a vector called
212 /// "from" into another vector called "to".
213 /// Input : from[3], to[3] which both must be *normalized* non-zero vectors
214 /// Output: mtx[3][3] -- a 3x3 matrix in column-major form
215 ///
216 /// Authors: Tomas Möller, John Hughes
217 /// "Efficiently Building a Matrix to Rotate One Vector to Another"
218 /// Journal of Graphics Tools, 4(4):1-4, 1999
219 
220 void REveTrans::SetupFromToVec(const REveVector& from, const REveVector& to)
221 {
222  static const float kFromToEpsilon = 0.000001f;
223 
224  ZeroTrans();
225 
226  Float_t e, f;
227  e = from.Dot(to);
228  f = (e < 0.0f) ? -e : e;
229 
230  if (f > 1.0f - kFromToEpsilon) /* "from" and "to"-vector almost parallel */
231  {
232  REveVector u, v; /* temporary storage vectors */
233  REveVector x; /* vector most nearly orthogonal to "from" */
234  Float_t c1, c2, c3; /* coefficients for later use */
235 
236  x.fX = (from.fX > 0.0f) ? from.fX : -from.fX;
237  x.fY = (from.fY > 0.0f) ? from.fY : -from.fY;
238  x.fZ = (from.fZ > 0.0f) ? from.fZ : -from.fZ;
239 
240  if (x.fX < x.fY)
241  {
242  if (x.fX < x.fZ) {
243  x.fX = 1.0f; x.fY = x.fZ = 0.0f;
244  } else {
245  x.fZ = 1.0f; x.fX = x.fY = 0.0f;
246  }
247  }
248  else
249  {
250  if (x.fY < x.fZ) {
251  x.fY = 1.0f; x.fX = x.fZ = 0.0f;
252  } else {
253  x.fZ = 1.0f; x.fX = x.fY = 0.0f;
254  }
255  }
256 
257  u.Sub(x, from);
258  v.Sub(x, to);
259 
260  c1 = 2.0f / u.Mag2();
261  c2 = 2.0f / v.Mag2();
262  c3 = c1 * c2 * u.Dot(v);
263 
264  for (int i = 0; i < 3; i++) {
265  for (int j = 0; j < 3; j++) {
266  CM(i, j) = - c1 * u[i] * u[j]
267  - c2 * v[i] * v[j]
268  + c3 * v[i] * u[j];
269  }
270  CM(i, i) += 1.0;
271  }
272  }
273  else /* the most common case, unless "from"="to", or "from"=-"to" */
274  {
275  REveVector v = from.Cross(to);
276 
277  Float_t h, hvx, hvz, hvxy, hvxz, hvyz;
278  h = 1.0f/(1.0f + e);
279  hvx = h * v.fX;
280  hvz = h * v.fZ;
281  hvxy = hvx * v.fY;
282  hvxz = hvx * v.fZ;
283  hvyz = hvz * v.fY;
284 
285  CM(0, 0) = e + hvx * v.fX;
286  CM(0, 1) = hvxy - v.fZ;
287  CM(0, 2) = hvxz + v.fY;
288 
289  CM(1, 0) = hvxy + v.fZ;
290  CM(1, 1) = e + h * v.fY * v.fY;
291  CM(1, 2) = hvyz - v.fX;
292 
293  CM(2, 0) = hvxz - v.fY;
294  CM(2, 1) = hvyz + v.fX;
295  CM(2, 2) = e + hvz * v.fZ;
296  }
297 }
298 
299 ////////////////////////////////////////////////////////////////////////////////
300 /// Multiply from left: this = t * this.
301 
302 void REveTrans::MultLeft(const REveTrans& t)
303 {
304  Double_t buf[4];
305  Double_t* col = fM;
306  for(int c=0; c<4; ++c, col+=4) {
307  const Double_t* row = t.fM;
308  for(int r=0; r<4; ++r, ++row)
309  buf[r] = row[0]*col[0] + row[4]*col[1] + row[8]*col[2] + row[12]*col[3];
310  col[0] = buf[0]; col[1] = buf[1]; col[2] = buf[2]; col[3] = buf[3];
311  }
312  fAsOK = kFALSE;
313 }
314 
315 ////////////////////////////////////////////////////////////////////////////////
316 /// Multiply from right: this = this * t.
317 
318 void REveTrans::MultRight(const REveTrans& t)
319 {
320  Double_t buf[4];
321  Double_t* row = fM;
322  for(int r=0; r<4; ++r, ++row) {
323  const Double_t* col = t.fM;
324  for(int c=0; c<4; ++c, col+=4)
325  buf[c] = row[0]*col[0] + row[4]*col[1] + row[8]*col[2] + row[12]*col[3];
326  row[0] = buf[0]; row[4] = buf[1]; row[8] = buf[2]; row[12] = buf[3];
327  }
328  fAsOK = kFALSE;
329 }
330 
331 ////////////////////////////////////////////////////////////////////////////////
332 /// Copy, multiply from right and return product.
333 /// Avoid unless necessary.
334 
335 REveTrans REveTrans::operator*(const REveTrans& t)
336 {
337  REveTrans b(*this);
338  b.MultRight(t);
339  return b;
340 }
341 
342 ////////////////////////////////////////////////////////////////////////////////
343 /// Transpose 3x3 rotation sub-matrix.
344 
345 void REveTrans::TransposeRotationPart()
346 {
347  Double_t x;
348  x = fM[F01]; fM[F01] = fM[F10]; fM[F10] = x;
349  x = fM[F02]; fM[F02] = fM[F20]; fM[F20] = x;
350  x = fM[F12]; fM[F12] = fM[F21]; fM[F21] = x;
351  fAsOK = kFALSE;
352 }
353 
354 ////////////////////////////////////////////////////////////////////////////////
355 /// Move in local-frame along axis with index ai.
356 
357 void REveTrans::MoveLF(Int_t ai, Double_t amount)
358 {
359  const Double_t *col = fM + 4*--ai;
360  fM[F03] += amount*col[0]; fM[F13] += amount*col[1]; fM[F23] += amount*col[2];
361 }
362 
363 ////////////////////////////////////////////////////////////////////////////////
364 /// General move in local-frame.
365 
366 void REveTrans::Move3LF(Double_t x, Double_t y, Double_t z)
367 {
368  fM[F03] += x*fM[0] + y*fM[4] + z*fM[8];
369  fM[F13] += x*fM[1] + y*fM[5] + z*fM[9];
370  fM[F23] += x*fM[2] + y*fM[6] + z*fM[10];
371 }
372 
373 ////////////////////////////////////////////////////////////////////////////////
374 /// Rotate in local frame. Does optimised version of MultRight.
375 
376 void REveTrans::RotateLF(Int_t i1, Int_t i2, Double_t amount)
377 {
378  if(i1 == i2) return;
379  // Algorithm: REveTrans a; a.SetupRotation(i1, i2, amount); MultRight(a);
380  // Optimized version:
381  const Double_t cos = TMath::Cos(amount), sin = TMath::Sin(amount);
382  Double_t b1, b2;
383  Double_t* row = fM;
384  --i1 <<= 2; --i2 <<= 2; // column major
385  for (int r=0; r<4; ++r, ++row) {
386  b1 = cos*row[i1] + sin*row[i2];
387  b2 = cos*row[i2] - sin*row[i1];
388  row[i1] = b1; row[i2] = b2;
389  }
390  fAsOK = kFALSE;
391 }
392 
393 ////////////////////////////////////////////////////////////////////////////////
394 /// Move in parent-frame along axis index ai.
395 
396 void REveTrans::MovePF(Int_t ai, Double_t amount)
397 {
398  fM[F03 + --ai] += amount;
399 }
400 
401 ////////////////////////////////////////////////////////////////////////////////
402 /// General move in parent-frame.
403 
404 void REveTrans::Move3PF(Double_t x, Double_t y, Double_t z)
405 {
406  fM[F03] += x;
407  fM[F13] += y;
408  fM[F23] += z;
409 }
410 
411 ////////////////////////////////////////////////////////////////////////////////
412 /// Rotate in parent frame. Does optimised version of MultLeft.
413 
414 void REveTrans::RotatePF(Int_t i1, Int_t i2, Double_t amount)
415 {
416  if(i1 == i2) return;
417  // Algorithm: REveTrans a; a.SetupRotation(i1, i2, amount); MultLeft(a);
418 
419  // Optimized version:
420  const Double_t cos = TMath::Cos(amount), sin = TMath::Sin(amount);
421  Double_t b1, b2;
422  Double_t* col = fM;
423  --i1; --i2;
424  for(int c=0; c<4; ++c, col+=4) {
425  b1 = cos*col[i1] - sin*col[i2];
426  b2 = cos*col[i2] + sin*col[i1];
427  col[i1] = b1; col[i2] = b2;
428  }
429  fAsOK = kFALSE;
430 }
431 
432 ////////////////////////////////////////////////////////////////////////////////
433 /// Move in a's coord-system along axis-index ai.
434 
435 void REveTrans::Move(const REveTrans& a, Int_t ai, Double_t amount)
436 {
437  const Double_t* vec = a.fM + 4*--ai;
438  fM[F03] += amount*vec[0];
439  fM[F13] += amount*vec[1];
440  fM[F23] += amount*vec[2];
441 }
442 
443 ////////////////////////////////////////////////////////////////////////////////
444 /// General move in a's coord-system.
445 
446 void REveTrans::Move3(const REveTrans& a, Double_t x, Double_t y, Double_t z)
447 {
448  const Double_t* m = a.fM;
449  fM[F03] += x*m[F00] + y*m[F01] + z*m[F02];
450  fM[F13] += x*m[F10] + y*m[F11] + z*m[F12];
451  fM[F23] += x*m[F20] + y*m[F21] + z*m[F22];
452 }
453 
454 ////////////////////////////////////////////////////////////////////////////////
455 /// Rotate in a's coord-system, rotating base vector with index i1
456 /// into i2.
457 
458 void REveTrans::Rotate(const REveTrans& a, Int_t i1, Int_t i2, Double_t amount)
459 {
460  if(i1 == i2) return;
461  REveTrans x(a);
462  x.Invert();
463  MultLeft(x);
464  RotatePF(i1, i2, amount);
465  MultLeft(a);
466  fAsOK = kFALSE;
467 }
468 
469 ////////////////////////////////////////////////////////////////////////////////
470 /// Set base-vector with index b.
471 
472 void REveTrans::SetBaseVec(Int_t b, Double_t x, Double_t y, Double_t z)
473 {
474  Double_t* col = fM + 4*--b;
475  col[0] = x; col[1] = y; col[2] = z;
476  fAsOK = kFALSE;
477 }
478 
479 ////////////////////////////////////////////////////////////////////////////////
480 /// Set base-vector with index b.
481 
482 void REveTrans::SetBaseVec(Int_t b, const TVector3& v)
483 {
484  Double_t* col = fM + 4*--b;
485  v.GetXYZ(col);
486  fAsOK = kFALSE;
487 }
488 
489 ////////////////////////////////////////////////////////////////////////////////
490 /// Get base-vector with index b.
491 
492 TVector3 REveTrans::GetBaseVec(Int_t b) const
493 {
494  return TVector3(&fM[4*--b]);
495 }
496 
497 void REveTrans::GetBaseVec(Int_t b, TVector3& v) const
498 {
499  // Get base-vector with index b.
500 
501  const Double_t* col = fM + 4*--b;
502  v.SetXYZ(col[0], col[1], col[2]);
503 }
504 
505 ////////////////////////////////////////////////////////////////////////////////
506 /// Set position (base-vec 4).
507 
508 void REveTrans::SetPos(Double_t x, Double_t y, Double_t z)
509 {
510  fM[F03] = x; fM[F13] = y; fM[F23] = z;
511 }
512 
513 void REveTrans::SetPos(Double_t* x)
514 {
515  // Set position (base-vec 4).
516  fM[F03] = x[0]; fM[F13] = x[1]; fM[F23] = x[2];
517 }
518 
519 void REveTrans::SetPos(Float_t* x)
520 {
521  // Set position (base-vec 4).
522  fM[F03] = x[0]; fM[F13] = x[1]; fM[F23] = x[2];
523 }
524 
525 void REveTrans::SetPos(const REveTrans& t)
526 {
527  // Set position (base-vec 4).
528  const Double_t* m = t.fM;
529  fM[F03] = m[F03]; fM[F13] = m[F13]; fM[F23] = m[F23];
530 }
531 
532 ////////////////////////////////////////////////////////////////////////////////
533 /// Get position (base-vec 4).
534 
535 void REveTrans::GetPos(Double_t& x, Double_t& y, Double_t& z) const
536 {
537  x = fM[F03]; y = fM[F13]; z = fM[F23];
538 }
539 
540 void REveTrans::GetPos(Double_t* x) const
541 {
542  // Get position (base-vec 4).
543  x[0] = fM[F03]; x[1] = fM[F13]; x[2] = fM[F23];
544 }
545 
546 void REveTrans::GetPos(Float_t* x) const
547 {
548  // Get position (base-vec 4).
549  x[0] = fM[F03]; x[1] = fM[F13]; x[2] = fM[F23];
550 }
551 
552 void REveTrans::GetPos(TVector3& v) const
553 {
554  // Get position (base-vec 4).
555  v.SetXYZ(fM[F03], fM[F13], fM[F23]);
556 }
557 
558 TVector3 REveTrans::GetPos() const
559 {
560  // Get position (base-vec 4).
561  return TVector3(fM[F03], fM[F13], fM[F23]);
562 }
563 
564 namespace
565 {
566 inline void clamp_angle(Float_t& a)
567 {
568  while(a < -TMath::TwoPi()) a += TMath::TwoPi();
569  while(a > TMath::TwoPi()) a -= TMath::TwoPi();
570 }
571 }
572 
573 void REveTrans::SetRotByAngles(Float_t a1, Float_t a2, Float_t a3)
574 {
575  // Sets Rotation part as given by angles:
576  // a1 around z, -a2 around y, a3 around x.
577 
578  clamp_angle(a1); clamp_angle(a2); clamp_angle(a3);
579 
580  Double_t a, b, c, d, e, f;
581  a = TMath::Cos(a3); b = TMath::Sin(a3);
582  c = TMath::Cos(a2); d = TMath::Sin(a2); // should be -sin(a2) for positive direction
583  e = TMath::Cos(a1); f = TMath::Sin(a1);
584  Double_t ad = a*d, bd = b*d;
585 
586  fM[F00] = c*e; fM[F01] = -bd*e - a*f; fM[F02] = -ad*e + b*f;
587  fM[F10] = c*f; fM[F11] = -bd*f + a*e; fM[F12] = -ad*f - b*e;
588  fM[F20] = d; fM[F21] = b*c; fM[F22] = a*c;
589 
590  fA1 = a1; fA2 = a2; fA3 = a3;
591  fAsOK = kTRUE;
592 }
593 
594 ////////////////////////////////////////////////////////////////////////////////
595 /// Sets Rotation part as given by angles a1, a1, a3 and pattern pat.
596 /// Pattern consists of "XxYyZz" characters.
597 /// eg: x means rotate about x axis, X means rotate in negative direction
598 /// xYz -> R_x(a3) * R_y(-a2) * R_z(a1); (standard Gled representation)
599 /// Note that angles and pattern elements have inverted order!
600 ///
601 /// Implements Eulerian/Cardanian angles in a uniform way.
602 
603 void REveTrans::SetRotByAnyAngles(Float_t a1, Float_t a2, Float_t a3,
604  const char* pat)
605 {
606  int n = strspn(pat, "XxYyZz"); if(n > 3) n = 3;
607  // Build Trans ... assign ...
608  Float_t a[] = { a3, a2, a1 };
609  UnitRot();
610  for(int i=0; i<n; i++) {
611  if(isupper(pat[i])) a[i] = -a[i];
612  switch(pat[i]) {
613  case 'x': case 'X': RotateLF(2, 3, a[i]); break;
614  case 'y': case 'Y': RotateLF(3, 1, a[i]); break;
615  case 'z': case 'Z': RotateLF(1, 2, a[i]); break;
616  }
617  }
618  fAsOK = kFALSE;
619 }
620 
621 ////////////////////////////////////////////////////////////////////////////////
622 /// Get Cardan rotation angles (pattern xYz above).
623 
624 void REveTrans::GetRotAngles(Float_t* x) const
625 {
626  if(!fAsOK) {
627  Double_t sx, sy, sz;
628  GetScale(sx, sy, sz);
629  Double_t d = fM[F20]/sx;
630  if(d>1) d=1; else if(d<-1) d=-1; // Fix numerical errors
631  fA2 = TMath::ASin(d);
632  Double_t cos = TMath::Cos(fA2);
633  if(TMath::Abs(cos) > 8.7e-6) {
634  fA1 = TMath::ATan2(fM[F10], fM[F00]);
635  fA3 = TMath::ATan2(fM[F21]/sy, fM[F22]/sz);
636  } else {
637  fA1 = TMath::ATan2(fM[F10]/sx, fM[F11]/sy);
638  fA3 = 0;
639  }
640  fAsOK = kTRUE;
641  }
642  x[0] = fA1; x[1] = fA2; x[2] = fA3;
643 }
644 
645 ////////////////////////////////////////////////////////////////////////////////
646 /// Scale matrix. Translation part untouched.
647 
648 void REveTrans::Scale(Double_t sx, Double_t sy, Double_t sz)
649 {
650  fM[F00] *= sx; fM[F10] *= sx; fM[F20] *= sx;
651  fM[F01] *= sy; fM[F11] *= sy; fM[F21] *= sy;
652  fM[F02] *= sz; fM[F12] *= sz; fM[F22] *= sz;
653 }
654 
655 ////////////////////////////////////////////////////////////////////////////////
656 /// Remove scaling, make all base vectors of unit length.
657 
658 Double_t REveTrans::Unscale()
659 {
660  Double_t sx, sy, sz;
661  Unscale(sx, sy, sz);
662  return (sx + sy + sz)/3;
663 }
664 
665 ////////////////////////////////////////////////////////////////////////////////
666 /// Remove scaling, make all base vectors of unit length.
667 
668 void REveTrans::Unscale(Double_t& sx, Double_t& sy, Double_t& sz)
669 {
670  GetScale(sx, sy, sz);
671  fM[F00] /= sx; fM[F10] /= sx; fM[F20] /= sx;
672  fM[F01] /= sy; fM[F11] /= sy; fM[F21] /= sy;
673  fM[F02] /= sz; fM[F12] /= sz; fM[F22] /= sz;
674 }
675 
676 ////////////////////////////////////////////////////////////////////////////////
677 /// Deduce scales from sizes of base vectors.
678 
679 void REveTrans::GetScale(Double_t& sx, Double_t& sy, Double_t& sz) const
680 {
681  sx = TMath::Sqrt( fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20] );
682  sy = TMath::Sqrt( fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21] );
683  sz = TMath::Sqrt( fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22] );
684 }
685 
686 ////////////////////////////////////////////////////////////////////////////////
687 /// Set scaling.
688 
689 void REveTrans::SetScale(Double_t sx, Double_t sy, Double_t sz)
690 {
691  sx /= TMath::Sqrt( fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20] );
692  sy /= TMath::Sqrt( fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21] );
693  sz /= TMath::Sqrt( fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22] );
694 
695  fM[F00] *= sx; fM[F10] *= sx; fM[F20] *= sx;
696  fM[F01] *= sy; fM[F11] *= sy; fM[F21] *= sy;
697  fM[F02] *= sz; fM[F12] *= sz; fM[F22] *= sz;
698 }
699 
700 ////////////////////////////////////////////////////////////////////////////////
701 /// Change x scaling.
702 
703 void REveTrans::SetScaleX(Double_t sx)
704 {
705  sx /= TMath::Sqrt( fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20] );
706  fM[F00] *= sx; fM[F10] *= sx; fM[F20] *= sx;
707 }
708 
709 ////////////////////////////////////////////////////////////////////////////////
710 /// Change y scaling.
711 
712 void REveTrans::SetScaleY(Double_t sy)
713 {
714  sy /= TMath::Sqrt( fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21] );
715  fM[F01] *= sy; fM[F11] *= sy; fM[F21] *= sy;
716 }
717 
718 ////////////////////////////////////////////////////////////////////////////////
719 /// Change z scaling.
720 
721 void REveTrans::SetScaleZ(Double_t sz)
722 {
723  sz /= TMath::Sqrt( fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22] );
724  fM[F02] *= sz; fM[F12] *= sz; fM[F22] *= sz;
725 }
726 
727 ////////////////////////////////////////////////////////////////////////////////
728 /// Multiply vector in-place.
729 
730 void REveTrans::MultiplyIP(TVector3& v, Double_t w) const
731 {
732  v.SetXYZ(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z() + fM[F03]*w,
733  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z() + fM[F13]*w,
734  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z() + fM[F23]*w);
735 }
736 
737 ////////////////////////////////////////////////////////////////////////////////
738 /// Multiply vector in-place.
739 
740 void REveTrans::MultiplyIP(Double_t* v, Double_t w) const
741 {
742  Double_t r[3] = { v[0], v[1], v[2] };
743  v[0] = fM[F00]*r[0] + fM[F01]*r[1] + fM[F02]*r[2] + fM[F03]*w;
744  v[1] = fM[F10]*r[0] + fM[F11]*r[1] + fM[F12]*r[2] + fM[F13]*w;
745  v[2] = fM[F20]*r[0] + fM[F21]*r[1] + fM[F22]*r[2] + fM[F23]*w;
746 }
747 
748 ////////////////////////////////////////////////////////////////////////////////
749 /// Multiply vector in-place.
750 
751 void REveTrans::MultiplyIP(Float_t* v, Double_t w) const
752 {
753  Double_t r[3] = { v[0], v[1], v[2] };
754  v[0] = fM[F00]*r[0] + fM[F01]*r[1] + fM[F02]*r[2] + fM[F03]*w;
755  v[1] = fM[F10]*r[0] + fM[F11]*r[1] + fM[F12]*r[2] + fM[F13]*w;
756  v[2] = fM[F20]*r[0] + fM[F21]*r[1] + fM[F22]*r[2] + fM[F23]*w;
757 }
758 
759 ////////////////////////////////////////////////////////////////////////////////
760 /// Multiply vector and return it.
761 
762 TVector3 REveTrans::Multiply(const TVector3& v, Double_t w) const
763 {
764  return TVector3(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z() + fM[F03]*w,
765  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z() + fM[F13]*w,
766  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z() + fM[F23]*w);
767 }
768 
769 ////////////////////////////////////////////////////////////////////////////////
770 /// Multiply vector and fill output array vout.
771 
772 void REveTrans::Multiply(const Double_t *vin, Double_t* vout, Double_t w) const
773 {
774  vout[0] = fM[F00]*vin[0] + fM[F01]*vin[1] + fM[F02]*vin[2] + fM[F03]*w;
775  vout[1] = fM[F10]*vin[0] + fM[F11]*vin[1] + fM[F12]*vin[1] + fM[F13]*w;
776  vout[2] = fM[F20]*vin[0] + fM[F21]*vin[1] + fM[F22]*vin[1] + fM[F23]*w;
777 }
778 
779 ////////////////////////////////////////////////////////////////////////////////
780 /// Rotate vector in-place. Translation is NOT applied.
781 
782 void REveTrans::RotateIP(TVector3& v) const
783 {
784  v.SetXYZ(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z(),
785  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z(),
786  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z());
787 }
788 
789 ////////////////////////////////////////////////////////////////////////////////
790 /// Rotate vector in-place. Translation is NOT applied.
791 
792 void REveTrans::RotateIP(Double_t* v) const
793 {
794  Double_t t[3] = { v[0], v[1], v[2] };
795 
796  v[0] = fM[F00]*t[0] + fM[F01]*t[1] + fM[F02]*t[2];
797  v[1] = fM[F10]*t[0] + fM[F11]*t[1] + fM[F12]*t[2];
798  v[2] = fM[F20]*t[0] + fM[F21]*t[1] + fM[F22]*t[2];
799 }
800 
801 ////////////////////////////////////////////////////////////////////////////////
802 /// Rotate vector in-place. Translation is NOT applied.
803 
804 void REveTrans::RotateIP(Float_t* v) const
805 {
806  Double_t t[3] = { v[0], v[1], v[2] };
807 
808  v[0] = fM[F00]*t[0] + fM[F01]*t[1] + fM[F02]*t[2];
809  v[1] = fM[F10]*t[0] + fM[F11]*t[1] + fM[F12]*t[2];
810  v[2] = fM[F20]*t[0] + fM[F21]*t[1] + fM[F22]*t[2];
811 }
812 
813 ////////////////////////////////////////////////////////////////////////////////
814 /// Rotate vector and return the rotated vector. Translation is NOT applied.
815 
816 TVector3 REveTrans::Rotate(const TVector3& v) const
817 {
818  return TVector3(fM[F00]*v.x() + fM[F01]*v.y() + fM[F02]*v.z(),
819  fM[F10]*v.x() + fM[F11]*v.y() + fM[F12]*v.z(),
820  fM[F20]*v.x() + fM[F21]*v.y() + fM[F22]*v.z());
821 }
822 
823 ////////////////////////////////////////////////////////////////////////////////
824 /// Norm 3-vector in column col.
825 
826 Double_t REveTrans::Norm3Column(Int_t col)
827 {
828  Double_t* c = fM + 4*--col;
829  const Double_t l = TMath::Sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
830  c[0] /= l; c[1] /= l; c[2] /= l;
831  return l;
832 }
833 
834 ////////////////////////////////////////////////////////////////////////////////
835 /// Orto-norm 3-vector in column col with respect to column ref.
836 
837 Double_t REveTrans::Orto3Column(Int_t col, Int_t ref)
838 {
839  Double_t* c = fM + 4*--col;
840  Double_t* rc = fM + 4*--ref;
841  const Double_t dp = c[0]*rc[0] + c[1]*rc[1] + c[2]*rc[2];
842  c[0] -= rc[0]*dp; c[1] -= rc[1]*dp; c[2] -= rc[2]*dp;
843  return dp;
844 }
845 
846 ////////////////////////////////////////////////////////////////////////////////
847 /// Orto-norm columns 1 to 3.
848 
849 void REveTrans::OrtoNorm3()
850 {
851  Norm3Column(1);
852  Orto3Column(2,1); Norm3Column(2);
853  fM[F02] = fM[F10]*fM[F21] - fM[F11]*fM[F20];
854  fM[F12] = fM[F20]*fM[F01] - fM[F21]*fM[F00];
855  fM[F22] = fM[F00]*fM[F11] - fM[F01]*fM[F10];
856  // Cross-product faster than the following.
857  // Orto3Column(3,1); Orto3Column(3,2); Norm3Column(3);
858 }
859 
860 ////////////////////////////////////////////////////////////////////////////////
861 /// Invert matrix.
862 /// Copied from ROOT's TMatrixFCramerInv.
863 
864 Double_t REveTrans::Invert()
865 {
866  static const REveException eh("REveTrans::Invert ");
867 
868  // Find all NECESSARY 2x2 dets: (18 of them)
869  const Double_t det2_12_01 = fM[F10]*fM[F21] - fM[F11]*fM[F20];
870  const Double_t det2_12_02 = fM[F10]*fM[F22] - fM[F12]*fM[F20];
871  const Double_t det2_12_03 = fM[F10]*fM[F23] - fM[F13]*fM[F20];
872  const Double_t det2_12_13 = fM[F11]*fM[F23] - fM[F13]*fM[F21];
873  const Double_t det2_12_23 = fM[F12]*fM[F23] - fM[F13]*fM[F22];
874  const Double_t det2_12_12 = fM[F11]*fM[F22] - fM[F12]*fM[F21];
875  const Double_t det2_13_01 = fM[F10]*fM[F31] - fM[F11]*fM[F30];
876  const Double_t det2_13_02 = fM[F10]*fM[F32] - fM[F12]*fM[F30];
877  const Double_t det2_13_03 = fM[F10]*fM[F33] - fM[F13]*fM[F30];
878  const Double_t det2_13_12 = fM[F11]*fM[F32] - fM[F12]*fM[F31];
879  const Double_t det2_13_13 = fM[F11]*fM[F33] - fM[F13]*fM[F31];
880  const Double_t det2_13_23 = fM[F12]*fM[F33] - fM[F13]*fM[F32];
881  const Double_t det2_23_01 = fM[F20]*fM[F31] - fM[F21]*fM[F30];
882  const Double_t det2_23_02 = fM[F20]*fM[F32] - fM[F22]*fM[F30];
883  const Double_t det2_23_03 = fM[F20]*fM[F33] - fM[F23]*fM[F30];
884  const Double_t det2_23_12 = fM[F21]*fM[F32] - fM[F22]*fM[F31];
885  const Double_t det2_23_13 = fM[F21]*fM[F33] - fM[F23]*fM[F31];
886  const Double_t det2_23_23 = fM[F22]*fM[F33] - fM[F23]*fM[F32];
887 
888  // Find all NECESSARY 3x3 dets: (16 of them)
889  const Double_t det3_012_012 = fM[F00]*det2_12_12 - fM[F01]*det2_12_02 + fM[F02]*det2_12_01;
890  const Double_t det3_012_013 = fM[F00]*det2_12_13 - fM[F01]*det2_12_03 + fM[F03]*det2_12_01;
891  const Double_t det3_012_023 = fM[F00]*det2_12_23 - fM[F02]*det2_12_03 + fM[F03]*det2_12_02;
892  const Double_t det3_012_123 = fM[F01]*det2_12_23 - fM[F02]*det2_12_13 + fM[F03]*det2_12_12;
893  const Double_t det3_013_012 = fM[F00]*det2_13_12 - fM[F01]*det2_13_02 + fM[F02]*det2_13_01;
894  const Double_t det3_013_013 = fM[F00]*det2_13_13 - fM[F01]*det2_13_03 + fM[F03]*det2_13_01;
895  const Double_t det3_013_023 = fM[F00]*det2_13_23 - fM[F02]*det2_13_03 + fM[F03]*det2_13_02;
896  const Double_t det3_013_123 = fM[F01]*det2_13_23 - fM[F02]*det2_13_13 + fM[F03]*det2_13_12;
897  const Double_t det3_023_012 = fM[F00]*det2_23_12 - fM[F01]*det2_23_02 + fM[F02]*det2_23_01;
898  const Double_t det3_023_013 = fM[F00]*det2_23_13 - fM[F01]*det2_23_03 + fM[F03]*det2_23_01;
899  const Double_t det3_023_023 = fM[F00]*det2_23_23 - fM[F02]*det2_23_03 + fM[F03]*det2_23_02;
900  const Double_t det3_023_123 = fM[F01]*det2_23_23 - fM[F02]*det2_23_13 + fM[F03]*det2_23_12;
901  const Double_t det3_123_012 = fM[F10]*det2_23_12 - fM[F11]*det2_23_02 + fM[F12]*det2_23_01;
902  const Double_t det3_123_013 = fM[F10]*det2_23_13 - fM[F11]*det2_23_03 + fM[F13]*det2_23_01;
903  const Double_t det3_123_023 = fM[F10]*det2_23_23 - fM[F12]*det2_23_03 + fM[F13]*det2_23_02;
904  const Double_t det3_123_123 = fM[F11]*det2_23_23 - fM[F12]*det2_23_13 + fM[F13]*det2_23_12;
905 
906  // Find the 4x4 det:
907  const Double_t det = fM[F00]*det3_123_123 - fM[F01]*det3_123_023 +
908  fM[F02]*det3_123_013 - fM[F03]*det3_123_012;
909 
910  if(det == 0) {
911  throw(eh + "matrix is singular.");
912  }
913 
914  const Double_t oneOverDet = 1.0/det;
915  const Double_t mn1OverDet = - oneOverDet;
916 
917  fM[F00] = det3_123_123 * oneOverDet;
918  fM[F01] = det3_023_123 * mn1OverDet;
919  fM[F02] = det3_013_123 * oneOverDet;
920  fM[F03] = det3_012_123 * mn1OverDet;
921 
922  fM[F10] = det3_123_023 * mn1OverDet;
923  fM[F11] = det3_023_023 * oneOverDet;
924  fM[F12] = det3_013_023 * mn1OverDet;
925  fM[F13] = det3_012_023 * oneOverDet;
926 
927  fM[F20] = det3_123_013 * oneOverDet;
928  fM[F21] = det3_023_013 * mn1OverDet;
929  fM[F22] = det3_013_013 * oneOverDet;
930  fM[F23] = det3_012_013 * mn1OverDet;
931 
932  fM[F30] = det3_123_012 * mn1OverDet;
933  fM[F31] = det3_023_012 * oneOverDet;
934  fM[F32] = det3_013_012 * mn1OverDet;
935  fM[F33] = det3_012_012 * oneOverDet;
936 
937  fAsOK = kFALSE;
938  return det;
939 }
940 
941 ////////////////////////////////////////////////////////////////////////////////
942 /// Stream an object of class REveTrans.
943 
944 void REveTrans::Streamer(TBuffer &R__b)
945 {
946  if (R__b.IsReading()) {
947  REveTrans::Class()->ReadBuffer(R__b, this);
948  fAsOK = kFALSE;
949  } else {
950  REveTrans::Class()->WriteBuffer(R__b, this);
951  }
952 }
953 
954 ////////////////////////////////////////////////////////////////////////////////
955 /// Print in reasonable format.
956 
957 void REveTrans::Print(Option_t* /*option*/) const
958 {
959  const Double_t* row = fM;
960  for(Int_t i=0; i<4; ++i, ++row)
961  printf("%8.3f %8.3f %8.3f | %8.3f\n", row[0], row[4], row[8], row[12]);
962 }
963 
964 #include <iomanip>
965 
966 ////////////////////////////////////////////////////////////////////////////////
967 /// Print to std::ostream.
968 
969 std::ostream& operator<<(std::ostream& s, const REveTrans& t)
970 {
971  s.setf(std::ios::fixed, std::ios::floatfield);
972  s.precision(3);
973  for(Int_t i=1; i<=4; i++)
974  for(Int_t j=1; j<=4; j++)
975  s << t(i,j) << ((j==4) ? "\n" : "\t");
976  return s;
977 }
978 
979 #include "TGeoMatrix.h"
980 #include "TBuffer3D.h"
981 
982 void REveTrans::SetFrom(Double_t* carr)
983 {
984  // Initialize from array.
985 
986  fUseTrans = kTRUE;
987  memcpy(fM, carr, 16*sizeof(Double_t));
988  fAsOK = kFALSE;
989 }
990 
991 ////////////////////////////////////////////////////////////////////////////////
992 /// Initialize from TGeoMatrix.
993 
994 void REveTrans::SetFrom(const TGeoMatrix& mat)
995 {
996  fUseTrans = kTRUE;
997  const Double_t *r = mat.GetRotationMatrix();
998  const Double_t *t = mat.GetTranslation();
999  Double_t *m = fM;
1000  if (mat.IsScale())
1001  {
1002  const Double_t *s = mat.GetScale();
1003  m[0] = r[0]*s[0]; m[1] = r[3]*s[0]; m[2] = r[6]*s[0]; m[3] = 0;
1004  m[4] = r[1]*s[1]; m[5] = r[4]*s[1]; m[6] = r[7]*s[1]; m[7] = 0;
1005  m[8] = r[2]*s[2]; m[9] = r[5]*s[2]; m[10] = r[8]*s[2]; m[11] = 0;
1006  m[12] = t[0]; m[13] = t[1]; m[14] = t[2]; m[15] = 1;
1007  }
1008  else
1009  {
1010  m[0] = r[0]; m[1] = r[3]; m[2] = r[6]; m[3] = 0;
1011  m[4] = r[1]; m[5] = r[4]; m[6] = r[7]; m[7] = 0;
1012  m[8] = r[2]; m[9] = r[5]; m[10] = r[8]; m[11] = 0;
1013  m[12] = t[0]; m[13] = t[1]; m[14] = t[2]; m[15] = 1;
1014  }
1015  fAsOK = kFALSE;
1016 }
1017 
1018 ////////////////////////////////////////////////////////////////////////////////
1019 /// Set TGeoHMatrix mat.
1020 
1021 void REveTrans::SetGeoHMatrix(TGeoHMatrix& mat)
1022 {
1023  Double_t *r = mat.GetRotationMatrix();
1024  Double_t *t = mat.GetTranslation();
1025  Double_t *s = mat.GetScale();
1026  if (fUseTrans)
1027  {
1028  mat.SetBit(TGeoMatrix::kGeoGenTrans);
1029  Double_t *m = fM;
1030  GetScale(s[0], s[1], s[2]);
1031  r[0] = m[0]/s[0]; r[3] = m[1]/s[0]; r[6] = m[2]/s[0]; m += 4;
1032  r[1] = m[0]/s[1]; r[4] = m[1]/s[1]; r[7] = m[2]/s[1]; m += 4;
1033  r[2] = m[0]/s[2]; r[5] = m[1]/s[2]; r[8] = m[2]/s[2]; m += 4;
1034  t[0] = m[0]; t[1] = m[1]; t[2] = m[2];
1035  }
1036  else
1037  {
1038  mat.ResetBit(TGeoMatrix::kGeoGenTrans);
1039  r[0] = 1; r[3] = 0; r[6] = 0;
1040  r[1] = 0; r[4] = 1; r[7] = 0;
1041  r[2] = 0; r[5] = 0; r[8] = 1;
1042  s[0] = s[1] = s[2] = 1;
1043  t[0] = t[1] = t[2] = 0;
1044  }
1045 }
1046 
1047 ////////////////////////////////////////////////////////////////////////////////
1048 /// Fill transformation part TBuffer3D core section.
1049 
1050 void REveTrans::SetBuffer3D(TBuffer3D& buff)
1051 {
1052  buff.fLocalFrame = fUseTrans;
1053  if (fUseTrans) {
1054  // In phys-shape ctor the rotation part is transposed, due to
1055  // TGeo's convention for rotation matrix. So we have to transpose
1056  // it here, also.
1057  Double_t *m = buff.fLocalMaster;
1058  m[0] = fM[0]; m[1] = fM[4]; m[2] = fM[8]; m[3] = fM[3];
1059  m[4] = fM[1]; m[5] = fM[5]; m[6] = fM[9]; m[7] = fM[7];
1060  m[8] = fM[2]; m[9] = fM[6]; m[10] = fM[10]; m[11] = fM[11];
1061  m[12] = fM[12]; m[13] = fM[13]; m[14] = fM[14]; m[15] = fM[15];
1062  // Otherwise this would do:
1063  // memcpy(buff.fLocalMaster, fM, 16*sizeof(Double_t));
1064  }
1065 }
1066 
1067 ////////////////////////////////////////////////////////////////////////////////
1068 /// Test if the transformation is a scale.
1069 /// To be used by ROOT TGLObject descendants that potentially need to
1070 /// use GL_NORMALIZE.
1071 /// The low/high limits are expected to be squares of actual limits.
1072 ///
1073 /// Ideally this should be done by the TGLViewer [but is not].
1074 
1075 Bool_t REveTrans::IsScale(Double_t low, Double_t high) const
1076 {
1077  if (!fUseTrans) return kFALSE;
1078  Double_t s;
1079  s = fM[F00]*fM[F00] + fM[F10]*fM[F10] + fM[F20]*fM[F20];
1080  if (s < low || s > high) return kTRUE;
1081  s = fM[F01]*fM[F01] + fM[F11]*fM[F11] + fM[F21]*fM[F21];
1082  if (s < low || s > high) return kTRUE;
1083  s = fM[F02]*fM[F02] + fM[F12]*fM[F12] + fM[F22]*fM[F22];
1084  if (s < low || s > high) return kTRUE;
1085  return kFALSE;
1086 }