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