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