Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/Boost.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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