CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // --------------------------------------------------------------------------- 00003 // 00004 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00005 // 00006 // This is the implementation of the HepBoost class. 00007 // 00008 00009 #ifdef GNUPRAGMA 00010 #pragma implementation 00011 #endif 00012 00013 #include "CLHEP/Vector/defs.h" 00014 #include "CLHEP/Vector/Boost.h" 00015 #include "CLHEP/Vector/Rotation.h" 00016 #include "CLHEP/Vector/LorentzRotation.h" 00017 #include "CLHEP/Vector/ZMxpv.h" 00018 00019 namespace CLHEP { 00020 00021 // ---------- Constructors and Assignment: 00022 00023 HepBoost & HepBoost::set (double bx, double by, double bz) { 00024 double bp2 = bx*bx + by*by + bz*bz; 00025 if (bp2 >= 1) { 00026 ZMthrowA (ZMxpvTachyonic( 00027 "Boost Vector supplied to set HepBoost represents speed >= c.")); 00028 } 00029 double gamma = 1.0 / sqrt(1.0 - bp2); 00030 double bgamma = gamma * gamma / (1.0 + gamma); 00031 rep_.xx_ = 1.0 + bgamma * bx * bx; 00032 rep_.yy_ = 1.0 + bgamma * by * by; 00033 rep_.zz_ = 1.0 + bgamma * bz * bz; 00034 rep_.xy_ = bgamma * bx * by; 00035 rep_.xz_ = bgamma * bx * bz; 00036 rep_.yz_ = bgamma * by * bz; 00037 rep_.xt_ = gamma * bx; 00038 rep_.yt_ = gamma * by; 00039 rep_.zt_ = gamma * bz; 00040 rep_.tt_ = gamma; 00041 return *this; 00042 } 00043 00044 HepBoost & HepBoost::set (const HepRep4x4Symmetric & m) { 00045 rep_ = m; 00046 return *this; 00047 } 00048 00049 HepBoost & HepBoost::set (Hep3Vector direction, double beta) { 00050 double length = direction.mag(); 00051 if (length <= 0) { // Nan-proofing 00052 ZMthrowA (ZMxpvZeroVector( 00053 "Direction supplied to set HepBoost is zero.")); 00054 set (0,0,0); 00055 return *this; 00056 } 00057 set(beta*direction.x()/length, 00058 beta*direction.y()/length, 00059 beta*direction.z()/length); 00060 return *this; 00061 } 00062 00063 HepBoost & HepBoost::set (const Hep3Vector & boost) { 00064 return set (boost.x(), boost.y(), boost.z()); 00065 } 00066 00067 // ---------- Accessors: 00068 00069 // ---------- Decomposition: 00070 00071 void HepBoost::decompose (HepRotation & rotation, HepBoost & boost) const { 00072 HepAxisAngle vdelta = HepAxisAngle(); 00073 rotation = HepRotation(vdelta); 00074 Hep3Vector beta = boostVector(); 00075 boost = HepBoost(beta); 00076 } 00077 00078 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const { 00079 rotation = HepAxisAngle(); 00080 boost = boostVector(); 00081 } 00082 00083 void HepBoost::decompose (HepBoost & boost, HepRotation & rotation) const { 00084 HepAxisAngle vdelta = HepAxisAngle(); 00085 rotation = HepRotation(vdelta); 00086 Hep3Vector beta = boostVector(); 00087 boost = HepBoost(beta); 00088 } 00089 00090 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const { 00091 rotation = HepAxisAngle(); 00092 boost = boostVector(); 00093 } 00094 00095 // ---------- Comparisons: 00096 00097 double HepBoost::distance2( const HepRotation & r ) const { 00098 double db2 = norm2(); 00099 double dr2 = r.norm2(); 00100 return (db2 + dr2); 00101 } 00102 00103 double HepBoost::distance2( const HepLorentzRotation & lt ) const { 00104 HepBoost b1; 00105 HepRotation r1; 00106 lt.decompose(b1,r1); 00107 double db2 = distance2(b1); 00108 double dr2 = r1.norm2(); 00109 return (db2 + dr2); 00110 } 00111 00112 double HepBoost::howNear ( const HepRotation & r ) const { 00113 return sqrt(distance2(r)); 00114 } 00115 00116 double HepBoost::howNear ( const HepLorentzRotation & lt ) const { 00117 return sqrt(distance2(lt)); 00118 } 00119 00120 bool HepBoost::isNear (const HepRotation & r, double epsilon) const { 00121 double db2 = norm2(); 00122 if (db2 > epsilon*epsilon) return false; 00123 double dr2 = r.norm2(); 00124 return (db2+dr2 <= epsilon*epsilon); 00125 } 00126 00127 bool HepBoost::isNear (const HepLorentzRotation & lt, 00128 double epsilon) const { 00129 HepBoost b1; 00130 HepRotation r1; 00131 double db2 = distance2(b1); 00132 lt.decompose(b1,r1); 00133 if (db2 > epsilon*epsilon) return false; 00134 double dr2 = r1.norm2(); 00135 return (db2 + dr2); 00136 } 00137 00138 // ---------- Properties: 00139 00140 double HepBoost::norm2() const { 00141 register double bgx = rep_.xt_; 00142 register double bgy = rep_.yt_; 00143 register double bgz = rep_.zt_; 00144 return bgx*bgx+bgy*bgy+bgz*bgz; 00145 } 00146 00147 void HepBoost::rectify() { 00148 // Assuming the representation of this is close to a true pure boost, 00149 // but may have drifted due to round-off error from many operations, 00150 // this forms an "exact" pure boost matrix for the LT again. 00151 00152 // The natural way to do this is to use the t column as a boost and set 00153 // based on that boost vector. 00154 00155 // There is perhaps danger that this boost vector will appear equal to or 00156 // greater than a unit vector; the best we can do for such a case is use 00157 // a boost in that direction but rescaled to just less than one. 00158 00159 // There is in principle no way that gamma could have become negative, 00160 // but if that happens, we ZMthrow and (if continuing) just rescale, which 00161 // will change the sign of the last column when computing the boost. 00162 00163 double gam = tt(); 00164 if (gam <= 0) { // 4/12/01 mf 00165 // ZMthrowA (ZMxpvTachyonic( 00166 ZMthrowC (ZMxpvTachyonic( 00167 "Attempt to rectify a boost with non-positive gamma.")); 00168 if (gam==0) return; // NaN-proofing 00169 } 00170 Hep3Vector boost (xt(), yt(), zt()); 00171 boost /= tt(); 00172 if ( boost.mag2() >= 1 ) { // NaN-proofing: 00173 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) ); // used to just check > 1 00174 } 00175 set ( boost ); 00176 } 00177 00178 // ---------- Application is all in .icc 00179 00180 // ---------- Operations in the group of 4-Rotations 00181 00182 HepLorentzRotation 00183 HepBoost::matrixMultiplication(const HepRep4x4 & m) const { 00184 HepRep4x4Symmetric r = rep4x4Symmetric(); 00185 return HepLorentzRotation( HepRep4x4 ( 00186 r.xx_*m.xx_ + r.xy_*m.yx_ + r.xz_*m.zx_ + r.xt_*m.tx_, 00187 r.xx_*m.xy_ + r.xy_*m.yy_ + r.xz_*m.zy_ + r.xt_*m.ty_, 00188 r.xx_*m.xz_ + r.xy_*m.yz_ + r.xz_*m.zz_ + r.xt_*m.tz_, 00189 r.xx_*m.xt_ + r.xy_*m.yt_ + r.xz_*m.zt_ + r.xt_*m.tt_, 00190 00191 r.xy_*m.xx_ + r.yy_*m.yx_ + r.yz_*m.zx_ + r.yt_*m.tx_, 00192 r.xy_*m.xy_ + r.yy_*m.yy_ + r.yz_*m.zy_ + r.yt_*m.ty_, 00193 r.xy_*m.xz_ + r.yy_*m.yz_ + r.yz_*m.zz_ + r.yt_*m.tz_, 00194 r.xy_*m.xt_ + r.yy_*m.yt_ + r.yz_*m.zt_ + r.yt_*m.tt_, 00195 00196 r.xz_*m.xx_ + r.yz_*m.yx_ + r.zz_*m.zx_ + r.zt_*m.tx_, 00197 r.xz_*m.xy_ + r.yz_*m.yy_ + r.zz_*m.zy_ + r.zt_*m.ty_, 00198 r.xz_*m.xz_ + r.yz_*m.yz_ + r.zz_*m.zz_ + r.zt_*m.tz_, 00199 r.xz_*m.xt_ + r.yz_*m.yt_ + r.zz_*m.zt_ + r.zt_*m.tt_, 00200 00201 r.xt_*m.xx_ + r.yt_*m.yx_ + r.zt_*m.zx_ + r.tt_*m.tx_, 00202 r.xt_*m.xy_ + r.yt_*m.yy_ + r.zt_*m.zy_ + r.tt_*m.ty_, 00203 r.xt_*m.xz_ + r.yt_*m.yz_ + r.zt_*m.zz_ + r.tt_*m.tz_, 00204 r.xt_*m.xt_ + r.yt_*m.yt_ + r.zt_*m.zt_ + r.tt_*m.tt_) ); 00205 } 00206 00207 HepLorentzRotation 00208 HepBoost::matrixMultiplication(const HepRep4x4Symmetric & m) const { 00209 HepRep4x4Symmetric r = rep4x4Symmetric(); 00210 return HepLorentzRotation( HepRep4x4 ( 00211 r.xx_*m.xx_ + r.xy_*m.xy_ + r.xz_*m.xz_ + r.xt_*m.xt_, 00212 r.xx_*m.xy_ + r.xy_*m.yy_ + r.xz_*m.yz_ + r.xt_*m.yt_, 00213 r.xx_*m.xz_ + r.xy_*m.yz_ + r.xz_*m.zz_ + r.xt_*m.zt_, 00214 r.xx_*m.xt_ + r.xy_*m.yt_ + r.xz_*m.zt_ + r.xt_*m.tt_, 00215 00216 r.xy_*m.xx_ + r.yy_*m.xy_ + r.yz_*m.xz_ + r.yt_*m.xt_, 00217 r.xy_*m.xy_ + r.yy_*m.yy_ + r.yz_*m.yz_ + r.yt_*m.yt_, 00218 r.xy_*m.xz_ + r.yy_*m.yz_ + r.yz_*m.zz_ + r.yt_*m.zt_, 00219 r.xy_*m.xt_ + r.yy_*m.yt_ + r.yz_*m.zt_ + r.yt_*m.tt_, 00220 00221 r.xz_*m.xx_ + r.yz_*m.xy_ + r.zz_*m.xz_ + r.zt_*m.xt_, 00222 r.xz_*m.xy_ + r.yz_*m.yy_ + r.zz_*m.yz_ + r.zt_*m.yt_, 00223 r.xz_*m.xz_ + r.yz_*m.yz_ + r.zz_*m.zz_ + r.zt_*m.zt_, 00224 r.xz_*m.xt_ + r.yz_*m.yt_ + r.zz_*m.zt_ + r.zt_*m.tt_, 00225 00226 r.xt_*m.xx_ + r.yt_*m.xy_ + r.zt_*m.xz_ + r.tt_*m.xt_, 00227 r.xt_*m.xy_ + r.yt_*m.yy_ + r.zt_*m.yz_ + r.tt_*m.yt_, 00228 r.xt_*m.xz_ + r.yt_*m.yz_ + r.zt_*m.zz_ + r.tt_*m.zt_, 00229 r.xt_*m.xt_ + r.yt_*m.yt_ + r.zt_*m.zt_ + r.tt_*m.tt_) ); 00230 } 00231 00232 HepLorentzRotation 00233 HepBoost::operator* (const HepLorentzRotation & lt) const { 00234 return matrixMultiplication(lt.rep4x4()); 00235 } 00236 00237 HepLorentzRotation 00238 HepBoost::operator* (const HepBoost & b) const { 00239 return matrixMultiplication(b.rep_); 00240 } 00241 00242 HepLorentzRotation 00243 HepBoost::operator* (const HepRotation & r) const { 00244 return matrixMultiplication(r.rep4x4()); 00245 } 00246 00247 // ---------- I/O: 00248 00249 std::ostream & HepBoost::print( std::ostream & os ) const { 00250 if ( rep_.tt_ <= 1 ) { 00251 os << "Lorentz Boost( IDENTITY )"; 00252 } else { 00253 double norm = boostVector().mag(); 00254 os << "\nLorentz Boost " << boostVector()/norm << 00255 "\n{beta = " << beta() << " gamma = " << gamma() << "}\n"; 00256 } 00257 return os; 00258 } 00259 00260 } // namespace CLHEP