CLHEP VERSION 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 those parts of the HepLorentzRotation class 00007 // which involve decomposition into Boost*Rotation. 00008 00009 #ifdef GNUPRAGMA 00010 #pragma implementation 00011 #endif 00012 00013 #include "CLHEP/Vector/defs.h" 00014 #include "CLHEP/Vector/LorentzRotation.h" 00015 00016 namespace CLHEP { 00017 00018 // ---------- Decomposition: 00019 00020 void HepLorentzRotation::decompose 00021 (HepBoost & bboost, HepRotation & rotation) const { 00022 00023 // The boost will be the pure boost based on column 4 of the transformation 00024 // matrix. Since the constructor takes the beta vector, and not beta*gamma, 00025 // we first divide through by gamma = the tt element. This of course can 00026 // never be zero since the last row has t**2 - v**2 = +1. 00027 00028 Hep3Vector betaVec ( xt(), yt(), zt() ); 00029 betaVec *= 1.0 / tt(); 00030 bboost.set( betaVec ); 00031 00032 // The rotation will be inverse of B times T. 00033 00034 HepBoost B( -betaVec ); 00035 HepLorentzRotation R( B * *this ); 00036 00037 HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(), 00038 R.yx(), R.yy(), R.yz(), 00039 R.zx(), R.zy(), R.zz() ); 00040 rotation.set( m1 ); 00041 rotation.rectify(); 00042 00043 return; 00044 00045 } 00046 00047 void HepLorentzRotation::decompose 00048 (Hep3Vector & bboost, HepAxisAngle & rotation) const { 00049 HepRotation r; 00050 HepBoost b; 00051 decompose(b,r); 00052 bboost = b.boostVector(); 00053 rotation = r.axisAngle(); 00054 return; 00055 } 00056 00057 void HepLorentzRotation::decompose 00058 (HepRotation & rotation, HepBoost & bboost) const { 00059 00060 // In this case the pure boost is based on row 4 of the matrix. 00061 00062 Hep3Vector betaVec( tx(), ty(), tz() ); 00063 betaVec *= 1.0 / tt(); 00064 bboost.set( betaVec ); 00065 00066 // The rotation will be T times the inverse of B. 00067 00068 HepBoost B( -betaVec ); 00069 HepLorentzRotation R( *this * B ); 00070 00071 HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(), 00072 R.yx(), R.yy(), R.yz(), 00073 R.zx(), R.zy(), R.zz() ); 00074 rotation.set( m1 ); 00075 rotation.rectify(); 00076 return; 00077 00078 } 00079 00080 void HepLorentzRotation::decompose 00081 (HepAxisAngle & rotation, Hep3Vector & bboost) const { 00082 HepRotation r; 00083 HepBoost b; 00084 decompose(r,b); 00085 rotation = r.axisAngle(); 00086 bboost = b.boostVector(); 00087 return; 00088 } 00089 00090 double HepLorentzRotation::distance2( const HepBoost & b ) const { 00091 HepBoost b1; 00092 HepRotation r1; 00093 decompose( b1, r1 ); 00094 double db2 = b1.distance2( b ); 00095 double dr2 = r1.norm2(); 00096 return ( db2 + dr2 ); 00097 } 00098 00099 double HepLorentzRotation::distance2( const HepRotation & r ) const { 00100 HepBoost b1; 00101 HepRotation r1; 00102 decompose( b1, r1 ); 00103 double db2 = b1.norm2( ); 00104 double dr2 = r1.distance2( r ); 00105 return ( db2 + dr2 ); 00106 } 00107 00108 double HepLorentzRotation::distance2( 00109 const HepLorentzRotation & lt ) const { 00110 HepBoost b1; 00111 HepRotation r1; 00112 decompose( b1, r1 ); 00113 HepBoost b2; 00114 HepRotation r2; 00115 lt.decompose (b2, r2); 00116 double db2 = b1.distance2( b2 ); 00117 double dr2 = r1.distance2( r2 ); 00118 return ( db2 + dr2 ); 00119 } 00120 00121 double HepLorentzRotation::howNear( const HepBoost & b ) const { 00122 return std::sqrt( distance2( b ) ); 00123 } 00124 double HepLorentzRotation::howNear( const HepRotation & r ) const { 00125 return std::sqrt( distance2( r ) ); 00126 } 00127 double HepLorentzRotation::howNear( const HepLorentzRotation & lt )const { 00128 return std::sqrt( distance2( lt ) ); 00129 } 00130 00131 bool HepLorentzRotation::isNear( 00132 const HepBoost & b, double epsilon ) const { 00133 HepBoost b1; 00134 HepRotation r1; 00135 decompose( b1, r1 ); 00136 double db2 = b1.distance2(b); 00137 if ( db2 > epsilon*epsilon ) { 00138 return false; // Saves the time-consuming Rotation::norm2 00139 } 00140 double dr2 = r1.norm2(); 00141 return ( (db2 + dr2) <= epsilon*epsilon ); 00142 } 00143 00144 bool HepLorentzRotation::isNear( 00145 const HepRotation & r, double epsilon ) const { 00146 HepBoost b1; 00147 HepRotation r1; 00148 decompose( b1, r1 ); 00149 double db2 = b1.norm2(); 00150 if ( db2 > epsilon*epsilon ) { 00151 return false; // Saves the time-consuming Rotation::distance2 00152 } 00153 double dr2 = r1.distance2(r); 00154 return ( (db2 + dr2) <= epsilon*epsilon ); 00155 } 00156 00157 bool HepLorentzRotation::isNear( 00158 const HepLorentzRotation & lt, double epsilon ) const { 00159 HepBoost b1; 00160 HepRotation r1; 00161 decompose( b1, r1 ); 00162 HepBoost b2; 00163 HepRotation r2; 00164 lt.decompose (b2, r2); 00165 double db2 = b1.distance2(b2); 00166 if ( db2 > epsilon*epsilon ) { 00167 return false; // Saves the time-consuming Rotation::distance2 00168 } 00169 double dr2 = r1.distance2(r2); 00170 return ( (db2 + dr2) <= epsilon*epsilon ); 00171 } 00172 00173 double HepLorentzRotation::norm2() const { 00174 HepBoost b; 00175 HepRotation r; 00176 decompose( b, r ); 00177 return b.norm2() + r.norm2(); 00178 } 00179 00180 void HepLorentzRotation::rectify() { 00181 00182 // Assuming the representation of this is close to a true LT, 00183 // but may have drifted due to round-off error from many operations, 00184 // this forms an "exact" orthosymplectic matrix for the LT again. 00185 00186 // There are several ways to do this, all equivalent to lowest order in 00187 // the corrected error. We choose to form an LT based on the inverse boost 00188 // extracted from row 4, and left-multiply by LT to form what would be 00189 // a rotation if the LT were kosher. We drop the possibly non-zero t 00190 // components of that, rectify that rotation and multiply back by the boost. 00191 00192 Hep3Vector beta (tx(), ty(), tz()); 00193 double gam = tt(); // NaN-proofing 00194 if ( gam <= 0 ) { 00195 ZMthrowA ( ZMxpvImproperTransformation ( 00196 "rectify() on a transformation with tt() <= 0 - will not help!" )); 00197 gam = 1; 00198 } 00199 beta *= 1.0/gam; 00200 HepLorentzRotation R = (*this) * HepBoost(-beta); 00201 00202 HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(), 00203 R.yx(), R.yy(), R.yz(), 00204 R.zx(), R.zy(), R.zz() ); 00205 00206 HepRotation Rgood (m1); 00207 Rgood.rectify(); 00208 00209 set ( Rgood, HepBoost(beta) ); 00210 } 00211 00212 } // namespace CLHEP