Geant4  10.01.p03
LorentzRotation.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id:$
3 // ---------------------------------------------------------------------------
4 //
5 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
6 //
7 // This is the implementation basic parts of the HepLorentzRotation class.
8 //
9 // Some ZOOM methods involving construction from columns and decomposition
10 // into boost*rotation are split off into LorentzRotationC and LorentzRotationD
11 
12 #ifdef GNUPRAGMA
13 #pragma implementation
14 #endif
15 
16 #include "CLHEP/Vector/LorentzRotation.h"
17 
18 #include <iostream>
19 #include <iomanip>
20 
21 namespace CLHEP {
22 
23 // ---------- Constructors and Assignment:
24 
25 
26 HepLorentzRotation & HepLorentzRotation::set
27  (double bx, double by, double bz) {
28  double bp2 = bx*bx + by*by + bz*bz;
29 // if (bp2 >= 1) {
30 // std::cerr << "HepLorentzRotation::set() - "
31 // << "Boost Vector supplied to set HepLorentzRotation represents speed >= c."
32 // << std::endl;
33 // }
34  double gamma = 1.0 / std::sqrt(1.0 - bp2);
35  double bgamma = gamma * gamma / (1.0 + gamma);
36  mxx = 1.0 + bgamma * bx * bx;
37  myy = 1.0 + bgamma * by * by;
38  mzz = 1.0 + bgamma * bz * bz;
39  mxy = myx = bgamma * bx * by;
40  mxz = mzx = bgamma * bx * bz;
41  myz = mzy = bgamma * by * bz;
42  mxt = mtx = gamma * bx;
43  myt = mty = gamma * by;
44  mzt = mtz = gamma * bz;
45  mtt = gamma;
46  return *this;
47 }
48 
49 HepLorentzRotation & HepLorentzRotation::set
50  (const HepBoost & B, const HepRotation & R) {
51  set (B.rep4x4());
52  *this = matrixMultiplication ( R.rep4x4() );
53  return *this;
54 }
55 
56 HepLorentzRotation & HepLorentzRotation::set
57  (const HepRotation & R, const HepBoost & B) {
58  set (R.rep4x4());
59  *this = matrixMultiplication ( B.rep4x4() );
60  return *this;
61 }
62 
63 // ---------- Accessors:
64 
65 // ------------ Subscripting:
66 
67 double HepLorentzRotation::operator () (int i, int j) const {
68  if (i == 0) {
69  if (j == 0) { return xx(); }
70  if (j == 1) { return xy(); }
71  if (j == 2) { return xz(); }
72  if (j == 3) { return xt(); }
73  } else if (i == 1) {
74  if (j == 0) { return yx(); }
75  if (j == 1) { return yy(); }
76  if (j == 2) { return yz(); }
77  if (j == 3) { return yt(); }
78  } else if (i == 2) {
79  if (j == 0) { return zx(); }
80  if (j == 1) { return zy(); }
81  if (j == 2) { return zz(); }
82  if (j == 3) { return zt(); }
83  } else if (i == 3) {
84  if (j == 0) { return tx(); }
85  if (j == 1) { return ty(); }
86  if (j == 2) { return tz(); }
87  if (j == 3) { return tt(); }
88  }
89  std::cerr << "HepLorentzRotation subscripting: bad indeces "
90  << "(" << i << "," << j << ")\n";
91  return 0.0;
92 }
93 
94 // ---------- Application:
95 
96 
97 // ---------- Comparison:
98 
99 int HepLorentzRotation::compare( const HepLorentzRotation & m1 ) const {
100  if (mtt<m1.mtt) return -1; else if (mtt>m1.mtt) return 1;
101  else if (mtz<m1.mtz) return -1; else if (mtz>m1.mtz) return 1;
102  else if (mty<m1.mty) return -1; else if (mty>m1.mty) return 1;
103  else if (mtx<m1.mtx) return -1; else if (mtx>m1.mtx) return 1;
104 
105  else if (mzt<m1.mzt) return -1; else if (mzt>m1.mzt) return 1;
106  else if (mzz<m1.mzz) return -1; else if (mzz>m1.mzz) return 1;
107  else if (mzy<m1.mzy) return -1; else if (mzy>m1.mzy) return 1;
108  else if (mzx<m1.mzx) return -1; else if (mzx>m1.mzx) return 1;
109 
110  else if (myt<m1.myt) return -1; else if (myt>m1.myt) return 1;
111  else if (myz<m1.myz) return -1; else if (myz>m1.myz) return 1;
112  else if (myy<m1.myy) return -1; else if (myy>m1.myy) return 1;
113  else if (myx<m1.myx) return -1; else if (myx>m1.myx) return 1;
114 
115  else if (mxt<m1.mxt) return -1; else if (mxt>m1.mxt) return 1;
116  else if (mxz<m1.mxz) return -1; else if (mxz>m1.mxz) return 1;
117  else if (mxy<m1.mxy) return -1; else if (mxy>m1.mxy) return 1;
118  else if (mxx<m1.mxx) return -1; else if (mxx>m1.mxx) return 1;
119 
120  else return 0;
121 }
122 
123 
124 // ---------- Operations in the group of 4-Rotations
125 
126 HepLorentzRotation
127 HepLorentzRotation::matrixMultiplication(const HepRep4x4 & m1) const {
128  return HepLorentzRotation(
129  mxx*m1.xx_ + mxy*m1.yx_ + mxz*m1.zx_ + mxt*m1.tx_,
130  mxx*m1.xy_ + mxy*m1.yy_ + mxz*m1.zy_ + mxt*m1.ty_,
131  mxx*m1.xz_ + mxy*m1.yz_ + mxz*m1.zz_ + mxt*m1.tz_,
132  mxx*m1.xt_ + mxy*m1.yt_ + mxz*m1.zt_ + mxt*m1.tt_,
133 
134  myx*m1.xx_ + myy*m1.yx_ + myz*m1.zx_ + myt*m1.tx_,
135  myx*m1.xy_ + myy*m1.yy_ + myz*m1.zy_ + myt*m1.ty_,
136  myx*m1.xz_ + myy*m1.yz_ + myz*m1.zz_ + myt*m1.tz_,
137  myx*m1.xt_ + myy*m1.yt_ + myz*m1.zt_ + myt*m1.tt_,
138 
139  mzx*m1.xx_ + mzy*m1.yx_ + mzz*m1.zx_ + mzt*m1.tx_,
140  mzx*m1.xy_ + mzy*m1.yy_ + mzz*m1.zy_ + mzt*m1.ty_,
141  mzx*m1.xz_ + mzy*m1.yz_ + mzz*m1.zz_ + mzt*m1.tz_,
142  mzx*m1.xt_ + mzy*m1.yt_ + mzz*m1.zt_ + mzt*m1.tt_,
143 
144  mtx*m1.xx_ + mty*m1.yx_ + mtz*m1.zx_ + mtt*m1.tx_,
145  mtx*m1.xy_ + mty*m1.yy_ + mtz*m1.zy_ + mtt*m1.ty_,
146  mtx*m1.xz_ + mty*m1.yz_ + mtz*m1.zz_ + mtt*m1.tz_,
147  mtx*m1.xt_ + mty*m1.yt_ + mtz*m1.zt_ + mtt*m1.tt_ );
148 }
149 
150 HepLorentzRotation & HepLorentzRotation::rotateX(double delta) {
151  double c1 = std::cos (delta);
152  double s1 = std::sin (delta);
153  HepLorentzVector rowy = row2();
154  HepLorentzVector rowz = row3();
155  HepLorentzVector r2 = c1 * rowy - s1 * rowz;
156  HepLorentzVector r3 = s1 * rowy + c1 * rowz;
157  myx = r2.x(); myy = r2.y(); myz = r2.z(); myt = r2.t();
158  mzx = r3.x(); mzy = r3.y(); mzz = r3.z(); mzt = r3.t();
159  return *this;
160 }
161 
162 HepLorentzRotation & HepLorentzRotation::rotateY(double delta) {
163  double c1 = std::cos (delta);
164  double s1 = std::sin (delta);
165  HepLorentzVector rowx = row1();
166  HepLorentzVector rowz = row3();
167  HepLorentzVector r1 = c1 * rowx + s1 * rowz;
168  HepLorentzVector r3 = -s1 * rowx + c1 * rowz;
169  mxx = r1.x(); mxy = r1.y(); mxz = r1.z(); mxt = r1.t();
170  mzx = r3.x(); mzy = r3.y(); mzz = r3.z(); mzt = r3.t();
171  return *this;
172 }
173 
174 HepLorentzRotation & HepLorentzRotation::rotateZ(double delta) {
175  double c1 = std::cos (delta);
176  double s1 = std::sin (delta);
177  HepLorentzVector rowx = row1();
178  HepLorentzVector rowy = row2();
179  HepLorentzVector r1 = c1 * rowx - s1 * rowy;
180  HepLorentzVector r2 = s1 * rowx + c1 * rowy;
181  mxx = r1.x(); mxy = r1.y(); mxz = r1.z(); mxt = r1.t();
182  myx = r2.x(); myy = r2.y(); myz = r2.z(); myt = r2.t();
183  return *this;
184 }
185 
186 HepLorentzRotation & HepLorentzRotation::boostX(double beta) {
187  double b2 = beta*beta;
188 // if (b2 >= 1) {
189 // std::cerr << "HepLorentzRotation::boostX() - "
190 // << "Beta supplied to HepLorentzRotation::boostX represents speed >= c."
191 // << std::endl;
192 // }
193  double g1 = 1.0/std::sqrt(1.0-b2);
194  double bg = beta*g1;
195  HepLorentzVector rowx = row1();
196  HepLorentzVector rowt = row4();
197  HepLorentzVector r1 = g1 * rowx + bg * rowt;
198  HepLorentzVector r4 = bg * rowx + g1 * rowt;
199  mxx = r1.x(); mxy = r1.y(); mxz = r1.z(); mxt = r1.t();
200  mtx = r4.x(); mty = r4.y(); mtz = r4.z(); mtt = r4.t();
201  return *this;
202 }
203 
204 HepLorentzRotation & HepLorentzRotation::boostY(double beta) {
205  double b2 = beta*beta;
206 // if (b2 >= 1) {
207 // std::cerr << "HepLorentzRotation::boostY() - "
208 // << "Beta supplied to HepLorentzRotation::boostY represents speed >= c."
209 // << std::endl;
210 // }
211  double g1 = 1.0/std::sqrt(1.0-b2);
212  double bg = beta*g1;
213  HepLorentzVector rowy = row2();
214  HepLorentzVector rowt = row4();
215  HepLorentzVector r2 = g1 * rowy + bg * rowt;
216  HepLorentzVector r4 = bg * rowy + g1 * rowt;
217  myx = r2.x(); myy = r2.y(); myz = r2.z(); myt = r2.t();
218  mtx = r4.x(); mty = r4.y(); mtz = r4.z(); mtt = r4.t();
219  return *this;
220 }
221 
222 HepLorentzRotation & HepLorentzRotation::boostZ(double beta) {
223  double b2 = beta*beta;
224 // if (b2 >= 1) {
225 // std::cerr << "HepLorentzRotation::boostZ() - "
226 // << "Beta supplied to HepLorentzRotation::boostZ represents speed >= c."
227 // << std::endl;
228 // }
229  double g1 = 1.0/std::sqrt(1.0-b2);
230  double bg = beta*g1;
231  HepLorentzVector rowz = row3();
232  HepLorentzVector rowt = row4();
233  HepLorentzVector r3 = g1 * rowz + bg * rowt;
234  HepLorentzVector r4 = bg * rowz + g1 * rowt;
235  mtx = r4.x(); mty = r4.y(); mtz = r4.z(); mtt = r4.t();
236  mzx = r3.x(); mzy = r3.y(); mzz = r3.z(); mzt = r3.t();
237  return *this;
238 }
239 
240 std::ostream & HepLorentzRotation::print( std::ostream & os ) const {
241  os << "\n [ ( " <<
242  std::setw(11) << std::setprecision(6) << xx() << " " <<
243  std::setw(11) << std::setprecision(6) << xy() << " " <<
244  std::setw(11) << std::setprecision(6) << xz() << " " <<
245  std::setw(11) << std::setprecision(6) << xt() << ")\n"
246  << " ( " <<
247  std::setw(11) << std::setprecision(6) << yx() << " " <<
248  std::setw(11) << std::setprecision(6) << yy() << " " <<
249  std::setw(11) << std::setprecision(6) << yz() << " " <<
250  std::setw(11) << std::setprecision(6) << yt() << ")\n"
251  << " ( " <<
252  std::setw(11) << std::setprecision(6) << zx() << " " <<
253  std::setw(11) << std::setprecision(6) << zy() << " " <<
254  std::setw(11) << std::setprecision(6) << zz() << " " <<
255  std::setw(11) << std::setprecision(6) << zt() << ")\n"
256  << " ( " <<
257  std::setw(11) << std::setprecision(6) << tx() << " " <<
258  std::setw(11) << std::setprecision(6) << ty() << " " <<
259  std::setw(11) << std::setprecision(6) << tz() << " " <<
260  std::setw(11) << std::setprecision(6) << tt() << ") ]\n";
261  return os;
262 }
263 
264 HepLorentzRotation operator* ( const HepRotation & r,
265  const HepLorentzRotation & lt) {
266  r.rep4x4();
267  lt.rep4x4();
268  return HepLorentzRotation( HepRep4x4(
269  r.xx()*lt.xx() + r.xy()*lt.yx() + r.xz()*lt.zx() + r.xt()*lt.tx(),
270  r.xx()*lt.xy() + r.xy()*lt.yy() + r.xz()*lt.zy() + r.xt()*lt.ty(),
271  r.xx()*lt.xz() + r.xy()*lt.yz() + r.xz()*lt.zz() + r.xt()*lt.tz(),
272  r.xx()*lt.xt() + r.xy()*lt.yt() + r.xz()*lt.zt() + r.xt()*lt.tt(),
273 
274  r.yx()*lt.xx() + r.yy()*lt.yx() + r.yz()*lt.zx() + r.yt()*lt.tx(),
275  r.yx()*lt.xy() + r.yy()*lt.yy() + r.yz()*lt.zy() + r.yt()*lt.ty(),
276  r.yx()*lt.xz() + r.yy()*lt.yz() + r.yz()*lt.zz() + r.yt()*lt.tz(),
277  r.yx()*lt.xt() + r.yy()*lt.yt() + r.yz()*lt.zt() + r.yt()*lt.tt(),
278 
279  r.zx()*lt.xx() + r.zy()*lt.yx() + r.zz()*lt.zx() + r.zt()*lt.tx(),
280  r.zx()*lt.xy() + r.zy()*lt.yy() + r.zz()*lt.zy() + r.zt()*lt.ty(),
281  r.zx()*lt.xz() + r.zy()*lt.yz() + r.zz()*lt.zz() + r.zt()*lt.tz(),
282  r.zx()*lt.xt() + r.zy()*lt.yt() + r.zz()*lt.zt() + r.zt()*lt.tt(),
283 
284  r.tx()*lt.xx() + r.ty()*lt.yx() + r.tz()*lt.zx() + r.tt()*lt.tx(),
285  r.tx()*lt.xy() + r.ty()*lt.yy() + r.tz()*lt.zy() + r.tt()*lt.ty(),
286  r.tx()*lt.xz() + r.ty()*lt.yz() + r.tz()*lt.zz() + r.tt()*lt.tz(),
287  r.tx()*lt.xt() + r.ty()*lt.yt() + r.tz()*lt.zt() + r.tt()*lt.tt() ) );
288 }
289 
290 
291 const HepLorentzRotation HepLorentzRotation::IDENTITY;
292 
293 } // namespace CLHEP
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
static const G4double b2
static const G4double c1
void print(const std::vector< T > &data)
Definition: DicomRun.hh:109