| 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 the HepBoost 6 // This is the implementation of the HepBoost class.
7 // 7 //
8 8
>> 9 #ifdef GNUPRAGMA
>> 10 #pragma implementation
>> 11 #endif
>> 12
9 #include "CLHEP/Vector/Boost.h" 13 #include "CLHEP/Vector/Boost.h"
10 #include "CLHEP/Vector/Rotation.h" 14 #include "CLHEP/Vector/Rotation.h"
11 #include "CLHEP/Vector/LorentzRotation.h" 15 #include "CLHEP/Vector/LorentzRotation.h"
12 <<
13 #include <cmath> <<
14 #include <iostream> <<
15 16
16 namespace CLHEP { 17 namespace CLHEP {
17 18
18 // ---------- Constructors and Assignment: 19 // ---------- Constructors and Assignment:
19 20
20 HepBoost & HepBoost::set (double bx, double by 21 HepBoost & HepBoost::set (double bx, double by, double bz) {
21 double bp2 = bx*bx + by*by + bz*bz; 22 double bp2 = bx*bx + by*by + bz*bz;
22 // if (bp2 >= 1) { 23 // if (bp2 >= 1) {
23 // std::cerr << "HepBoost::set() - " 24 // std::cerr << "HepBoost::set() - "
24 // << "Boost Vector supplied to set HepBo 25 // << "Boost Vector supplied to set HepBoost represents speed >= c." << std::endl;
25 // } 26 // }
26 double ggamma = 1.0 / std::sqrt(1.0 - bp2); 27 double ggamma = 1.0 / std::sqrt(1.0 - bp2);
27 double bgamma = ggamma * ggamma / (1.0 + gga 28 double bgamma = ggamma * ggamma / (1.0 + ggamma);
28 rep_.xx_ = 1.0 + bgamma * bx * bx; 29 rep_.xx_ = 1.0 + bgamma * bx * bx;
29 rep_.yy_ = 1.0 + bgamma * by * by; 30 rep_.yy_ = 1.0 + bgamma * by * by;
30 rep_.zz_ = 1.0 + bgamma * bz * bz; 31 rep_.zz_ = 1.0 + bgamma * bz * bz;
31 rep_.xy_ = bgamma * bx * by; 32 rep_.xy_ = bgamma * bx * by;
32 rep_.xz_ = bgamma * bx * bz; 33 rep_.xz_ = bgamma * bx * bz;
33 rep_.yz_ = bgamma * by * bz; 34 rep_.yz_ = bgamma * by * bz;
34 rep_.xt_ = ggamma * bx; 35 rep_.xt_ = ggamma * bx;
35 rep_.yt_ = ggamma * by; 36 rep_.yt_ = ggamma * by;
36 rep_.zt_ = ggamma * bz; 37 rep_.zt_ = ggamma * bz;
37 rep_.tt_ = ggamma; 38 rep_.tt_ = ggamma;
38 return *this; 39 return *this;
39 } 40 }
40 41
41 HepBoost & HepBoost::set (const HepRep4x4Symme 42 HepBoost & HepBoost::set (const HepRep4x4Symmetric & m1) {
42 rep_ = m1; 43 rep_ = m1;
43 return *this; 44 return *this;
44 } 45 }
45 46
46 HepBoost & HepBoost::set (Hep3Vector ddirectio 47 HepBoost & HepBoost::set (Hep3Vector ddirection, double bbeta) {
47 double length = ddirection.mag(); 48 double length = ddirection.mag();
48 if (length <= 0) { // Nan-proofing 49 if (length <= 0) { // Nan-proofing
49 std::cerr << "HepBoost::set() - " 50 std::cerr << "HepBoost::set() - "
50 << "Direction supplied to set HepBoost i 51 << "Direction supplied to set HepBoost is zero." << std::endl;
51 set (0,0,0); 52 set (0,0,0);
52 return *this; 53 return *this;
53 } 54 }
54 set(bbeta*ddirection.x()/length, 55 set(bbeta*ddirection.x()/length,
55 bbeta*ddirection.y()/length, 56 bbeta*ddirection.y()/length,
56 bbeta*ddirection.z()/length); 57 bbeta*ddirection.z()/length);
57 return *this; 58 return *this;
58 } 59 }
59 60
60 HepBoost & HepBoost::set (const Hep3Vector & b 61 HepBoost & HepBoost::set (const Hep3Vector & bboost) {
61 return set (bboost.x(), bboost.y(), bboost.z 62 return set (bboost.x(), bboost.y(), bboost.z());
62 } 63 }
63 64
64 // ---------- Accessors: 65 // ---------- Accessors:
65 66
66 // ---------- Decomposition: 67 // ---------- Decomposition:
67 68
68 void HepBoost::decompose (HepRotation & rotati 69 void HepBoost::decompose (HepRotation & rotation, HepBoost & boost) const {
69 HepAxisAngle vdelta = HepAxisAngle(); 70 HepAxisAngle vdelta = HepAxisAngle();
70 rotation = HepRotation(vdelta); 71 rotation = HepRotation(vdelta);
71 Hep3Vector bbeta = boostVector(); 72 Hep3Vector bbeta = boostVector();
72 boost = HepBoost(bbeta); 73 boost = HepBoost(bbeta);
73 } 74 }
74 75
75 void HepBoost::decompose (HepAxisAngle & rotat 76 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const {
76 rotation = HepAxisAngle(); 77 rotation = HepAxisAngle();
77 boost = boostVector(); 78 boost = boostVector();
78 } 79 }
79 80
80 void HepBoost::decompose (HepBoost & boost, He 81 void HepBoost::decompose (HepBoost & boost, HepRotation & rotation) const {
81 HepAxisAngle vdelta = HepAxisAngle(); 82 HepAxisAngle vdelta = HepAxisAngle();
82 rotation = HepRotation(vdelta); 83 rotation = HepRotation(vdelta);
83 Hep3Vector bbeta = boostVector(); 84 Hep3Vector bbeta = boostVector();
84 boost = HepBoost(bbeta); 85 boost = HepBoost(bbeta);
85 } 86 }
86 87
87 void HepBoost::decompose (Hep3Vector & boost, 88 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const {
88 rotation = HepAxisAngle(); 89 rotation = HepAxisAngle();
89 boost = boostVector(); 90 boost = boostVector();
90 } 91 }
91 92
92 // ---------- Comparisons: 93 // ---------- Comparisons:
93 94
94 double HepBoost::distance2( const HepRotation 95 double HepBoost::distance2( const HepRotation & r ) const {
95 double db2 = norm2(); 96 double db2 = norm2();
96 double dr2 = r.norm2(); 97 double dr2 = r.norm2();
97 return (db2 + dr2); 98 return (db2 + dr2);
98 } 99 }
99 100
100 double HepBoost::distance2( const HepLorentzRo 101 double HepBoost::distance2( const HepLorentzRotation & lt ) const {
101 HepBoost b1; 102 HepBoost b1;
102 HepRotation r1; 103 HepRotation r1;
103 lt.decompose(b1,r1); 104 lt.decompose(b1,r1);
104 double db2 = distance2(b1); 105 double db2 = distance2(b1);
105 double dr2 = r1.norm2(); 106 double dr2 = r1.norm2();
106 return (db2 + dr2); 107 return (db2 + dr2);
107 } 108 }
108 109
109 double HepBoost::howNear ( const HepRotation & 110 double HepBoost::howNear ( const HepRotation & r ) const {
110 return std::sqrt(distance2(r)); 111 return std::sqrt(distance2(r));
111 } 112 }
112 113
113 double HepBoost::howNear ( const HepLorentzRot 114 double HepBoost::howNear ( const HepLorentzRotation & lt ) const {
114 return std::sqrt(distance2(lt)); 115 return std::sqrt(distance2(lt));
115 } 116 }
116 117
117 bool HepBoost::isNear (const HepRotation & r, 118 bool HepBoost::isNear (const HepRotation & r, double epsilon) const {
118 double db2 = norm2(); 119 double db2 = norm2();
119 if (db2 > epsilon*epsilon) return false; 120 if (db2 > epsilon*epsilon) return false;
120 double dr2 = r.norm2(); 121 double dr2 = r.norm2();
121 return (db2+dr2 <= epsilon*epsilon); 122 return (db2+dr2 <= epsilon*epsilon);
122 } 123 }
123 124
124 bool HepBoost::isNear (const HepLorentzRotatio 125 bool HepBoost::isNear (const HepLorentzRotation & lt,
125 double epsilon) const { 126 double epsilon) const {
126 HepBoost b1; 127 HepBoost b1;
127 HepRotation r1; 128 HepRotation r1;
128 double db2 = distance2(b1); 129 double db2 = distance2(b1);
129 lt.decompose(b1,r1); 130 lt.decompose(b1,r1);
130 if (db2 > epsilon*epsilon) return false; 131 if (db2 > epsilon*epsilon) return false;
131 double dr2 = r1.norm2(); 132 double dr2 = r1.norm2();
132 return (db2 + dr2); 133 return (db2 + dr2);
133 } 134 }
134 135
135 // ---------- Properties: 136 // ---------- Properties:
136 137
137 double HepBoost::norm2() const { 138 double HepBoost::norm2() const {
138 double bgx = rep_.xt_; 139 double bgx = rep_.xt_;
139 double bgy = rep_.yt_; 140 double bgy = rep_.yt_;
140 double bgz = rep_.zt_; 141 double bgz = rep_.zt_;
141 return bgx*bgx+bgy*bgy+bgz*bgz; 142 return bgx*bgx+bgy*bgy+bgz*bgz;
142 } 143 }
143 144
144 void HepBoost::rectify() { 145 void HepBoost::rectify() {
145 // Assuming the representation of this is cl 146 // Assuming the representation of this is close to a true pure boost,
146 // but may have drifted due to round-off err 147 // but may have drifted due to round-off error from many operations,
147 // this forms an "exact" pure boost matrix f 148 // this forms an "exact" pure boost matrix for the LT again.
148 149
149 // The natural way to do this is to use the 150 // The natural way to do this is to use the t column as a boost and set
150 // based on that boost vector. 151 // based on that boost vector.
151 152
152 // There is perhaps danger that this boost v 153 // There is perhaps danger that this boost vector will appear equal to or
153 // greater than a unit vector; the best we c 154 // greater than a unit vector; the best we can do for such a case is use
154 // a boost in that direction but rescaled to 155 // a boost in that direction but rescaled to just less than one.
155 156
156 // There is in principle no way that gamma c 157 // There is in principle no way that gamma could have become negative,
157 // but if that happens, we ZMthrow and (if c 158 // but if that happens, we ZMthrow and (if continuing) just rescale, which
158 // will change the sign of the last column w 159 // will change the sign of the last column when computing the boost.
159 160
160 double gam = tt(); 161 double gam = tt();
161 if (gam <= 0) { // 4/12/01 mf 162 if (gam <= 0) { // 4/12/01 mf
162 std::cerr << "HepBoost::rectify() - " 163 std::cerr << "HepBoost::rectify() - "
163 << "Attempt to rectify a boost with non- 164 << "Attempt to rectify a boost with non-positive gamma." << std::endl;
164 if (gam==0) return; // NaN-proof 165 if (gam==0) return; // NaN-proofing
165 } 166 }
166 Hep3Vector boost (xt(), yt(), zt()); 167 Hep3Vector boost (xt(), yt(), zt());
167 boost /= tt(); 168 boost /= tt();
168 if ( boost.mag2() >= 1 ) { // NaN-p 169 if ( boost.mag2() >= 1 ) { // NaN-proofing:
169 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) 170 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) ); // used to just check > 1
170 } 171 }
171 set ( boost ); 172 set ( boost );
172 } 173 }
173 174
174 // ---------- Application is all in .icc 175 // ---------- Application is all in .icc
175 176
176 // ---------- Operations in the group of 4-Rot 177 // ---------- Operations in the group of 4-Rotations
177 178
178 HepLorentzRotation 179 HepLorentzRotation
179 HepBoost::matrixMultiplication(const HepRep4x4 180 HepBoost::matrixMultiplication(const HepRep4x4 & m1) const {
180 HepRep4x4Symmetric r = rep4x4Symmetric(); 181 HepRep4x4Symmetric r = rep4x4Symmetric();
181 return HepLorentzRotation( HepRep4x4 ( 182 return HepLorentzRotation( HepRep4x4 (
182 r.xx_*m1.xx_ + r.xy_*m1.yx_ + r.xz_*m1.zx_ 183 r.xx_*m1.xx_ + r.xy_*m1.yx_ + r.xz_*m1.zx_ + r.xt_*m1.tx_,
183 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.zy_ 184 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.zy_ + r.xt_*m1.ty_,
184 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ 185 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.tz_,
185 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ 186 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
186 187
187 r.xy_*m1.xx_ + r.yy_*m1.yx_ + r.yz_*m1.zx_ 188 r.xy_*m1.xx_ + r.yy_*m1.yx_ + r.yz_*m1.zx_ + r.yt_*m1.tx_,
188 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.zy_ 189 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.zy_ + r.yt_*m1.ty_,
189 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ 190 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.tz_,
190 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ 191 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
191 192
192 r.xz_*m1.xx_ + r.yz_*m1.yx_ + r.zz_*m1.zx_ 193 r.xz_*m1.xx_ + r.yz_*m1.yx_ + r.zz_*m1.zx_ + r.zt_*m1.tx_,
193 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.zy_ 194 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.zy_ + r.zt_*m1.ty_,
194 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ 195 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.tz_,
195 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ 196 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
196 197
197 r.xt_*m1.xx_ + r.yt_*m1.yx_ + r.zt_*m1.zx_ 198 r.xt_*m1.xx_ + r.yt_*m1.yx_ + r.zt_*m1.zx_ + r.tt_*m1.tx_,
198 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.zy_ 199 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.zy_ + r.tt_*m1.ty_,
199 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ 200 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.tz_,
200 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ 201 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
201 } 202 }
202 203
203 HepLorentzRotation 204 HepLorentzRotation
204 HepBoost::matrixMultiplication(const HepRep4x4 205 HepBoost::matrixMultiplication(const HepRep4x4Symmetric & m1) const {
205 HepRep4x4Symmetric r = rep4x4Symmetric(); 206 HepRep4x4Symmetric r = rep4x4Symmetric();
206 return HepLorentzRotation( HepRep4x4 ( 207 return HepLorentzRotation( HepRep4x4 (
207 r.xx_*m1.xx_ + r.xy_*m1.xy_ + r.xz_*m1.xz_ 208 r.xx_*m1.xx_ + r.xy_*m1.xy_ + r.xz_*m1.xz_ + r.xt_*m1.xt_,
208 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.yz_ 209 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.yz_ + r.xt_*m1.yt_,
209 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ 210 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.zt_,
210 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ 211 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
211 212
212 r.xy_*m1.xx_ + r.yy_*m1.xy_ + r.yz_*m1.xz_ 213 r.xy_*m1.xx_ + r.yy_*m1.xy_ + r.yz_*m1.xz_ + r.yt_*m1.xt_,
213 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.yz_ 214 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.yz_ + r.yt_*m1.yt_,
214 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ 215 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.zt_,
215 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ 216 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
216 217
217 r.xz_*m1.xx_ + r.yz_*m1.xy_ + r.zz_*m1.xz_ 218 r.xz_*m1.xx_ + r.yz_*m1.xy_ + r.zz_*m1.xz_ + r.zt_*m1.xt_,
218 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.yz_ 219 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.yz_ + r.zt_*m1.yt_,
219 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ 220 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.zt_,
220 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ 221 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
221 222
222 r.xt_*m1.xx_ + r.yt_*m1.xy_ + r.zt_*m1.xz_ 223 r.xt_*m1.xx_ + r.yt_*m1.xy_ + r.zt_*m1.xz_ + r.tt_*m1.xt_,
223 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.yz_ 224 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.yz_ + r.tt_*m1.yt_,
224 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ 225 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.zt_,
225 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ 226 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
226 } 227 }
227 228
228 HepLorentzRotation 229 HepLorentzRotation
229 HepBoost::operator* (const HepLorentzRotation 230 HepBoost::operator* (const HepLorentzRotation & lt) const {
230 return matrixMultiplication(lt.rep4x4()); 231 return matrixMultiplication(lt.rep4x4());
231 } 232 }
232 233
233 HepLorentzRotation 234 HepLorentzRotation
234 HepBoost::operator* (const HepBoost & b) const 235 HepBoost::operator* (const HepBoost & b) const {
235 return matrixMultiplication(b.rep_); 236 return matrixMultiplication(b.rep_);
236 } 237 }
237 238
238 HepLorentzRotation 239 HepLorentzRotation
239 HepBoost::operator* (const HepRotation & r) co 240 HepBoost::operator* (const HepRotation & r) const {
240 return matrixMultiplication(r.rep4x4()); 241 return matrixMultiplication(r.rep4x4());
241 } 242 }
242 243
243 // ---------- I/O: 244 // ---------- I/O:
244 245
245 std::ostream & HepBoost::print( std::ostream & 246 std::ostream & HepBoost::print( std::ostream & os ) const {
246 if ( rep_.tt_ <= 1 ) { 247 if ( rep_.tt_ <= 1 ) {
247 os << "Lorentz Boost( IDENTITY )"; 248 os << "Lorentz Boost( IDENTITY )";
248 } else { 249 } else {
249 double norm = boostVector().mag(); 250 double norm = boostVector().mag();
250 os << "\nLorentz Boost " << boostVector()/ 251 os << "\nLorentz Boost " << boostVector()/norm <<
251 "\n{beta = " << beta() << " gamma = 252 "\n{beta = " << beta() << " gamma = " << gamma() << "}\n";
252 } 253 }
253 return os; 254 return os;
254 } 255 }
255 256
256 } // namespace CLHEP 257 } // namespace CLHEP
257 258