source: trunk/CLHEP/src/LorentzRotation.cc@ 19

Last change on this file since 19 was 4, checked in by Pavel Demin, 16 years ago

first commit
File size: 10.3 KB
Line 
1// -*- C++ -*-
2// $Id: LorentzRotation.cc,v 1.1 2008-06-04 14:15:05 demin Exp $
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/defs.h"
17#include "CLHEP/Vector/LorentzRotation.h"
18#include "CLHEP/Vector/ZMxpv.h"
19
20#include <iostream>
21#include <iomanip>
22
23namespace CLHEP {
24
25// ---------- Constructors and Assignment:
26
27
28HepLorentzRotation & HepLorentzRotation::set
29 (double bx, double by, double bz) {
30 double bp2 = bx*bx + by*by + bz*bz;
31 if (bp2 >= 1) {
32 ZMthrowA (ZMxpvTachyonic(
33 "Boost Vector supplied to set HepLorentzRotation represents speed >= c."));
34 }
35 double gamma = 1.0 / sqrt(1.0 - bp2);
36 double bgamma = gamma * gamma / (1.0 + gamma);
37 mxx = 1.0 + bgamma * bx * bx;
38 myy = 1.0 + bgamma * by * by;
39 mzz = 1.0 + bgamma * bz * bz;
40 mxy = myx = bgamma * bx * by;
41 mxz = mzx = bgamma * bx * bz;
42 myz = mzy = bgamma * by * bz;
43 mxt = mtx = gamma * bx;
44 myt = mty = gamma * by;
45 mzt = mtz = gamma * bz;
46 mtt = gamma;
47 return *this;
48}
49
50HepLorentzRotation & HepLorentzRotation::set
51 (const HepBoost & B, const HepRotation & R) {
52 set (B.rep4x4());
53 *this = matrixMultiplication ( R.rep4x4() );
54 return *this;
55}
56
57HepLorentzRotation & HepLorentzRotation::set
58 (const HepRotation & R, const HepBoost & B) {
59 set (R.rep4x4());
60 *this = matrixMultiplication ( B.rep4x4() );
61 return *this;
62}
63
64// ---------- Accessors:
65
66// ------------ Subscripting:
67
68double HepLorentzRotation::operator () (int i, int j) const {
69 if (i == 0) {
70 if (j == 0) { return xx(); }
71 if (j == 1) { return xy(); }
72 if (j == 2) { return xz(); }
73 if (j == 3) { return xt(); }
74 } else if (i == 1) {
75 if (j == 0) { return yx(); }
76 if (j == 1) { return yy(); }
77 if (j == 2) { return yz(); }
78 if (j == 3) { return yt(); }
79 } else if (i == 2) {
80 if (j == 0) { return zx(); }
81 if (j == 1) { return zy(); }
82 if (j == 2) { return zz(); }
83 if (j == 3) { return zt(); }
84 } else if (i == 3) {
85 if (j == 0) { return tx(); }
86 if (j == 1) { return ty(); }
87 if (j == 2) { return tz(); }
88 if (j == 3) { return tt(); }
89 }
90 std::cerr << "HepLorentzRotation subscripting: bad indeces "
91 << "(" << i << "," << j << ")\n";
92 return 0.0;
93}
94
95// ---------- Application:
96
97
98// ---------- Comparison:
99
100int HepLorentzRotation::compare( const HepLorentzRotation & m ) const {
101 if (mtt<m.mtt) return -1; else if (mtt>m.mtt) return 1;
102 else if (mtz<m.mtz) return -1; else if (mtz>m.mtz) return 1;
103 else if (mty<m.mty) return -1; else if (mty>m.mty) return 1;
104 else if (mtx<m.mtx) return -1; else if (mtx>m.mtx) return 1;
105
106 else if (mzt<m.mzt) return -1; else if (mzt>m.mzt) return 1;
107 else if (mzz<m.mzz) return -1; else if (mzz>m.mzz) return 1;
108 else if (mzy<m.mzy) return -1; else if (mzy>m.mzy) return 1;
109 else if (mzx<m.mzx) return -1; else if (mzx>m.mzx) return 1;
110
111 else if (myt<m.myt) return -1; else if (myt>m.myt) return 1;
112 else if (myz<m.myz) return -1; else if (myz>m.myz) return 1;
113 else if (myy<m.myy) return -1; else if (myy>m.myy) return 1;
114 else if (myx<m.myx) return -1; else if (myx>m.myx) return 1;
115
116 else if (mxt<m.mxt) return -1; else if (mxt>m.mxt) return 1;
117 else if (mxz<m.mxz) return -1; else if (mxz>m.mxz) return 1;
118 else if (mxy<m.mxy) return -1; else if (mxy>m.mxy) return 1;
119 else if (mxx<m.mxx) return -1; else if (mxx>m.mxx) return 1;
120
121 else return 0;
122}
123
124
125// ---------- Operations in the group of 4-Rotations
126
127HepLorentzRotation
128HepLorentzRotation::matrixMultiplication(const HepRep4x4 & m) const {
129 return HepLorentzRotation(
130 mxx*m.xx_ + mxy*m.yx_ + mxz*m.zx_ + mxt*m.tx_,
131 mxx*m.xy_ + mxy*m.yy_ + mxz*m.zy_ + mxt*m.ty_,
132 mxx*m.xz_ + mxy*m.yz_ + mxz*m.zz_ + mxt*m.tz_,
133 mxx*m.xt_ + mxy*m.yt_ + mxz*m.zt_ + mxt*m.tt_,
134
135 myx*m.xx_ + myy*m.yx_ + myz*m.zx_ + myt*m.tx_,
136 myx*m.xy_ + myy*m.yy_ + myz*m.zy_ + myt*m.ty_,
137 myx*m.xz_ + myy*m.yz_ + myz*m.zz_ + myt*m.tz_,
138 myx*m.xt_ + myy*m.yt_ + myz*m.zt_ + myt*m.tt_,
139
140 mzx*m.xx_ + mzy*m.yx_ + mzz*m.zx_ + mzt*m.tx_,
141 mzx*m.xy_ + mzy*m.yy_ + mzz*m.zy_ + mzt*m.ty_,
142 mzx*m.xz_ + mzy*m.yz_ + mzz*m.zz_ + mzt*m.tz_,
143 mzx*m.xt_ + mzy*m.yt_ + mzz*m.zt_ + mzt*m.tt_,
144
145 mtx*m.xx_ + mty*m.yx_ + mtz*m.zx_ + mtt*m.tx_,
146 mtx*m.xy_ + mty*m.yy_ + mtz*m.zy_ + mtt*m.ty_,
147 mtx*m.xz_ + mty*m.yz_ + mtz*m.zz_ + mtt*m.tz_,
148 mtx*m.xt_ + mty*m.yt_ + mtz*m.zt_ + mtt*m.tt_ );
149}
150
151HepLorentzRotation & HepLorentzRotation::rotateX(double delta) {
152 double c = cos (delta);
153 double s = sin (delta);
154 HepLorentzVector rowy = row2();
155 HepLorentzVector rowz = row3();
156 HepLorentzVector r2 = c * rowy - s * rowz;
157 HepLorentzVector r3 = s * rowy + c * rowz;
158 myx = r2.x(); myy = r2.y(); myz = r2.z(); myt = r2.t();
159 mzx = r3.x(); mzy = r3.y(); mzz = r3.z(); mzt = r3.t();
160 return *this;
161}
162
163HepLorentzRotation & HepLorentzRotation::rotateY(double delta) {
164 double c = cos (delta);
165 double s = sin (delta);
166 HepLorentzVector rowx = row1();
167 HepLorentzVector rowz = row3();
168 HepLorentzVector r1 = c * rowx + s * rowz;
169 HepLorentzVector r3 = -s * rowx + c * rowz;
170 mxx = r1.x(); mxy = r1.y(); mxz = r1.z(); mxt = r1.t();
171 mzx = r3.x(); mzy = r3.y(); mzz = r3.z(); mzt = r3.t();
172 return *this;
173}
174
175HepLorentzRotation & HepLorentzRotation::rotateZ(double delta) {
176 double c = cos (delta);
177 double s = sin (delta);
178 HepLorentzVector rowx = row1();
179 HepLorentzVector rowy = row2();
180 HepLorentzVector r1 = c * rowx - s * rowy;
181 HepLorentzVector r2 = s * rowx + c * rowy;
182 mxx = r1.x(); mxy = r1.y(); mxz = r1.z(); mxt = r1.t();
183 myx = r2.x(); myy = r2.y(); myz = r2.z(); myt = r2.t();
184 return *this;
185}
186
187HepLorentzRotation & HepLorentzRotation::boostX(double beta) {
188 double b2 = beta*beta;
189 if (b2 >= 1) {
190 ZMthrowA (ZMxpvTachyonic(
191 "Beta supplied to HepLorentzRotation::boostX represents speed >= c."));
192 }
193 double g = 1.0/sqrt(1.0-b2);
194 double bg = beta*g;
195 HepLorentzVector rowx = row1();
196 HepLorentzVector rowt = row4();
197 HepLorentzVector r1 = g * rowx + bg * rowt;
198 HepLorentzVector r4 = bg * rowx + g * 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
204HepLorentzRotation & HepLorentzRotation::boostY(double beta) {
205 double b2 = beta*beta;
206 if (b2 >= 1) {
207 ZMthrowA (ZMxpvTachyonic(
208 "Beta supplied to HepLorentzRotation::boostY represents speed >= c."));
209 }
210 double g = 1.0/sqrt(1.0-b2);
211 double bg = beta*g;
212 HepLorentzVector rowy = row2();
213 HepLorentzVector rowt = row4();
214 HepLorentzVector r2 = g * rowy + bg * rowt;
215 HepLorentzVector r4 = bg * rowy + g * rowt;
216 myx = r2.x(); myy = r2.y(); myz = r2.z(); myt = r2.t();
217 mtx = r4.x(); mty = r4.y(); mtz = r4.z(); mtt = r4.t();
218 return *this;
219}
220
221HepLorentzRotation & HepLorentzRotation::boostZ(double beta) {
222 double b2 = beta*beta;
223 if (b2 >= 1) {
224 ZMthrowA (ZMxpvTachyonic(
225 "Beta supplied to HepLorentzRotation::boostZ represents speed >= c."));
226 }
227 double g = 1.0/sqrt(1.0-b2);
228 double bg = beta*g;
229 HepLorentzVector rowz = row3();
230 HepLorentzVector rowt = row4();
231 HepLorentzVector r3 = g * rowz + bg * rowt;
232 HepLorentzVector r4 = bg * rowz + g * rowt;
233 mtx = r4.x(); mty = r4.y(); mtz = r4.z(); mtt = r4.t();
234 mzx = r3.x(); mzy = r3.y(); mzz = r3.z(); mzt = r3.t();
235 return *this;
236}
237
238std::ostream & HepLorentzRotation::print( std::ostream & os ) const {
239// using std::setw;
240// using std::setprecision;
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
264HepLorentzRotation 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
291const HepLorentzRotation HepLorentzRotation::IDENTITY;
292
293} // namespace CLHEP
Note: See TracBrowser for help on using the repository browser.