Geant4  10.01.p02
LorentzRotationD.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 // This is the implementation of those parts of the HepLorentzRotation class
7 // which involve decomposition into Boost*Rotation.
8 
9 #ifdef GNUPRAGMA
10 #pragma implementation
11 #endif
12 
13 #include "CLHEP/Vector/LorentzRotation.h"
14 
15 namespace CLHEP {
16 
17 // ---------- Decomposition:
18 
19 void HepLorentzRotation::decompose
20  (HepBoost & bboost, HepRotation & rotation) const {
21 
22  // The boost will be the pure boost based on column 4 of the transformation
23  // matrix. Since the constructor takes the beta vector, and not beta*gamma,
24  // we first divide through by gamma = the tt element. This of course can
25  // never be zero since the last row has t**2 - v**2 = +1.
26 
27  Hep3Vector betaVec ( xt(), yt(), zt() );
28  betaVec *= 1.0 / tt();
29  bboost.set( betaVec );
30 
31  // The rotation will be inverse of B times T.
32 
33  HepBoost B( -betaVec );
34  HepLorentzRotation R( B * *this );
35 
36  HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
37  R.yx(), R.yy(), R.yz(),
38  R.zx(), R.zy(), R.zz() );
39  rotation.set( m1 );
40  rotation.rectify();
41 
42  return;
43 
44 }
45 
46 void HepLorentzRotation::decompose
47  (Hep3Vector & bboost, HepAxisAngle & rotation) const {
48  HepRotation r;
49  HepBoost b;
50  decompose(b,r);
51  bboost = b.boostVector();
52  rotation = r.axisAngle();
53  return;
54 }
55 
56 void HepLorentzRotation::decompose
57  (HepRotation & rotation, HepBoost & bboost) const {
58 
59  // In this case the pure boost is based on row 4 of the matrix.
60 
61  Hep3Vector betaVec( tx(), ty(), tz() );
62  betaVec *= 1.0 / tt();
63  bboost.set( betaVec );
64 
65  // The rotation will be T times the inverse of B.
66 
67  HepBoost B( -betaVec );
68  HepLorentzRotation R( *this * B );
69 
70  HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
71  R.yx(), R.yy(), R.yz(),
72  R.zx(), R.zy(), R.zz() );
73  rotation.set( m1 );
74  rotation.rectify();
75  return;
76 
77 }
78 
79 void HepLorentzRotation::decompose
80  (HepAxisAngle & rotation, Hep3Vector & bboost) const {
81  HepRotation r;
82  HepBoost b;
83  decompose(r,b);
84  rotation = r.axisAngle();
85  bboost = b.boostVector();
86  return;
87 }
88 
89 double HepLorentzRotation::distance2( const HepBoost & b ) const {
90  HepBoost b1;
91  HepRotation r1;
92  decompose( b1, r1 );
93  double db2 = b1.distance2( b );
94  double dr2 = r1.norm2();
95  return ( db2 + dr2 );
96 }
97 
98 double HepLorentzRotation::distance2( const HepRotation & r ) const {
99  HepBoost b1;
100  HepRotation r1;
101  decompose( b1, r1 );
102  double db2 = b1.norm2( );
103  double dr2 = r1.distance2( r );
104  return ( db2 + dr2 );
105 }
106 
107 double HepLorentzRotation::distance2(
108  const HepLorentzRotation & lt ) const {
109  HepBoost b1;
110  HepRotation r1;
111  decompose( b1, r1 );
112  HepBoost b2;
113  HepRotation r2;
114  lt.decompose (b2, r2);
115  double db2 = b1.distance2( b2 );
116  double dr2 = r1.distance2( r2 );
117  return ( db2 + dr2 );
118 }
119 
120 double HepLorentzRotation::howNear( const HepBoost & b ) const {
121  return std::sqrt( distance2( b ) );
122 }
123 double HepLorentzRotation::howNear( const HepRotation & r ) const {
124  return std::sqrt( distance2( r ) );
125 }
126 double HepLorentzRotation::howNear( const HepLorentzRotation & lt )const {
127  return std::sqrt( distance2( lt ) );
128 }
129 
130 bool HepLorentzRotation::isNear(
131  const HepBoost & b, double epsilon ) const {
132  HepBoost b1;
133  HepRotation r1;
134  decompose( b1, r1 );
135  double db2 = b1.distance2(b);
136  if ( db2 > epsilon*epsilon ) {
137  return false; // Saves the time-consuming Rotation::norm2
138  }
139  double dr2 = r1.norm2();
140  return ( (db2 + dr2) <= epsilon*epsilon );
141 }
142 
143 bool HepLorentzRotation::isNear(
144  const HepRotation & r, double epsilon ) const {
145  HepBoost b1;
146  HepRotation r1;
147  decompose( b1, r1 );
148  double db2 = b1.norm2();
149  if ( db2 > epsilon*epsilon ) {
150  return false; // Saves the time-consuming Rotation::distance2
151  }
152  double dr2 = r1.distance2(r);
153  return ( (db2 + dr2) <= epsilon*epsilon );
154 }
155 
156 bool HepLorentzRotation::isNear(
157  const HepLorentzRotation & lt, double epsilon ) const {
158  HepBoost b1;
159  HepRotation r1;
160  decompose( b1, r1 );
161  HepBoost b2;
162  HepRotation r2;
163  lt.decompose (b2, r2);
164  double db2 = b1.distance2(b2);
165  if ( db2 > epsilon*epsilon ) {
166  return false; // Saves the time-consuming Rotation::distance2
167  }
168  double dr2 = r1.distance2(r2);
169  return ( (db2 + dr2) <= epsilon*epsilon );
170 }
171 
172 double HepLorentzRotation::norm2() const {
173  HepBoost b;
174  HepRotation r;
175  decompose( b, r );
176  return b.norm2() + r.norm2();
177 }
178 
179 void HepLorentzRotation::rectify() {
180 
181  // Assuming the representation of this is close to a true LT,
182  // but may have drifted due to round-off error from many operations,
183  // this forms an "exact" orthosymplectic matrix for the LT again.
184 
185  // There are several ways to do this, all equivalent to lowest order in
186  // the corrected error. We choose to form an LT based on the inverse boost
187  // extracted from row 4, and left-multiply by LT to form what would be
188  // a rotation if the LT were kosher. We drop the possibly non-zero t
189  // components of that, rectify that rotation and multiply back by the boost.
190 
191  Hep3Vector beta (tx(), ty(), tz());
192  double gam = tt(); // NaN-proofing
193  if ( gam <= 0 ) {
194  std::cerr << "HepLorentzRotation::rectify() - "
195  << "rectify() on a transformation with tt() <= 0 - will not help!"
196  << std::endl;
197  gam = 1;
198  }
199  beta *= 1.0/gam;
200  HepLorentzRotation R = (*this) * HepBoost(-beta);
201 
202  HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
203  R.yx(), R.yy(), R.yz(),
204  R.zx(), R.zy(), R.zz() );
205 
206  HepRotation Rgood (m1);
207  Rgood.rectify();
208 
209  set ( Rgood, HepBoost(beta) );
210 }
211 
212 } // namespace CLHEP
static const G4double b2
static const G4double b1