CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

Boost.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 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 
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 
45  rep_ = m1;
46  return *this;
47 }
48 
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 
63 HepBoost & HepBoost::set (const Hep3Vector & bboost) {
64  return set (bboost.x(), bboost.y(), bboost.z());
65 }
66 
67 // ---------- Accessors:
68 
69 // ---------- Decomposition:
70 
71 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 
78 void HepBoost::decompose (HepAxisAngle & rotation, Hep3Vector & boost) const {
79  rotation = HepAxisAngle();
80  boost = boostVector();
81 }
82 
83 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 
90 void HepBoost::decompose (Hep3Vector & boost, HepAxisAngle & rotation) const {
91  rotation = HepAxisAngle();
92  boost = boostVector();
93 }
94 
95 // ---------- Comparisons:
96 
97 double HepBoost::distance2( const HepRotation & r ) const {
98  double db2 = norm2();
99  double dr2 = r.norm2();
100  return (db2 + dr2);
101 }
102 
103 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 
112 double HepBoost::howNear ( const HepRotation & r ) const {
113  return std::sqrt(distance2(r));
114 }
115 
116 double HepBoost::howNear ( const HepLorentzRotation & lt ) const {
117  return std::sqrt(distance2(lt));
118 }
119 
120 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 
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 
140 double HepBoost::norm2() const {
141  register double bgx = rep_.xt_;
142  register double bgy = rep_.yt_;
143  register double bgz = rep_.zt_;
144  return bgx*bgx+bgy*bgy+bgz*bgz;
145 }
146 
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 
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 
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 
234  return matrixMultiplication(lt.rep4x4());
235 }
236 
239  return matrixMultiplication(b.rep_);
240 }
241 
244  return matrixMultiplication(r.rep4x4());
245 }
246 
247 // ---------- I/O:
248 
249 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