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