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