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