DD4hep  1.28.0
Detector Description Toolkit for High Energy Physics
MatrixHelpers.cpp
Go to the documentation of this file.
1 //==========================================================================
2 // AIDA Detector description implementation
3 //--------------------------------------------------------------------------
4 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
5 // All rights reserved.
6 //
7 // For the licensing terms see $DD4hepINSTALL/LICENSE.
8 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
9 //
10 // Author : M.Frank
11 //
12 //==========================================================================
13 
14 // Framework include files
15 #include <DD4hep/MatrixHelpers.h>
16 #include <DD4hep/DD4hepUnits.h>
17 
18 #ifdef DD4HEP_USE_GEANT4_UNITS
19 #define MM_2_CM 1.0
20 #else
21 #define MM_2_CM 0.1
22 #endif
23 
24 // ROOT includes
25 #include <TGeoMatrix.h>
26 
27 TGeoIdentity* dd4hep::detail::matrix::_identity() {
28  return gGeoIdentity;
29 }
30 
31 ROOT::Math::XYZVector dd4hep::detail::matrix::_scale(const TGeoMatrix* matrix) {
32  if ( matrix->IsScale() ) {
33  const Double_t* mat_scale = matrix->GetScale();
34  return ROOT::Math::XYZVector(mat_scale[0],mat_scale[1],mat_scale[2]);
35  }
36  return ROOT::Math::XYZVector(1e0,1e0,1e0);
37 }
38 
39 ROOT::Math::XYZVector dd4hep::detail::matrix::_scale(const TGeoMatrix& matrix) {
40  return _scale(&matrix);
41 }
42 
43 dd4hep::Position dd4hep::detail::matrix::_translation(const TGeoMatrix* matrix) {
44  if ( matrix->IsTranslation() ) {
45  const Double_t* trans = matrix->GetTranslation();
46  return Position(trans[0]*MM_2_CM,trans[1]*MM_2_CM,trans[2]*MM_2_CM);
47  }
48  return Position(0e0,0e0,0e0);
49 }
50 
51 dd4hep::Position dd4hep::detail::matrix::_translation(const TGeoMatrix& matrix) {
52  return _translation(&matrix);
53 }
54 
55 TGeoTranslation* dd4hep::detail::matrix::_translation(const Position& pos) {
56  return new TGeoTranslation("", pos.X(), pos.Y(), pos.Z());
57 }
58 
59 TGeoRotation* dd4hep::detail::matrix::_rotationZYX(const RotationZYX& rot) {
60  return new TGeoRotation("", rot.Phi() * RAD_2_DEGREE, rot.Theta() * RAD_2_DEGREE, rot.Psi() * RAD_2_DEGREE);
61 }
62 
63 TGeoRotation* dd4hep::detail::matrix::_rotation3D(const Rotation3D& rot3D) {
64  EulerAngles rot(rot3D);
65  return new TGeoRotation("", rot.Phi() * RAD_2_DEGREE, rot.Theta() * RAD_2_DEGREE, rot.Psi() * RAD_2_DEGREE);
66 }
67 
69 dd4hep::Rotation3D dd4hep::detail::matrix::_rotation3D(const TGeoMatrix* matrix) {
70  if ( matrix->IsRotation() ) {
71  const Double_t* rot = matrix->GetRotationMatrix();
72  return Rotation3D(rot[0],rot[1],rot[2],
73  rot[3],rot[4],rot[5],
74  rot[6],rot[7],rot[8]);
75  }
76  return Rotation3D(0e0,0e0,0e0,
77  0e0,0e0,0e0,
78  0e0,0e0,0e0);
79 }
80 
82 dd4hep::Rotation3D dd4hep::detail::matrix::_rotation3D(const TGeoMatrix& matrix) {
83  return _rotation3D(&matrix);
84 }
85 
87 TGeoHMatrix& dd4hep::detail::matrix::_transform(TGeoHMatrix& tr, const RotationZYX& rot) {
88  tr.RotateZ(rot.Phi() * RAD_2_DEGREE);
89  tr.RotateY(rot.Theta() * RAD_2_DEGREE);
90  tr.RotateX(rot.Psi() * RAD_2_DEGREE);
91  return tr;
92 }
93 
95 TGeoHMatrix& dd4hep::detail::matrix::_transform(TGeoHMatrix& tr, const Position& pos) {
96  double t[3];
97  pos.GetCoordinates(t);
98  tr.SetDx(t[0]);
99  tr.SetDy(t[1]);
100  tr.SetDz(t[2]);
101  return tr;
102 }
103 
105 TGeoHMatrix& dd4hep::detail::matrix::_transform(TGeoHMatrix& tr, const Rotation3D& rot) {
106  Double_t* r = tr.GetRotationMatrix();
107  rot.GetComponents(r);
108  return tr;
109 }
110 
112 TGeoHMatrix& dd4hep::detail::matrix::_transform(TGeoHMatrix& tr, const Transform3D& trans) {
113  Position pos;
114  RotationZYX rot;
115  trans.GetDecomposition(rot, pos);
116  return _transform(tr, pos, rot);
117 }
118 
120 TGeoHMatrix& dd4hep::detail::matrix::_transform(TGeoHMatrix& tr, const Position& pos, const RotationZYX& rot) {
121  return _transform(_transform(tr, rot), pos);
122 }
123 
125 TGeoHMatrix* dd4hep::detail::matrix::_transform(const Position& pos) {
126  return &_transform(*(new TGeoHMatrix()), pos);
127 }
128 
130 TGeoHMatrix* dd4hep::detail::matrix::_transform(const RotationZYX& rot) {
131  return &_transform(*(new TGeoHMatrix()), rot);
132 }
133 
135 TGeoHMatrix* dd4hep::detail::matrix::_transform(const Rotation3D& rot) {
136  return &_transform(*(new TGeoHMatrix()), rot);
137 }
138 
140 TGeoHMatrix* dd4hep::detail::matrix::_transform(const Transform3D& trans) {
141  return &_transform(*(new TGeoHMatrix()), trans);
142 }
143 
145 TGeoHMatrix* dd4hep::detail::matrix::_transform(const Position& pos, const RotationZYX& rot) {
146  return &_transform(*(new TGeoHMatrix()), pos, rot);
147 }
148 
150 dd4hep::Transform3D dd4hep::detail::matrix::_transform(const TGeoMatrix* matrix) {
151  const Double_t* t = matrix->GetTranslation();
152  if ( matrix->IsRotation() ) {
153  const Double_t* rot = matrix->GetRotationMatrix();
154  return Transform3D(rot[0],rot[1],rot[2],t[0]*MM_2_CM,
155  rot[3],rot[4],rot[5],t[1]*MM_2_CM,
156  rot[6],rot[7],rot[8],t[2]*MM_2_CM);
157  }
158  return Transform3D(0e0,0e0,0e0,t[0]*MM_2_CM,
159  0e0,0e0,0e0,t[1]*MM_2_CM,
160  0e0,0e0,0e0,t[2]*MM_2_CM);
161 }
162 
163 // add another implementation that takes const reference
164 dd4hep::Transform3D dd4hep::detail::matrix::_transform(const TGeoMatrix& matrix) {
165  const Double_t* t = matrix.GetTranslation();
166  if ( matrix.IsRotation() ) {
167  const Double_t* r = matrix.GetRotationMatrix();
168  return Transform3D(r[0],r[1],r[2],t[0]*MM_2_CM,
169  r[3],r[4],r[5],t[1]*MM_2_CM,
170  r[6],r[7],r[8],t[2]*MM_2_CM);
171  }
172  return Transform3D(0e0,0e0,0e0,t[0]*MM_2_CM,
173  0e0,0e0,0e0,t[1]*MM_2_CM,
174  0e0,0e0,0e0,t[2]*MM_2_CM);
175 }
176 
178 void dd4hep::detail::matrix::_transform(const TGeoMatrix& matrix, Transform3D& tr) {
179  tr = _transform(matrix);
180 }
181 
183 void dd4hep::detail::matrix::_transform(const TGeoMatrix* matrix, Transform3D& tr) {
184  if ( matrix ) {
185  _transform(*matrix, tr);
186  return;
187  }
188  tr = Transform3D();
189 }
190 
191 dd4hep::XYZAngles dd4hep::detail::matrix::_xyzAngles(const TGeoMatrix* matrix) {
192  return matrix->IsRotation() ? _xyzAngles(matrix->GetRotationMatrix()) : XYZAngles(0,0,0);
193 }
194 
195 dd4hep::XYZAngles dd4hep::detail::matrix::_xyzAngles(const double* rot) {
196  Double_t cosb = std::sqrt(rot[0]*rot[0] + rot[1]*rot[1]);
197  if (cosb > 0.00001) {
198  return XYZAngles(atan2(rot[5], rot[8]), atan2(-rot[2], cosb), atan2(rot[1], rot[0]));
199  }
200  return XYZAngles(atan2(-rot[7], rot[4]),atan2(-rot[2], cosb),0);
201 }
202 
203 void dd4hep::detail::matrix::_decompose(const TGeoMatrix& trafo, Position& pos, Rotation3D& rot) {
204  _decompose(_transform(trafo), pos, rot);
205 }
206 
207 void dd4hep::detail::matrix::_decompose(const TGeoMatrix& trafo, Position& pos, RotationZYX& rot) {
208  _decompose(_transform(trafo), pos, rot);
209 }
210 
211 void dd4hep::detail::matrix::_decompose(const TGeoMatrix& trafo, Position& pos, XYZAngles& rot) {
212  _decompose(_transform(trafo), pos, rot);
213 }
214 
215 void dd4hep::detail::matrix::_decompose(const Transform3D& trafo, Position& pos, Rotation3D& rot) {
216  trafo.GetDecomposition(rot, pos);
217 }
218 
219 void dd4hep::detail::matrix::_decompose(const Rotation3D& rot, Position& x, Position& y, Position& z) {
220  rot.GetComponents(x,y,z);
221 }
222 
223 void dd4hep::detail::matrix::_decompose(const Transform3D& trafo, Translation3D& pos, RotationZYX& rot) {
224  trafo.GetDecomposition(rot, pos);
225 }
226 
227 void dd4hep::detail::matrix::_decompose(const Transform3D& trafo, Translation3D& pos, XYZAngles& rot) {
228  EulerAngles r;
229  trafo.GetDecomposition(r,pos);
230  rot.SetXYZ(r.Psi(),r.Theta(),r.Phi());
231 }
232 
233 void dd4hep::detail::matrix::_decompose(const Transform3D& trafo, Position& pos, RotationZYX& rot) {
234  trafo.GetDecomposition(rot,pos);
235 }
236 
237 void dd4hep::detail::matrix::_decompose(const Transform3D& trafo, Position& pos, XYZAngles& rot) {
238  EulerAngles r;
239  trafo.GetDecomposition(r,pos);
240  rot.SetXYZ(r.Psi(),r.Theta(),r.Phi());
241 }
242 
244 double dd4hep::detail::matrix::determinant(const TGeoMatrix* matrix) {
245  return (matrix) ? determinant(*matrix) : 0e0;
246 }
247 
249 double dd4hep::detail::matrix::determinant(const TGeoMatrix& matrix) {
250  const double* r = matrix.GetRotationMatrix();
251  if ( r ) {
252  double det =
253  r[0]*r[4]*r[8] + r[3]*r[7]*r[2] + r[6]*r[1]*r[5] -
254  r[2]*r[4]*r[6] - r[5]*r[7]*r[0] - r[8]*r[1]*r[3];
255  return det;
256  }
257  return 0.0;
258 }
259 
261 double dd4hep::detail::matrix::determinant(const Transform3D& tr) {
262  Position p;
263  Rotation3D r;
264  tr.GetDecomposition(r, p);
265  return determinant(r);
266 }
267 
269 double dd4hep::detail::matrix::determinant(const Rotation3D& rot) {
270  Position x, y, z;
271  rot.GetComponents(x,y,z);
272  double det = (x.Cross(y)).Dot(z);
273  return det;
274 }
275 
277 int dd4hep::detail::matrix::_matrixEqual(const TGeoMatrix& left, const TGeoMatrix& right) {
278  static constexpr double epsilon = 1e-12;
279  int result = MATRICES_EQUAL;
280  const Double_t* t1 = left.GetTranslation();
281  const Double_t* t2 = right.GetTranslation();
282  for(int i=0; i<3; ++i) {
283  if ( std::fabs(t1[i]-t2[i]) > epsilon ) {
284  result = MATRICES_DIFFER_TRANSLATION;
285  break;
286  }
287  }
288  const Double_t* r1 = left.GetRotationMatrix();
289  const Double_t* r2 = right.GetRotationMatrix();
290  for(int i=0; i<9; ++i) {
291  if ( std::fabs(r1[i]-r2[i]) > epsilon ) {
292  result |= MATRICES_DIFFER_ROTATION;
293  break;
294  }
295  }
296  return result;
297 }
298 
299 
dd4hep::detail::matrix::_translation
TGeoTranslation * _translation(const Position &pos)
Convert a Position object to a TGeoTranslation.
Definition: MatrixHelpers.cpp:55
dd4hep::detail::matrix::_decompose
void _decompose(const TGeoMatrix &matrix, Position &pos, Rotation3D &rot)
Decompose a generic ROOT Matrix into a translation (Position) and a Rotation3D.
Definition: MatrixHelpers.cpp:203
dd4hep::detail::matrix::_rotationZYX
TGeoRotation * _rotationZYX(const RotationZYX &rot)
Convert a RotationZYX object to a newly created TGeoRotation.
Definition: MatrixHelpers.cpp:59
MatrixHelpers.h
dd4hep::XYZAngles
ROOT::Math::XYZVector XYZAngles
Definition: Objects.h:83
dd4hep::detail::matrix::_scale
ROOT::Math::XYZVector _scale(const TGeoMatrix *matrix)
Extract the scale part of a TGeoMatrix.
Definition: MatrixHelpers.cpp:31
dd4hep::detail::matrix::_identity
TGeoIdentity * _identity()
Access the TGeo identity transformation.
Definition: MatrixHelpers.cpp:27
dd4hep::detail::matrix::_transform
TGeoHMatrix * _transform(const Transform3D &trans)
Convert a Transform3D object to a newly created TGeoHMatrix.
Definition: MatrixHelpers.cpp:140
dd4hep::detail::matrix::_rotation3D
TGeoRotation * _rotation3D(const Rotation3D &rot)
Convert a Rotation3D object to a newly createdTGeoRotation.
Definition: MatrixHelpers.cpp:63
dd4hep::Rotation3D
ROOT::Math::Rotation3D Rotation3D
Definition: Objects.h:113
epsilon
const double epsilon
Definition: test_cellid_position_converter.cpp:42
dd4hep::Translation3D
ROOT::Math::Translation3D Translation3D
Definition: Objects.h:119
dd4hep::detail::matrix::MATRICES_EQUAL
@ MATRICES_EQUAL
Definition: MatrixHelpers.h:121
dd4hep::detail::matrix::MATRICES_DIFFER_ROTATION
@ MATRICES_DIFFER_ROTATION
Definition: MatrixHelpers.h:123
dd4hep::detail::matrix::_xyzAngles
XYZAngles _xyzAngles(const double *matrix)
Convert a 3x3 rotation matrix to XYZAngles.
Definition: MatrixHelpers.cpp:195
dd4hep::Transform3D
ROOT::Math::Transform3D Transform3D
Definition: Objects.h:117
dd4hep::Position
ROOT::Math::XYZVector Position
Definition: Objects.h:81
dd4hep::EulerAngles
ROOT::Math::EulerAngles EulerAngles
Definition: Objects.h:115
dd4hep::detail::matrix::MATRICES_DIFFER_TRANSLATION
@ MATRICES_DIFFER_TRANSLATION
Definition: MatrixHelpers.h:122
det
DetElement::Object * det
Definition: AlignmentsCalculator.cpp:66
dd4hep::detail::matrix::_matrixEqual
int _matrixEqual(const TGeoMatrix &left, const TGeoMatrix &right)
Check matrices for equality.
Definition: MatrixHelpers.cpp:277
MM_2_CM
#define MM_2_CM
Definition: MatrixHelpers.cpp:21
dd4hep::RotationZYX
ROOT::Math::RotationZYX RotationZYX
Definition: Objects.h:105
dd4hep::detail::matrix::determinant
double determinant(const TGeoMatrix *matrix)
Access determinant of rotation component (0 if no rotation present)
Definition: MatrixHelpers.cpp:244
DD4hepUnits.h
RAD_2_DEGREE
#define RAD_2_DEGREE
Definition: Handle.h:26
Generated on Thu May 16 2024 15:36:32 for DD4hep by   doxygen 1.8.18