Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoMatrix.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 25/10/01
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, 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 /** \class TGeoMatrix
13 \ingroup Geometry_classes
14 
15 Geometrical transformation package.
16 
17  All geometrical transformations handled by the modeller are provided as a
18 built-in package. This was designed to minimize memory requirements and
19 optimize performance of point/vector master-to-local and local-to-master
20 computation. We need to have in mind that a transformation in TGeo has 2
21 major use-cases. The first one is for defining the placement of a volume
22 with respect to its container reference frame. This frame will be called
23 'master' and the frame of the positioned volume - 'local'. If T is a
24 transformation used for positioning volume daughters, then:
25 
26 ~~~ {.cpp}
27  MASTER = T * LOCAL
28 ~~~
29 
30  Therefore a local-to-master conversion will be performed by using T, while
31 a master-to-local by using its inverse. The second use case is the computation
32 of the global transformation of a given object in the geometry. Since the
33 geometry is built as 'volumes-inside-volumes', this global transformation
34 represent the pile-up of all local transformations in the corresponding
35 branch. The conversion from the global reference frame and the given object
36 is also called master-to-local, but it is handled by the manager class.
37  A general homogenous transformation is defined as a 4x4 matrix embedding
38 a rotation, a translation and a scale. The advantage of this description
39 is that each basic transformation can be represented as a homogenous matrix,
40 composition being performed as simple matrix multiplication.
41 
42  Rotation: Inverse rotation:
43 
44 ~~~ {.cpp}
45  r11 r12 r13 0 r11 r21 r31 0
46  r21 r22 r23 0 r12 r22 r32 0
47  r31 r32 r33 0 r13 r23 r33 0
48  0 0 0 1 0 0 0 1
49 ~~~
50 
51  Translation: Inverse translation:
52 
53 ~~~ {.cpp}
54  1 0 0 tx 1 0 0 -tx
55  0 1 0 ty 0 1 0 -ty
56  0 0 1 tz 0 0 1 -tz
57  0 0 0 1 0 0 0 1
58 ~~~
59 
60  Scale: Inverse scale:
61 
62 ~~~ {.cpp}
63  sx 0 0 0 1/sx 0 0 0
64  0 sy 0 0 0 1/sy 0 0
65  0 0 sz 0 0 0 1/sz 0
66  0 0 0 1 0 0 0 1
67 ~~~
68 
69  where:
70  - `rij` are the 3x3 rotation matrix components,
71  - `tx`, `ty`, `tz` are the translation components
72  - `sx`, `sy`, `sz` are arbitrary scale constants on each axis,
73 
74  The disadvantage in using this approach is that computation for 4x4 matrices
75 is expensive. Even combining two translation would become a multiplication
76 of their corresponding matrices, which is quite an undesired effect. On the
77 other hand, it is not a good idea to store a translation as a block of 16
78 numbers. We have therefore chosen to implement each basic transformation type
79 as a class deriving from the same basic abstract class and handling its specific
80 data and point/vector transformation algorithms.
81 
82 \image html geom_transf.jpg
83 
84 ### The base class TGeoMatrix defines abstract metods for:
85 
86 #### translation, rotation and scale getters. Every derived class stores only
87  its specific data, e.g. a translation stores an array of 3 doubles and a
88  rotation an array of 9. However, asking which is the rotation array of a
89  TGeoTranslation through the base TGeoMatrix interface is a legal operation.
90  The answer in this case is a pointer to a global constant array representing
91  an identity rotation.
92 
93 ~~~ {.cpp}
94  Double_t *TGeoMatrix::GetTranslation()
95  Double_t *TGeoMatrix::GetRotation()
96  Double_t *TGeoMatrix::GetScale()
97 ~~~
98 
99 #### MasterToLocal() and LocalToMaster() point and vector transformations :
100 
101 ~~~ {.cpp}
102  void TGeoMatrix::MasterToLocal(const Double_t *master, Double_t *local)
103  void TGeoMatrix::LocalToMaster(const Double_t *local, Double_t *master)
104  void TGeoMatrix::MasterToLocalVect(const Double_t *master, Double_t *local)
105  void TGeoMatrix::LocalToMasterVect(const Double_t *local, Double_t *master)
106 ~~~
107 
108  These allow correct conversion also for reflections.
109 
110 #### Transformation type getters :
111 
112 ~~~ {.cpp}
113  Bool_t TGeoMatrix::IsIdentity()
114  Bool_t TGeoMatrix::IsTranslation()
115  Bool_t TGeoMatrix::IsRotation()
116  Bool_t TGeoMatrix::IsScale()
117  Bool_t TGeoMatrix::IsCombi() (translation + rotation)
118  Bool_t TGeoMatrix::IsGeneral() (translation + rotation + scale)
119 ~~~
120 
121  Combinations of basic transformations are represented by specific classes
122 deriving from TGeoMatrix. In order to define a matrix as a combination of several
123 others, a special class TGeoHMatrix is provided. Here is an example of matrix
124 creation :
125 
126 ### Matrix creation example:
127 
128 ~~~ {.cpp}
129  root[0] TGeoRotation r1,r2;
130  r1.SetAngles(90,0,30); // rotation defined by Euler angles
131  r2.SetAngles(90,90,90,180,0,0); // rotation defined by GEANT3 angles
132  TGeoTranslation t1(-10,10,0);
133  TGeoTranslation t2(10,-10,5);
134  TGeoCombiTrans c1(t1,r1);
135  TGeoCombiTrans c2(t2,r2);
136  TGeoHMatrix h = c1 * c2; // composition is done via TGeoHMatrix class
137  root[7] TGeoHMatrix *ph = new TGeoHMatrix(hm); // this is the one we want to
138  // use for positioning a volume
139  root[8] ph->Print();
140  ...
141  pVolume->AddNode(pVolDaughter,id,ph) // now ph is owned by the manager
142 ~~~
143 
144 ### Rule for matrix creation:
145  Unless explicitly used for positioning nodes (TGeoVolume::AddNode()) all
146 matrices deletion have to be managed by users. Matrices passed to geometry
147 have to be created by using new() operator and their deletion is done by
148 TGeoManager class.
149 
150 ### Available geometrical transformations
151 
152 #### TGeoTranslation
153 Represent a (dx,dy,dz) translation. Data members:
154  Double_t fTranslation[3]. Translations can be added/subtracted.
155 
156 ~~~ {.cpp}
157  TGeoTranslation t1;
158  t1->SetTranslation(-5,10,4);
159  TGeoTranslation *t2 = new TGeoTranslation(4,3,10);
160  t2->Subtract(&t1);
161 ~~~
162 
163 #### Rotations
164  Represent a pure rotation. Data members: Double_t fRotationMatrix[3*3].
165  Rotations can be defined either by Euler angles, either, by GEANT3 angles :
166 
167 ~~~ {.cpp}
168  TGeoRotation *r1 = new TGeoRotation();
169  r1->SetAngles(phi, theta, psi); // all angles in degrees
170 ~~~
171 
172  This represent the composition of : first a rotation about Z axis with
173  angle phi, then a rotation with theta about the rotated X axis, and
174  finally a rotation with psi about the new Z axis.
175 
176 ~~~ {.cpp}
177  r1->SetAngles(th1,phi1, th2,phi2, th3,phi3)
178 ~~~
179 
180  This is a rotation defined in GEANT3 style. Theta and phi are the spherical
181  angles of each axis of the rotated coordinate system with respect to the
182  initial one. This construction allows definition of malformed rotations,
183  e.g. not orthogonal. A check is performed and an error message is issued
184  in this case.
185 
186  Specific utilities : determinant, inverse.
187 
188 #### Scale transformations
189  Represent a scale shrinking/enlargement. Data
190  members :Double_t fScale[3]. Not fully implemented yet.
191 
192 #### Combined transformations
193 Represent a rotation followed by a translation.
194 Data members: Double_t fTranslation[3], TGeoRotation *fRotation.
195 
196 ~~~ {.cpp}
197  TGeoRotation *rot = new TGeoRotation("rot",10,20,30);
198  TGeoTranslation trans;
199  ...
200  TGeoCombiTrans *c1 = new TGeoCombiTrans(trans, rot);
201  TGeoCombiTrans *c2 = new TGeoCombiTrans("somename",10,20,30,rot)
202 ~~~
203 
204 
205 #### TGeoGenTrans
206 Combined transformations including a scale. Not implemented.
207 
208 #### TGeoIdentity
209 A generic singleton matrix representing a identity transformation
210  NOTE: identified by the global variable gGeoIdentity.
211 */
212 
213 #include "Riostream.h"
214 #include "TObjArray.h"
215 
216 #include "TGeoManager.h"
217 #include "TGeoMatrix.h"
218 #include "TMath.h"
219 
220 TGeoIdentity *gGeoIdentity = 0;
221 const Int_t kN3 = 3*sizeof(Double_t);
222 const Int_t kN9 = 9*sizeof(Double_t);
223 
224 // statics and globals
225 
226 ClassImp(TGeoMatrix);
227 
228 ////////////////////////////////////////////////////////////////////////////////
229 /// dummy constructor
230 
231 TGeoMatrix::TGeoMatrix()
232 {
233  ResetBit(kGeoMatrixBits);
234 }
235 
236 ////////////////////////////////////////////////////////////////////////////////
237 /// copy constructor
238 
239 TGeoMatrix::TGeoMatrix(const TGeoMatrix &other)
240  :TNamed(other)
241 {
242  ResetBit(kGeoRegistered);
243 }
244 
245 ////////////////////////////////////////////////////////////////////////////////
246 /// Constructor
247 
248 TGeoMatrix::TGeoMatrix(const char *name)
249  :TNamed(name, "")
250 {
251 }
252 
253 ////////////////////////////////////////////////////////////////////////////////
254 /// Destructor
255 
256 TGeoMatrix::~TGeoMatrix()
257 {
258  if (IsRegistered() && gGeoManager) {
259  if (!gGeoManager->IsCleaning()) {
260  gGeoManager->GetListOfMatrices()->Remove(this);
261  Warning("dtor", "Registered matrix %s was removed", GetName());
262  }
263  }
264 }
265 
266 ////////////////////////////////////////////////////////////////////////////////
267 /// Returns true if no rotation or the rotation is about Z axis
268 
269 Bool_t TGeoMatrix::IsRotAboutZ() const
270 {
271  if (IsIdentity()) return kTRUE;
272  const Double_t *rot = GetRotationMatrix();
273  if (TMath::Abs(rot[6])>1E-9) return kFALSE;
274  if (TMath::Abs(rot[7])>1E-9) return kFALSE;
275  if ((1.-TMath::Abs(rot[8]))>1E-9) return kFALSE;
276  return kTRUE;
277 }
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 /// Get total size in bytes of this
281 
282 Int_t TGeoMatrix::GetByteCount() const
283 {
284  Int_t count = 4+28+strlen(GetName())+strlen(GetTitle()); // fId + TNamed
285  if (IsTranslation()) count += 12;
286  if (IsScale()) count += 12;
287  if (IsCombi() || IsGeneral()) count += 4 + 36;
288  return count;
289 }
290 
291 ////////////////////////////////////////////////////////////////////////////////
292 /// Provide a pointer name containing uid.
293 
294 char *TGeoMatrix::GetPointerName() const
295 {
296  static TString name;
297  name = TString::Format("pMatrix%d", GetUniqueID());
298  return (char*)name.Data();
299 }
300 
301 ////////////////////////////////////////////////////////////////////////////////
302 /// The homogenous matrix associated with the transformation is used for
303 /// piling up's and visualization. A homogenous matrix is a 4*4 array
304 /// containing the translation, the rotation and the scale components
305 /// ~~~ {.cpp}
306 /// | R00*sx R01 R02 dx |
307 /// | R10 R11*sy R12 dy |
308 /// | R20 R21 R22*sz dz |
309 /// | 0 0 0 1 |
310 /// ~~~
311 /// where Rij is the rotation matrix, (sx, sy, sz) is the scale
312 /// transformation and (dx, dy, dz) is the translation.
313 
314 void TGeoMatrix::GetHomogenousMatrix(Double_t *hmat) const
315 {
316  Double_t *hmatrix = hmat;
317  const Double_t *mat = GetRotationMatrix();
318  for (Int_t i=0; i<3; i++) {
319  memcpy(hmatrix, mat, kN3);
320  mat += 3;
321  hmatrix += 3;
322  *hmatrix = 0.0;
323  hmatrix++;
324  }
325  memcpy(hmatrix, GetTranslation(), kN3);
326  hmatrix = hmat;
327  if (IsScale()) {
328  for (Int_t i=0; i<3; i++) {
329  *hmatrix *= GetScale()[i];
330  hmatrix += 5;
331  }
332  }
333  hmatrix[15] = 1.;
334 }
335 
336 ////////////////////////////////////////////////////////////////////////////////
337 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
338 
339 void TGeoMatrix::LocalToMaster(const Double_t *local, Double_t *master) const
340 {
341  if (IsIdentity()) {
342  memcpy(master, local, kN3);
343  return;
344  }
345  Int_t i;
346  const Double_t *tr = GetTranslation();
347  if (!IsRotation()) {
348  for (i=0; i<3; i++) master[i] = tr[i] + local[i];
349  return;
350  }
351  const Double_t *rot = GetRotationMatrix();
352  for (i=0; i<3; i++) {
353  master[i] = tr[i]
354  + local[0]*rot[3*i]
355  + local[1]*rot[3*i+1]
356  + local[2]*rot[3*i+2];
357  }
358 }
359 
360 ////////////////////////////////////////////////////////////////////////////////
361 /// convert a vector by multiplying its column vector (x, y, z, 1) to matrix inverse
362 
363 void TGeoMatrix::LocalToMasterVect(const Double_t *local, Double_t *master) const
364 {
365  if (!IsRotation()) {
366  memcpy(master, local, kN3);
367  return;
368  }
369  const Double_t *rot = GetRotationMatrix();
370  for (Int_t i=0; i<3; i++) {
371  master[i] = local[0]*rot[3*i]
372  + local[1]*rot[3*i+1]
373  + local[2]*rot[3*i+2];
374  }
375 }
376 
377 ////////////////////////////////////////////////////////////////////////////////
378 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
379 
380 void TGeoMatrix::LocalToMasterBomb(const Double_t *local, Double_t *master) const
381 {
382  if (IsIdentity()) {
383  memcpy(master, local, kN3);
384  return;
385  }
386  Int_t i;
387  const Double_t *tr = GetTranslation();
388  Double_t bombtr[3] = {0.,0.,0.};
389  gGeoManager->BombTranslation(tr, &bombtr[0]);
390  if (!IsRotation()) {
391  for (i=0; i<3; i++) master[i] = bombtr[i] + local[i];
392  return;
393  }
394  const Double_t *rot = GetRotationMatrix();
395  for (i=0; i<3; i++) {
396  master[i] = bombtr[i]
397  + local[0]*rot[3*i]
398  + local[1]*rot[3*i+1]
399  + local[2]*rot[3*i+2];
400  }
401 }
402 
403 ////////////////////////////////////////////////////////////////////////////////
404 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix
405 
406 void TGeoMatrix::MasterToLocal(const Double_t *master, Double_t *local) const
407 {
408  if (IsIdentity()) {
409  memcpy(local, master, kN3);
410  return;
411  }
412  const Double_t *tr = GetTranslation();
413  Double_t mt0 = master[0]-tr[0];
414  Double_t mt1 = master[1]-tr[1];
415  Double_t mt2 = master[2]-tr[2];
416  if (!IsRotation()) {
417  local[0] = mt0;
418  local[1] = mt1;
419  local[2] = mt2;
420  return;
421  }
422  const Double_t *rot = GetRotationMatrix();
423  local[0] = mt0*rot[0] + mt1*rot[3] + mt2*rot[6];
424  local[1] = mt0*rot[1] + mt1*rot[4] + mt2*rot[7];
425  local[2] = mt0*rot[2] + mt1*rot[5] + mt2*rot[8];
426 }
427 
428 ////////////////////////////////////////////////////////////////////////////////
429 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix
430 
431 void TGeoMatrix::MasterToLocalVect(const Double_t *master, Double_t *local) const
432 {
433  if (!IsRotation()) {
434  memcpy(local, master, kN3);
435  return;
436  }
437  const Double_t *rot = GetRotationMatrix();
438  for (Int_t i=0; i<3; i++) {
439  local[i] = master[0]*rot[i]
440  + master[1]*rot[i+3]
441  + master[2]*rot[i+6];
442  }
443 }
444 
445 ////////////////////////////////////////////////////////////////////////////////
446 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix
447 
448 void TGeoMatrix::MasterToLocalBomb(const Double_t *master, Double_t *local) const
449 {
450  if (IsIdentity()) {
451  memcpy(local, master, kN3);
452  return;
453  }
454  const Double_t *tr = GetTranslation();
455  Double_t bombtr[3] = {0.,0.,0.};
456  Int_t i;
457  gGeoManager->UnbombTranslation(tr, &bombtr[0]);
458  if (!IsRotation()) {
459  for (i=0; i<3; i++) local[i] = master[i]-bombtr[i];
460  return;
461  }
462  const Double_t *rot = GetRotationMatrix();
463  for (i=0; i<3; i++) {
464  local[i] = (master[0]-bombtr[0])*rot[i]
465  + (master[1]-bombtr[1])*rot[i+3]
466  + (master[2]-bombtr[2])*rot[i+6];
467  }
468 }
469 
470 ////////////////////////////////////////////////////////////////////////////////
471 /// Normalize a vector.
472 
473 void TGeoMatrix::Normalize(Double_t *vect)
474 {
475  Double_t normfactor = vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2];
476  if (normfactor <= 1E-10) return;
477  normfactor = 1./TMath::Sqrt(normfactor);
478  vect[0] *= normfactor;
479  vect[1] *= normfactor;
480  vect[2] *= normfactor;
481 }
482 
483 ////////////////////////////////////////////////////////////////////////////////
484 /// print the matrix in 4x4 format
485 
486 void TGeoMatrix::Print(Option_t *) const
487 {
488  const Double_t *rot = GetRotationMatrix();
489  const Double_t *tr = GetTranslation();
490  printf("matrix %s - tr=%d rot=%d refl=%d scl=%d shr=%d reg=%d own=%d\n", GetName(),(Int_t)IsTranslation(),
491  (Int_t)IsRotation(), (Int_t)IsReflection(), (Int_t)IsScale(), (Int_t)IsShared(), (Int_t)IsRegistered(),
492  (Int_t)IsOwned());
493  printf("%10.6f%12.6f%12.6f Tx = %10.6f\n", rot[0], rot[1], rot[2], tr[0]);
494  printf("%10.6f%12.6f%12.6f Ty = %10.6f\n", rot[3], rot[4], rot[5], tr[1]);
495  printf("%10.6f%12.6f%12.6f Tz = %10.6f\n", rot[6], rot[7], rot[8], tr[2]);
496  if (IsScale()) {
497  const Double_t *scl = GetScale();
498  printf("Sx=%10.6fSy=%12.6fSz=%12.6f\n", scl[0], scl[1], scl[2]);
499  }
500 }
501 
502 ////////////////////////////////////////////////////////////////////////////////
503 /// Multiply by a reflection respect to YZ.
504 
505 void TGeoMatrix::ReflectX(Bool_t, Bool_t)
506 {
507 }
508 
509 ////////////////////////////////////////////////////////////////////////////////
510 /// Multiply by a reflection respect to ZX.
511 
512 void TGeoMatrix::ReflectY(Bool_t, Bool_t)
513 {
514 }
515 
516 ////////////////////////////////////////////////////////////////////////////////
517 /// Multiply by a reflection respect to XY.
518 
519 void TGeoMatrix::ReflectZ(Bool_t, Bool_t)
520 {
521 }
522 
523 ////////////////////////////////////////////////////////////////////////////////
524 /// Register the matrix in the current manager, which will become the owner.
525 
526 void TGeoMatrix::RegisterYourself()
527 {
528  if (!gGeoManager) {
529  Warning("RegisterYourself", "cannot register without geometry");
530  return;
531  }
532  if (!IsRegistered()) {
533  gGeoManager->RegisterMatrix(this);
534  SetBit(kGeoRegistered);
535  }
536 }
537 
538 ////////////////////////////////////////////////////////////////////////////////
539 /// If no name was supplied in the ctor, the type of transformation is checked.
540 /// A letter will be prepended to the name :
541 /// - t - translation
542 /// - r - rotation
543 /// - s - scale
544 /// - c - combi (translation + rotation)
545 /// - g - general (tr+rot+scale)
546 /// The index of the transformation in gGeoManager list of transformations will
547 /// be appended.
548 
549 void TGeoMatrix::SetDefaultName()
550 {
551  if (!gGeoManager) return;
552  if (strlen(GetName())) return;
553  char type = 'n';
554  if (IsTranslation()) type = 't';
555  if (IsRotation()) type = 'r';
556  if (IsScale()) type = 's';
557  if (IsCombi()) type = 'c';
558  if (IsGeneral()) type = 'g';
559  TObjArray *matrices = gGeoManager->GetListOfMatrices();
560  Int_t index = 0;
561  if (matrices) index =matrices->GetEntriesFast() - 1;
562  TString name = TString::Format("%c%d", type, index);
563  SetName(name);
564 }
565 
566 /** \class TGeoTranslation
567 \ingroup Geometry_classes
568 
569 Class describing translations. A translation is
570 basically an array of 3 doubles matching the positions 12, 13
571 and 14 in the homogenous matrix description.
572 */
573 
574 ClassImp(TGeoTranslation);
575 
576 ////////////////////////////////////////////////////////////////////////////////
577 /// Default constructor
578 
579 TGeoTranslation::TGeoTranslation()
580 {
581  for (Int_t i=0; i<3; i++) fTranslation[i] = 0;
582 }
583 
584 ////////////////////////////////////////////////////////////////////////////////
585 /// Copy ctor.
586 
587 TGeoTranslation::TGeoTranslation(const TGeoTranslation &other)
588  :TGeoMatrix(other)
589 {
590  SetTranslation(other);
591 }
592 
593 ////////////////////////////////////////////////////////////////////////////////
594 /// Ctor. based on a general matrix
595 
596 TGeoTranslation::TGeoTranslation(const TGeoMatrix &other)
597  :TGeoMatrix(other)
598 {
599  ResetBit(kGeoRotation);
600  ResetBit(kGeoScale);
601  SetTranslation(other);
602 }
603 
604 ////////////////////////////////////////////////////////////////////////////////
605 /// Default constructor defining the translation
606 
607 TGeoTranslation::TGeoTranslation(Double_t dx, Double_t dy, Double_t dz)
608  :TGeoMatrix("")
609 {
610  if (dx || dy || dz) SetBit(kGeoTranslation);
611  SetTranslation(dx, dy, dz);
612 }
613 
614 ////////////////////////////////////////////////////////////////////////////////
615 /// Default constructor defining the translation
616 
617 TGeoTranslation::TGeoTranslation(const char *name, Double_t dx, Double_t dy, Double_t dz)
618  :TGeoMatrix(name)
619 {
620  if (dx || dy || dz) SetBit(kGeoTranslation);
621  SetTranslation(dx, dy, dz);
622 }
623 
624 ////////////////////////////////////////////////////////////////////////////////
625 /// Assignment from a general matrix
626 
627 TGeoTranslation& TGeoTranslation::operator = (const TGeoMatrix &matrix)
628 {
629  if (&matrix == this) return *this;
630  Bool_t registered = TestBit(kGeoRegistered);
631  TNamed::operator=(matrix);
632  SetTranslation(matrix);
633  SetBit(kGeoRegistered,registered);
634  ResetBit(kGeoRotation);
635  ResetBit(kGeoScale);
636  return *this;
637 }
638 
639 ////////////////////////////////////////////////////////////////////////////////
640 /// Translation composition
641 
642 TGeoTranslation &TGeoTranslation::operator*=(const TGeoTranslation &right)
643 {
644  const Double_t *tr = right.GetTranslation();
645  fTranslation[0] += tr[0];
646  fTranslation[1] += tr[1];
647  fTranslation[2] += tr[2];
648  if (!IsTranslation()) SetBit(kGeoTranslation, right.IsTranslation());
649  return *this;
650 }
651 
652 TGeoTranslation TGeoTranslation::operator*(const TGeoTranslation &right) const
653 {
654  TGeoTranslation t = *this;
655  t *= right;
656  return t;
657 }
658 
659 TGeoHMatrix TGeoTranslation::operator*(const TGeoMatrix &right) const
660 {
661  TGeoHMatrix t = *this;
662  t *= right;
663  return t;
664 }
665 
666 ////////////////////////////////////////////////////////////////////////////////
667 /// Is-equal operator
668 
669 Bool_t TGeoTranslation::operator ==(const TGeoTranslation &other) const
670 {
671  if (&other == this) return kTRUE;
672  const Double_t *tr = GetTranslation();
673  const Double_t *otr = other.GetTranslation();
674  for (auto i=0; i<3; i++)
675  if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
676  return kTRUE;
677 }
678 
679 ////////////////////////////////////////////////////////////////////////////////
680 /// Return a temporary inverse of this.
681 
682 TGeoHMatrix TGeoTranslation::Inverse() const
683 {
684  TGeoHMatrix h;
685  h = *this;
686  Double_t tr[3];
687  tr[0] = -fTranslation[0];
688  tr[1] = -fTranslation[1];
689  tr[2] = -fTranslation[2];
690  h.SetTranslation(tr);
691  return h;
692 }
693 
694 ////////////////////////////////////////////////////////////////////////////////
695 /// Adding a translation to this one
696 
697 void TGeoTranslation::Add(const TGeoTranslation *other)
698 {
699  const Double_t *trans = other->GetTranslation();
700  for (Int_t i=0; i<3; i++)
701  fTranslation[i] += trans[i];
702 }
703 
704 ////////////////////////////////////////////////////////////////////////////////
705 /// Make a clone of this matrix.
706 
707 TGeoMatrix *TGeoTranslation::MakeClone() const
708 {
709  TGeoMatrix *matrix = new TGeoTranslation(*this);
710  return matrix;
711 }
712 
713 ////////////////////////////////////////////////////////////////////////////////
714 /// Rotate about X axis of the master frame with angle expressed in degrees.
715 
716 void TGeoTranslation::RotateX(Double_t /*angle*/)
717 {
718  Warning("RotateX", "Not implemented. Use TGeoCombiTrans instead");
719 }
720 
721 ////////////////////////////////////////////////////////////////////////////////
722 /// Rotate about Y axis of the master frame with angle expressed in degrees.
723 
724 void TGeoTranslation::RotateY(Double_t /*angle*/)
725 {
726  Warning("RotateY", "Not implemented. Use TGeoCombiTrans instead");
727 }
728 
729 ////////////////////////////////////////////////////////////////////////////////
730 /// Rotate about Z axis of the master frame with angle expressed in degrees.
731 
732 void TGeoTranslation::RotateZ(Double_t /*angle*/)
733 {
734  Warning("RotateZ", "Not implemented. Use TGeoCombiTrans instead");
735 }
736 
737 ////////////////////////////////////////////////////////////////////////////////
738 /// Subtracting a translation from this one
739 
740 void TGeoTranslation::Subtract(const TGeoTranslation *other)
741 {
742  const Double_t *trans = other->GetTranslation();
743  for (Int_t i=0; i<3; i++)
744  fTranslation[i] -= trans[i];
745 }
746 
747 ////////////////////////////////////////////////////////////////////////////////
748 /// Set translation components
749 
750 void TGeoTranslation::SetTranslation(Double_t dx, Double_t dy, Double_t dz)
751 {
752  fTranslation[0] = dx;
753  fTranslation[1] = dy;
754  fTranslation[2] = dz;
755  if (dx || dy || dz) SetBit(kGeoTranslation);
756  else ResetBit(kGeoTranslation);
757 }
758 
759 ////////////////////////////////////////////////////////////////////////////////
760 /// Set translation components
761 
762 void TGeoTranslation::SetTranslation(const TGeoMatrix &other)
763 {
764  SetBit(kGeoTranslation, other.IsTranslation());
765  const Double_t *transl = other.GetTranslation();
766  memcpy(fTranslation, transl, kN3);
767 }
768 
769 ////////////////////////////////////////////////////////////////////////////////
770 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
771 
772 void TGeoTranslation::LocalToMaster(const Double_t *local, Double_t *master) const
773 {
774  const Double_t *tr = GetTranslation();
775  for (Int_t i=0; i<3; i++)
776  master[i] = tr[i] + local[i];
777 }
778 
779 ////////////////////////////////////////////////////////////////////////////////
780 /// convert a vector to MARS
781 
782 void TGeoTranslation::LocalToMasterVect(const Double_t *local, Double_t *master) const
783 {
784  memcpy(master, local, kN3);
785 }
786 
787 ////////////////////////////////////////////////////////////////////////////////
788 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
789 
790 void TGeoTranslation::LocalToMasterBomb(const Double_t *local, Double_t *master) const
791 {
792  const Double_t *tr = GetTranslation();
793  Double_t bombtr[3] = {0.,0.,0.};
794  gGeoManager->BombTranslation(tr, &bombtr[0]);
795  for (Int_t i=0; i<3; i++)
796  master[i] = bombtr[i] + local[i];
797 }
798 
799 ////////////////////////////////////////////////////////////////////////////////
800 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix
801 
802 void TGeoTranslation::MasterToLocal(const Double_t *master, Double_t *local) const
803 {
804  const Double_t *tr = GetTranslation();
805  for (Int_t i=0; i<3; i++)
806  local[i] = master[i]-tr[i];
807 }
808 
809 ////////////////////////////////////////////////////////////////////////////////
810 /// convert a vector from MARS to local
811 
812 void TGeoTranslation::MasterToLocalVect(const Double_t *master, Double_t *local) const
813 {
814  memcpy(local, master, kN3);
815 }
816 
817 ////////////////////////////////////////////////////////////////////////////////
818 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix
819 
820 void TGeoTranslation::MasterToLocalBomb(const Double_t *master, Double_t *local) const
821 {
822  const Double_t *tr = GetTranslation();
823  Double_t bombtr[3] = {0.,0.,0.};
824  gGeoManager->UnbombTranslation(tr, &bombtr[0]);
825  for (Int_t i=0; i<3; i++)
826  local[i] = master[i]-bombtr[i];
827 }
828 
829 ////////////////////////////////////////////////////////////////////////////////
830 /// Save a primitive as a C++ statement(s) on output stream "out".
831 
832 void TGeoTranslation::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
833 {
834  if (TestBit(kGeoSavePrimitive)) return;
835  out << " // Translation: " << GetName() << std::endl;
836  out << " dx = " << fTranslation[0] << ";" << std::endl;
837  out << " dy = " << fTranslation[1] << ";" << std::endl;
838  out << " dz = " << fTranslation[2] << ";" << std::endl;
839  out << " TGeoTranslation *" << GetPointerName() << " = new TGeoTranslation(\"" << GetName() << "\",dx,dy,dz);" << std::endl;
840  TObject::SetBit(kGeoSavePrimitive);
841 }
842 
843 /** \class TGeoRotation
844 \ingroup Geometry_classes
845 Class describing rotations. A rotation is a 3*3 array
846 Column vectors has to be orthogonal unit vectors.
847 */
848 
849 ClassImp(TGeoRotation);
850 
851 ////////////////////////////////////////////////////////////////////////////////
852 /// Default constructor.
853 
854 TGeoRotation::TGeoRotation()
855 {
856  for (Int_t i=0; i<9; i++) {
857  if (i%4) fRotationMatrix[i] = 0;
858  else fRotationMatrix[i] = 1.0;
859  }
860 }
861 
862 ////////////////////////////////////////////////////////////////////////////////
863 /// Copy ctor.
864 
865 TGeoRotation::TGeoRotation(const TGeoRotation &other)
866  :TGeoMatrix(other)
867 {
868  SetRotation(other);
869 }
870 
871 ////////////////////////////////////////////////////////////////////////////////
872 /// Copy ctor.
873 
874 TGeoRotation::TGeoRotation(const TGeoMatrix &other)
875  :TGeoMatrix(other)
876 {
877  ResetBit(kGeoTranslation);
878  ResetBit(kGeoScale);
879  SetRotation(other);
880 }
881 
882 ////////////////////////////////////////////////////////////////////////////////
883 /// Named rotation constructor
884 
885 TGeoRotation::TGeoRotation(const char *name)
886  :TGeoMatrix(name)
887 {
888  for (Int_t i=0; i<9; i++) {
889  if (i%4) fRotationMatrix[i] = 0;
890  else fRotationMatrix[i] = 1.0;
891  }
892 }
893 
894 ////////////////////////////////////////////////////////////////////////////////
895 /// Default rotation constructor with Euler angles. Phi is the rotation angle about
896 /// Z axis and is done first, theta is the rotation about new X and is done
897 /// second, psi is the rotation angle about new Z and is done third. All angles are in
898 /// degrees.
899 
900 TGeoRotation::TGeoRotation(const char *name, Double_t phi, Double_t theta, Double_t psi)
901  :TGeoMatrix(name)
902 {
903  SetAngles(phi, theta, psi);
904 }
905 
906 ////////////////////////////////////////////////////////////////////////////////
907 /// Rotation constructor a la GEANT3. Angles theta(i), phi(i) are the polar and azimuthal
908 /// angles of the (i) axis of the rotated system with respect to the initial non-rotated
909 /// system.
910 /// Example : the identity matrix (no rotation) is composed by
911 /// theta1=90, phi1=0, theta2=90, phi2=90, theta3=0, phi3=0
912 /// SetBit(kGeoRotation);
913 
914 TGeoRotation::TGeoRotation(const char *name, Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2,
915  Double_t theta3, Double_t phi3)
916  :TGeoMatrix(name)
917 {
918  SetAngles(theta1, phi1, theta2, phi2, theta3, phi3);
919 }
920 
921 ////////////////////////////////////////////////////////////////////////////////
922 /// Assignment from a general matrix
923 
924 TGeoRotation& TGeoRotation::operator = (const TGeoMatrix &other)
925 {
926  if (&other == this) return *this;
927  Bool_t registered = TestBit(kGeoRegistered);
928  TNamed::operator=(other);
929  SetRotation(other);
930  SetBit(kGeoRegistered,registered);
931  ResetBit(kGeoTranslation);
932  ResetBit(kGeoScale);
933  return *this;
934 }
935 
936 ////////////////////////////////////////////////////////////////////////////////
937 /// Composition
938 
939 TGeoRotation &TGeoRotation::operator*=(const TGeoRotation &right)
940 {
941  if (!right.IsRotation()) return *this;
942  MultiplyBy(&right, true);
943  return *this;
944 }
945 
946 TGeoRotation TGeoRotation::operator*(const TGeoRotation &right) const
947 {
948  TGeoRotation r = *this;
949  r *= right;
950  return r;
951 }
952 
953 TGeoHMatrix TGeoRotation::operator*(const TGeoMatrix &right) const
954 {
955  TGeoHMatrix t = *this;
956  t *= right;
957  return t;
958 }
959 
960 ////////////////////////////////////////////////////////////////////////////////
961 /// Is-equal operator
962 
963 Bool_t TGeoRotation::operator ==(const TGeoRotation &other) const
964 {
965  if (&other == this) return kTRUE;
966  const Double_t *rot = GetRotationMatrix();
967  const Double_t *orot = other.GetRotationMatrix();
968  for (auto i=0; i<9; i++)
969  if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
970  return kTRUE;
971 }
972 
973 ////////////////////////////////////////////////////////////////////////////////
974 /// Return a temporary inverse of this.
975 
976 TGeoHMatrix TGeoRotation::Inverse() const
977 {
978  TGeoHMatrix h;
979  h = *this;
980  Double_t newrot[9];
981  newrot[0] = fRotationMatrix[0];
982  newrot[1] = fRotationMatrix[3];
983  newrot[2] = fRotationMatrix[6];
984  newrot[3] = fRotationMatrix[1];
985  newrot[4] = fRotationMatrix[4];
986  newrot[5] = fRotationMatrix[7];
987  newrot[6] = fRotationMatrix[2];
988  newrot[7] = fRotationMatrix[5];
989  newrot[8] = fRotationMatrix[8];
990  h.SetRotation(newrot);
991  return h;
992 }
993 
994 ////////////////////////////////////////////////////////////////////////////////
995 /// Perform orthogonality test for rotation.
996 
997 Bool_t TGeoRotation::IsValid() const
998 {
999  const Double_t *r = fRotationMatrix;
1000  Double_t cij;
1001  for (Int_t i=0; i<2; i++) {
1002  for (Int_t j=i+1; j<3; j++) {
1003  // check columns
1004  cij = TMath::Abs(r[i]*r[j]+r[i+3]*r[j+3]+r[i+6]*r[j+6]);
1005  if (cij>1E-4) return kFALSE;
1006  // check rows
1007  cij = TMath::Abs(r[3*i]*r[3*j]+r[3*i+1]*r[3*j+1]+r[3*i+2]*r[3*j+2]);
1008  if (cij>1E-4) return kFALSE;
1009  }
1010  }
1011  return kTRUE;
1012 }
1013 
1014 ////////////////////////////////////////////////////////////////////////////////
1015 /// reset data members
1016 
1017 void TGeoRotation::Clear(Option_t *)
1018 {
1019  memcpy(fRotationMatrix,kIdentityMatrix,kN9);
1020  ResetBit(kGeoRotation);
1021 }
1022 
1023 ////////////////////////////////////////////////////////////////////////////////
1024 /// Perform a rotation about Z having the sine/cosine of the rotation angle.
1025 
1026 void TGeoRotation::FastRotZ(const Double_t *sincos)
1027 {
1028  fRotationMatrix[0] = sincos[1];
1029  fRotationMatrix[1] = -sincos[0];
1030  fRotationMatrix[3] = sincos[0];
1031  fRotationMatrix[4] = sincos[1];
1032  SetBit(kGeoRotation);
1033 }
1034 
1035 ////////////////////////////////////////////////////////////////////////////////
1036 /// Returns rotation angle about Z axis in degrees. If the rotation is a pure
1037 /// rotation about Z, fixX parameter does not matter, otherwise its meaning is:
1038 /// - fixX = true : result is the phi angle of the projection of the rotated X axis in the un-rotated XY
1039 /// - fixX = false : result is the phi angle of the projection of the rotated Y axis - 90 degrees
1040 
1041 Double_t TGeoRotation::GetPhiRotation(Bool_t fixX) const
1042 {
1043  Double_t phi;
1044  if (fixX) phi = 180.*TMath::ATan2(-fRotationMatrix[1],fRotationMatrix[4])/TMath::Pi();
1045  else phi = 180.*TMath::ATan2(fRotationMatrix[3], fRotationMatrix[0])/TMath::Pi();
1046  return phi;
1047 }
1048 
1049 ////////////////////////////////////////////////////////////////////////////////
1050 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
1051 
1052 void TGeoRotation::LocalToMaster(const Double_t *local, Double_t *master) const
1053 {
1054  const Double_t *rot = GetRotationMatrix();
1055  for (Int_t i=0; i<3; i++) {
1056  master[i] = local[0]*rot[3*i]
1057  + local[1]*rot[3*i+1]
1058  + local[2]*rot[3*i+2];
1059  }
1060 }
1061 
1062 ////////////////////////////////////////////////////////////////////////////////
1063 /// convert a point by multiplying its column vector (x, y, z, 1) to matrix
1064 
1065 void TGeoRotation::MasterToLocal(const Double_t *master, Double_t *local) const
1066 {
1067  const Double_t *rot = GetRotationMatrix();
1068  for (Int_t i=0; i<3; i++) {
1069  local[i] = master[0]*rot[i]
1070  + master[1]*rot[i+3]
1071  + master[2]*rot[i+6];
1072  }
1073 }
1074 
1075 ////////////////////////////////////////////////////////////////////////////////
1076 /// Make a clone of this matrix.
1077 
1078 TGeoMatrix *TGeoRotation::MakeClone() const
1079 {
1080  TGeoMatrix *matrix = new TGeoRotation(*this);
1081  return matrix;
1082 }
1083 
1084 ////////////////////////////////////////////////////////////////////////////////
1085 /// Rotate about X axis of the master frame with angle expressed in degrees.
1086 
1087 void TGeoRotation::RotateX(Double_t angle)
1088 {
1089  SetBit(kGeoRotation);
1090  Double_t phi = angle*TMath::DegToRad();
1091  Double_t c = TMath::Cos(phi);
1092  Double_t s = TMath::Sin(phi);
1093  Double_t v[9];
1094  v[0] = fRotationMatrix[0];
1095  v[1] = fRotationMatrix[1];
1096  v[2] = fRotationMatrix[2];
1097  v[3] = c*fRotationMatrix[3]-s*fRotationMatrix[6];
1098  v[4] = c*fRotationMatrix[4]-s*fRotationMatrix[7];
1099  v[5] = c*fRotationMatrix[5]-s*fRotationMatrix[8];
1100  v[6] = s*fRotationMatrix[3]+c*fRotationMatrix[6];
1101  v[7] = s*fRotationMatrix[4]+c*fRotationMatrix[7];
1102  v[8] = s*fRotationMatrix[5]+c*fRotationMatrix[8];
1103 
1104  memcpy(fRotationMatrix, v, kN9);
1105 }
1106 
1107 ////////////////////////////////////////////////////////////////////////////////
1108 /// Rotate about Y axis of the master frame with angle expressed in degrees.
1109 
1110 void TGeoRotation::RotateY(Double_t angle)
1111 {
1112  SetBit(kGeoRotation);
1113  Double_t phi = angle*TMath::DegToRad();
1114  Double_t c = TMath::Cos(phi);
1115  Double_t s = TMath::Sin(phi);
1116  Double_t v[9];
1117  v[0] = c*fRotationMatrix[0]+s*fRotationMatrix[6];
1118  v[1] = c*fRotationMatrix[1]+s*fRotationMatrix[7];
1119  v[2] = c*fRotationMatrix[2]+s*fRotationMatrix[8];
1120  v[3] = fRotationMatrix[3];
1121  v[4] = fRotationMatrix[4];
1122  v[5] = fRotationMatrix[5];
1123  v[6] = -s*fRotationMatrix[0]+c*fRotationMatrix[6];
1124  v[7] = -s*fRotationMatrix[1]+c*fRotationMatrix[7];
1125  v[8] = -s*fRotationMatrix[2]+c*fRotationMatrix[8];
1126 
1127  memcpy(fRotationMatrix, v, kN9);
1128 }
1129 
1130 ////////////////////////////////////////////////////////////////////////////////
1131 /// Rotate about Z axis of the master frame with angle expressed in degrees.
1132 
1133 void TGeoRotation::RotateZ(Double_t angle)
1134 {
1135  SetBit(kGeoRotation);
1136  Double_t phi = angle*TMath::DegToRad();
1137  Double_t c = TMath::Cos(phi);
1138  Double_t s = TMath::Sin(phi);
1139  Double_t v[9];
1140  v[0] = c*fRotationMatrix[0]-s*fRotationMatrix[3];
1141  v[1] = c*fRotationMatrix[1]-s*fRotationMatrix[4];
1142  v[2] = c*fRotationMatrix[2]-s*fRotationMatrix[5];
1143  v[3] = s*fRotationMatrix[0]+c*fRotationMatrix[3];
1144  v[4] = s*fRotationMatrix[1]+c*fRotationMatrix[4];
1145  v[5] = s*fRotationMatrix[2]+c*fRotationMatrix[5];
1146  v[6] = fRotationMatrix[6];
1147  v[7] = fRotationMatrix[7];
1148  v[8] = fRotationMatrix[8];
1149 
1150  memcpy(&fRotationMatrix[0],v,kN9);
1151 }
1152 
1153 ////////////////////////////////////////////////////////////////////////////////
1154 /// Multiply by a reflection respect to YZ.
1155 
1156 void TGeoRotation::ReflectX(Bool_t leftside, Bool_t)
1157 {
1158  if (leftside) {
1159  fRotationMatrix[0]=-fRotationMatrix[0];
1160  fRotationMatrix[1]=-fRotationMatrix[1];
1161  fRotationMatrix[2]=-fRotationMatrix[2];
1162  } else {
1163  fRotationMatrix[0]=-fRotationMatrix[0];
1164  fRotationMatrix[3]=-fRotationMatrix[3];
1165  fRotationMatrix[6]=-fRotationMatrix[6];
1166  }
1167  SetBit(kGeoRotation);
1168  SetBit(kGeoReflection, !IsReflection());
1169 }
1170 
1171 ////////////////////////////////////////////////////////////////////////////////
1172 /// Multiply by a reflection respect to ZX.
1173 
1174 void TGeoRotation::ReflectY(Bool_t leftside, Bool_t)
1175 {
1176  if (leftside) {
1177  fRotationMatrix[3]=-fRotationMatrix[3];
1178  fRotationMatrix[4]=-fRotationMatrix[4];
1179  fRotationMatrix[5]=-fRotationMatrix[5];
1180  } else {
1181  fRotationMatrix[1]=-fRotationMatrix[1];
1182  fRotationMatrix[4]=-fRotationMatrix[4];
1183  fRotationMatrix[7]=-fRotationMatrix[7];
1184  }
1185  SetBit(kGeoRotation);
1186  SetBit(kGeoReflection, !IsReflection());
1187 }
1188 
1189 ////////////////////////////////////////////////////////////////////////////////
1190 /// Multiply by a reflection respect to XY.
1191 
1192 void TGeoRotation::ReflectZ(Bool_t leftside, Bool_t)
1193 {
1194  if (leftside) {
1195  fRotationMatrix[6]=-fRotationMatrix[6];
1196  fRotationMatrix[7]=-fRotationMatrix[7];
1197  fRotationMatrix[8]=-fRotationMatrix[8];
1198  } else {
1199  fRotationMatrix[2]=-fRotationMatrix[2];
1200  fRotationMatrix[5]=-fRotationMatrix[5];
1201  fRotationMatrix[8]=-fRotationMatrix[8];
1202  }
1203  SetBit(kGeoRotation);
1204  SetBit(kGeoReflection, !IsReflection());
1205 }
1206 
1207 ////////////////////////////////////////////////////////////////////////////////
1208 /// Save a primitive as a C++ statement(s) on output stream "out".
1209 
1210 void TGeoRotation::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1211 {
1212  if (TestBit(kGeoSavePrimitive)) return;
1213  out << " // Rotation: " << GetName() << std::endl;
1214  Double_t th1,ph1,th2,ph2,th3,ph3;
1215  GetAngles(th1,ph1,th2,ph2,th3,ph3);
1216  out << " thx = " << th1 << "; phx = " << ph1 << ";" << std::endl;
1217  out << " thy = " << th2 << "; phy = " << ph2 << ";" << std::endl;
1218  out << " thz = " << th3 << "; phz = " << ph3 << ";" << std::endl;
1219  out << " TGeoRotation *" << GetPointerName() << " = new TGeoRotation(\"" << GetName() << "\",thx,phx,thy,phy,thz,phz);" << std::endl;
1220  TObject::SetBit(kGeoSavePrimitive);
1221 }
1222 
1223 ////////////////////////////////////////////////////////////////////////////////
1224 /// Copy rotation elements from other rotation matrix.
1225 
1226 void TGeoRotation::SetRotation(const TGeoMatrix &other)
1227 {
1228  SetBit(kGeoRotation, other.IsRotation());
1229  SetMatrix(other.GetRotationMatrix());
1230 }
1231 
1232 ////////////////////////////////////////////////////////////////////////////////
1233 /// Set matrix elements according to Euler angles. Phi is the rotation angle about
1234 /// Z axis and is done first, theta is the rotation about new X and is done
1235 /// second, psi is the rotation angle about new Z and is done third. All angles are in
1236 /// degrees.
1237 
1238 void TGeoRotation::SetAngles(Double_t phi, Double_t theta, Double_t psi)
1239 {
1240  Double_t degrad = TMath::Pi()/180.;
1241  Double_t sinphi = TMath::Sin(degrad*phi);
1242  Double_t cosphi = TMath::Cos(degrad*phi);
1243  Double_t sinthe = TMath::Sin(degrad*theta);
1244  Double_t costhe = TMath::Cos(degrad*theta);
1245  Double_t sinpsi = TMath::Sin(degrad*psi);
1246  Double_t cospsi = TMath::Cos(degrad*psi);
1247 
1248  fRotationMatrix[0] = cospsi*cosphi - costhe*sinphi*sinpsi;
1249  fRotationMatrix[1] = -sinpsi*cosphi - costhe*sinphi*cospsi;
1250  fRotationMatrix[2] = sinthe*sinphi;
1251  fRotationMatrix[3] = cospsi*sinphi + costhe*cosphi*sinpsi;
1252  fRotationMatrix[4] = -sinpsi*sinphi + costhe*cosphi*cospsi;
1253  fRotationMatrix[5] = -sinthe*cosphi;
1254  fRotationMatrix[6] = sinpsi*sinthe;
1255  fRotationMatrix[7] = cospsi*sinthe;
1256  fRotationMatrix[8] = costhe;
1257 
1258  if (!IsValid()) Error("SetAngles", "invalid rotation (Euler angles : phi=%f theta=%f psi=%f)",phi,theta,psi);
1259  CheckMatrix();
1260 }
1261 
1262 ////////////////////////////////////////////////////////////////////////////////
1263 /// Set matrix elements in the GEANT3 way
1264 
1265 void TGeoRotation::SetAngles(Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2,
1266  Double_t theta3, Double_t phi3)
1267 {
1268  Double_t degrad = TMath::Pi()/180.;
1269  fRotationMatrix[0] = TMath::Cos(degrad*phi1)*TMath::Sin(degrad*theta1);
1270  fRotationMatrix[3] = TMath::Sin(degrad*phi1)*TMath::Sin(degrad*theta1);
1271  fRotationMatrix[6] = TMath::Cos(degrad*theta1);
1272  fRotationMatrix[1] = TMath::Cos(degrad*phi2)*TMath::Sin(degrad*theta2);
1273  fRotationMatrix[4] = TMath::Sin(degrad*phi2)*TMath::Sin(degrad*theta2);
1274  fRotationMatrix[7] = TMath::Cos(degrad*theta2);
1275  fRotationMatrix[2] = TMath::Cos(degrad*phi3)*TMath::Sin(degrad*theta3);
1276  fRotationMatrix[5] = TMath::Sin(degrad*phi3)*TMath::Sin(degrad*theta3);
1277  fRotationMatrix[8] = TMath::Cos(degrad*theta3);
1278  // do the trick to eliminate most of the floating point errors
1279  for (Int_t i=0; i<9; i++) {
1280  if (TMath::Abs(fRotationMatrix[i])<1E-15) fRotationMatrix[i] = 0;
1281  if (TMath::Abs(fRotationMatrix[i]-1)<1E-15) fRotationMatrix[i] = 1;
1282  if (TMath::Abs(fRotationMatrix[i]+1)<1E-15) fRotationMatrix[i] = -1;
1283  }
1284  if (!IsValid()) Error("SetAngles", "invalid rotation (G3 angles, th1=%f phi1=%f, th2=%f ph2=%f, th3=%f phi3=%f)",
1285  theta1,phi1,theta2,phi2,theta3,phi3);
1286  CheckMatrix();
1287 }
1288 
1289 ////////////////////////////////////////////////////////////////////////////////
1290 /// Retrieve rotation angles
1291 
1292 void TGeoRotation::GetAngles(Double_t &theta1, Double_t &phi1, Double_t &theta2, Double_t &phi2,
1293  Double_t &theta3, Double_t &phi3) const
1294 {
1295  Double_t raddeg = 180./TMath::Pi();
1296  theta1 = raddeg*TMath::ACos(fRotationMatrix[6]);
1297  theta2 = raddeg*TMath::ACos(fRotationMatrix[7]);
1298  theta3 = raddeg*TMath::ACos(fRotationMatrix[8]);
1299  if (TMath::Abs(fRotationMatrix[0])<1E-6 && TMath::Abs(fRotationMatrix[3])<1E-6) phi1=0.;
1300  else phi1 = raddeg*TMath::ATan2(fRotationMatrix[3],fRotationMatrix[0]);
1301  if (phi1<0) phi1+=360.;
1302  if (TMath::Abs(fRotationMatrix[1])<1E-6 && TMath::Abs(fRotationMatrix[4])<1E-6) phi2=0.;
1303  else phi2 = raddeg*TMath::ATan2(fRotationMatrix[4],fRotationMatrix[1]);
1304  if (phi2<0) phi2+=360.;
1305  if (TMath::Abs(fRotationMatrix[2])<1E-6 && TMath::Abs(fRotationMatrix[5])<1E-6) phi3=0.;
1306  else phi3 = raddeg*TMath::ATan2(fRotationMatrix[5],fRotationMatrix[2]);
1307  if (phi3<0) phi3+=360.;
1308 }
1309 
1310 ////////////////////////////////////////////////////////////////////////////////
1311 /// Retrieve Euler angles.
1312 
1313 void TGeoRotation::GetAngles(Double_t &phi, Double_t &theta, Double_t &psi) const
1314 {
1315  const Double_t *m = fRotationMatrix;
1316  // Check if theta is 0 or 180.
1317  if (TMath::Abs(1.-TMath::Abs(m[8]))<1.e-9) {
1318  theta = TMath::ACos(m[8])*TMath::RadToDeg();
1319  phi = TMath::ATan2(-m[8]*m[1],m[0])*TMath::RadToDeg();
1320  psi = 0.; // convention, phi+psi matters
1321  return;
1322  }
1323  // sin(theta) != 0
1324  phi = TMath::ATan2(m[2],-m[5]);
1325  Double_t sphi = TMath::Sin(phi);
1326  if (TMath::Abs(sphi)<1.e-9) theta = -TMath::ASin(m[5]/TMath::Cos(phi))*TMath::RadToDeg();
1327  else theta = TMath::ASin(m[2]/sphi)*TMath::RadToDeg();
1328  phi *= TMath::RadToDeg();
1329  psi = TMath::ATan2(m[6],m[7])*TMath::RadToDeg();
1330 }
1331 
1332 ////////////////////////////////////////////////////////////////////////////////
1333 /// computes determinant of the rotation matrix
1334 
1335 Double_t TGeoRotation::Determinant() const
1336 {
1337  Double_t
1338  det = fRotationMatrix[0]*fRotationMatrix[4]*fRotationMatrix[8] +
1339  fRotationMatrix[3]*fRotationMatrix[7]*fRotationMatrix[2] +
1340  fRotationMatrix[6]*fRotationMatrix[1]*fRotationMatrix[5] -
1341  fRotationMatrix[2]*fRotationMatrix[4]*fRotationMatrix[6] -
1342  fRotationMatrix[5]*fRotationMatrix[7]*fRotationMatrix[0] -
1343  fRotationMatrix[8]*fRotationMatrix[1]*fRotationMatrix[3];
1344  return det;
1345 }
1346 
1347 ////////////////////////////////////////////////////////////////////////////////
1348 /// performes an orthogonality check and finds if the matrix is a reflection
1349 /// Warning("CheckMatrix", "orthogonality check not performed yet");
1350 
1351 void TGeoRotation::CheckMatrix()
1352 {
1353  if (Determinant() < 0) SetBit(kGeoReflection);
1354  Double_t dd = fRotationMatrix[0] + fRotationMatrix[4] + fRotationMatrix[8] - 3.;
1355  if (TMath::Abs(dd) < 1.E-12) ResetBit(kGeoRotation);
1356  else SetBit(kGeoRotation);
1357 }
1358 
1359 ////////////////////////////////////////////////////////////////////////////////
1360 /// Get the inverse rotation matrix (which is simply the transpose)
1361 
1362 void TGeoRotation::GetInverse(Double_t *invmat) const
1363 {
1364  if (!invmat) {
1365  Error("GetInverse", "no place to store the inverse matrix");
1366  return;
1367  }
1368  for (Int_t i=0; i<3; i++) {
1369  for (Int_t j=0; j<3; j++) {
1370  invmat[3*i+j] = fRotationMatrix[3*j+i];
1371  }
1372  }
1373 }
1374 
1375 ////////////////////////////////////////////////////////////////////////////////
1376 /// Multiply this rotation with the one specified by ROT.
1377 /// - after=TRUE (default): THIS*ROT
1378 /// - after=FALSE : ROT*THIS
1379 
1380 void TGeoRotation::MultiplyBy(const TGeoRotation *rot, Bool_t after)
1381 {
1382  const Double_t *matleft, *matright;
1383  SetBit(kGeoRotation);
1384  Double_t newmat[9] = {0};
1385  if (after) {
1386  matleft = &fRotationMatrix[0];
1387  matright = rot->GetRotationMatrix();
1388  } else {
1389  matleft = rot->GetRotationMatrix();
1390  matright = &fRotationMatrix[0];
1391  }
1392  for (Int_t i=0; i<3; i++) {
1393  for (Int_t j=0; j<3; j++) {
1394  for (Int_t k=0; k<3; k++) {
1395  newmat[3*i+j] += matleft[3*i+k] * matright[3*k+j];
1396  }
1397  }
1398  }
1399  memcpy(&fRotationMatrix[0], &newmat[0], kN9);
1400 }
1401 
1402 /** \class TGeoScale
1403 \ingroup Geometry_classes
1404 Class describing scale transformations. A scale is an
1405 array of 3 doubles (sx, sy, sz) multiplying elements 0, 5 and 10
1406 of the homogenous matrix. A scale is normalized : sx*sy*sz = 1
1407 */
1408 
1409 ClassImp(TGeoScale);
1410 
1411 ////////////////////////////////////////////////////////////////////////////////
1412 /// default constructor
1413 
1414 TGeoScale::TGeoScale()
1415 {
1416  SetBit(kGeoScale);
1417  for (Int_t i=0; i<3; i++) fScale[i] = 1.;
1418 }
1419 
1420 ////////////////////////////////////////////////////////////////////////////////
1421 /// Copy constructor
1422 
1423 TGeoScale::TGeoScale(const TGeoScale &other)
1424  :TGeoMatrix(other)
1425 {
1426  SetScale(other);
1427 }
1428 
1429 ////////////////////////////////////////////////////////////////////////////////
1430 /// Ctor. based on a general matrix
1431 
1432 TGeoScale::TGeoScale(const TGeoMatrix &other)
1433  :TGeoMatrix(other)
1434 {
1435  ResetBit(kGeoTranslation);
1436  ResetBit(kGeoRotation);
1437  SetScale(other);
1438 }
1439 
1440 ////////////////////////////////////////////////////////////////////////////////
1441 /// default constructor
1442 
1443 TGeoScale::TGeoScale(Double_t sx, Double_t sy, Double_t sz)
1444  :TGeoMatrix("")
1445 {
1446  SetBit(kGeoScale);
1447  SetScale(sx, sy, sz);
1448 }
1449 
1450 ////////////////////////////////////////////////////////////////////////////////
1451 /// default constructor
1452 
1453 TGeoScale::TGeoScale(const char *name, Double_t sx, Double_t sy, Double_t sz)
1454  :TGeoMatrix(name)
1455 {
1456  SetBit(kGeoScale);
1457  SetScale(sx, sy, sz);
1458 }
1459 
1460 ////////////////////////////////////////////////////////////////////////////////
1461 /// destructor
1462 
1463 TGeoScale::~TGeoScale()
1464 {
1465 }
1466 
1467 ////////////////////////////////////////////////////////////////////////////////
1468 /// Assignment from a general matrix
1469 
1470 TGeoScale& TGeoScale::operator = (const TGeoMatrix &matrix)
1471 {
1472  if (&matrix == this) return *this;
1473  Bool_t registered = TestBit(kGeoRegistered);
1474  TNamed::operator=(matrix);
1475  SetScale(matrix);
1476  SetBit(kGeoRegistered,registered);
1477  ResetBit(kGeoTranslation);
1478  ResetBit(kGeoRotation);
1479  return *this;
1480 }
1481 
1482 ////////////////////////////////////////////////////////////////////////////////
1483 /// Scale composition
1484 
1485 TGeoScale &TGeoScale::operator*=(const TGeoScale &right)
1486 {
1487  const Double_t *scl = right.GetScale();
1488  fScale[0] *= scl[0];
1489  fScale[1] *= scl[1];
1490  fScale[2] *= scl[2];
1491  SetBit(kGeoReflection, fScale[0] * fScale[1] * fScale[2] < 0);
1492  if (!IsScale()) SetBit(kGeoScale, right.IsScale());
1493  return *this;
1494 }
1495 
1496 TGeoScale TGeoScale::operator*(const TGeoScale &right) const
1497 {
1498  TGeoScale s = *this;
1499  s *= right;
1500  return s;
1501 }
1502 
1503 TGeoHMatrix TGeoScale::operator*(const TGeoMatrix &right) const
1504 {
1505  TGeoHMatrix t = *this;
1506  t *= right;
1507  return t;
1508 }
1509 
1510 ////////////////////////////////////////////////////////////////////////////////
1511 /// Is-equal operator
1512 
1513 Bool_t TGeoScale::operator ==(const TGeoScale &other) const
1514 {
1515  if (&other == this) return kTRUE;
1516  const Double_t *scl = GetScale();
1517  const Double_t *oscl = other.GetScale();
1518  for (auto i=0; i<3; i++)
1519  if (TMath::Abs(scl[i]-oscl[i])>1.E-10) return kFALSE;
1520  return kTRUE;
1521 }
1522 
1523 ////////////////////////////////////////////////////////////////////////////////
1524 /// Return a temporary inverse of this.
1525 
1526 TGeoHMatrix TGeoScale::Inverse() const
1527 {
1528  TGeoHMatrix h;
1529  h = *this;
1530  Double_t scale[3];
1531  scale[0] = 1./fScale[0];
1532  scale[1] = 1./fScale[1];
1533  scale[2] = 1./fScale[2];
1534  h.SetScale(scale);
1535  return h;
1536 }
1537 
1538 ////////////////////////////////////////////////////////////////////////////////
1539 /// scale setter
1540 
1541 void TGeoScale::SetScale(Double_t sx, Double_t sy, Double_t sz)
1542 {
1543  if (TMath::Abs(sx*sy*sz) < 1.E-10) {
1544  Error("SetScale", "Invalid scale %f, %f, %f for transformation %s",sx,sy,sx,GetName());
1545  return;
1546  }
1547  fScale[0] = sx;
1548  fScale[1] = sy;
1549  fScale[2] = sz;
1550  if (sx*sy*sz<0) SetBit(kGeoReflection);
1551  else SetBit(kGeoReflection, kFALSE);
1552 }
1553 
1554 ////////////////////////////////////////////////////////////////////////////////
1555 /// Set scale from other transformation
1556 
1557 void TGeoScale::SetScale(const TGeoMatrix &other)
1558 {
1559  SetBit(kGeoScale, other.IsScale());
1560  SetBit(kGeoReflection, other.IsReflection());
1561  memcpy(fScale, other.GetScale(), kN3);
1562 }
1563 
1564 ////////////////////////////////////////////////////////////////////////////////
1565 /// Convert a local point to the master frame.
1566 
1567 void TGeoScale::LocalToMaster(const Double_t *local, Double_t *master) const
1568 {
1569  master[0] = local[0]*fScale[0];
1570  master[1] = local[1]*fScale[1];
1571  master[2] = local[2]*fScale[2];
1572 }
1573 
1574 ////////////////////////////////////////////////////////////////////////////////
1575 /// Convert the local distance along unit vector DIR to master frame. If DIR
1576 /// is not specified perform a conversion such as the returned distance is the
1577 /// the minimum for all possible directions.
1578 
1579 Double_t TGeoScale::LocalToMaster(Double_t dist, const Double_t *dir) const
1580 {
1581  Double_t scale;
1582  if (!dir) {
1583  scale = TMath::Abs(fScale[0]);
1584  if (TMath::Abs(fScale[1])<scale) scale = TMath::Abs(fScale[1]);
1585  if (TMath::Abs(fScale[2])<scale) scale = TMath::Abs(fScale[2]);
1586  } else {
1587  scale = fScale[0]*fScale[0]*dir[0]*dir[0] +
1588  fScale[1]*fScale[1]*dir[1]*dir[1] +
1589  fScale[2]*fScale[2]*dir[2]*dir[2];
1590  scale = TMath::Sqrt(scale);
1591  }
1592  return scale*dist;
1593 }
1594 
1595 ////////////////////////////////////////////////////////////////////////////////
1596 /// Make a clone of this matrix.
1597 
1598 TGeoMatrix *TGeoScale::MakeClone() const
1599 {
1600  TGeoMatrix *matrix = new TGeoScale(*this);
1601  return matrix;
1602 }
1603 
1604 ////////////////////////////////////////////////////////////////////////////////
1605 /// Convert a global point to local frame.
1606 
1607 void TGeoScale::MasterToLocal(const Double_t *master, Double_t *local) const
1608 {
1609  local[0] = master[0]/fScale[0];
1610  local[1] = master[1]/fScale[1];
1611  local[2] = master[2]/fScale[2];
1612 }
1613 
1614 ////////////////////////////////////////////////////////////////////////////////
1615 /// Convert the distance along unit vector DIR to local frame. If DIR
1616 /// is not specified perform a conversion such as the returned distance is the
1617 /// the minimum for all possible directions.
1618 
1619 Double_t TGeoScale::MasterToLocal(Double_t dist, const Double_t *dir) const
1620 {
1621  Double_t scale;
1622  if (!dir) {
1623  scale = TMath::Abs(fScale[0]);
1624  if (TMath::Abs(fScale[1])>scale) scale = TMath::Abs(fScale[1]);
1625  if (TMath::Abs(fScale[2])>scale) scale = TMath::Abs(fScale[2]);
1626  scale = 1./scale;
1627  } else {
1628  scale = (dir[0]*dir[0])/(fScale[0]*fScale[0]) +
1629  (dir[1]*dir[1])/(fScale[1]*fScale[1]) +
1630  (dir[2]*dir[2])/(fScale[2]*fScale[2]);
1631  scale = TMath::Sqrt(scale);
1632  }
1633  return scale*dist;
1634 }
1635 
1636 /** \class TGeoCombiTrans
1637 \ingroup Geometry_classes
1638 Class describing rotation + translation. Most frequently used in the description
1639 of TGeoNode 's
1640 */
1641 
1642 ClassImp(TGeoCombiTrans);
1643 
1644 ////////////////////////////////////////////////////////////////////////////////
1645 /// dummy ctor
1646 
1647 TGeoCombiTrans::TGeoCombiTrans()
1648 {
1649  for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
1650  fRotation = 0;
1651 }
1652 
1653 ////////////////////////////////////////////////////////////////////////////////
1654 /// Copy ctor from generic matrix.
1655 
1656 TGeoCombiTrans::TGeoCombiTrans(const TGeoMatrix &other)
1657  :TGeoMatrix(other)
1658 {
1659  ResetBit(kGeoScale);
1660  if (other.IsTranslation()) {
1661  SetBit(kGeoTranslation);
1662  memcpy(fTranslation,other.GetTranslation(),kN3);
1663  } else {
1664  for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
1665  }
1666  if (other.IsRotation()) {
1667  SetBit(kGeoRotation);
1668  SetBit(kGeoMatrixOwned);
1669  fRotation = new TGeoRotation(other);
1670  }
1671  else fRotation = 0;
1672 }
1673 
1674 ////////////////////////////////////////////////////////////////////////////////
1675 /// Constructor from a translation and a rotation.
1676 
1677 TGeoCombiTrans::TGeoCombiTrans(const TGeoTranslation &tr, const TGeoRotation &rot)
1678 {
1679  if (tr.IsTranslation()) {
1680  SetBit(kGeoTranslation);
1681  const Double_t *trans = tr.GetTranslation();
1682  memcpy(fTranslation, trans, kN3);
1683  } else {
1684  for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
1685  }
1686  if (rot.IsRotation()) {
1687  SetBit(kGeoRotation);
1688  SetBit(kGeoMatrixOwned);
1689  fRotation = new TGeoRotation(rot);
1690  SetBit(kGeoReflection, rot.TestBit(kGeoReflection));
1691  }
1692  else fRotation = 0;
1693 }
1694 
1695 ////////////////////////////////////////////////////////////////////////////////
1696 /// Named ctor.
1697 
1698 TGeoCombiTrans::TGeoCombiTrans(const char *name)
1699  :TGeoMatrix(name)
1700 {
1701  for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
1702  fRotation = 0;
1703 }
1704 
1705 ////////////////////////////////////////////////////////////////////////////////
1706 /// Constructor from a translation specified by X,Y,Z and a pointer to a rotation. The rotation will not be owned by this.
1707 
1708 TGeoCombiTrans::TGeoCombiTrans(Double_t dx, Double_t dy, Double_t dz, TGeoRotation *rot)
1709  :TGeoMatrix("")
1710 {
1711  SetTranslation(dx, dy, dz);
1712  fRotation = 0;
1713  SetRotation(rot);
1714 }
1715 
1716 ////////////////////////////////////////////////////////////////////////////////
1717 /// Named ctor
1718 
1719 TGeoCombiTrans::TGeoCombiTrans(const char *name, Double_t dx, Double_t dy, Double_t dz, TGeoRotation *rot)
1720  :TGeoMatrix(name)
1721 {
1722  SetTranslation(dx, dy, dz);
1723  fRotation = 0;
1724  SetRotation(rot);
1725 }
1726 
1727 ////////////////////////////////////////////////////////////////////////////////
1728 /// Assignment operator with generic matrix.
1729 
1730 TGeoCombiTrans &TGeoCombiTrans::operator=(const TGeoMatrix &matrix)
1731 {
1732  if (&matrix == this) return *this;
1733  Bool_t registered = TestBit(kGeoRegistered);
1734  Clear();
1735  TNamed::operator=(matrix);
1736 
1737  if (matrix.IsTranslation()) {
1738  memcpy(fTranslation,matrix.GetTranslation(),kN3);
1739  }
1740  if (matrix.IsRotation()) {
1741  if (!fRotation) {
1742  fRotation = new TGeoRotation();
1743  SetBit(kGeoMatrixOwned);
1744  } else {
1745  if (!TestBit(kGeoMatrixOwned)) {
1746  fRotation = new TGeoRotation();
1747  SetBit(kGeoMatrixOwned);
1748  }
1749  }
1750  fRotation->SetMatrix(matrix.GetRotationMatrix());
1751  fRotation->SetBit(kGeoReflection, matrix.TestBit(kGeoReflection));
1752  fRotation->SetBit(kGeoRotation);
1753  } else {
1754  if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
1755  ResetBit(kGeoMatrixOwned);
1756  fRotation = 0;
1757  }
1758  SetBit(kGeoRegistered,registered);
1759  ResetBit(kGeoScale);
1760  return *this;
1761 }
1762 
1763 ////////////////////////////////////////////////////////////////////////////////
1764 /// Is-equal operator
1765 
1766 Bool_t TGeoCombiTrans::operator==(const TGeoMatrix &other) const
1767 {
1768  if (&other == this) return kTRUE;
1769  const Double_t *tr = GetTranslation();
1770  const Double_t *otr = other.GetTranslation();
1771  for (auto i=0; i<3; i++) if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
1772  const Double_t *rot = GetRotationMatrix();
1773  const Double_t *orot = other.GetRotationMatrix();
1774  for (auto i=0; i<9; i++) if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
1775  return kTRUE;
1776 }
1777 
1778 ////////////////////////////////////////////////////////////////////////////////
1779 /// Composition
1780 
1781 TGeoCombiTrans &TGeoCombiTrans::operator*=(const TGeoMatrix &right)
1782 {
1783  Multiply(&right);
1784  return *this;
1785 }
1786 
1787 TGeoCombiTrans TGeoCombiTrans::operator*(const TGeoMatrix &right) const
1788 {
1789  TGeoHMatrix h = *this;
1790  h *= right;
1791  return h;
1792 }
1793 
1794 ////////////////////////////////////////////////////////////////////////////////
1795 /// destructor
1796 
1797 TGeoCombiTrans::~TGeoCombiTrans()
1798 {
1799  if (fRotation) {
1800  if(TestBit(TGeoMatrix::kGeoMatrixOwned) && !fRotation->IsRegistered()) delete fRotation;
1801  }
1802 }
1803 
1804 ////////////////////////////////////////////////////////////////////////////////
1805 /// Reset translation/rotation to identity
1806 
1807 void TGeoCombiTrans::Clear(Option_t *)
1808 {
1809  if (IsTranslation()) {
1810  ResetBit(kGeoTranslation);
1811  memset(fTranslation, 0, kN3);
1812  }
1813  if (fRotation) {
1814  if (TestBit(kGeoMatrixOwned)) delete fRotation;
1815  fRotation = 0;
1816  }
1817  ResetBit(kGeoRotation);
1818  ResetBit(kGeoTranslation);
1819  ResetBit(kGeoMatrixOwned);
1820 }
1821 
1822 ////////////////////////////////////////////////////////////////////////////////
1823 /// Return a temporary inverse of this.
1824 
1825 TGeoHMatrix TGeoCombiTrans::Inverse() const
1826 {
1827  TGeoHMatrix h;
1828  h = *this;
1829  Bool_t is_tr = IsTranslation();
1830  Bool_t is_rot = IsRotation();
1831  Double_t tr[3];
1832  Double_t newrot[9];
1833  const Double_t *rot = GetRotationMatrix();
1834  tr[0] = -fTranslation[0]*rot[0] - fTranslation[1]*rot[3] - fTranslation[2]*rot[6];
1835  tr[1] = -fTranslation[0]*rot[1] - fTranslation[1]*rot[4] - fTranslation[2]*rot[7];
1836  tr[2] = -fTranslation[0]*rot[2] - fTranslation[1]*rot[5] - fTranslation[2]*rot[8];
1837  h.SetTranslation(tr);
1838  newrot[0] = rot[0];
1839  newrot[1] = rot[3];
1840  newrot[2] = rot[6];
1841  newrot[3] = rot[1];
1842  newrot[4] = rot[4];
1843  newrot[5] = rot[7];
1844  newrot[6] = rot[2];
1845  newrot[7] = rot[5];
1846  newrot[8] = rot[8];
1847  h.SetRotation(newrot);
1848  h.SetBit(kGeoTranslation,is_tr);
1849  h.SetBit(kGeoRotation,is_rot);
1850  return h;
1851 }
1852 
1853 ////////////////////////////////////////////////////////////////////////////////
1854 /// Make a clone of this matrix.
1855 
1856 TGeoMatrix *TGeoCombiTrans::MakeClone() const
1857 {
1858  TGeoMatrix *matrix = new TGeoCombiTrans(*this);
1859  return matrix;
1860 }
1861 
1862 ////////////////////////////////////////////////////////////////////////////////
1863 /// multiply to the right with an other transformation
1864 /// if right is identity matrix, just return
1865 
1866 void TGeoCombiTrans::Multiply(const TGeoMatrix *right)
1867 {
1868  if (right->IsIdentity()) return;
1869  TGeoHMatrix h = *this;
1870  h.Multiply(right);
1871  operator=(h);
1872 }
1873 
1874 ////////////////////////////////////////////////////////////////////////////////
1875 /// Register the matrix in the current manager, which will become the owner.
1876 
1877 void TGeoCombiTrans::RegisterYourself()
1878 {
1879  TGeoMatrix::RegisterYourself();
1880  if (fRotation && fRotation->IsRotation()) fRotation->RegisterYourself();
1881 }
1882 
1883 ////////////////////////////////////////////////////////////////////////////////
1884 /// Rotate about X axis with angle expressed in degrees.
1885 
1886 void TGeoCombiTrans::RotateX(Double_t angle)
1887 {
1888  if (!fRotation || !TestBit(kGeoMatrixOwned)) {
1889  if (fRotation) fRotation = new TGeoRotation(*fRotation);
1890  else fRotation = new TGeoRotation();
1891  SetBit(kGeoMatrixOwned);
1892  }
1893  SetBit(kGeoRotation);
1894  const Double_t *rot = fRotation->GetRotationMatrix();
1895  Double_t phi = angle*TMath::DegToRad();
1896  Double_t c = TMath::Cos(phi);
1897  Double_t s = TMath::Sin(phi);
1898  Double_t v[9];
1899  v[0] = rot[0];
1900  v[1] = rot[1];
1901  v[2] = rot[2];
1902  v[3] = c*rot[3]-s*rot[6];
1903  v[4] = c*rot[4]-s*rot[7];
1904  v[5] = c*rot[5]-s*rot[8];
1905  v[6] = s*rot[3]+c*rot[6];
1906  v[7] = s*rot[4]+c*rot[7];
1907  v[8] = s*rot[5]+c*rot[8];
1908  fRotation->SetMatrix(v);
1909  fRotation->SetBit(kGeoRotation);
1910  if (!IsTranslation()) return;
1911  v[0] = fTranslation[0];
1912  v[1] = c*fTranslation[1]-s*fTranslation[2];
1913  v[2] = s*fTranslation[1]+c*fTranslation[2];
1914  memcpy(fTranslation,v,kN3);
1915 }
1916 
1917 ////////////////////////////////////////////////////////////////////////////////
1918 /// Rotate about Y axis with angle expressed in degrees.
1919 
1920 void TGeoCombiTrans::RotateY(Double_t angle)
1921 {
1922  if (!fRotation || !TestBit(kGeoMatrixOwned)) {
1923  if (fRotation) fRotation = new TGeoRotation(*fRotation);
1924  else fRotation = new TGeoRotation();
1925  SetBit(kGeoMatrixOwned);
1926  }
1927  SetBit(kGeoRotation);
1928  const Double_t *rot = fRotation->GetRotationMatrix();
1929  Double_t phi = angle*TMath::DegToRad();
1930  Double_t c = TMath::Cos(phi);
1931  Double_t s = TMath::Sin(phi);
1932  Double_t v[9];
1933  v[0] = c*rot[0]+s*rot[6];
1934  v[1] = c*rot[1]+s*rot[7];
1935  v[2] = c*rot[2]+s*rot[8];
1936  v[3] = rot[3];
1937  v[4] = rot[4];
1938  v[5] = rot[5];
1939  v[6] = -s*rot[0]+c*rot[6];
1940  v[7] = -s*rot[1]+c*rot[7];
1941  v[8] = -s*rot[2]+c*rot[8];
1942  fRotation->SetMatrix(v);
1943  fRotation->SetBit(kGeoRotation);
1944  if (!IsTranslation()) return;
1945  v[0] = c*fTranslation[0]+s*fTranslation[2];
1946  v[1] = fTranslation[1];
1947  v[2] = -s*fTranslation[0]+c*fTranslation[2];
1948  memcpy(fTranslation,v,kN3);
1949 }
1950 
1951 ////////////////////////////////////////////////////////////////////////////////
1952 /// Rotate about Z axis with angle expressed in degrees.
1953 
1954 void TGeoCombiTrans::RotateZ(Double_t angle)
1955 {
1956  if (!fRotation || !TestBit(kGeoMatrixOwned)) {
1957  if (fRotation) fRotation = new TGeoRotation(*fRotation);
1958  else fRotation = new TGeoRotation();
1959  SetBit(kGeoMatrixOwned);
1960  }
1961  SetBit(kGeoRotation);
1962  const Double_t *rot = fRotation->GetRotationMatrix();
1963  Double_t phi = angle*TMath::DegToRad();
1964  Double_t c = TMath::Cos(phi);
1965  Double_t s = TMath::Sin(phi);
1966  Double_t v[9];
1967  v[0] = c*rot[0]-s*rot[3];
1968  v[1] = c*rot[1]-s*rot[4];
1969  v[2] = c*rot[2]-s*rot[5];
1970  v[3] = s*rot[0]+c*rot[3];
1971  v[4] = s*rot[1]+c*rot[4];
1972  v[5] = s*rot[2]+c*rot[5];
1973  v[6] = rot[6];
1974  v[7] = rot[7];
1975  v[8] = rot[8];
1976  fRotation->SetMatrix(v);
1977  fRotation->SetBit(kGeoRotation);
1978  if (!IsTranslation()) return;
1979  v[0] = c*fTranslation[0]-s*fTranslation[1];
1980  v[1] = s*fTranslation[0]+c*fTranslation[1];
1981  v[2] = fTranslation[2];
1982  memcpy(fTranslation,v,kN3);
1983 }
1984 
1985 ////////////////////////////////////////////////////////////////////////////////
1986 /// Multiply by a reflection respect to YZ.
1987 
1988 void TGeoCombiTrans::ReflectX(Bool_t leftside, Bool_t rotonly)
1989 {
1990  if (leftside && !rotonly) fTranslation[0] = -fTranslation[0];
1991  if (!fRotation || !TestBit(kGeoMatrixOwned)) {
1992  if (fRotation) fRotation = new TGeoRotation(*fRotation);
1993  else fRotation = new TGeoRotation();
1994  SetBit(kGeoMatrixOwned);
1995  }
1996  SetBit(kGeoRotation);
1997  fRotation->ReflectX(leftside);
1998  SetBit(kGeoReflection, !IsReflection());
1999 }
2000 
2001 ////////////////////////////////////////////////////////////////////////////////
2002 /// Multiply by a reflection respect to ZX.
2003 
2004 void TGeoCombiTrans::ReflectY(Bool_t leftside, Bool_t rotonly)
2005 {
2006  if (leftside && !rotonly) fTranslation[1] = -fTranslation[1];
2007  if (!fRotation || !TestBit(kGeoMatrixOwned)) {
2008  if (fRotation) fRotation = new TGeoRotation(*fRotation);
2009  else fRotation = new TGeoRotation();
2010  SetBit(kGeoMatrixOwned);
2011  }
2012  SetBit(kGeoRotation);
2013  fRotation->ReflectY(leftside);
2014  SetBit(kGeoReflection, !IsReflection());
2015 }
2016 
2017 ////////////////////////////////////////////////////////////////////////////////
2018 /// Multiply by a reflection respect to XY.
2019 
2020 void TGeoCombiTrans::ReflectZ(Bool_t leftside, Bool_t rotonly)
2021 {
2022  if (leftside && !rotonly) fTranslation[2] = -fTranslation[2];
2023  if (!fRotation || !TestBit(kGeoMatrixOwned)) {
2024  if (fRotation) fRotation = new TGeoRotation(*fRotation);
2025  else fRotation = new TGeoRotation();
2026  SetBit(kGeoMatrixOwned);
2027  }
2028  SetBit(kGeoRotation);
2029  fRotation->ReflectZ(leftside);
2030  SetBit(kGeoReflection, !IsReflection());
2031 }
2032 
2033 ////////////////////////////////////////////////////////////////////////////////
2034 /// Save a primitive as a C++ statement(s) on output stream "out".
2035 
2036 void TGeoCombiTrans::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
2037 {
2038  if (TestBit(kGeoSavePrimitive)) return;
2039  out << " // Combi transformation: " << GetName() << std::endl;
2040  out << " dx = " << fTranslation[0] << ";" << std::endl;
2041  out << " dy = " << fTranslation[1] << ";" << std::endl;
2042  out << " dz = " << fTranslation[2] << ";" << std::endl;
2043  if (fRotation && fRotation->IsRotation()) {
2044  fRotation->SavePrimitive(out,option);
2045  out << " " << GetPointerName() << " = new TGeoCombiTrans(\"" << GetName() << "\", dx,dy,dz,";
2046  out << fRotation->GetPointerName() << ");" << std::endl;
2047  } else {
2048  out << " " << GetPointerName() << " = new TGeoCombiTrans(\"" << GetName() << "\");" << std::endl;
2049  out << " " << GetPointerName() << "->SetTranslation(dx,dy,dz);" << std::endl;
2050  }
2051  TObject::SetBit(kGeoSavePrimitive);
2052 }
2053 
2054 ////////////////////////////////////////////////////////////////////////////////
2055 /// Assign a foreign rotation to the combi. The rotation is NOT owned by this.
2056 
2057 void TGeoCombiTrans::SetRotation(const TGeoRotation *rot)
2058 {
2059  if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
2060  fRotation = 0;
2061  ResetBit(TGeoMatrix::kGeoMatrixOwned);
2062  ResetBit(kGeoRotation);
2063  ResetBit(kGeoReflection);
2064  if (!rot) return;
2065  if (!rot->IsRotation()) return;
2066 
2067  SetBit(kGeoRotation);
2068  SetBit(kGeoReflection, rot->TestBit(kGeoReflection));
2069  TGeoRotation *rr = (TGeoRotation*)rot;
2070  fRotation = rr;
2071 }
2072 
2073 ////////////////////////////////////////////////////////////////////////////////
2074 /// Copy the rotation from another one.
2075 
2076 void TGeoCombiTrans::SetRotation(const TGeoRotation &rot)
2077 {
2078  if (fRotation && TestBit(kGeoMatrixOwned)) delete fRotation;
2079  fRotation = 0;
2080  if (!rot.IsRotation()) {
2081  ResetBit(kGeoRotation);
2082  ResetBit(kGeoReflection);
2083  ResetBit(TGeoMatrix::kGeoMatrixOwned);
2084  return;
2085  }
2086 
2087  SetBit(kGeoRotation);
2088  SetBit(kGeoReflection, rot.TestBit(kGeoReflection));
2089  fRotation = new TGeoRotation(rot);
2090  SetBit(kGeoMatrixOwned);
2091 }
2092 
2093 ////////////////////////////////////////////////////////////////////////////////
2094 /// copy the translation component
2095 
2096 void TGeoCombiTrans::SetTranslation(const TGeoTranslation &tr)
2097 {
2098  if (tr.IsTranslation()) {
2099  SetBit(kGeoTranslation);
2100  const Double_t *trans = tr.GetTranslation();
2101  memcpy(fTranslation, trans, kN3);
2102  } else {
2103  if (!IsTranslation()) return;
2104  memset(fTranslation, 0, kN3);
2105  ResetBit(kGeoTranslation);
2106  }
2107 }
2108 
2109 ////////////////////////////////////////////////////////////////////////////////
2110 /// set the translation component
2111 
2112 void TGeoCombiTrans::SetTranslation(Double_t dx, Double_t dy, Double_t dz)
2113 {
2114  fTranslation[0] = dx;
2115  fTranslation[1] = dy;
2116  fTranslation[2] = dz;
2117  if (fTranslation[0] || fTranslation[1] || fTranslation[2]) SetBit(kGeoTranslation);
2118  else ResetBit(kGeoTranslation);
2119 }
2120 
2121 ////////////////////////////////////////////////////////////////////////////////
2122 /// set the translation component
2123 
2124 void TGeoCombiTrans::SetTranslation(Double_t *vect)
2125 {
2126  fTranslation[0] = vect[0];
2127  fTranslation[1] = vect[1];
2128  fTranslation[2] = vect[2];
2129  if (fTranslation[0] || fTranslation[1] || fTranslation[2]) SetBit(kGeoTranslation);
2130  else ResetBit(kGeoTranslation);
2131 }
2132 
2133 ////////////////////////////////////////////////////////////////////////////////
2134 /// get the rotation array
2135 
2136 const Double_t *TGeoCombiTrans::GetRotationMatrix() const
2137 {
2138  if (!fRotation) return kIdentityMatrix;
2139  return fRotation->GetRotationMatrix();
2140 }
2141 
2142 /** \class TGeoGenTrans
2143 \ingroup Geometry_classes
2144 Most general transformation, holding a translation, a rotation and a scale
2145 */
2146 
2147 ClassImp(TGeoGenTrans);
2148 
2149 ////////////////////////////////////////////////////////////////////////////////
2150 /// dummy ctor
2151 
2152 TGeoGenTrans::TGeoGenTrans()
2153 {
2154  SetBit(kGeoGenTrans);
2155  for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
2156  for (Int_t j=0; j<3; j++) fScale[j] = 1.0;
2157  fRotation = 0;
2158 }
2159 
2160 ////////////////////////////////////////////////////////////////////////////////
2161 /// constructor
2162 
2163 TGeoGenTrans::TGeoGenTrans(const char *name)
2164  :TGeoCombiTrans(name)
2165 {
2166  SetBit(kGeoGenTrans);
2167  for (Int_t i=0; i<3; i++) fTranslation[i] = 0.0;
2168  for (Int_t j=0; j<3; j++) fScale[j] = 1.0;
2169  fRotation = 0;
2170 }
2171 
2172 ////////////////////////////////////////////////////////////////////////////////
2173 /// constructor
2174 
2175 TGeoGenTrans::TGeoGenTrans(Double_t dx, Double_t dy, Double_t dz,
2176  Double_t sx, Double_t sy, Double_t sz, TGeoRotation *rot)
2177  :TGeoCombiTrans("")
2178 {
2179  SetBit(kGeoGenTrans);
2180  SetTranslation(dx, dy, dz);
2181  SetScale(sx, sy, sz);
2182  SetRotation(rot);
2183 }
2184 
2185 ////////////////////////////////////////////////////////////////////////////////
2186 /// constructor
2187 
2188 TGeoGenTrans::TGeoGenTrans(const char *name, Double_t dx, Double_t dy, Double_t dz,
2189  Double_t sx, Double_t sy, Double_t sz, TGeoRotation *rot)
2190  :TGeoCombiTrans(name)
2191 {
2192  SetBit(kGeoGenTrans);
2193  SetTranslation(dx, dy, dz);
2194  SetScale(sx, sy, sz);
2195  SetRotation(rot);
2196 }
2197 
2198 ////////////////////////////////////////////////////////////////////////////////
2199 /// destructor
2200 
2201 TGeoGenTrans::~TGeoGenTrans()
2202 {
2203 }
2204 
2205 ////////////////////////////////////////////////////////////////////////////////
2206 /// clear the fields of this transformation
2207 
2208 void TGeoGenTrans::Clear(Option_t *)
2209 {
2210  memset(&fTranslation[0], 0, kN3);
2211  memset(&fScale[0], 0, kN3);
2212  if (fRotation) fRotation->Clear();
2213 }
2214 
2215 ////////////////////////////////////////////////////////////////////////////////
2216 /// set the scale
2217 
2218 void TGeoGenTrans::SetScale(Double_t sx, Double_t sy, Double_t sz)
2219 {
2220  if (sx<1.E-5 || sy<1.E-5 || sz<1.E-5) {
2221  Error("ctor", "Invalid scale");
2222  return;
2223  }
2224  fScale[0] = sx;
2225  fScale[1] = sy;
2226  fScale[2] = sz;
2227 }
2228 
2229 ////////////////////////////////////////////////////////////////////////////////
2230 /// Return a temporary inverse of this.
2231 
2232 TGeoHMatrix TGeoGenTrans::Inverse() const
2233 {
2234  TGeoHMatrix h = *this;
2235  return h;
2236 }
2237 
2238 ////////////////////////////////////////////////////////////////////////////////
2239 /// A scale transformation should be normalized by sx*sy*sz factor
2240 
2241 Bool_t TGeoGenTrans::Normalize()
2242 {
2243  Double_t normfactor = fScale[0]*fScale[1]*fScale[2];
2244  if (normfactor <= 1E-5) return kFALSE;
2245  for (Int_t i=0; i<3; i++)
2246  fScale[i] /= normfactor;
2247  return kTRUE;
2248 }
2249 
2250 /** \class TGeoIdentity
2251 \ingroup Geometry_classes
2252 An identity transformation. It holds no data member
2253 and returns pointers to static null translation and identity
2254 transformations for rotation and scale
2255 */
2256 
2257 ClassImp(TGeoIdentity);
2258 
2259 ////////////////////////////////////////////////////////////////////////////////
2260 /// dummy ctor
2261 
2262 TGeoIdentity::TGeoIdentity()
2263 {
2264  if (!gGeoIdentity) gGeoIdentity = this;
2265  RegisterYourself();
2266 }
2267 
2268 ////////////////////////////////////////////////////////////////////////////////
2269 /// constructor
2270 
2271 TGeoIdentity::TGeoIdentity(const char *name)
2272  :TGeoMatrix(name)
2273 {
2274  if (!gGeoIdentity) gGeoIdentity = this;
2275  RegisterYourself();
2276 }
2277 
2278 ////////////////////////////////////////////////////////////////////////////////
2279 /// Return a temporary inverse of this.
2280 
2281 TGeoHMatrix TGeoIdentity::Inverse() const
2282 {
2283  TGeoHMatrix h = *gGeoIdentity;
2284  return h;
2285 }
2286 
2287 /** \class TGeoHMatrix
2288 \ingroup Geometry_classes
2289 
2290 Matrix class used for computing global transformations
2291 Should NOT be used for node definition. An instance of this class
2292 is generally used to pile-up local transformations starting from
2293 the top level physical node, down to the current node.
2294 */
2295 
2296 ClassImp(TGeoHMatrix);
2297 
2298 ////////////////////////////////////////////////////////////////////////////////
2299 /// dummy ctor
2300 
2301 TGeoHMatrix::TGeoHMatrix()
2302 {
2303  memset(&fTranslation[0], 0, kN3);
2304  memcpy(fRotationMatrix,kIdentityMatrix,kN9);
2305  memcpy(fScale,kUnitScale,kN3);
2306 }
2307 
2308 ////////////////////////////////////////////////////////////////////////////////
2309 /// constructor
2310 
2311 TGeoHMatrix::TGeoHMatrix(const char* name)
2312  :TGeoMatrix(name)
2313 {
2314  memset(&fTranslation[0], 0, kN3);
2315  memcpy(fRotationMatrix,kIdentityMatrix,kN9);
2316  memcpy(fScale,kUnitScale,kN3);
2317 }
2318 
2319 ////////////////////////////////////////////////////////////////////////////////
2320 /// assignment
2321 
2322 TGeoHMatrix::TGeoHMatrix(const TGeoMatrix &matrix)
2323  :TGeoMatrix(matrix)
2324 {
2325  memset(&fTranslation[0], 0, kN3);
2326  memcpy(fRotationMatrix,kIdentityMatrix,kN9);
2327  memcpy(fScale,kUnitScale,kN3);
2328  if (matrix.IsIdentity()) return;
2329  if (matrix.IsTranslation())
2330  SetTranslation(matrix.GetTranslation());
2331  if (matrix.IsRotation())
2332  memcpy(fRotationMatrix,matrix.GetRotationMatrix(),kN9);
2333  if (matrix.IsScale())
2334  memcpy(fScale,matrix.GetScale(),kN3);
2335 }
2336 
2337 ////////////////////////////////////////////////////////////////////////////////
2338 /// destructor
2339 
2340 TGeoHMatrix::~TGeoHMatrix()
2341 {
2342 }
2343 
2344 ////////////////////////////////////////////////////////////////////////////////
2345 /// assignment
2346 
2347 TGeoHMatrix &TGeoHMatrix::operator=(const TGeoMatrix *matrix)
2348 {
2349  return TGeoHMatrix::operator=(*matrix);
2350 }
2351 
2352 ////////////////////////////////////////////////////////////////////////////////
2353 /// assignment
2354 
2355 TGeoHMatrix &TGeoHMatrix::operator=(const TGeoMatrix &matrix)
2356 {
2357  if (&matrix == this) return *this;
2358  Clear();
2359  Bool_t registered = TestBit(kGeoRegistered);
2360  TNamed::operator=(matrix);
2361  if (matrix.IsIdentity()) return *this;
2362  if (matrix.IsTranslation())
2363  memcpy(fTranslation,matrix.GetTranslation(),kN3);
2364  if (matrix.IsRotation())
2365  memcpy(fRotationMatrix,matrix.GetRotationMatrix(),kN9);
2366  if (matrix.IsScale())
2367  memcpy(fScale,matrix.GetScale(),kN3);
2368  SetBit(kGeoRegistered,registered);
2369  return *this;
2370 }
2371 
2372 ////////////////////////////////////////////////////////////////////////////////
2373 /// Composition
2374 
2375 TGeoHMatrix &TGeoHMatrix::operator*=(const TGeoMatrix &right)
2376 {
2377  Multiply(&right);
2378  return *this;
2379 }
2380 
2381 TGeoHMatrix TGeoHMatrix::operator*(const TGeoMatrix &right) const
2382 {
2383  TGeoHMatrix h = *this;
2384  h *= right;
2385  return h;
2386 }
2387 
2388 ////////////////////////////////////////////////////////////////////////////////
2389 /// Is-equal operator
2390 
2391 Bool_t TGeoHMatrix::operator==(const TGeoMatrix &other) const
2392 {
2393  if (&other == this) return kTRUE;
2394  const Double_t *tr = GetTranslation();
2395  const Double_t *otr = other.GetTranslation();
2396  for (auto i=0; i<3; i++) if (TMath::Abs(tr[i]-otr[i])>1.E-10) return kFALSE;
2397  const Double_t *rot = GetRotationMatrix();
2398  const Double_t *orot = other.GetRotationMatrix();
2399  for (auto i=0; i<9; i++) if (TMath::Abs(rot[i]-orot[i])>1.E-10) return kFALSE;
2400  const Double_t *scl = GetScale();
2401  const Double_t *oscl = other.GetScale();
2402  for (auto i=0; i<3; i++) if (TMath::Abs(scl[i]-oscl[i])>1.E-10) return kFALSE;
2403  return kTRUE;
2404 }
2405 
2406 ////////////////////////////////////////////////////////////////////////////////
2407 /// Fast copy method.
2408 
2409 void TGeoHMatrix::CopyFrom(const TGeoMatrix *other)
2410 {
2411  SetBit(kGeoTranslation, other->IsTranslation());
2412  SetBit(kGeoRotation, other->IsRotation());
2413  SetBit(kGeoReflection, other->IsReflection());
2414  memcpy(fTranslation,other->GetTranslation(),kN3);
2415  memcpy(fRotationMatrix,other->GetRotationMatrix(),kN9);
2416 }
2417 
2418 ////////////////////////////////////////////////////////////////////////////////
2419 /// clear the data for this matrix
2420 
2421 void TGeoHMatrix::Clear(Option_t *)
2422 {
2423  SetBit(kGeoReflection, kFALSE);
2424  if (IsIdentity()) return;
2425  ResetBit(kGeoTranslation);
2426  ResetBit(kGeoRotation);
2427  ResetBit(kGeoScale);
2428  memcpy(fTranslation,kNullVector,kN3);
2429  memcpy(fRotationMatrix,kIdentityMatrix,kN9);
2430  memcpy(fScale,kUnitScale,kN3);
2431 }
2432 
2433 ////////////////////////////////////////////////////////////////////////////////
2434 /// Make a clone of this matrix.
2435 
2436 TGeoMatrix *TGeoHMatrix::MakeClone() const
2437 {
2438  TGeoMatrix *matrix = new TGeoHMatrix(*this);
2439  return matrix;
2440 }
2441 
2442 ////////////////////////////////////////////////////////////////////////////////
2443 /// Perform a rotation about Z having the sine/cosine of the rotation angle.
2444 
2445 void TGeoHMatrix::FastRotZ(const Double_t *sincos)
2446 {
2447  fRotationMatrix[0] = sincos[1];
2448  fRotationMatrix[1] = -sincos[0];
2449  fRotationMatrix[3] = sincos[0];
2450  fRotationMatrix[4] = sincos[1];
2451  SetBit(kGeoRotation);
2452 }
2453 
2454 ////////////////////////////////////////////////////////////////////////////////
2455 /// Return a temporary inverse of this.
2456 
2457 TGeoHMatrix TGeoHMatrix::Inverse() const
2458 {
2459  TGeoHMatrix h;
2460  h = *this;
2461  if (IsTranslation()) {
2462  Double_t tr[3];
2463  tr[0] = -fTranslation[0]*fRotationMatrix[0] - fTranslation[1]*fRotationMatrix[3] - fTranslation[2]*fRotationMatrix[6];
2464  tr[1] = -fTranslation[0]*fRotationMatrix[1] - fTranslation[1]*fRotationMatrix[4] - fTranslation[2]*fRotationMatrix[7];
2465  tr[2] = -fTranslation[0]*fRotationMatrix[2] - fTranslation[1]*fRotationMatrix[5] - fTranslation[2]*fRotationMatrix[8];
2466  h.SetTranslation(tr);
2467  }
2468  if (IsRotation()) {
2469  Double_t newrot[9];
2470  newrot[0] = fRotationMatrix[0];
2471  newrot[1] = fRotationMatrix[3];
2472  newrot[2] = fRotationMatrix[6];
2473  newrot[3] = fRotationMatrix[1];
2474  newrot[4] = fRotationMatrix[4];
2475  newrot[5] = fRotationMatrix[7];
2476  newrot[6] = fRotationMatrix[2];
2477  newrot[7] = fRotationMatrix[5];
2478  newrot[8] = fRotationMatrix[8];
2479  h.SetRotation(newrot);
2480  }
2481  if (IsScale()) {
2482  Double_t sc[3];
2483  sc[0] = 1./fScale[0];
2484  sc[1] = 1./fScale[1];
2485  sc[2] = 1./fScale[2];
2486  h.SetScale(sc);
2487  }
2488  return h;
2489 }
2490 
2491 ////////////////////////////////////////////////////////////////////////////////
2492 /// computes determinant of the rotation matrix
2493 
2494 Double_t TGeoHMatrix::Determinant() const
2495 {
2496  Double_t
2497  det = fRotationMatrix[0]*fRotationMatrix[4]*fRotationMatrix[8] +
2498  fRotationMatrix[3]*fRotationMatrix[7]*fRotationMatrix[2] +
2499  fRotationMatrix[6]*fRotationMatrix[1]*fRotationMatrix[5] -
2500  fRotationMatrix[2]*fRotationMatrix[4]*fRotationMatrix[6] -
2501  fRotationMatrix[5]*fRotationMatrix[7]*fRotationMatrix[0] -
2502  fRotationMatrix[8]*fRotationMatrix[1]*fRotationMatrix[3];
2503  return det;
2504 }
2505 
2506 ////////////////////////////////////////////////////////////////////////////////
2507 /// multiply to the right with an other transformation
2508 /// if right is identity matrix, just return
2509 
2510 void TGeoHMatrix::Multiply(const TGeoMatrix *right)
2511 {
2512  if (right->IsIdentity()) return;
2513  const Double_t *r_tra = right->GetTranslation();
2514  const Double_t *r_rot = right->GetRotationMatrix();
2515  const Double_t *r_scl = right->GetScale();
2516  if (IsIdentity()) {
2517  if (right->IsRotation()) {
2518  SetBit(kGeoRotation);
2519  memcpy(fRotationMatrix,r_rot,kN9);
2520  if (right->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
2521  }
2522  if (right->IsScale()) {
2523  SetBit(kGeoScale);
2524  memcpy(fScale,r_scl,kN3);
2525  }
2526  if (right->IsTranslation()) {
2527  SetBit(kGeoTranslation);
2528  memcpy(fTranslation,r_tra,kN3);
2529  }
2530  return;
2531  }
2532  Int_t i, j;
2533  Double_t new_rot[9];
2534 
2535  if (right->IsRotation()) {
2536  SetBit(kGeoRotation);
2537  if (right->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
2538  }
2539  if (right->IsScale()) SetBit(kGeoScale);
2540  if (right->IsTranslation()) SetBit(kGeoTranslation);
2541 
2542  // new translation
2543  if (IsTranslation()) {
2544  for (i=0; i<3; i++) {
2545  fTranslation[i] += fRotationMatrix[3*i]*r_tra[0]
2546  + fRotationMatrix[3*i+1]*r_tra[1]
2547  + fRotationMatrix[3*i+2]*r_tra[2];
2548  }
2549  }
2550  if (IsRotation()) {
2551  // new rotation
2552  for (i=0; i<3; i++) {
2553  for (j=0; j<3; j++) {
2554  new_rot[3*i+j] = fRotationMatrix[3*i]*r_rot[j] +
2555  fRotationMatrix[3*i+1]*r_rot[3+j] +
2556  fRotationMatrix[3*i+2]*r_rot[6+j];
2557  }
2558  }
2559  memcpy(fRotationMatrix,new_rot,kN9);
2560  }
2561  // new scale
2562  if (IsScale()) {
2563  for (i=0; i<3; i++) fScale[i] *= r_scl[i];
2564  }
2565 }
2566 
2567 ////////////////////////////////////////////////////////////////////////////////
2568 /// multiply to the left with an other transformation
2569 /// if right is identity matrix, just return
2570 
2571 void TGeoHMatrix::MultiplyLeft(const TGeoMatrix *left)
2572 {
2573  if (left == gGeoIdentity) return;
2574  const Double_t *l_tra = left->GetTranslation();
2575  const Double_t *l_rot = left->GetRotationMatrix();
2576  const Double_t *l_scl = left->GetScale();
2577  if (IsIdentity()) {
2578  if (left->IsRotation()) {
2579  if (left->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
2580  SetBit(kGeoRotation);
2581  memcpy(fRotationMatrix,l_rot,kN9);
2582  }
2583  if (left->IsScale()) {
2584  SetBit(kGeoScale);
2585  memcpy(fScale,l_scl,kN3);
2586  }
2587  if (left->IsTranslation()) {
2588  SetBit(kGeoTranslation);
2589  memcpy(fTranslation,l_tra,kN3);
2590  }
2591  return;
2592  }
2593  Int_t i, j;
2594  Double_t new_tra[3];
2595  Double_t new_rot[9];
2596 
2597  if (left->IsRotation()) {
2598  SetBit(kGeoRotation);
2599  if (left->IsReflection()) SetBit(kGeoReflection, !TestBit(kGeoReflection));
2600  }
2601  if (left->IsScale()) SetBit(kGeoScale);
2602  if (left->IsTranslation()) SetBit(kGeoTranslation);
2603 
2604  // new translation
2605  if (IsTranslation()) {
2606  for (i=0; i<3; i++) {
2607  new_tra[i] = l_tra[i]
2608  + l_rot[3*i]* fTranslation[0]
2609  + l_rot[3*i+1]*fTranslation[1]
2610  + l_rot[3*i+2]*fTranslation[2];
2611  }
2612  memcpy(fTranslation,new_tra,kN3);
2613  }
2614  if (IsRotation()) {
2615  // new rotation
2616  for (i=0; i<3; i++) {
2617  for (j=0; j<3; j++) {
2618  new_rot[3*i+j] = l_rot[3*i]*fRotationMatrix[j] +
2619  l_rot[3*i+1]*fRotationMatrix[3+j] +
2620  l_rot[3*i+2]*fRotationMatrix[6+j];
2621  }
2622  }
2623  memcpy(fRotationMatrix,new_rot,kN9);
2624  }
2625  // new scale
2626  if (IsScale()) {
2627  for (i=0; i<3; i++) fScale[i] *= l_scl[i];
2628  }
2629 }
2630 
2631 ////////////////////////////////////////////////////////////////////////////////
2632 /// Rotate about X axis with angle expressed in degrees.
2633 
2634 void TGeoHMatrix::RotateX(Double_t angle)
2635 {
2636  SetBit(kGeoRotation);
2637  Double_t phi = angle*TMath::DegToRad();
2638  Double_t c = TMath::Cos(phi);
2639  Double_t s = TMath::Sin(phi);
2640  Double_t v[9];
2641  v[0] = fRotationMatrix[0];
2642  v[1] = fRotationMatrix[1];
2643  v[2] = fRotationMatrix[2];
2644  v[3] = c*fRotationMatrix[3]-s*fRotationMatrix[6];
2645  v[4] = c*fRotationMatrix[4]-s*fRotationMatrix[7];
2646  v[5] = c*fRotationMatrix[5]-s*fRotationMatrix[8];
2647  v[6] = s*fRotationMatrix[3]+c*fRotationMatrix[6];
2648  v[7] = s*fRotationMatrix[4]+c*fRotationMatrix[7];
2649  v[8] = s*fRotationMatrix[5]+c*fRotationMatrix[8];
2650  memcpy(fRotationMatrix, v, kN9);
2651 
2652  v[0] = fTranslation[0];
2653  v[1] = c*fTranslation[1]-s*fTranslation[2];
2654  v[2] = s*fTranslation[1]+c*fTranslation[2];
2655  memcpy(fTranslation,v,kN3);
2656 }
2657 
2658 ////////////////////////////////////////////////////////////////////////////////
2659 /// Rotate about Y axis with angle expressed in degrees.
2660 
2661 void TGeoHMatrix::RotateY(Double_t angle)
2662 {
2663  SetBit(kGeoRotation);
2664  Double_t phi = angle*TMath::DegToRad();
2665  Double_t c = TMath::Cos(phi);
2666  Double_t s = TMath::Sin(phi);
2667  Double_t v[9];
2668  v[0] = c*fRotationMatrix[0]+s*fRotationMatrix[6];
2669  v[1] = c*fRotationMatrix[1]+s*fRotationMatrix[7];
2670  v[2] = c*fRotationMatrix[2]+s*fRotationMatrix[8];
2671  v[3] = fRotationMatrix[3];
2672  v[4] = fRotationMatrix[4];
2673  v[5] = fRotationMatrix[5];
2674  v[6] = -s*fRotationMatrix[0]+c*fRotationMatrix[6];
2675  v[7] = -s*fRotationMatrix[1]+c*fRotationMatrix[7];
2676  v[8] = -s*fRotationMatrix[2]+c*fRotationMatrix[8];
2677  memcpy(fRotationMatrix, v, kN9);
2678 
2679  v[0] = c*fTranslation[0]+s*fTranslation[2];
2680  v[1] = fTranslation[1];
2681  v[2] = -s*fTranslation[0]+c*fTranslation[2];
2682  memcpy(fTranslation,v,kN3);
2683 }
2684 
2685 ////////////////////////////////////////////////////////////////////////////////
2686 /// Rotate about Z axis with angle expressed in degrees.
2687 
2688 void TGeoHMatrix::RotateZ(Double_t angle)
2689 {
2690  SetBit(kGeoRotation);
2691  Double_t phi = angle*TMath::DegToRad();
2692  Double_t c = TMath::Cos(phi);
2693  Double_t s = TMath::Sin(phi);
2694  Double_t v[9];
2695  v[0] = c*fRotationMatrix[0]-s*fRotationMatrix[3];
2696  v[1] = c*fRotationMatrix[1]-s*fRotationMatrix[4];
2697  v[2] = c*fRotationMatrix[2]-s*fRotationMatrix[5];
2698  v[3] = s*fRotationMatrix[0]+c*fRotationMatrix[3];
2699  v[4] = s*fRotationMatrix[1]+c*fRotationMatrix[4];
2700  v[5] = s*fRotationMatrix[2]+c*fRotationMatrix[5];
2701  v[6] = fRotationMatrix[6];
2702  v[7] = fRotationMatrix[7];
2703  v[8] = fRotationMatrix[8];
2704  memcpy(&fRotationMatrix[0],v,kN9);
2705 
2706  v[0] = c*fTranslation[0]-s*fTranslation[1];
2707  v[1] = s*fTranslation[0]+c*fTranslation[1];
2708  v[2] = fTranslation[2];
2709  memcpy(fTranslation,v,kN3);
2710 }
2711 
2712 ////////////////////////////////////////////////////////////////////////////////
2713 /// Multiply by a reflection respect to YZ.
2714 
2715 void TGeoHMatrix::ReflectX(Bool_t leftside, Bool_t rotonly)
2716 {
2717  if (leftside && !rotonly) fTranslation[0] = -fTranslation[0];
2718  if (leftside) {
2719  fRotationMatrix[0]=-fRotationMatrix[0];
2720  fRotationMatrix[1]=-fRotationMatrix[1];
2721  fRotationMatrix[2]=-fRotationMatrix[2];
2722  } else {
2723  fRotationMatrix[0]=-fRotationMatrix[0];
2724  fRotationMatrix[3]=-fRotationMatrix[3];
2725  fRotationMatrix[6]=-fRotationMatrix[6];
2726  }
2727  SetBit(kGeoRotation);
2728  SetBit(kGeoReflection, !IsReflection());
2729 }
2730 
2731 ////////////////////////////////////////////////////////////////////////////////
2732 /// Multiply by a reflection respect to ZX.
2733 
2734 void TGeoHMatrix::ReflectY(Bool_t leftside, Bool_t rotonly)
2735 {
2736  if (leftside && !rotonly) fTranslation[1] = -fTranslation[1];
2737  if (leftside) {
2738  fRotationMatrix[3]=-fRotationMatrix[3];
2739  fRotationMatrix[4]=-fRotationMatrix[4];
2740  fRotationMatrix[5]=-fRotationMatrix[5];
2741  } else {
2742  fRotationMatrix[1]=-fRotationMatrix[1];
2743  fRotationMatrix[4]=-fRotationMatrix[4];
2744  fRotationMatrix[7]=-fRotationMatrix[7];
2745  }
2746  SetBit(kGeoRotation);
2747  SetBit(kGeoReflection, !IsReflection());
2748 }
2749 
2750 ////////////////////////////////////////////////////////////////////////////////
2751 /// Multiply by a reflection respect to XY.
2752 
2753 void TGeoHMatrix::ReflectZ(Bool_t leftside, Bool_t rotonly)
2754 {
2755  if (leftside && !rotonly) fTranslation[2] = -fTranslation[2];
2756  if (leftside) {
2757  fRotationMatrix[6]=-fRotationMatrix[6];
2758  fRotationMatrix[7]=-fRotationMatrix[7];
2759  fRotationMatrix[8]=-fRotationMatrix[8];
2760  } else {
2761  fRotationMatrix[2]=-fRotationMatrix[2];
2762  fRotationMatrix[5]=-fRotationMatrix[5];
2763  fRotationMatrix[8]=-fRotationMatrix[8];
2764  }
2765  SetBit(kGeoRotation);
2766  SetBit(kGeoReflection, !IsReflection());
2767 }
2768 
2769 ////////////////////////////////////////////////////////////////////////////////
2770 /// Save a primitive as a C++ statement(s) on output stream "out".
2771 
2772 void TGeoHMatrix::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
2773 {
2774  if (TestBit(kGeoSavePrimitive)) return;
2775  const Double_t *tr = fTranslation;
2776  const Double_t *rot = fRotationMatrix;
2777  out << " // HMatrix: " << GetName() << std::endl;
2778  out << " tr[0] = " << tr[0] << "; " << "tr[1] = " << tr[1] << "; " << "tr[2] = " << tr[2] << ";" << std::endl;
2779  out << " rot[0] =" << rot[0] << "; " << "rot[1] = " << rot[1] << "; " << "rot[2] = " << rot[2] << ";" << std::endl;
2780  out << " rot[3] =" << rot[3] << "; " << "rot[4] = " << rot[4] << "; " << "rot[5] = " << rot[5] << ";" << std::endl;
2781  out << " rot[6] =" << rot[6] << "; " << "rot[7] = " << rot[7] << "; " << "rot[8] = " << rot[8] << ";" << std::endl;
2782  char *name = GetPointerName();
2783  out << " TGeoHMatrix *" << name << " = new TGeoHMatrix(\"" << GetName() << "\");" << std::endl;
2784  out << " " << name << "->SetTranslation(tr);" << std::endl;
2785  out << " " << name << "->SetRotation(rot);" << std::endl;
2786  if (IsTranslation()) out << " " << name << "->SetBit(TGeoMatrix::kGeoTranslation);" << std::endl;
2787  if (IsRotation()) out << " " << name << "->SetBit(TGeoMatrix::kGeoRotation);" << std::endl;
2788  if (IsReflection()) out << " " << name << "->SetBit(TGeoMatrix::kGeoReflection);" << std::endl;
2789  TObject::SetBit(kGeoSavePrimitive);
2790 }