Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/LorentzRotationD.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 those parts of the HepLorentzRotation class
  7 // which involve decomposition into Boost*Rotation.
  8 
  9 #include "CLHEP/Vector/LorentzRotation.h"
 10 
 11 #include <cmath>
 12 #include <iostream>
 13 
 14 namespace CLHEP  {
 15 
 16 // ----------  Decomposition:
 17 
 18 void HepLorentzRotation::decompose 
 19   (HepBoost & bboost, HepRotation & rotation) const {
 20 
 21   // The boost will be the pure boost based on column 4 of the transformation
 22   // matrix.  Since the constructor takes the beta vector, and not beta*gamma,
 23   // we first divide through by gamma = the tt element.  This of course can
 24   // never be zero since the last row has t**2 - v**2 = +1.
 25 
 26   Hep3Vector betaVec ( xt(), yt(), zt() );
 27   betaVec *= 1.0 / tt();
 28   bboost.set( betaVec );
 29 
 30   // The rotation will be inverse of B times T.
 31 
 32   HepBoost B( -betaVec );
 33   HepLorentzRotation R( B * *this );
 34 
 35   HepRep3x3 m1  ( R.xx(), R.xy(), R.xz(),
 36                   R.yx(), R.yy(), R.yz(),
 37                   R.zx(), R.zy(), R.zz() );
 38   rotation.set( m1 );
 39   rotation.rectify();
 40   
 41   return;
 42 
 43 }
 44 
 45 void HepLorentzRotation::decompose 
 46   (Hep3Vector & bboost, HepAxisAngle & rotation) const {
 47   HepRotation r;
 48   HepBoost b;
 49   decompose(b,r);
 50   bboost = b.boostVector();
 51   rotation = r.axisAngle();
 52   return;
 53 }
 54 
 55 void HepLorentzRotation::decompose 
 56   (HepRotation & rotation, HepBoost & bboost) const {
 57 
 58   // In this case the pure boost is based on row 4 of the matrix.  
 59 
 60   Hep3Vector betaVec( tx(), ty(), tz() );
 61   betaVec *= 1.0 / tt();
 62   bboost.set( betaVec );
 63 
 64   // The rotation will be T times the inverse of B.
 65 
 66   HepBoost B( -betaVec );
 67   HepLorentzRotation R( *this * B );
 68 
 69   HepRep3x3 m1 ( R.xx(), R.xy(), R.xz(),
 70                  R.yx(), R.yy(), R.yz(),
 71                  R.zx(), R.zy(), R.zz() );
 72   rotation.set( m1 );
 73   rotation.rectify();
 74   return;
 75 
 76 }
 77 
 78 void HepLorentzRotation::decompose 
 79   (HepAxisAngle & rotation, Hep3Vector & bboost) const {
 80   HepRotation r;
 81   HepBoost b;
 82   decompose(r,b);
 83   rotation = r.axisAngle();
 84   bboost = b.boostVector();
 85   return;
 86 }
 87 
 88 double HepLorentzRotation::distance2( const HepBoost & b ) const {
 89   HepBoost    b1;
 90   HepRotation r1; 
 91   decompose( b1, r1 );
 92   double db2 = b1.distance2( b );
 93   double dr2 = r1.norm2(); 
 94   return ( db2 + dr2 );
 95 }
 96 
 97 double HepLorentzRotation::distance2( const HepRotation & r ) const {
 98   HepBoost    b1;
 99   HepRotation r1; 
100   decompose( b1, r1 );
101   double db2 = b1.norm2( );
102   double dr2 = r1.distance2( r ); 
103   return ( db2 + dr2 );
104 }
105 
106 double HepLorentzRotation::distance2( 
107            const HepLorentzRotation & lt  ) const {
108   HepBoost    b1;
109   HepRotation r1; 
110   decompose( b1, r1 );
111   HepBoost    b2;
112   HepRotation r2; 
113   lt.decompose (b2, r2);
114   double db2 = b1.distance2( b2 );
115   double dr2 = r1.distance2( r2 ); 
116   return ( db2 + dr2 );
117 }
118 
119 double HepLorentzRotation::howNear( const HepBoost & b ) const {
120   return std::sqrt( distance2( b ) );
121 }
122 double HepLorentzRotation::howNear( const HepRotation & r ) const {
123   return std::sqrt( distance2( r ) );
124 }
125 double HepLorentzRotation::howNear( const HepLorentzRotation & lt )const {
126   return std::sqrt( distance2( lt ) );
127 }
128 
129 bool HepLorentzRotation::isNear(
130     const HepBoost & b, double epsilon ) const {
131   HepBoost    b1;
132   HepRotation r1; 
133   decompose( b1, r1 );
134   double db2 = b1.distance2(b);
135   if ( db2 > epsilon*epsilon ) {
136      return false;       // Saves the time-consuming Rotation::norm2
137   }
138   double dr2 = r1.norm2();
139   return ( (db2 + dr2) <= epsilon*epsilon );
140 }
141 
142 bool HepLorentzRotation::isNear(
143     const HepRotation & r, double epsilon ) const {
144   HepBoost    b1;
145   HepRotation r1; 
146   decompose( b1, r1 );
147   double db2 = b1.norm2();
148   if ( db2 > epsilon*epsilon ) {
149      return false;       // Saves the time-consuming Rotation::distance2
150   }
151   double dr2 = r1.distance2(r);
152   return ( (db2 + dr2) <= epsilon*epsilon );
153 }
154 
155 bool HepLorentzRotation::isNear(
156     const HepLorentzRotation & lt, double epsilon ) const {
157   HepBoost    b1;
158   HepRotation r1; 
159   decompose( b1, r1 );
160   HepBoost    b2;
161   HepRotation r2; 
162   lt.decompose (b2, r2);
163   double db2 = b1.distance2(b2);
164   if ( db2 > epsilon*epsilon ) {
165      return false;       // Saves the time-consuming Rotation::distance2
166   }
167   double dr2 = r1.distance2(r2);
168   return ( (db2 + dr2) <= epsilon*epsilon );
169 }
170 
171 double HepLorentzRotation::norm2() const {
172   HepBoost    b;
173   HepRotation r;
174   decompose( b, r );
175   return b.norm2() + r.norm2();
176 }
177 
178 void HepLorentzRotation::rectify() {
179   
180   // Assuming the representation of this is close to a true LT,
181   // but may have drifted due to round-off error from many operations,
182   // this forms an "exact" orthosymplectic matrix for the LT again.
183  
184   // There are several ways to do this, all equivalent to lowest order in
185   // the corrected error.  We choose to form an LT based on the inverse boost
186   // extracted from row 4, and left-multiply by LT to form what would be
187   // a rotation if the LT were kosher.  We drop the possibly non-zero t
188   // components of that, rectify that rotation and multiply back by the boost.
189                     
190   Hep3Vector beta (tx(), ty(), tz());
191   double gam = tt();      // NaN-proofing
192   if ( gam <= 0 ) {
193     std::cerr << "HepLorentzRotation::rectify() - "
194   << "rectify() on a transformation with tt() <= 0 - will not help!"
195         << std::endl;
196     gam = 1;
197   }
198   beta *= 1.0/gam;
199   HepLorentzRotation R = (*this) * HepBoost(-beta);
200 
201   HepRep3x3  m1 ( R.xx(), R.xy(), R.xz(),
202                   R.yx(), R.yy(), R.yz(),
203                   R.zx(), R.zy(), R.zz() );
204 
205   HepRotation Rgood (m1);
206   Rgood.rectify();
207 
208   set ( Rgood, HepBoost(beta) );
209 }
210 
211 }  // namespace CLHEP
212