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