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 that part of the HepLorentzRotation class 00007 // which is concerned with setting or constructing the transformation based 00008 // on 4 supplied columns or rows. 00009 00010 #ifdef GNUPRAGMA 00011 #pragma implementation 00012 #endif 00013 00014 #include "CLHEP/Vector/defs.h" 00015 #include "CLHEP/Vector/LorentzRotation.h" 00016 #include "CLHEP/Vector/LorentzVector.h" 00017 #include "CLHEP/Vector/ZMxpv.h" 00018 00019 #include <cmath> 00020 00021 namespace CLHEP { 00022 00023 // ---------- Constructors and Assignment: 00024 00025 HepLorentzRotation & HepLorentzRotation::set (const HepLorentzVector & ccol1, 00026 const HepLorentzVector & ccol2, 00027 const HepLorentzVector & ccol3, 00028 const HepLorentzVector & ccol4) { 00029 // First, test that the four cols do represent something close to a 00030 // true LT: 00031 00032 ZMpvMetric_t savedMetric = HepLorentzVector::setMetric (TimePositive); 00033 00034 if ( ccol4.getT() < 0 ) { 00035 ZMthrowC (ZMxpvImproperTransformation( 00036 "column 4 supplied to define transformation has negative T component")); 00037 *this = HepLorentzRotation(); 00038 return *this; 00039 } 00040 00041 double u1u1 = ccol1.dot(ccol1); 00042 double f11 = std::fabs(u1u1 + 1.0); 00043 if ( f11 > Hep4RotationInterface::tolerance ) { 00044 ZMthrowC (ZMxpvNotSymplectic( 00045 "column 1 supplied for HepLorentzRotation has w*w != -1")); 00046 } 00047 double u2u2 = ccol2.dot(ccol2); 00048 double f22 = std::fabs(u2u2 + 1.0); 00049 if ( f22 > Hep4RotationInterface::tolerance ) { 00050 ZMthrowC (ZMxpvNotSymplectic( 00051 "column 2 supplied for HepLorentzRotation has w*w != -1")); 00052 } 00053 double u3u3 = ccol3.dot(ccol3); 00054 double f33 = std::fabs(u3u3 + 1.0); 00055 if ( f33 > Hep4RotationInterface::tolerance ) { 00056 ZMthrowC (ZMxpvNotSymplectic( 00057 "column 3 supplied for HepLorentzRotation has w*w != -1")); 00058 } 00059 double u4u4 = ccol4.dot(ccol4); 00060 double f44 = std::fabs(u4u4 - 1.0); 00061 if ( f44 > Hep4RotationInterface::tolerance ) { 00062 ZMthrowC (ZMxpvNotSymplectic( 00063 "column 4 supplied for HepLorentzRotation has w*w != +1")); 00064 } 00065 00066 double u1u2 = ccol1.dot(ccol2); 00067 double f12 = std::fabs(u1u2); 00068 if ( f12 > Hep4RotationInterface::tolerance ) { 00069 ZMthrowC (ZMxpvNotOrthogonal( 00070 "columns 1 and 2 supplied for HepLorentzRotation have non-zero dot")); 00071 } 00072 double u1u3 = ccol1.dot(ccol3); 00073 double f13 = std::fabs(u1u3); 00074 00075 if ( f13 > Hep4RotationInterface::tolerance ) { 00076 ZMthrowC (ZMxpvNotOrthogonal( 00077 "columns 1 and 3 supplied for HepLorentzRotation have non-zero dot")); 00078 } 00079 double u1u4 = ccol1.dot(ccol4); 00080 double f14 = std::fabs(u1u4); 00081 if ( f14 > Hep4RotationInterface::tolerance ) { 00082 ZMthrowC (ZMxpvNotOrthogonal( 00083 "columns 1 and 4 supplied for HepLorentzRotation have non-zero dot")); 00084 } 00085 double u2u3 = ccol2.dot(ccol3); 00086 double f23 = std::fabs(u2u3); 00087 if ( f23 > Hep4RotationInterface::tolerance ) { 00088 ZMthrowC (ZMxpvNotOrthogonal( 00089 "columns 2 and 3 supplied for HepLorentzRotation have non-zero dot")); 00090 } 00091 double u2u4 = ccol2.dot(ccol4); 00092 double f24 = std::fabs(u2u4); 00093 if ( f24 > Hep4RotationInterface::tolerance ) { 00094 ZMthrowC (ZMxpvNotOrthogonal( 00095 "columns 2 and 4 supplied for HepLorentzRotation have non-zero dot")); 00096 } 00097 double u3u4 = ccol3.dot(ccol4); 00098 double f34 = std::fabs(u3u4); 00099 if ( f34 > Hep4RotationInterface::tolerance ) { 00100 ZMthrowC (ZMxpvNotOrthogonal( 00101 "columns 3 and 4 supplied for HepLorentzRotation have non-zero dot")); 00102 } 00103 00104 // Our strategy will be to order the cols, then do gram-schmidt on them 00105 // (that is, remove the components of col d that make it non-orthogonal to 00106 // col c, normalize that, then remove the components of b that make it 00107 // non-orthogonal to d and to c, normalize that, etc. 00108 00109 // Because col4, the time col, is most likely to be computed directly, we 00110 // will start from there and work left-ward. 00111 00112 HepLorentzVector a, b, c, d; 00113 bool isLorentzTransformation = true; 00114 double norm; 00115 00116 d = ccol4; 00117 norm = d.dot(d); 00118 if (norm <= 0.0) { 00119 isLorentzTransformation = false; 00120 if (norm == 0.0) { 00121 d = T_HAT4; // Moot, but let's keep going... 00122 00123 norm = 1.0; 00124 } 00125 } 00126 d /= norm; 00127 00128 c = ccol3 - ccol3.dot(d) * d; 00129 norm = -c.dot(c); 00130 if (norm <= 0.0) { 00131 isLorentzTransformation = false; 00132 if (norm == 0.0) { 00133 c = Z_HAT4; // Moot 00134 norm = 1.0; 00135 } 00136 } 00137 c /= norm; 00138 00139 b = ccol2 + ccol2.dot(c) * c - ccol2.dot(d) * d; 00140 norm = -b.dot(b); 00141 if (norm <= 0.0) { 00142 isLorentzTransformation = false; 00143 if (norm == 0.0) { 00144 b = Y_HAT4; // Moot 00145 norm = 1.0; 00146 } 00147 } 00148 b /= norm; 00149 00150 a = ccol1 + ccol1.dot(b) * b + ccol1.dot(c) * c - ccol1.dot(d) * d; 00151 norm = -a.dot(a); 00152 if (norm <= 0.0) { 00153 isLorentzTransformation = false; 00154 if (norm == 0.0) { 00155 a = X_HAT4; // Moot 00156 norm = 1.0; 00157 } 00158 } 00159 a /= norm; 00160 00161 if ( !isLorentzTransformation ) { 00162 ZMthrowC (ZMxpvImproperTransformation( 00163 "cols 1-4 supplied to define transformation form either \n" 00164 " a boosted reflection or a tachyonic transformation -- \n" 00165 " transformation will be set to Identity ")); 00166 00167 00168 *this = HepLorentzRotation(); 00169 } 00170 00171 if ( isLorentzTransformation ) { 00172 mxx = a.x(); myx = a.y(); mzx = a.z(); mtx = a.t(); 00173 mxy = b.x(); myy = b.y(); mzy = b.z(); mty = b.t(); 00174 mxz = c.x(); myz = c.y(); mzz = c.z(); mtz = c.t(); 00175 mxt = d.x(); myt = d.y(); mzt = d.z(); mtt = d.t(); 00176 } 00177 00178 HepLorentzVector::setMetric (savedMetric); 00179 return *this; 00180 00181 } // set ( col1, col2, col3, col4 ) 00182 00183 HepLorentzRotation & HepLorentzRotation::setRows 00184 (const HepLorentzVector & rrow1, 00185 const HepLorentzVector & rrow2, 00186 const HepLorentzVector & rrow3, 00187 const HepLorentzVector & rrow4) { 00188 // Set based on using those rows as columns: 00189 set (rrow1, rrow2, rrow3, rrow4); 00190 // Now transpose in place: 00191 register double q1, q2, q3; 00192 q1 = mxy; q2 = mxz; q3 = mxt; 00193 mxy = myx; mxz = mzx; mxt = mtx; 00194 myx = q1; mzx = q2; mtx = q3; 00195 q1 = myz; q2 = myt; q3 = mzt; 00196 myz = mzy; myt = mty; mzt = mtz; 00197 mzy = q1; mty = q2; mtz = q3; 00198 return *this; 00199 } // LorentzTransformation::setRows(row1 ... row4) 00200 00201 HepLorentzRotation::HepLorentzRotation ( const HepLorentzVector & ccol1, 00202 const HepLorentzVector & ccol2, 00203 const HepLorentzVector & ccol3, 00204 const HepLorentzVector & ccol4 ) 00205 { 00206 set ( ccol1, ccol2, ccol3, ccol4 ); 00207 } 00208 00209 } // namespace CLHEP 00210