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