Geant4 Cross Reference |
1 // -*- C++ -*- 1 // -*- C++ -*- 2 // ------------------------------------------- 2 // --------------------------------------------------------------------------- 3 // 3 // 4 // This file is a part of the CLHEP - a Class 4 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 5 // 5 // 6 // This is the implementation of the HepBoostX 6 // This is the implementation of the HepBoostX class. 7 // 7 // 8 8 >> 9 #ifdef GNUPRAGMA >> 10 #pragma implementation >> 11 #endif >> 12 9 #include "CLHEP/Vector/BoostX.h" 13 #include "CLHEP/Vector/BoostX.h" 10 #include "CLHEP/Vector/Boost.h" 14 #include "CLHEP/Vector/Boost.h" 11 #include "CLHEP/Vector/Rotation.h" 15 #include "CLHEP/Vector/Rotation.h" 12 #include "CLHEP/Vector/LorentzRotation.h" 16 #include "CLHEP/Vector/LorentzRotation.h" 13 << 14 #include <cmath> << 15 #include <iostream> << 16 17 17 namespace CLHEP { 18 namespace CLHEP { 18 19 19 20 20 // ---------- Constructors and Assignment: 21 // ---------- Constructors and Assignment: 21 22 22 HepBoostX & HepBoostX::set (double bbeta) { 23 HepBoostX & HepBoostX::set (double bbeta) { 23 double b2 = bbeta*bbeta; 24 double b2 = bbeta*bbeta; 24 if (b2 >= 1) { 25 if (b2 >= 1) { 25 std::cerr << "HepBoostX::set() - " 26 std::cerr << "HepBoostX::set() - " 26 << "Beta supplied to set HepBoostX repre 27 << "Beta supplied to set HepBoostX represents speed >= c." << std::endl; 27 beta_ = 1.0 - 1.0E-8; // NaN-proofing 28 beta_ = 1.0 - 1.0E-8; // NaN-proofing 28 gamma_ = 1.0 / std::sqrt(1.0 - b2); 29 gamma_ = 1.0 / std::sqrt(1.0 - b2); 29 return *this; 30 return *this; 30 } 31 } 31 beta_ = bbeta; 32 beta_ = bbeta; 32 gamma_ = 1.0 / std::sqrt(1.0 - b2); 33 gamma_ = 1.0 / std::sqrt(1.0 - b2); 33 return *this; 34 return *this; 34 } 35 } 35 36 36 // ---------- Accessors: 37 // ---------- Accessors: 37 38 38 HepRep4x4 HepBoostX::rep4x4() const { 39 HepRep4x4 HepBoostX::rep4x4() const { 39 double bg = beta_*gamma_; 40 double bg = beta_*gamma_; 40 return HepRep4x4( gamma_, 0, 0, bg, 41 return HepRep4x4( gamma_, 0, 0, bg, 41 0, 1, 0, 0, 42 0, 1, 0, 0, 42 0, 0, 1, 0, 43 0, 0, 1, 0, 43 bg, 0, 0, gamma_ 44 bg, 0, 0, gamma_ ); 44 } 45 } 45 46 46 HepRep4x4Symmetric HepBoostX::rep4x4Symmetric( 47 HepRep4x4Symmetric HepBoostX::rep4x4Symmetric() const { 47 double bg = beta_*gamma_; 48 double bg = beta_*gamma_; 48 return HepRep4x4Symmetric( gamma_, 0, 0 49 return HepRep4x4Symmetric( gamma_, 0, 0, bg, 49 1, 0, 50 1, 0, 0, 50 1, 51 1, 0, 51 gam 52 gamma_ ); 52 } 53 } 53 54 54 // ---------- Decomposition: 55 // ---------- Decomposition: 55 56 56 void HepBoostX::decompose (HepRotation & rotat 57 void HepBoostX::decompose (HepRotation & rotation, HepBoost & boost) const { 57 HepAxisAngle vdelta = HepAxisAngle(); 58 HepAxisAngle vdelta = HepAxisAngle(); 58 rotation = HepRotation(vdelta); 59 rotation = HepRotation(vdelta); 59 Hep3Vector bbeta = boostVector(); 60 Hep3Vector bbeta = boostVector(); 60 boost = HepBoost(bbeta); 61 boost = HepBoost(bbeta); 61 } 62 } 62 63 63 void HepBoostX::decompose (HepAxisAngle & rota 64 void HepBoostX::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const { 64 rotation = HepAxisAngle(); 65 rotation = HepAxisAngle(); 65 boost = boostVector(); 66 boost = boostVector(); 66 } 67 } 67 68 68 void HepBoostX::decompose (HepBoost & boost, H 69 void HepBoostX::decompose (HepBoost & boost, HepRotation & rotation) const { 69 HepAxisAngle vdelta = HepAxisAngle(); 70 HepAxisAngle vdelta = HepAxisAngle(); 70 rotation = HepRotation(vdelta); 71 rotation = HepRotation(vdelta); 71 Hep3Vector bbeta = boostVector(); 72 Hep3Vector bbeta = boostVector(); 72 boost = HepBoost(bbeta); 73 boost = HepBoost(bbeta); 73 } 74 } 74 75 75 void HepBoostX::decompose (Hep3Vector & boost, 76 void HepBoostX::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const { 76 rotation = HepAxisAngle(); 77 rotation = HepAxisAngle(); 77 boost = boostVector(); 78 boost = boostVector(); 78 } 79 } 79 80 80 // ---------- Comparisons: 81 // ---------- Comparisons: 81 82 82 double HepBoostX::distance2( const HepBoost & 83 double HepBoostX::distance2( const HepBoost & b ) const { 83 return b.distance2(*this); 84 return b.distance2(*this); 84 } 85 } 85 86 86 double HepBoostX::distance2( const HepRotation 87 double HepBoostX::distance2( const HepRotation & r ) const { 87 double db2 = norm2(); 88 double db2 = norm2(); 88 double dr2 = r.norm2(); 89 double dr2 = r.norm2(); 89 return (db2 + dr2); 90 return (db2 + dr2); 90 } 91 } 91 92 92 double HepBoostX::distance2( const HepLorentzR 93 double HepBoostX::distance2( const HepLorentzRotation & lt ) const { 93 HepBoost b1; 94 HepBoost b1; 94 HepRotation r1; 95 HepRotation r1; 95 lt.decompose(b1,r1); 96 lt.decompose(b1,r1); 96 double db2 = distance2(b1); 97 double db2 = distance2(b1); 97 double dr2 = r1.norm2(); 98 double dr2 = r1.norm2(); 98 return (db2 + dr2); 99 return (db2 + dr2); 99 } 100 } 100 101 101 bool HepBoostX::isNear (const HepRotation & r, 102 bool HepBoostX::isNear (const HepRotation & r, double epsilon) const { 102 double db2 = norm2(); 103 double db2 = norm2(); 103 if (db2 > epsilon*epsilon) return false; 104 if (db2 > epsilon*epsilon) return false; 104 double dr2 = r.norm2(); 105 double dr2 = r.norm2(); 105 return (db2+dr2 <= epsilon*epsilon); 106 return (db2+dr2 <= epsilon*epsilon); 106 } 107 } 107 108 108 bool HepBoostX::isNear ( const HepLorentzRotat 109 bool HepBoostX::isNear ( const HepLorentzRotation & lt, 109 double epsilon) const { 110 double epsilon) const { 110 HepBoost b1; 111 HepBoost b1; 111 HepRotation r1; 112 HepRotation r1; 112 double db2 = distance2(b1); 113 double db2 = distance2(b1); 113 lt.decompose(b1,r1); 114 lt.decompose(b1,r1); 114 if (db2 > epsilon*epsilon) return false; 115 if (db2 > epsilon*epsilon) return false; 115 double dr2 = r1.norm2(); 116 double dr2 = r1.norm2(); 116 return (db2 + dr2); 117 return (db2 + dr2); 117 } 118 } 118 119 119 // ---------- Properties: 120 // ---------- Properties: 120 121 121 void HepBoostX::rectify() { 122 void HepBoostX::rectify() { 122 // Assuming the representation of this is cl 123 // Assuming the representation of this is close to a true pure boost, 123 // but may have drifted due to round-off err 124 // but may have drifted due to round-off error from many operations, 124 // this forms an "exact" pure boostX matrix 125 // this forms an "exact" pure boostX matrix for again. 125 126 126 double b2 = beta_*beta_; 127 double b2 = beta_*beta_; 127 if (b2 >= 1) { 128 if (b2 >= 1) { 128 beta_ = 1.0 - 1.0e-8; // NaN-proofing 129 beta_ = 1.0 - 1.0e-8; // NaN-proofing 129 b2 = beta_*beta_; 130 b2 = beta_*beta_; 130 } 131 } 131 gamma_ = 1.0 / std::sqrt(1.0 - b2); 132 gamma_ = 1.0 / std::sqrt(1.0 - b2); 132 } 133 } 133 134 134 // ---------- Application: 135 // ---------- Application: 135 136 136 // ---------- Operations in the group of 4-Rot 137 // ---------- Operations in the group of 4-Rotations 137 138 138 HepBoostX HepBoostX::operator * (const HepBoos 139 HepBoostX HepBoostX::operator * (const HepBoostX & b) const { 139 return HepBoostX ( (beta()+b.beta()) / (1+be 140 return HepBoostX ( (beta()+b.beta()) / (1+beta()*b.beta()) ); 140 } 141 } 141 HepLorentzRotation HepBoostX::operator * (cons 142 HepLorentzRotation HepBoostX::operator * (const HepBoost & b) const { 142 HepLorentzRotation me (*this); 143 HepLorentzRotation me (*this); 143 return me*b; 144 return me*b; 144 } 145 } 145 HepLorentzRotation HepBoostX::operator * (cons 146 HepLorentzRotation HepBoostX::operator * (const HepRotation & r) const { 146 HepLorentzRotation me (*this); 147 HepLorentzRotation me (*this); 147 return me*r; 148 return me*r; 148 } 149 } 149 HepLorentzRotation HepBoostX::operator * (cons 150 HepLorentzRotation HepBoostX::operator * (const HepLorentzRotation & lt) const { 150 HepLorentzRotation me (*this); 151 HepLorentzRotation me (*this); 151 return me*lt; 152 return me*lt; 152 } 153 } 153 154 154 // ---------- I/O: 155 // ---------- I/O: 155 156 156 std::ostream & HepBoostX::print( std::ostream 157 std::ostream & HepBoostX::print( std::ostream & os ) const { 157 os << "Boost in X direction (beta = " << bet 158 os << "Boost in X direction (beta = " << beta_ 158 << ", gamma = " << gamma_ << ") "; 159 << ", gamma = " << gamma_ << ") "; 159 return os; 160 return os; 160 } 161 } 161 162 162 } // namespace CLHEP 163 } // namespace CLHEP 163 164