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 that part of 6 // This is the implementation of that part of the HepLorentzRotation class 7 // which is concerned with setting or construc 7 // which is concerned with setting or constructing the transformation based 8 // on 4 supplied columns or rows. 8 // on 4 supplied columns or rows. 9 9 >> 10 #ifdef GNUPRAGMA >> 11 #pragma implementation >> 12 #endif >> 13 10 #include "CLHEP/Vector/LorentzRotation.h" 14 #include "CLHEP/Vector/LorentzRotation.h" 11 #include "CLHEP/Vector/LorentzVector.h" 15 #include "CLHEP/Vector/LorentzVector.h" 12 16 13 #include <cmath> 17 #include <cmath> 14 #include <iostream> << 15 18 16 namespace CLHEP { 19 namespace CLHEP { 17 20 18 // ---------- Constructors and Assignment: 21 // ---------- Constructors and Assignment: 19 22 20 HepLorentzRotation & HepLorentzRotation::set ( 23 HepLorentzRotation & HepLorentzRotation::set (const HepLorentzVector & ccol1, 21 const HepLor 24 const HepLorentzVector & ccol2, 22 const HepLorentzVe 25 const HepLorentzVector & ccol3, 23 const HepLorentzVe 26 const HepLorentzVector & ccol4) { 24 // First, test that the four cols do represe 27 // First, test that the four cols do represent something close to a 25 // true LT: 28 // true LT: 26 29 27 ZMpvMetric_t savedMetric = HepLorentzVector: 30 ZMpvMetric_t savedMetric = HepLorentzVector::setMetric (TimePositive); 28 31 29 if ( ccol4.getT() < 0 ) { 32 if ( ccol4.getT() < 0 ) { 30 std::cerr << "HepLorentzRotation::set() - 33 std::cerr << "HepLorentzRotation::set() - " 31 << "column 4 supplied to define transfor 34 << "column 4 supplied to define transformation has negative T component" 32 << std::endl; 35 << std::endl; 33 *this = HepLorentzRotation(); 36 *this = HepLorentzRotation(); 34 return *this; 37 return *this; 35 } 38 } 36 /* 39 /* 37 double u1u1 = ccol1.dot(ccol1); 40 double u1u1 = ccol1.dot(ccol1); 38 double f11 = std::fabs(u1u1 + 1.0); 41 double f11 = std::fabs(u1u1 + 1.0); 39 if ( f11 > Hep4RotationInterface::tolerance 42 if ( f11 > Hep4RotationInterface::tolerance ) { 40 std::cerr << "HepLorentzRotation::set() - 43 std::cerr << "HepLorentzRotation::set() - " 41 << "column 1 supplied for HepLorentzRota 44 << "column 1 supplied for HepLorentzRotation has w*w != -1" << std::endl; 42 } 45 } 43 double u2u2 = ccol2.dot(ccol2); 46 double u2u2 = ccol2.dot(ccol2); 44 double f22 = std::fabs(u2u2 + 1.0); 47 double f22 = std::fabs(u2u2 + 1.0); 45 if ( f22 > Hep4RotationInterface::tolerance 48 if ( f22 > Hep4RotationInterface::tolerance ) { 46 std::cerr << "HepLorentzRotation::set() - 49 std::cerr << "HepLorentzRotation::set() - " 47 << "column 2 supplied for HepLorentzRota 50 << "column 2 supplied for HepLorentzRotation has w*w != -1" << std::endl; 48 } 51 } 49 double u3u3 = ccol3.dot(ccol3); 52 double u3u3 = ccol3.dot(ccol3); 50 double f33 = std::fabs(u3u3 + 1.0); 53 double f33 = std::fabs(u3u3 + 1.0); 51 if ( f33 > Hep4RotationInterface::tolerance 54 if ( f33 > Hep4RotationInterface::tolerance ) { 52 std::cerr << "HepLorentzRotation::set() - 55 std::cerr << "HepLorentzRotation::set() - " 53 << "column 3 supplied for HepLorentzRota 56 << "column 3 supplied for HepLorentzRotation has w*w != -1" << std::endl; 54 } 57 } 55 double u4u4 = ccol4.dot(ccol4); 58 double u4u4 = ccol4.dot(ccol4); 56 double f44 = std::fabs(u4u4 - 1.0); 59 double f44 = std::fabs(u4u4 - 1.0); 57 if ( f44 > Hep4RotationInterface::tolerance 60 if ( f44 > Hep4RotationInterface::tolerance ) { 58 std::cerr << "HepLorentzRotation::set() - 61 std::cerr << "HepLorentzRotation::set() - " 59 << "column 4 supplied for HepLorentzRota 62 << "column 4 supplied for HepLorentzRotation has w*w != +1" << std::endl; 60 } 63 } 61 64 62 double u1u2 = ccol1.dot(ccol2); 65 double u1u2 = ccol1.dot(ccol2); 63 double f12 = std::fabs(u1u2); 66 double f12 = std::fabs(u1u2); 64 if ( f12 > Hep4RotationInterface::tolerance 67 if ( f12 > Hep4RotationInterface::tolerance ) { 65 std::cerr << "HepLorentzRotation::set() - 68 std::cerr << "HepLorentzRotation::set() - " 66 << "columns 1 and 2 supplied for HepLore 69 << "columns 1 and 2 supplied for HepLorentzRotation have non-zero dot" << std::endl; 67 } 70 } 68 double u1u3 = ccol1.dot(ccol3); 71 double u1u3 = ccol1.dot(ccol3); 69 double f13 = std::fabs(u1u3); 72 double f13 = std::fabs(u1u3); 70 73 71 if ( f13 > Hep4RotationInterface::tolerance 74 if ( f13 > Hep4RotationInterface::tolerance ) { 72 std::cerr << "HepLorentzRotation::set() - 75 std::cerr << "HepLorentzRotation::set() - " 73 << "columns 1 and 3 supplied for HepLore 76 << "columns 1 and 3 supplied for HepLorentzRotation have non-zero dot" << std::endl; 74 } 77 } 75 double u1u4 = ccol1.dot(ccol4); 78 double u1u4 = ccol1.dot(ccol4); 76 double f14 = std::fabs(u1u4); 79 double f14 = std::fabs(u1u4); 77 if ( f14 > Hep4RotationInterface::tolerance 80 if ( f14 > Hep4RotationInterface::tolerance ) { 78 std::cerr << "HepLorentzRotation::set() - 81 std::cerr << "HepLorentzRotation::set() - " 79 << "columns 1 and 4 supplied for HepLore 82 << "columns 1 and 4 supplied for HepLorentzRotation have non-zero dot" << std::endl; 80 } 83 } 81 double u2u3 = ccol2.dot(ccol3); 84 double u2u3 = ccol2.dot(ccol3); 82 double f23 = std::fabs(u2u3); 85 double f23 = std::fabs(u2u3); 83 if ( f23 > Hep4RotationInterface::tolerance 86 if ( f23 > Hep4RotationInterface::tolerance ) { 84 std::cerr << "HepLorentzRotation::set() - 87 std::cerr << "HepLorentzRotation::set() - " 85 << "columns 2 and 3 supplied for HepLore 88 << "columns 2 and 3 supplied for HepLorentzRotation have non-zero dot" << std::endl; 86 } 89 } 87 double u2u4 = ccol2.dot(ccol4); 90 double u2u4 = ccol2.dot(ccol4); 88 double f24 = std::fabs(u2u4); 91 double f24 = std::fabs(u2u4); 89 if ( f24 > Hep4RotationInterface::tolerance 92 if ( f24 > Hep4RotationInterface::tolerance ) { 90 std::cerr << "HepLorentzRotation::set() - 93 std::cerr << "HepLorentzRotation::set() - " 91 << "columns 2 and 4 supplied for HepLore 94 << "columns 2 and 4 supplied for HepLorentzRotation have non-zero dot" << std::endl; 92 } 95 } 93 double u3u4 = ccol3.dot(ccol4); 96 double u3u4 = ccol3.dot(ccol4); 94 double f34 = std::fabs(u3u4); 97 double f34 = std::fabs(u3u4); 95 if ( f34 > Hep4RotationInterface::tolerance 98 if ( f34 > Hep4RotationInterface::tolerance ) { 96 std::cerr << "HepLorentzRotation::set() - 99 std::cerr << "HepLorentzRotation::set() - " 97 << "columns 3 and 4 supplied for HepLore 100 << "columns 3 and 4 supplied for HepLorentzRotation have non-zero dot" << std::endl; 98 } 101 } 99 */ 102 */ 100 // Our strategy will be to order the cols, t 103 // Our strategy will be to order the cols, then do gram-schmidt on them 101 // (that is, remove the components of col d 104 // (that is, remove the components of col d that make it non-orthogonal to 102 // col c, normalize that, then remove the co 105 // col c, normalize that, then remove the components of b that make it 103 // non-orthogonal to d and to c, normalize t 106 // non-orthogonal to d and to c, normalize that, etc. 104 107 105 // Because col4, the time col, is most likel 108 // Because col4, the time col, is most likely to be computed directly, we 106 // will start from there and work left-ward. 109 // will start from there and work left-ward. 107 110 108 HepLorentzVector a, b, c, d; 111 HepLorentzVector a, b, c, d; 109 bool isLorentzTransformation = true; 112 bool isLorentzTransformation = true; 110 double norm; 113 double norm; 111 114 112 d = ccol4; 115 d = ccol4; 113 norm = d.dot(d); 116 norm = d.dot(d); 114 if (norm <= 0.0) { 117 if (norm <= 0.0) { 115 isLorentzTransformation = false; 118 isLorentzTransformation = false; 116 if (norm == 0.0) { 119 if (norm == 0.0) { 117 d = T_HAT4; // Moot, but let's kee 120 d = T_HAT4; // Moot, but let's keep going... 118 norm = 1.0; 121 norm = 1.0; 119 } 122 } 120 } 123 } 121 d /= norm; 124 d /= norm; 122 125 123 c = ccol3 - ccol3.dot(d) * d; 126 c = ccol3 - ccol3.dot(d) * d; 124 norm = -c.dot(c); 127 norm = -c.dot(c); 125 if (norm <= 0.0) { 128 if (norm <= 0.0) { 126 isLorentzTransformation = false; 129 isLorentzTransformation = false; 127 if (norm == 0.0) { 130 if (norm == 0.0) { 128 c = Z_HAT4; // Moot 131 c = Z_HAT4; // Moot 129 norm = 1.0; 132 norm = 1.0; 130 } 133 } 131 } 134 } 132 c /= norm; 135 c /= norm; 133 136 134 b = ccol2 + ccol2.dot(c) * c - ccol2.dot(d) 137 b = ccol2 + ccol2.dot(c) * c - ccol2.dot(d) * d; 135 norm = -b.dot(b); 138 norm = -b.dot(b); 136 if (norm <= 0.0) { 139 if (norm <= 0.0) { 137 isLorentzTransformation = false; 140 isLorentzTransformation = false; 138 if (norm == 0.0) { 141 if (norm == 0.0) { 139 b = Y_HAT4; // Moot 142 b = Y_HAT4; // Moot 140 norm = 1.0; 143 norm = 1.0; 141 } 144 } 142 } 145 } 143 b /= norm; 146 b /= norm; 144 147 145 a = ccol1 + ccol1.dot(b) * b + ccol1.dot(c) 148 a = ccol1 + ccol1.dot(b) * b + ccol1.dot(c) * c - ccol1.dot(d) * d; 146 norm = -a.dot(a); 149 norm = -a.dot(a); 147 if (norm <= 0.0) { 150 if (norm <= 0.0) { 148 isLorentzTransformation = false; 151 isLorentzTransformation = false; 149 if (norm == 0.0) { 152 if (norm == 0.0) { 150 a = X_HAT4; // Moot 153 a = X_HAT4; // Moot 151 norm = 1.0; 154 norm = 1.0; 152 } 155 } 153 } 156 } 154 a /= norm; 157 a /= norm; 155 158 156 if ( !isLorentzTransformation ) { 159 if ( !isLorentzTransformation ) { 157 std::cerr << "HepLorentzRotation::set() 160 std::cerr << "HepLorentzRotation::set() - " 158 << "cols 1-4 supplied to define transf 161 << "cols 1-4 supplied to define transformation form either \n" 159 << " a boosted reflection or a t 162 << " a boosted reflection or a tachyonic transformation -- \n" 160 << " transformation will be set 163 << " transformation will be set to Identity " << std::endl; 161 164 162 *this = HepLorentzRotation(); 165 *this = HepLorentzRotation(); 163 } 166 } 164 167 165 if ( isLorentzTransformation ) { 168 if ( isLorentzTransformation ) { 166 mxx = a.x(); myx = a.y(); mzx = a.z(); mtx 169 mxx = a.x(); myx = a.y(); mzx = a.z(); mtx = a.t(); 167 mxy = b.x(); myy = b.y(); mzy = b.z(); mty 170 mxy = b.x(); myy = b.y(); mzy = b.z(); mty = b.t(); 168 mxz = c.x(); myz = c.y(); mzz = c.z(); mtz 171 mxz = c.x(); myz = c.y(); mzz = c.z(); mtz = c.t(); 169 mxt = d.x(); myt = d.y(); mzt = d.z(); mtt 172 mxt = d.x(); myt = d.y(); mzt = d.z(); mtt = d.t(); 170 } 173 } 171 174 172 HepLorentzVector::setMetric (savedMetric); 175 HepLorentzVector::setMetric (savedMetric); 173 return *this; 176 return *this; 174 177 175 } // set ( col1, col2, col3, col4 ) 178 } // set ( col1, col2, col3, col4 ) 176 179 177 HepLorentzRotation & HepLorentzRotation::setRo 180 HepLorentzRotation & HepLorentzRotation::setRows 178 (const HepLorentzVector & rrow1, 181 (const HepLorentzVector & rrow1, 179 const HepLorentzVector & rrow2, 182 const HepLorentzVector & rrow2, 180 const HepLorentzVector & rrow3, 183 const HepLorentzVector & rrow3, 181 const HepLorentzVector & rrow4) { 184 const HepLorentzVector & rrow4) { 182 // Set based on using those rows as columns: 185 // Set based on using those rows as columns: 183 set (rrow1, rrow2, rrow3, rrow4); 186 set (rrow1, rrow2, rrow3, rrow4); 184 // Now transpose in place: 187 // Now transpose in place: 185 double q1, q2, q3; 188 double q1, q2, q3; 186 q1 = mxy; q2 = mxz; q3 = mxt; 189 q1 = mxy; q2 = mxz; q3 = mxt; 187 mxy = myx; mxz = mzx; mxt = mtx; 190 mxy = myx; mxz = mzx; mxt = mtx; 188 myx = q1; mzx = q2; mtx = q3; 191 myx = q1; mzx = q2; mtx = q3; 189 q1 = myz; q2 = myt; q3 = mzt; 192 q1 = myz; q2 = myt; q3 = mzt; 190 myz = mzy; myt = mty; mzt = mtz; 193 myz = mzy; myt = mty; mzt = mtz; 191 mzy = q1; mty = q2; mtz = q3; 194 mzy = q1; mty = q2; mtz = q3; 192 return *this; 195 return *this; 193 } // LorentzTransformation::setRows(row1 ... r 196 } // LorentzTransformation::setRows(row1 ... row4) 194 197 195 HepLorentzRotation::HepLorentzRotation ( const 198 HepLorentzRotation::HepLorentzRotation ( const HepLorentzVector & ccol1, 196 const HepLorentzV 199 const HepLorentzVector & ccol2, 197 const HepLorentzV 200 const HepLorentzVector & ccol3, 198 const HepLorentzV 201 const HepLorentzVector & ccol4 ) 199 { 202 { 200 set ( ccol1, ccol2, ccol3, ccol4 ); 203 set ( ccol1, ccol2, ccol3, ccol4 ); 201 } 204 } 202 205 203 } // namespace CLHEP 206 } // namespace CLHEP 204 207 205 208