10 #pragma implementation
13 #include "CLHEP/Vector/Boost.h"
14 #include "CLHEP/Vector/Rotation.h"
15 #include "CLHEP/Vector/LorentzRotation.h"
21 HepBoost & HepBoost::set (
double bx,
double by,
double bz) {
22 double bp2 = bx*bx + by*by + bz*bz;
27 double ggamma = 1.0 / std::sqrt(1.0 - bp2);
28 double bgamma = ggamma * ggamma / (1.0 + ggamma);
29 rep_.xx_ = 1.0 + bgamma * bx * bx;
30 rep_.yy_ = 1.0 + bgamma * by * by;
31 rep_.zz_ = 1.0 + bgamma * bz * bz;
32 rep_.xy_ = bgamma * bx * by;
33 rep_.xz_ = bgamma * bx * bz;
34 rep_.yz_ = bgamma * by * bz;
35 rep_.xt_ = ggamma * bx;
36 rep_.yt_ = ggamma * by;
37 rep_.zt_ = ggamma * bz;
42 HepBoost & HepBoost::set (
const HepRep4x4Symmetric & m1) {
47 HepBoost & HepBoost::set (Hep3Vector ddirection,
double bbeta) {
48 double length = ddirection.mag();
50 std::cerr <<
"HepBoost::set() - "
51 <<
"Direction supplied to set HepBoost is zero." << std::endl;
55 set(bbeta*ddirection.x()/length,
56 bbeta*ddirection.y()/length,
57 bbeta*ddirection.z()/length);
61 HepBoost & HepBoost::set (
const Hep3Vector & bboost) {
62 return set (bboost.x(), bboost.y(), bboost.z());
69 void HepBoost::decompose (HepRotation & rotation, HepBoost & boost)
const {
70 HepAxisAngle vdelta = HepAxisAngle();
71 rotation = HepRotation(vdelta);
72 Hep3Vector bbeta = boostVector();
73 boost = HepBoost(bbeta);
76 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost)
const {
77 rotation = HepAxisAngle();
78 boost = boostVector();
81 void HepBoost::decompose (HepBoost & boost, HepRotation & rotation)
const {
82 HepAxisAngle vdelta = HepAxisAngle();
83 rotation = HepRotation(vdelta);
84 Hep3Vector bbeta = boostVector();
85 boost = HepBoost(bbeta);
88 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation)
const {
89 rotation = HepAxisAngle();
90 boost = boostVector();
95 double HepBoost::distance2(
const HepRotation & r )
const {
97 double dr2 = r.norm2();
101 double HepBoost::distance2(
const HepLorentzRotation & lt )
const {
105 double db2 = distance2(b1);
106 double dr2 = r1.norm2();
110 double HepBoost::howNear (
const HepRotation & r )
const {
111 return std::sqrt(distance2(r));
114 double HepBoost::howNear (
const HepLorentzRotation & lt )
const {
115 return std::sqrt(distance2(lt));
118 bool HepBoost::isNear (
const HepRotation & r,
double epsilon)
const {
119 double db2 = norm2();
120 if (db2 > epsilon*epsilon)
return false;
121 double dr2 = r.norm2();
122 return (db2+dr2 <= epsilon*epsilon);
125 bool HepBoost::isNear (
const HepLorentzRotation & lt,
126 double epsilon)
const {
129 double db2 = distance2(b1);
131 if (db2 > epsilon*epsilon)
return false;
132 double dr2 = r1.norm2();
138 double HepBoost::norm2()
const {
139 double bgx = rep_.xt_;
140 double bgy = rep_.yt_;
141 double bgz = rep_.zt_;
142 return bgx*bgx+bgy*bgy+bgz*bgz;
145 void HepBoost::rectify() {
163 std::cerr <<
"HepBoost::rectify() - "
164 <<
"Attempt to rectify a boost with non-positive gamma." << std::endl;
167 Hep3Vector boost (xt(), yt(), zt());
169 if ( boost.mag2() >= 1 ) {
170 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) );
180 HepBoost::matrixMultiplication(
const HepRep4x4 & m1)
const {
181 HepRep4x4Symmetric r = rep4x4Symmetric();
182 return HepLorentzRotation( HepRep4x4 (
183 r.xx_*m1.xx_ + r.xy_*m1.yx_ + r.xz_*m1.zx_ + r.xt_*m1.tx_,
184 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.zy_ + r.xt_*m1.ty_,
185 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.tz_,
186 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
188 r.xy_*m1.xx_ + r.yy_*m1.yx_ + r.yz_*m1.zx_ + r.yt_*m1.tx_,
189 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.zy_ + r.yt_*m1.ty_,
190 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.tz_,
191 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
193 r.xz_*m1.xx_ + r.yz_*m1.yx_ + r.zz_*m1.zx_ + r.zt_*m1.tx_,
194 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.zy_ + r.zt_*m1.ty_,
195 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.tz_,
196 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
198 r.xt_*m1.xx_ + r.yt_*m1.yx_ + r.zt_*m1.zx_ + r.tt_*m1.tx_,
199 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.zy_ + r.tt_*m1.ty_,
200 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.tz_,
201 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
205 HepBoost::matrixMultiplication(
const HepRep4x4Symmetric & m1)
const {
206 HepRep4x4Symmetric r = rep4x4Symmetric();
207 return HepLorentzRotation( HepRep4x4 (
208 r.xx_*m1.xx_ + r.xy_*m1.xy_ + r.xz_*m1.xz_ + r.xt_*m1.xt_,
209 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.yz_ + r.xt_*m1.yt_,
210 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.zt_,
211 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
213 r.xy_*m1.xx_ + r.yy_*m1.xy_ + r.yz_*m1.xz_ + r.yt_*m1.xt_,
214 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.yz_ + r.yt_*m1.yt_,
215 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.zt_,
216 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
218 r.xz_*m1.xx_ + r.yz_*m1.xy_ + r.zz_*m1.xz_ + r.zt_*m1.xt_,
219 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.yz_ + r.zt_*m1.yt_,
220 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.zt_,
221 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
223 r.xt_*m1.xx_ + r.yt_*m1.xy_ + r.zt_*m1.xz_ + r.tt_*m1.xt_,
224 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.yz_ + r.tt_*m1.yt_,
225 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.zt_,
226 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
231 return matrixMultiplication(lt.rep4x4());
236 return matrixMultiplication(b.rep_);
241 return matrixMultiplication(r.rep4x4());
247 if ( rep_.tt_ <= 1 ) {
248 os <<
"Lorentz Boost( IDENTITY )";
250 double norm = boostVector().mag();
251 os <<
"\nLorentz Boost " << boostVector()/norm <<
252 "\n{beta = " << beta() <<
" gamma = " << gamma() <<
"}\n";
G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
void print(const std::vector< T > &data)