00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifdef GNUPRAGMA
00010 #pragma implementation
00011 #endif
00012
00013 #include "CLHEP/Vector/Boost.h"
00014 #include "CLHEP/Vector/Rotation.h"
00015 #include "CLHEP/Vector/LorentzRotation.h"
00016
00017 namespace CLHEP {
00018
00019
00020
00021 HepBoost & HepBoost::set (double bx, double by, double bz) {
00022 double bp2 = bx*bx + by*by + bz*bz;
00023
00024
00025
00026
00027 double ggamma = 1.0 / std::sqrt(1.0 - bp2);
00028 double bgamma = ggamma * ggamma / (1.0 + ggamma);
00029 rep_.xx_ = 1.0 + bgamma * bx * bx;
00030 rep_.yy_ = 1.0 + bgamma * by * by;
00031 rep_.zz_ = 1.0 + bgamma * bz * bz;
00032 rep_.xy_ = bgamma * bx * by;
00033 rep_.xz_ = bgamma * bx * bz;
00034 rep_.yz_ = bgamma * by * bz;
00035 rep_.xt_ = ggamma * bx;
00036 rep_.yt_ = ggamma * by;
00037 rep_.zt_ = ggamma * bz;
00038 rep_.tt_ = ggamma;
00039 return *this;
00040 }
00041
00042 HepBoost & HepBoost::set (const HepRep4x4Symmetric & m1) {
00043 rep_ = m1;
00044 return *this;
00045 }
00046
00047 HepBoost & HepBoost::set (Hep3Vector ddirection, double bbeta) {
00048 double length = ddirection.mag();
00049 if (length <= 0) {
00050 std::cerr << "HepBoost::set() - "
00051 << "Direction supplied to set HepBoost is zero." << std::endl;
00052 set (0,0,0);
00053 return *this;
00054 }
00055 set(bbeta*ddirection.x()/length,
00056 bbeta*ddirection.y()/length,
00057 bbeta*ddirection.z()/length);
00058 return *this;
00059 }
00060
00061 HepBoost & HepBoost::set (const Hep3Vector & bboost) {
00062 return set (bboost.x(), bboost.y(), bboost.z());
00063 }
00064
00065
00066
00067
00068
00069 void HepBoost::decompose (HepRotation & rotation, HepBoost & boost) const {
00070 HepAxisAngle vdelta = HepAxisAngle();
00071 rotation = HepRotation(vdelta);
00072 Hep3Vector bbeta = boostVector();
00073 boost = HepBoost(bbeta);
00074 }
00075
00076 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const {
00077 rotation = HepAxisAngle();
00078 boost = boostVector();
00079 }
00080
00081 void HepBoost::decompose (HepBoost & boost, HepRotation & rotation) const {
00082 HepAxisAngle vdelta = HepAxisAngle();
00083 rotation = HepRotation(vdelta);
00084 Hep3Vector bbeta = boostVector();
00085 boost = HepBoost(bbeta);
00086 }
00087
00088 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const {
00089 rotation = HepAxisAngle();
00090 boost = boostVector();
00091 }
00092
00093
00094
00095 double HepBoost::distance2( const HepRotation & r ) const {
00096 double db2 = norm2();
00097 double dr2 = r.norm2();
00098 return (db2 + dr2);
00099 }
00100
00101 double HepBoost::distance2( const HepLorentzRotation & lt ) const {
00102 HepBoost b1;
00103 HepRotation r1;
00104 lt.decompose(b1,r1);
00105 double db2 = distance2(b1);
00106 double dr2 = r1.norm2();
00107 return (db2 + dr2);
00108 }
00109
00110 double HepBoost::howNear ( const HepRotation & r ) const {
00111 return std::sqrt(distance2(r));
00112 }
00113
00114 double HepBoost::howNear ( const HepLorentzRotation & lt ) const {
00115 return std::sqrt(distance2(lt));
00116 }
00117
00118 bool HepBoost::isNear (const HepRotation & r, double epsilon) const {
00119 double db2 = norm2();
00120 if (db2 > epsilon*epsilon) return false;
00121 double dr2 = r.norm2();
00122 return (db2+dr2 <= epsilon*epsilon);
00123 }
00124
00125 bool HepBoost::isNear (const HepLorentzRotation & lt,
00126 double epsilon) const {
00127 HepBoost b1;
00128 HepRotation r1;
00129 double db2 = distance2(b1);
00130 lt.decompose(b1,r1);
00131 if (db2 > epsilon*epsilon) return false;
00132 double dr2 = r1.norm2();
00133 return (db2 + dr2);
00134 }
00135
00136
00137
00138 double HepBoost::norm2() const {
00139 register double bgx = rep_.xt_;
00140 register double bgy = rep_.yt_;
00141 register double bgz = rep_.zt_;
00142 return bgx*bgx+bgy*bgy+bgz*bgz;
00143 }
00144
00145 void HepBoost::rectify() {
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161 double gam = tt();
00162 if (gam <= 0) {
00163 std::cerr << "HepBoost::rectify() - "
00164 << "Attempt to rectify a boost with non-positive gamma." << std::endl;
00165 if (gam==0) return;
00166 }
00167 Hep3Vector boost (xt(), yt(), zt());
00168 boost /= tt();
00169 if ( boost.mag2() >= 1 ) {
00170 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) );
00171 }
00172 set ( boost );
00173 }
00174
00175
00176
00177
00178
00179 HepLorentzRotation
00180 HepBoost::matrixMultiplication(const HepRep4x4 & m1) const {
00181 HepRep4x4Symmetric r = rep4x4Symmetric();
00182 return HepLorentzRotation( HepRep4x4 (
00183 r.xx_*m1.xx_ + r.xy_*m1.yx_ + r.xz_*m1.zx_ + r.xt_*m1.tx_,
00184 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.zy_ + r.xt_*m1.ty_,
00185 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.tz_,
00186 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
00187
00188 r.xy_*m1.xx_ + r.yy_*m1.yx_ + r.yz_*m1.zx_ + r.yt_*m1.tx_,
00189 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.zy_ + r.yt_*m1.ty_,
00190 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.tz_,
00191 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
00192
00193 r.xz_*m1.xx_ + r.yz_*m1.yx_ + r.zz_*m1.zx_ + r.zt_*m1.tx_,
00194 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.zy_ + r.zt_*m1.ty_,
00195 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.tz_,
00196 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
00197
00198 r.xt_*m1.xx_ + r.yt_*m1.yx_ + r.zt_*m1.zx_ + r.tt_*m1.tx_,
00199 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.zy_ + r.tt_*m1.ty_,
00200 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.tz_,
00201 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
00202 }
00203
00204 HepLorentzRotation
00205 HepBoost::matrixMultiplication(const HepRep4x4Symmetric & m1) const {
00206 HepRep4x4Symmetric r = rep4x4Symmetric();
00207 return HepLorentzRotation( HepRep4x4 (
00208 r.xx_*m1.xx_ + r.xy_*m1.xy_ + r.xz_*m1.xz_ + r.xt_*m1.xt_,
00209 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.yz_ + r.xt_*m1.yt_,
00210 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.zt_,
00211 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
00212
00213 r.xy_*m1.xx_ + r.yy_*m1.xy_ + r.yz_*m1.xz_ + r.yt_*m1.xt_,
00214 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.yz_ + r.yt_*m1.yt_,
00215 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.zt_,
00216 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
00217
00218 r.xz_*m1.xx_ + r.yz_*m1.xy_ + r.zz_*m1.xz_ + r.zt_*m1.xt_,
00219 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.yz_ + r.zt_*m1.yt_,
00220 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.zt_,
00221 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
00222
00223 r.xt_*m1.xx_ + r.yt_*m1.xy_ + r.zt_*m1.xz_ + r.tt_*m1.xt_,
00224 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.yz_ + r.tt_*m1.yt_,
00225 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.zt_,
00226 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
00227 }
00228
00229 HepLorentzRotation
00230 HepBoost::operator* (const HepLorentzRotation & lt) const {
00231 return matrixMultiplication(lt.rep4x4());
00232 }
00233
00234 HepLorentzRotation
00235 HepBoost::operator* (const HepBoost & b) const {
00236 return matrixMultiplication(b.rep_);
00237 }
00238
00239 HepLorentzRotation
00240 HepBoost::operator* (const HepRotation & r) const {
00241 return matrixMultiplication(r.rep4x4());
00242 }
00243
00244
00245
00246 std::ostream & HepBoost::print( std::ostream & os ) const {
00247 if ( rep_.tt_ <= 1 ) {
00248 os << "Lorentz Boost( IDENTITY )";
00249 } else {
00250 double norm = boostVector().mag();
00251 os << "\nLorentz Boost " << boostVector()/norm <<
00252 "\n{beta = " << beta() << " gamma = " << gamma() << "}\n";
00253 }
00254 return os;
00255 }
00256
00257 }