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 ]

Diff markup

Differences between /externals/clhep/src/LorentzRotationD.cc (Version 11.3.0) and /externals/clhep/src/LorentzRotationD.cc (Version 10.7.p4)


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