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 the HepBoost class.
7 //
8
9 #ifdef GNUPRAGMA
10 #pragma implementation
11 #endif
12
13 #include "CLHEP/Vector/defs.h"
14 #include "CLHEP/Vector/Boost.h"
15 #include "CLHEP/Vector/Rotation.h"
16 #include "CLHEP/Vector/LorentzRotation.h"
17 #include "CLHEP/Vector/ZMxpv.h"
18
19 namespace CLHEP {
20
21 // ---------- Constructors and Assignment:
22
set(double bx,double by,double bz)23 HepBoost & HepBoost::set (double bx, double by, double bz) {
24 double bp2 = bx*bx + by*by + bz*bz;
25 if (bp2 >= 1) {
26 ZMthrowA (ZMxpvTachyonic(
27 "Boost Vector supplied to set HepBoost represents speed >= c."));
28 }
29 double ggamma = 1.0 / std::sqrt(1.0 - bp2);
30 double bgamma = ggamma * ggamma / (1.0 + ggamma);
31 rep_.xx_ = 1.0 + bgamma * bx * bx;
32 rep_.yy_ = 1.0 + bgamma * by * by;
33 rep_.zz_ = 1.0 + bgamma * bz * bz;
34 rep_.xy_ = bgamma * bx * by;
35 rep_.xz_ = bgamma * bx * bz;
36 rep_.yz_ = bgamma * by * bz;
37 rep_.xt_ = ggamma * bx;
38 rep_.yt_ = ggamma * by;
39 rep_.zt_ = ggamma * bz;
40 rep_.tt_ = ggamma;
41 return *this;
42 }
43
set(const HepRep4x4Symmetric & m1)44 HepBoost & HepBoost::set (const HepRep4x4Symmetric & m1) {
45 rep_ = m1;
46 return *this;
47 }
48
set(Hep3Vector ddirection,double bbeta)49 HepBoost & HepBoost::set (Hep3Vector ddirection, double bbeta) {
50 double length = ddirection.mag();
51 if (length <= 0) { // Nan-proofing
52 ZMthrowA (ZMxpvZeroVector(
53 "Direction supplied to set HepBoost is zero."));
54 set (0,0,0);
55 return *this;
56 }
57 set(bbeta*ddirection.x()/length,
58 bbeta*ddirection.y()/length,
59 bbeta*ddirection.z()/length);
60 return *this;
61 }
62
set(const Hep3Vector & bboost)63 HepBoost & HepBoost::set (const Hep3Vector & bboost) {
64 return set (bboost.x(), bboost.y(), bboost.z());
65 }
66
67 // ---------- Accessors:
68
69 // ---------- Decomposition:
70
decompose(HepRotation & rotation,HepBoost & boost) const71 void HepBoost::decompose (HepRotation & rotation, HepBoost & boost) const {
72 HepAxisAngle vdelta = HepAxisAngle();
73 rotation = HepRotation(vdelta);
74 Hep3Vector bbeta = boostVector();
75 boost = HepBoost(bbeta);
76 }
77
decompose(HepAxisAngle & rotation,Hep3Vector & boost) const78 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const {
79 rotation = HepAxisAngle();
80 boost = boostVector();
81 }
82
decompose(HepBoost & boost,HepRotation & rotation) const83 void HepBoost::decompose (HepBoost & boost, HepRotation & rotation) const {
84 HepAxisAngle vdelta = HepAxisAngle();
85 rotation = HepRotation(vdelta);
86 Hep3Vector bbeta = boostVector();
87 boost = HepBoost(bbeta);
88 }
89
decompose(Hep3Vector & boost,HepAxisAngle & rotation) const90 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const {
91 rotation = HepAxisAngle();
92 boost = boostVector();
93 }
94
95 // ---------- Comparisons:
96
distance2(const HepRotation & r) const97 double HepBoost::distance2( const HepRotation & r ) const {
98 double db2 = norm2();
99 double dr2 = r.norm2();
100 return (db2 + dr2);
101 }
102
distance2(const HepLorentzRotation & lt) const103 double HepBoost::distance2( const HepLorentzRotation & lt ) const {
104 HepBoost b1;
105 HepRotation r1;
106 lt.decompose(b1,r1);
107 double db2 = distance2(b1);
108 double dr2 = r1.norm2();
109 return (db2 + dr2);
110 }
111
howNear(const HepRotation & r) const112 double HepBoost::howNear ( const HepRotation & r ) const {
113 return std::sqrt(distance2(r));
114 }
115
howNear(const HepLorentzRotation & lt) const116 double HepBoost::howNear ( const HepLorentzRotation & lt ) const {
117 return std::sqrt(distance2(lt));
118 }
119
isNear(const HepRotation & r,double epsilon) const120 bool HepBoost::isNear (const HepRotation & r, double epsilon) const {
121 double db2 = norm2();
122 if (db2 > epsilon*epsilon) return false;
123 double dr2 = r.norm2();
124 return (db2+dr2 <= epsilon*epsilon);
125 }
126
isNear(const HepLorentzRotation & lt,double epsilon) const127 bool HepBoost::isNear (const HepLorentzRotation & lt,
128 double epsilon) const {
129 HepBoost b1;
130 HepRotation r1;
131 double db2 = distance2(b1);
132 lt.decompose(b1,r1);
133 if (db2 > epsilon*epsilon) return false;
134 double dr2 = r1.norm2();
135 return (db2 + dr2);
136 }
137
138 // ---------- Properties:
139
norm2() const140 double HepBoost::norm2() const {
141 double bgx = rep_.xt_;
142 double bgy = rep_.yt_;
143 double bgz = rep_.zt_;
144 return bgx*bgx+bgy*bgy+bgz*bgz;
145 }
146
rectify()147 void HepBoost::rectify() {
148 // Assuming the representation of this is close to a true pure boost,
149 // but may have drifted due to round-off error from many operations,
150 // this forms an "exact" pure boost matrix for the LT again.
151
152 // The natural way to do this is to use the t column as a boost and set
153 // based on that boost vector.
154
155 // There is perhaps danger that this boost vector will appear equal to or
156 // greater than a unit vector; the best we can do for such a case is use
157 // a boost in that direction but rescaled to just less than one.
158
159 // There is in principle no way that gamma could have become negative,
160 // but if that happens, we ZMthrow and (if continuing) just rescale, which
161 // will change the sign of the last column when computing the boost.
162
163 double gam = tt();
164 if (gam <= 0) { // 4/12/01 mf
165 // ZMthrowA (ZMxpvTachyonic(
166 ZMthrowC (ZMxpvTachyonic(
167 "Attempt to rectify a boost with non-positive gamma."));
168 if (gam==0) return; // NaN-proofing
169 }
170 Hep3Vector boost (xt(), yt(), zt());
171 boost /= tt();
172 if ( boost.mag2() >= 1 ) { // NaN-proofing:
173 boost /= ( boost.mag() * ( 1.0 + 1.0e-16 ) ); // used to just check > 1
174 }
175 set ( boost );
176 }
177
178 // ---------- Application is all in .icc
179
180 // ---------- Operations in the group of 4-Rotations
181
182 HepLorentzRotation
matrixMultiplication(const HepRep4x4 & m1) const183 HepBoost::matrixMultiplication(const HepRep4x4 & m1) const {
184 HepRep4x4Symmetric r = rep4x4Symmetric();
185 return HepLorentzRotation( HepRep4x4 (
186 r.xx_*m1.xx_ + r.xy_*m1.yx_ + r.xz_*m1.zx_ + r.xt_*m1.tx_,
187 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.zy_ + r.xt_*m1.ty_,
188 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.tz_,
189 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
190
191 r.xy_*m1.xx_ + r.yy_*m1.yx_ + r.yz_*m1.zx_ + r.yt_*m1.tx_,
192 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.zy_ + r.yt_*m1.ty_,
193 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.tz_,
194 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
195
196 r.xz_*m1.xx_ + r.yz_*m1.yx_ + r.zz_*m1.zx_ + r.zt_*m1.tx_,
197 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.zy_ + r.zt_*m1.ty_,
198 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.tz_,
199 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
200
201 r.xt_*m1.xx_ + r.yt_*m1.yx_ + r.zt_*m1.zx_ + r.tt_*m1.tx_,
202 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.zy_ + r.tt_*m1.ty_,
203 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.tz_,
204 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
205 }
206
207 HepLorentzRotation
matrixMultiplication(const HepRep4x4Symmetric & m1) const208 HepBoost::matrixMultiplication(const HepRep4x4Symmetric & m1) const {
209 HepRep4x4Symmetric r = rep4x4Symmetric();
210 return HepLorentzRotation( HepRep4x4 (
211 r.xx_*m1.xx_ + r.xy_*m1.xy_ + r.xz_*m1.xz_ + r.xt_*m1.xt_,
212 r.xx_*m1.xy_ + r.xy_*m1.yy_ + r.xz_*m1.yz_ + r.xt_*m1.yt_,
213 r.xx_*m1.xz_ + r.xy_*m1.yz_ + r.xz_*m1.zz_ + r.xt_*m1.zt_,
214 r.xx_*m1.xt_ + r.xy_*m1.yt_ + r.xz_*m1.zt_ + r.xt_*m1.tt_,
215
216 r.xy_*m1.xx_ + r.yy_*m1.xy_ + r.yz_*m1.xz_ + r.yt_*m1.xt_,
217 r.xy_*m1.xy_ + r.yy_*m1.yy_ + r.yz_*m1.yz_ + r.yt_*m1.yt_,
218 r.xy_*m1.xz_ + r.yy_*m1.yz_ + r.yz_*m1.zz_ + r.yt_*m1.zt_,
219 r.xy_*m1.xt_ + r.yy_*m1.yt_ + r.yz_*m1.zt_ + r.yt_*m1.tt_,
220
221 r.xz_*m1.xx_ + r.yz_*m1.xy_ + r.zz_*m1.xz_ + r.zt_*m1.xt_,
222 r.xz_*m1.xy_ + r.yz_*m1.yy_ + r.zz_*m1.yz_ + r.zt_*m1.yt_,
223 r.xz_*m1.xz_ + r.yz_*m1.yz_ + r.zz_*m1.zz_ + r.zt_*m1.zt_,
224 r.xz_*m1.xt_ + r.yz_*m1.yt_ + r.zz_*m1.zt_ + r.zt_*m1.tt_,
225
226 r.xt_*m1.xx_ + r.yt_*m1.xy_ + r.zt_*m1.xz_ + r.tt_*m1.xt_,
227 r.xt_*m1.xy_ + r.yt_*m1.yy_ + r.zt_*m1.yz_ + r.tt_*m1.yt_,
228 r.xt_*m1.xz_ + r.yt_*m1.yz_ + r.zt_*m1.zz_ + r.tt_*m1.zt_,
229 r.xt_*m1.xt_ + r.yt_*m1.yt_ + r.zt_*m1.zt_ + r.tt_*m1.tt_) );
230 }
231
232 HepLorentzRotation
operator *(const HepLorentzRotation & lt) const233 HepBoost::operator* (const HepLorentzRotation & lt) const {
234 return matrixMultiplication(lt.rep4x4());
235 }
236
237 HepLorentzRotation
operator *(const HepBoost & b) const238 HepBoost::operator* (const HepBoost & b) const {
239 return matrixMultiplication(b.rep_);
240 }
241
242 HepLorentzRotation
operator *(const HepRotation & r) const243 HepBoost::operator* (const HepRotation & r) const {
244 return matrixMultiplication(r.rep4x4());
245 }
246
247 // ---------- I/O:
248
print(std::ostream & os) const249 std::ostream & HepBoost::print( std::ostream & os ) const {
250 if ( rep_.tt_ <= 1 ) {
251 os << "Lorentz Boost( IDENTITY )";
252 } else {
253 double norm = boostVector().mag();
254 os << "\nLorentz Boost " << boostVector()/norm <<
255 "\n{beta = " << beta() << " gamma = " << gamma() << "}\n";
256 }
257 return os;
258 }
259
260 } // namespace CLHEP
261