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