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