Geant4 Cross Reference |
1 // -*- C++ -*- 2 // --------------------------------------------------------------------------- 3 // 4 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 5 // 6 // This is the implementation of the HepBoost class. 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, double bz) { 21 double bp2 = bx*bx + by*by + bz*bz; 22 // if (bp2 >= 1) { 23 // std::cerr << "HepBoost::set() - " 24 // << "Boost Vector supplied to set HepBoost represents speed >= c." << std::endl; 25 // } 26 double ggamma = 1.0 / std::sqrt(1.0 - bp2); 27 double bgamma = ggamma * ggamma / (1.0 + ggamma); 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 HepRep4x4Symmetric & m1) { 42 rep_ = m1; 43 return *this; 44 } 45 46 HepBoost & HepBoost::set (Hep3Vector ddirection, double bbeta) { 47 double length = ddirection.mag(); 48 if (length <= 0) { // Nan-proofing 49 std::cerr << "HepBoost::set() - " 50 << "Direction supplied to set HepBoost is zero." << std::endl; 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 & bboost) { 61 return set (bboost.x(), bboost.y(), bboost.z()); 62 } 63 64 // ---------- Accessors: 65 66 // ---------- Decomposition: 67 68 void HepBoost::decompose (HepRotation & rotation, HepBoost & boost) const { 69 HepAxisAngle vdelta = HepAxisAngle(); 70 rotation = HepRotation(vdelta); 71 Hep3Vector bbeta = boostVector(); 72 boost = HepBoost(bbeta); 73 } 74 75 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const { 76 rotation = HepAxisAngle(); 77 boost = boostVector(); 78 } 79 80 void HepBoost::decompose (HepBoost & boost, HepRotation & rotation) const { 81 HepAxisAngle vdelta = HepAxisAngle(); 82 rotation = HepRotation(vdelta); 83 Hep3Vector bbeta = boostVector(); 84 boost = HepBoost(bbeta); 85 } 86 87 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const { 88 rotation = HepAxisAngle(); 89 boost = boostVector(); 90 } 91 92 // ---------- Comparisons: 93 94 double HepBoost::distance2( const HepRotation & r ) const { 95 double db2 = norm2(); 96 double dr2 = r.norm2(); 97 return (db2 + dr2); 98 } 99 100 double HepBoost::distance2( const HepLorentzRotation & lt ) const { 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 & r ) const { 110 return std::sqrt(distance2(r)); 111 } 112 113 double HepBoost::howNear ( const HepLorentzRotation & lt ) const { 114 return std::sqrt(distance2(lt)); 115 } 116 117 bool HepBoost::isNear (const HepRotation & r, double epsilon) const { 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 HepLorentzRotation & lt, 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 close to a true pure boost, 146 // but may have drifted due to round-off error from many operations, 147 // this forms an "exact" pure boost matrix for the LT again. 148 149 // The natural way to do this is to use the t column as a boost and set 150 // based on that boost vector. 151 152 // There is perhaps danger that this boost vector will appear equal to or 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 just less than one. 155 156 // There is in principle no way that gamma could have become negative, 157 // but if that happens, we ZMthrow and (if continuing) just rescale, which 158 // will change the sign of the last column when computing the boost. 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-positive gamma." << std::endl; 164 if (gam==0) return; // NaN-proofing 165 } 166 Hep3Vector boost (xt(), yt(), zt()); 167 boost /= tt(); 168 if ( boost.mag2() >= 1 ) { // NaN-proofing: 169 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) ); // used to just check > 1 170 } 171 set ( boost ); 172 } 173 174 // ---------- Application is all in .icc 175 176 // ---------- Operations in the group of 4-Rotations 177 178 HepLorentzRotation 179 HepBoost::matrixMultiplication(const HepRep4x4 & m1) const { 180 HepRep4x4Symmetric r = rep4x4Symmetric(); 181 return HepLorentzRotation( HepRep4x4 ( 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_ + r.xt_*m1.ty_, 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_ + r.xt_*m1.tt_, 186 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_ + r.yt_*m1.ty_, 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_ + r.yt_*m1.tt_, 191 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_ + r.zt_*m1.ty_, 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_ + r.zt_*m1.tt_, 196 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_ + r.tt_*m1.ty_, 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_ + r.tt_*m1.tt_) ); 201 } 202 203 HepLorentzRotation 204 HepBoost::matrixMultiplication(const HepRep4x4Symmetric & m1) const { 205 HepRep4x4Symmetric r = rep4x4Symmetric(); 206 return HepLorentzRotation( HepRep4x4 ( 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_ + r.xt_*m1.yt_, 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_ + r.xt_*m1.tt_, 211 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_ + r.yt_*m1.yt_, 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_ + r.yt_*m1.tt_, 216 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_ + r.zt_*m1.yt_, 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_ + r.zt_*m1.tt_, 221 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_ + r.tt_*m1.yt_, 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_ + r.tt_*m1.tt_) ); 226 } 227 228 HepLorentzRotation 229 HepBoost::operator* (const HepLorentzRotation & lt) const { 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) const { 240 return matrixMultiplication(r.rep4x4()); 241 } 242 243 // ---------- I/O: 244 245 std::ostream & HepBoost::print( std::ostream & os ) const { 246 if ( rep_.tt_ <= 1 ) { 247 os << "Lorentz Boost( IDENTITY )"; 248 } else { 249 double norm = boostVector().mag(); 250 os << "\nLorentz Boost " << boostVector()/norm << 251 "\n{beta = " << beta() << " gamma = " << gamma() << "}\n"; 252 } 253 return os; 254 } 255 256 } // namespace CLHEP 257