| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4AffineTransformation Inline implementation
27 //
28 // --------------------------------------------------------------------
29
30 inline G4AffineTransform::G4AffineTransform()
31 : rxx(1),rxy(0),rxz(0),
32 ryx(0),ryy(1),ryz(0),
33 rzx(0),rzy(0),rzz(1),
34 tx(0),ty(0),tz(0)
35 {}
36
37 inline G4AffineTransform::G4AffineTransform(const G4ThreeVector& tlate)
38 : rxx(1),rxy(0),rxz(0),
39 ryx(0),ryy(1),ryz(0),
40 rzx(0),rzy(0),rzz(1),
41 tx(tlate.x()),ty(tlate.y()),tz(tlate.z())
42 {}
43
44 inline G4AffineTransform::G4AffineTransform(const G4RotationMatrix& rot)
45 : rxx(rot.xx()),rxy(rot.xy()),rxz(rot.xz()),
46 ryx(rot.yx()),ryy(rot.yy()),ryz(rot.yz()),
47 rzx(rot.zx()),rzy(rot.zy()),rzz(rot.zz()),
48 tx(0),ty(0),tz(0)
49 {}
50
51 inline G4AffineTransform::G4AffineTransform( const G4RotationMatrix& rot,
52 const G4ThreeVector& tlate )
53 : rxx(rot.xx()),rxy(rot.xy()),rxz(rot.xz()),
54 ryx(rot.yx()),ryy(rot.yy()),ryz(rot.yz()),
55 rzx(rot.zx()),rzy(rot.zy()),rzz(rot.zz()),
56 tx(tlate.x()),ty(tlate.y()),tz(tlate.z())
57 {}
58
59 inline G4AffineTransform::G4AffineTransform( const G4RotationMatrix* rot,
60 const G4ThreeVector& tlate )
61 : tx(tlate.x()),ty(tlate.y()),tz(tlate.z())
62 {
63 if (rot != nullptr)
64 {
65 rxx=rot->xx();rxy=rot->xy();rxz=rot->xz();
66 ryx=rot->yx();ryy=rot->yy();ryz=rot->yz();
67 rzx=rot->zx();rzy=rot->zy();rzz=rot->zz();
68 }
69 else
70 {
71 rxx=1; rxy=0; rxz=0;
72 ryx=0; ryy=1; ryz=0;
73 rzx=0; rzy=0; rzz=1;
74 }
75 }
76
77 inline
78 G4AffineTransform::
79 G4AffineTransform( const G4double prxx,const G4double prxy,const G4double prxz,
80 const G4double pryx,const G4double pryy,const G4double pryz,
81 const G4double przx,const G4double przy,const G4double przz,
82 const G4double ptx,const G4double pty,const G4double ptz )
83 : rxx(prxx),rxy(prxy),rxz(prxz),
84 ryx(pryx),ryy(pryy),ryz(pryz),
85 rzx(przx),rzy(przy),rzz(przz),
86 tx(ptx),ty(pty),tz(ptz)
87 {}
88
89 inline G4AffineTransform&
90 G4AffineTransform::operator = (const G4AffineTransform& rhs)
91 {
92 if (this == &rhs) { return *this; }
93 rxx = rhs.rxx; rxy = rhs.rxy; rxz = rhs.rxz;
94 ryx = rhs.ryx; ryy = rhs.ryy; ryz = rhs.ryz;
95 rzx = rhs.rzx; rzy = rhs.rzy; rzz = rhs.rzz;
96 tx = rhs.tx; ty = rhs.ty; tz = rhs.tz;
97 return *this;
98 }
99
100 inline G4AffineTransform
101 G4AffineTransform::operator * (const G4AffineTransform& tf) const
102 {
103 return G4AffineTransform(rxx*tf.rxx+rxy*tf.ryx+rxz*tf.rzx,
104 rxx*tf.rxy+rxy*tf.ryy+rxz*tf.rzy,
105 rxx*tf.rxz+rxy*tf.ryz+rxz*tf.rzz,
106
107 ryx*tf.rxx+ryy*tf.ryx+ryz*tf.rzx,
108 ryx*tf.rxy+ryy*tf.ryy+ryz*tf.rzy,
109 ryx*tf.rxz+ryy*tf.ryz+ryz*tf.rzz,
110
111 rzx*tf.rxx+rzy*tf.ryx+rzz*tf.rzx,
112 rzx*tf.rxy+rzy*tf.ryy+rzz*tf.rzy,
113 rzx*tf.rxz+rzy*tf.ryz+rzz*tf.rzz,
114
115 tx*tf.rxx+ty*tf.ryx+tz*tf.rzx+tf.tx,
116 tx*tf.rxy+ty*tf.ryy+tz*tf.rzy+tf.ty,
117 tx*tf.rxz+ty*tf.ryz+tz*tf.rzz+tf.tz);
118 }
119
120 inline G4AffineTransform&
121 G4AffineTransform::operator *= (const G4AffineTransform& tf)
122 {
123 // Use temporaries for `in place' compound transform computation
124 //
125 G4double nrxx=rxx*tf.rxx+rxy*tf.ryx+rxz*tf.rzx;
126 G4double nrxy=rxx*tf.rxy+rxy*tf.ryy+rxz*tf.rzy;
127 G4double nrxz=rxx*tf.rxz+rxy*tf.ryz+rxz*tf.rzz;
128
129 G4double nryx=ryx*tf.rxx+ryy*tf.ryx+ryz*tf.rzx;
130 G4double nryy=ryx*tf.rxy+ryy*tf.ryy+ryz*tf.rzy;
131 G4double nryz=ryx*tf.rxz+ryy*tf.ryz+ryz*tf.rzz;
132
133 G4double nrzx=rzx*tf.rxx+rzy*tf.ryx+rzz*tf.rzx;
134 G4double nrzy=rzx*tf.rxy+rzy*tf.ryy+rzz*tf.rzy;
135 G4double nrzz=rzx*tf.rxz+rzy*tf.ryz+rzz*tf.rzz;
136
137 G4double ntx=tx*tf.rxx+ty*tf.ryx+tz*tf.rzx+tf.tx;
138 G4double nty=tx*tf.rxy+ty*tf.ryy+tz*tf.rzy+tf.ty;
139 G4double ntz=tx*tf.rxz+ty*tf.ryz+tz*tf.rzz+tf.tz;
140
141 tx=ntx; ty=nty; tz=ntz;
142 rxx=nrxx; rxy=nrxy; rxz=nrxz;
143 ryx=nryx; ryy=nryy; ryz=nryz;
144 rzx=nrzx; rzy=nrzy; rzz=nrzz;
145
146 return *this;
147 }
148
149 inline G4AffineTransform&
150 G4AffineTransform::Product( const G4AffineTransform& tf1,
151 const G4AffineTransform& tf2 )
152 {
153 rxx = tf1.rxx*tf2.rxx + tf1.rxy*tf2.ryx + tf1.rxz*tf2.rzx;
154 rxy = tf1.rxx*tf2.rxy + tf1.rxy*tf2.ryy + tf1.rxz*tf2.rzy;
155 rxz = tf1.rxx*tf2.rxz + tf1.rxy*tf2.ryz + tf1.rxz*tf2.rzz;
156
157 ryx = tf1.ryx*tf2.rxx + tf1.ryy*tf2.ryx + tf1.ryz*tf2.rzx;
158 ryy = tf1.ryx*tf2.rxy + tf1.ryy*tf2.ryy + tf1.ryz*tf2.rzy;
159 ryz = tf1.ryx*tf2.rxz + tf1.ryy*tf2.ryz + tf1.ryz*tf2.rzz;
160
161 rzx = tf1.rzx*tf2.rxx + tf1.rzy*tf2.ryx + tf1.rzz*tf2.rzx;
162 rzy = tf1.rzx*tf2.rxy + tf1.rzy*tf2.ryy + tf1.rzz*tf2.rzy;
163 rzz = tf1.rzx*tf2.rxz + tf1.rzy*tf2.ryz + tf1.rzz*tf2.rzz;
164
165 tx = tf1.tx*tf2.rxx + tf1.ty*tf2.ryx + tf1.tz*tf2.rzx + tf2.tx;
166 ty = tf1.tx*tf2.rxy + tf1.ty*tf2.ryy + tf1.tz*tf2.rzy + tf2.ty;
167 tz = tf1.tx*tf2.rxz + tf1.ty*tf2.ryz + tf1.tz*tf2.rzz + tf2.tz;
168
169 return *this;
170 }
171
172 inline G4AffineTransform&
173 G4AffineTransform::InverseProduct( const G4AffineTransform& tf1,
174 const G4AffineTransform& tf2 )
175 {
176 if ((tf2.rxx + tf2.ryy + tf2.rzz) == 3.)
177 {
178 rxx = tf1.rxx;
179 rxy = tf1.rxy;
180 rxz = tf1.rxz;
181
182 ryx = tf1.ryx;
183 ryy = tf1.ryy;
184 ryz = tf1.ryz;
185
186 rzx = tf1.rzx;
187 rzy = tf1.rzy;
188 rzz = tf1.rzz;
189
190 tx = tf1.tx - tf2.tx;
191 ty = tf1.ty - tf2.ty;
192 tz = tf1.tz - tf2.tz;
193 }
194 else
195 {
196 G4double tf1rxx = tf1.rxx, tf1rxy = tf1.rxy, tf1rxz = tf1.rxz;
197 rxx = tf1rxx*tf2.rxx + tf1rxy*tf2.rxy + tf1rxz*tf2.rxz;
198 rxy = tf1rxx*tf2.ryx + tf1rxy*tf2.ryy + tf1rxz*tf2.ryz;
199 rxz = tf1rxx*tf2.rzx + tf1rxy*tf2.rzy + tf1rxz*tf2.rzz;
200
201 G4double tf1ryx = tf1.ryx, tf1ryy = tf1.ryy, tf1ryz = tf1.ryz;
202 ryx = tf1ryx*tf2.rxx + tf1ryy*tf2.rxy + tf1ryz*tf2.rxz;
203 ryy = tf1ryx*tf2.ryx + tf1ryy*tf2.ryy + tf1ryz*tf2.ryz;
204 ryz = tf1ryx*tf2.rzx + tf1ryy*tf2.rzy + tf1ryz*tf2.rzz;
205
206 G4double tf1rzx = tf1.rzx, tf1rzy = tf1.rzy, tf1rzz = tf1.rzz;
207 rzx = tf1rzx*tf2.rxx + tf1rzy*tf2.rxy + tf1rzz*tf2.rxz;
208 rzy = tf1rzx*tf2.ryx + tf1rzy*tf2.ryy + tf1rzz*tf2.ryz;
209 rzz = tf1rzx*tf2.rzx + tf1rzy*tf2.rzy + tf1rzz*tf2.rzz;
210
211 G4double tf1_2tx = tf1.tx - tf2.tx;
212 G4double tf1_2ty = tf1.ty - tf2.ty;
213 G4double tf1_2tz = tf1.tz - tf2.tz;
214 tx = tf1_2tx*tf2.rxx + tf1_2ty*tf2.rxy + tf1_2tz*tf2.rxz;
215 ty = tf1_2tx*tf2.ryx + tf1_2ty*tf2.ryy + tf1_2tz*tf2.ryz;
216 tz = tf1_2tx*tf2.rzx + tf1_2ty*tf2.rzy + tf1_2tz*tf2.rzz;
217 }
218 return *this;
219 }
220
221 inline G4ThreeVector
222 G4AffineTransform::TransformPoint(const G4ThreeVector& vec) const
223 {
224 G4double vecx = vec.x(), vecy = vec.y(), vecz = vec.z();
225 return { vecx*rxx + vecy*ryx + vecz*rzx + tx,
226 vecx*rxy + vecy*ryy + vecz*rzy + ty,
227 vecx*rxz + vecy*ryz + vecz*rzz + tz };
228 }
229
230 inline G4ThreeVector
231 G4AffineTransform::InverseTransformPoint(const G4ThreeVector& vec) const
232 {
233 G4double vecx = vec.x()-tx, vecy = vec.y()-ty, vecz = vec.z()-tz;
234 return { vecx*rxx + vecy*rxy + vecz*rxz,
235 vecx*ryx + vecy*ryy + vecz*ryz,
236 vecx*rzx + vecy*rzy + vecz*rzz };
237 }
238
239 inline G4ThreeVector
240 G4AffineTransform::TransformAxis(const G4ThreeVector& axis) const
241 {
242 G4double axisx = axis.x(), axisy = axis.y(), axisz = axis.z();
243 return { axisx*rxx + axisy*ryx + axisz*rzx,
244 axisx*rxy + axisy*ryy + axisz*rzy,
245 axisx*rxz + axisy*ryz + axisz*rzz };
246 }
247
248 inline G4ThreeVector
249 G4AffineTransform::InverseTransformAxis(const G4ThreeVector& axis) const
250 {
251 G4double axisx = axis.x(), axisy = axis.y(), axisz = axis.z();
252 return { axisx*rxx + axisy*rxy + axisz*rxz,
253 axisx*ryx + axisy*ryy + axisz*ryz,
254 axisx*rzx + axisy*rzy + axisz*rzz };
255 }
256
257 inline
258 void G4AffineTransform::ApplyPointTransform(G4ThreeVector& vec) const
259 {
260 G4double vecx = vec.x(), vecy = vec.y(), vecz = vec.z();
261 vec.setX( vecx*rxx + vecy*ryx + vecz*rzx + tx );
262 vec.setY( vecx*rxy + vecy*ryy + vecz*rzy + ty );
263 vec.setZ( vecx*rxz + vecy*ryz + vecz*rzz + tz );
264 }
265
266 inline
267 void G4AffineTransform::ApplyAxisTransform(G4ThreeVector& axis) const
268 {
269 G4double axisx = axis.x(), axisy = axis.y(), axisz = axis.z();
270 axis.setX( axisx*rxx + axisy*ryx + axisz*rzx );
271 axis.setY( axisx*rxy + axisy*ryy + axisz*rzy );
272 axis.setZ( axisx*rxz + axisy*ryz + axisz*rzz );
273 }
274
275 inline
276 G4AffineTransform G4AffineTransform::Inverse() const
277 {
278 G4double ttx = -tx, tty = -ty, ttz = -tz;
279 return G4AffineTransform( rxx, ryx, rzx,
280 rxy, ryy, rzy,
281 rxz, ryz, rzz,
282 ttx*rxx + tty*rxy + ttz*rxz,
283 ttx*ryx + tty*ryy + ttz*ryz,
284 ttx*rzx + tty*rzy + ttz*rzz );
285 }
286
287 inline
288 G4AffineTransform& G4AffineTransform::Invert()
289 {
290 G4double ttx = -tx, tty = -ty, ttz = -tz;
291 tx = ttx*rxx + tty*rxy + ttz*rxz;
292 ty = ttx*ryx + tty*ryy + ttz*ryz;
293 tz = ttx*rzx + tty*rzy + ttz*rzz;
294
295 G4double tmp1=ryx; ryx=rxy; rxy=tmp1;
296 G4double tmp2=rzx; rzx=rxz; rxz=tmp2;
297 G4double tmp3=rzy; rzy=ryz; ryz=tmp3;
298
299 return *this;
300 }
301
302 inline
303 G4AffineTransform& G4AffineTransform::operator +=(const G4ThreeVector& tlate)
304 {
305 tx += tlate.x();
306 ty += tlate.y();
307 tz += tlate.z();
308
309 return *this;
310 }
311
312 inline
313 G4AffineTransform& G4AffineTransform::operator -=(const G4ThreeVector& tlate)
314 {
315 tx -= tlate.x();
316 ty -= tlate.y();
317 tz -= tlate.z();
318
319 return *this;
320 }
321
322 inline
323 G4bool G4AffineTransform::operator == (const G4AffineTransform& tf) const
324 {
325 return tx==tf.tx&&ty==tf.ty&&tz==tf.tz&&
326 rxx==tf.rxx&&rxy==tf.rxy&&rxz==tf.rxz&&
327 ryx==tf.ryx&&ryy==tf.ryy&&ryz==tf.ryz&&
328 rzx==tf.rzx&&rzy==tf.rzy&&rzz==tf.rzz;
329 }
330
331 inline
332 G4bool G4AffineTransform::operator != (const G4AffineTransform& tf) const
333 {
334 return tx!=tf.tx||ty!=tf.ty||tz!=tf.tz||
335 rxx!=tf.rxx||rxy!=tf.rxy||rxz!=tf.rxz||
336 ryx!=tf.ryx||ryy!=tf.ryy||ryz!=tf.ryz||
337 rzx!=tf.rzx||rzy!=tf.rzy||rzz!=tf.rzz;
338 }
339
340 inline
341 G4double G4AffineTransform::operator [] (const G4int n) const
342 {
343 G4double v = 0.0;
344 switch(n)
345 {
346 case 0:
347 v=rxx;
348 break;
349 case 1:
350 v=rxy;
351 break;
352 case 2:
353 v=rxz;
354 break;
355 case 4:
356 v=ryx;
357 break;
358 case 5:
359 v=ryy;
360 break;
361 case 6:
362 v=ryz;
363 break;
364 case 8:
365 v=rzx;
366 break;
367 case 9:
368 v=rzy;
369 break;
370 case 10:
371 v=rzz;
372 break;
373 case 12:
374 v=tx;
375 break;
376 case 13:
377 v=ty;
378 break;
379 case 14:
380 v=tz;
381 break;
382 case 3:
383 case 7:
384 case 11:
385 break;
386 case 15:
387 v=1.0;
388 break;
389 }
390 return v;
391 }
392
393 inline
394 G4bool G4AffineTransform::IsRotated() const
395 {
396 return !(rxx==1.0 && ryy==1.0 && rzz==1.0);
397 }
398
399 inline
400 G4bool G4AffineTransform::IsTranslated() const
401 {
402 return (tx != 0.0) || (ty != 0.0) || (tz != 0.0);
403 }
404
405 inline G4RotationMatrix G4AffineTransform::NetRotation() const
406 {
407 return G4Rep3x3(rxx,rxy,rxz,
408 ryx,ryy,ryz,
409 rzx,rzy,rzz);
410 }
411
412 inline G4RotationMatrix G4AffineTransform::InverseNetRotation() const
413 {
414 return G4Rep3x3(rxx,ryx,rzx,
415 rxy,ryy,rzy,
416 rxz,ryz,rzz);
417 }
418
419 inline
420 G4ThreeVector G4AffineTransform::NetTranslation() const
421 {
422 return {tx,ty,tz};
423 }
424
425 inline
426 G4ThreeVector G4AffineTransform::InverseNetTranslation() const
427 {
428 G4double ttx = -tx, tty = -ty, ttz = -tz;
429 G4double invtx = ttx*rxx + tty*rxy + ttz*rxz;
430 G4double invty = ttx*ryx + tty*ryy + ttz*ryz;
431 G4double invtz = ttx*rzx + tty*rzy + ttz*rzz;
432 return {invtx,invty,invtz};
433 }
434
435 inline
436 void G4AffineTransform::SetNetRotation(const G4RotationMatrix& rot)
437 {
438 rxx=rot.xx();
439 rxy=rot.xy();
440 rxz=rot.xz();
441 ryx=rot.yx();
442 ryy=rot.yy();
443 ryz=rot.yz();
444 rzx=rot.zx();
445 rzy=rot.zy();
446 rzz=rot.zz();
447 }
448
449 inline
450 void G4AffineTransform::SetNetTranslation(const G4ThreeVector& tlate)
451 {
452 tx=tlate.x();
453 ty=tlate.y();
454 tz=tlate.z();
455 }
456
457 inline
458 G4AffineTransform::operator G4Transform3D () const
459 {
460 return G4Transform3D(NetRotation().inverse(),NetTranslation());
461 }
462
463 inline
464 std::ostream& operator << (std::ostream& os, const G4AffineTransform& transf)
465 {
466 std::streamsize oldPrec = os.precision(6);
467 G4double DeviationTolerance = 1.0e-05;
468
469 G4double diagDeviation = 0.0;
470 G4double offdDeviationUL = 0.0;
471 G4double offdDeviationDR = 0.0;
472 G4double offdDeviation = 0.0;
473
474 os << " Transformation: " << G4endl;
475
476 // double a = std::max ( 1, 2, 3 ) ;
477
478 G4bool UnitTr = ! transf.IsRotated();
479 diagDeviation = std::max( std::fabs( transf[0] - 1.0 ) , // |rxx - 1|
480 std::fabs( transf[5] - 1.0 ) ); // |ryy - 1|
481 diagDeviation = std::max( diagDeviation,
482 std::fabs( transf[10] - 1.0 ) ); // |rzz - 1|
483
484 offdDeviationUL = std::max( std::fabs( transf[1] ) , // |rxy|
485 std::fabs( transf[2] ) ); // |rxz|
486 offdDeviationUL = std::max( offdDeviationUL,
487 std::fabs( transf[4] ) ); // |ryx|
488
489 offdDeviationDR = std::max( std::fabs( transf[6] ) , // |ryz|
490 std::fabs( transf[8] ) ); // |rzx|
491 offdDeviationDR = std::max( offdDeviationDR,
492 std::fabs( transf[9] ) ); // |rzy|
493 offdDeviation = std::max( offdDeviationUL, offdDeviationDR );
494
495 if( UnitTr || std::max(diagDeviation, offdDeviation) < DeviationTolerance )
496 {
497 os << " UNIT Rotation " << G4endl;
498 }
499 else
500 {
501 os << "rx/x,y,z: "
502 << transf[0] << " " << transf[1] << " " << transf[2] << G4endl
503 << "ry/x,y,z: "
504 << transf[4] << " " << transf[5] << " " << transf[6] << G4endl
505 << "rz/x,y,z: "
506 << transf[8] << " " << transf[9] << " " << transf[10] << G4endl;
507 }
508
509 os << "tr/x,y,z: " << transf[12] << " " << transf[13] << " " << transf[14]
510 << G4endl;
511
512 os.precision(oldPrec);
513
514 return os;
515 }
516