Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/LorentzRotationC.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/LorentzRotationC.cc (Version 11.3.0) and /externals/clhep/src/LorentzRotationC.cc (Version 11.0.p2)


  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 that part of       6 // This is the implementation of that part of the HepLorentzRotation class
  7 // which is concerned with setting or construc      7 // which is concerned with setting or constructing the transformation based 
  8 // on 4 supplied columns or rows.                   8 // on 4 supplied columns or rows.
  9                                                     9 
 10 #include "CLHEP/Vector/LorentzRotation.h"          10 #include "CLHEP/Vector/LorentzRotation.h"
 11 #include "CLHEP/Vector/LorentzVector.h"            11 #include "CLHEP/Vector/LorentzVector.h"
 12                                                    12 
 13 #include <cmath>                                   13 #include <cmath>
 14 #include <iostream>                                14 #include <iostream>
 15                                                    15 
 16 namespace CLHEP  {                                 16 namespace CLHEP  {
 17                                                    17 
 18 // ----------  Constructors and Assignment:        18 // ----------  Constructors and Assignment:
 19                                                    19 
 20 HepLorentzRotation & HepLorentzRotation::set (     20 HepLorentzRotation & HepLorentzRotation::set (const HepLorentzVector & ccol1,
 21                                   const HepLor     21                                   const HepLorentzVector & ccol2,
 22                             const HepLorentzVe     22                             const HepLorentzVector & ccol3,
 23                             const HepLorentzVe     23                             const HepLorentzVector & ccol4) {
 24   // First, test that the four cols do represe     24   // First, test that the four cols do represent something close to a
 25   // true LT:                                      25   // true LT:
 26                                                    26 
 27   ZMpvMetric_t savedMetric = HepLorentzVector:     27   ZMpvMetric_t savedMetric = HepLorentzVector::setMetric (TimePositive);
 28                                                    28 
 29   if ( ccol4.getT() < 0 ) {                        29   if ( ccol4.getT() < 0 ) {
 30     std::cerr << "HepLorentzRotation::set() -      30     std::cerr << "HepLorentzRotation::set() - "
 31       << "column 4 supplied to define transfor     31       << "column 4 supplied to define transformation has negative T component"
 32       << std::endl;                                32       << std::endl;
 33     *this = HepLorentzRotation();                  33     *this = HepLorentzRotation();
 34     return *this;                                  34     return *this;
 35   }                                                35   }
 36 /*                                                 36 /*
 37   double u1u1 = ccol1.dot(ccol1);                  37   double u1u1 = ccol1.dot(ccol1);
 38   double f11  = std::fabs(u1u1 + 1.0);             38   double f11  = std::fabs(u1u1 + 1.0);
 39   if ( f11 > Hep4RotationInterface::tolerance      39   if ( f11 > Hep4RotationInterface::tolerance ) {
 40     std::cerr << "HepLorentzRotation::set() -      40     std::cerr << "HepLorentzRotation::set() - "
 41       << "column 1 supplied for HepLorentzRota     41       << "column 1 supplied for HepLorentzRotation has w*w != -1" << std::endl;
 42   }                                                42   }
 43   double u2u2 = ccol2.dot(ccol2);                  43   double u2u2 = ccol2.dot(ccol2);
 44   double f22  = std::fabs(u2u2 + 1.0);             44   double f22  = std::fabs(u2u2 + 1.0);
 45   if ( f22 > Hep4RotationInterface::tolerance      45   if ( f22 > Hep4RotationInterface::tolerance ) {
 46     std::cerr << "HepLorentzRotation::set() -      46     std::cerr << "HepLorentzRotation::set() - "
 47       << "column 2 supplied for HepLorentzRota     47       << "column 2 supplied for HepLorentzRotation has w*w != -1" << std::endl;
 48   }                                                48   }
 49   double u3u3 = ccol3.dot(ccol3);                  49   double u3u3 = ccol3.dot(ccol3);
 50   double f33  = std::fabs(u3u3 + 1.0);             50   double f33  = std::fabs(u3u3 + 1.0);
 51   if ( f33 > Hep4RotationInterface::tolerance      51   if ( f33 > Hep4RotationInterface::tolerance ) {
 52     std::cerr << "HepLorentzRotation::set() -      52     std::cerr << "HepLorentzRotation::set() - "
 53       << "column 3 supplied for HepLorentzRota     53       << "column 3 supplied for HepLorentzRotation has w*w != -1" << std::endl;
 54   }                                                54   }
 55   double u4u4 = ccol4.dot(ccol4);                  55   double u4u4 = ccol4.dot(ccol4);
 56   double f44  = std::fabs(u4u4 - 1.0);             56   double f44  = std::fabs(u4u4 - 1.0);
 57   if ( f44 > Hep4RotationInterface::tolerance      57   if ( f44 > Hep4RotationInterface::tolerance ) {
 58     std::cerr << "HepLorentzRotation::set() -      58     std::cerr << "HepLorentzRotation::set() - "
 59       << "column 4 supplied for HepLorentzRota     59       << "column 4 supplied for HepLorentzRotation has w*w != +1" << std::endl;
 60   }                                                60   }
 61                                                    61 
 62   double u1u2 = ccol1.dot(ccol2);                  62   double u1u2 = ccol1.dot(ccol2);
 63   double f12  = std::fabs(u1u2);                   63   double f12  = std::fabs(u1u2);
 64   if ( f12 > Hep4RotationInterface::tolerance      64   if ( f12 > Hep4RotationInterface::tolerance ) {
 65     std::cerr << "HepLorentzRotation::set() -      65     std::cerr << "HepLorentzRotation::set() - "
 66       << "columns 1 and 2 supplied for HepLore     66       << "columns 1 and 2 supplied for HepLorentzRotation have non-zero dot" << std::endl;
 67   }                                                67   }
 68   double u1u3 = ccol1.dot(ccol3);                  68   double u1u3 = ccol1.dot(ccol3);
 69   double f13  = std::fabs(u1u3);                   69   double f13  = std::fabs(u1u3);
 70                                                    70 
 71   if ( f13 > Hep4RotationInterface::tolerance      71   if ( f13 > Hep4RotationInterface::tolerance ) {
 72     std::cerr << "HepLorentzRotation::set() -      72     std::cerr << "HepLorentzRotation::set() - "
 73       << "columns 1 and 3 supplied for HepLore     73       << "columns 1 and 3 supplied for HepLorentzRotation have non-zero dot" << std::endl;
 74   }                                                74   }
 75   double u1u4 = ccol1.dot(ccol4);                  75   double u1u4 = ccol1.dot(ccol4);
 76   double f14  = std::fabs(u1u4);                   76   double f14  = std::fabs(u1u4);
 77   if ( f14 > Hep4RotationInterface::tolerance      77   if ( f14 > Hep4RotationInterface::tolerance ) {
 78     std::cerr << "HepLorentzRotation::set() -      78     std::cerr << "HepLorentzRotation::set() - "
 79       << "columns 1 and 4 supplied for HepLore     79       << "columns 1 and 4 supplied for HepLorentzRotation have non-zero dot" << std::endl;
 80   }                                                80   }
 81   double u2u3 = ccol2.dot(ccol3);                  81   double u2u3 = ccol2.dot(ccol3);
 82   double f23  = std::fabs(u2u3);                   82   double f23  = std::fabs(u2u3);
 83   if ( f23 > Hep4RotationInterface::tolerance      83   if ( f23 > Hep4RotationInterface::tolerance ) {
 84     std::cerr << "HepLorentzRotation::set() -      84     std::cerr << "HepLorentzRotation::set() - "
 85       << "columns 2 and 3 supplied for HepLore     85       << "columns 2 and 3 supplied for HepLorentzRotation have non-zero dot" << std::endl;
 86   }                                                86   }
 87   double u2u4 = ccol2.dot(ccol4);                  87   double u2u4 = ccol2.dot(ccol4);
 88   double f24  = std::fabs(u2u4);                   88   double f24  = std::fabs(u2u4);
 89   if ( f24 > Hep4RotationInterface::tolerance      89   if ( f24 > Hep4RotationInterface::tolerance ) {
 90     std::cerr << "HepLorentzRotation::set() -      90     std::cerr << "HepLorentzRotation::set() - "
 91       << "columns 2 and 4 supplied for HepLore     91       << "columns 2 and 4 supplied for HepLorentzRotation have non-zero dot" << std::endl;
 92   }                                                92   }
 93   double u3u4 = ccol3.dot(ccol4);                  93   double u3u4 = ccol3.dot(ccol4);
 94   double f34  = std::fabs(u3u4);                   94   double f34  = std::fabs(u3u4);
 95   if ( f34 > Hep4RotationInterface::tolerance      95   if ( f34 > Hep4RotationInterface::tolerance ) {
 96     std::cerr << "HepLorentzRotation::set() -      96     std::cerr << "HepLorentzRotation::set() - "
 97       << "columns 3 and 4 supplied for HepLore     97       << "columns 3 and 4 supplied for HepLorentzRotation have non-zero dot" << std::endl;
 98   }                                                98   }
 99 */                                                 99 */
100   // Our strategy will be to order the cols, t    100   // Our strategy will be to order the cols, then do gram-schmidt on them
101   // (that is, remove the components of col d     101   // (that is, remove the components of col d that make it non-orthogonal to
102   // col c, normalize that, then remove the co    102   // col c, normalize that, then remove the components of b that make it
103   // non-orthogonal to d and to c, normalize t    103   // non-orthogonal to d and to c, normalize that, etc.
104                                                   104 
105   // Because col4, the time col, is most likel    105   // Because col4, the time col, is most likely to be computed directly, we
106   // will start from there and work left-ward.    106   // will start from there and work left-ward.
107                                                   107 
108   HepLorentzVector a, b, c, d;                    108   HepLorentzVector a, b, c, d;
109   bool isLorentzTransformation = true;            109   bool isLorentzTransformation = true;
110   double norm;                                    110   double norm;
111                                                   111 
112   d = ccol4;                                      112   d = ccol4;
113   norm = d.dot(d);                                113   norm = d.dot(d);
114   if (norm <= 0.0) {                              114   if (norm <= 0.0) {
115     isLorentzTransformation = false;              115     isLorentzTransformation = false;
116     if (norm == 0.0) {                            116     if (norm == 0.0) {
117       d = T_HAT4;       // Moot, but let's kee    117       d = T_HAT4;       // Moot, but let's keep going...
118       norm = 1.0;                                 118       norm = 1.0;
119     }                                             119     }
120   }                                               120   }
121   d /= norm;                                      121   d /= norm;
122                                                   122 
123   c = ccol3 - ccol3.dot(d) * d;                   123   c = ccol3 - ccol3.dot(d) * d;
124   norm = -c.dot(c);                               124   norm = -c.dot(c);
125   if (norm <= 0.0) {                              125   if (norm <= 0.0) {
126     isLorentzTransformation = false;              126     isLorentzTransformation = false;
127     if (norm == 0.0) {                            127     if (norm == 0.0) {
128       c = Z_HAT4;       // Moot                   128       c = Z_HAT4;       // Moot
129       norm = 1.0;                                 129       norm = 1.0;
130     }                                             130     }
131   }                                               131   }
132   c /= norm;                                      132   c /= norm;
133                                                   133 
134   b = ccol2 + ccol2.dot(c) * c - ccol2.dot(d)     134   b = ccol2 + ccol2.dot(c) * c - ccol2.dot(d) * d;
135   norm = -b.dot(b);                               135   norm = -b.dot(b);
136   if (norm <= 0.0) {                              136   if (norm <= 0.0) {
137     isLorentzTransformation = false;              137     isLorentzTransformation = false;
138     if (norm == 0.0) {                            138     if (norm == 0.0) {
139       b = Y_HAT4;       // Moot                   139       b = Y_HAT4;       // Moot
140       norm = 1.0;                                 140       norm = 1.0;
141     }                                             141     }
142   }                                               142   }
143   b /= norm;                                      143   b /= norm;
144                                                   144 
145   a = ccol1 + ccol1.dot(b) * b + ccol1.dot(c)     145   a = ccol1 + ccol1.dot(b) * b + ccol1.dot(c) * c - ccol1.dot(d) * d;
146   norm = -a.dot(a);                               146   norm = -a.dot(a);
147   if (norm <= 0.0) {                              147   if (norm <= 0.0) {
148     isLorentzTransformation = false;              148     isLorentzTransformation = false;
149     if (norm == 0.0) {                            149     if (norm == 0.0) {
150       a = X_HAT4;       // Moot                   150       a = X_HAT4;       // Moot
151       norm = 1.0;                                 151       norm = 1.0;
152     }                                             152     }
153   }                                               153   }
154   a /= norm;                                      154   a /= norm;
155                                                   155 
156   if ( !isLorentzTransformation ) {               156   if ( !isLorentzTransformation ) {
157       std::cerr << "HepLorentzRotation::set()     157       std::cerr << "HepLorentzRotation::set() - "
158         << "cols 1-4 supplied to define transf    158         << "cols 1-4 supplied to define transformation form either \n"
159         << "       a boosted reflection or a t    159         << "       a boosted reflection or a tachyonic transformation -- \n"
160         << "       transformation will be set     160         << "       transformation will be set to Identity " << std::endl;
161                                                   161 
162     *this = HepLorentzRotation();                 162     *this = HepLorentzRotation();
163   }                                               163   }
164                                                   164 
165   if ( isLorentzTransformation ) {                165   if ( isLorentzTransformation ) {
166     mxx = a.x(); myx = a.y(); mzx = a.z(); mtx    166     mxx = a.x(); myx = a.y(); mzx = a.z(); mtx = a.t();
167     mxy = b.x(); myy = b.y(); mzy = b.z(); mty    167     mxy = b.x(); myy = b.y(); mzy = b.z(); mty = b.t();
168     mxz = c.x(); myz = c.y(); mzz = c.z(); mtz    168     mxz = c.x(); myz = c.y(); mzz = c.z(); mtz = c.t();
169     mxt = d.x(); myt = d.y(); mzt = d.z(); mtt    169     mxt = d.x(); myt = d.y(); mzt = d.z(); mtt = d.t();
170   }                                               170   }
171                                                   171 
172   HepLorentzVector::setMetric (savedMetric);      172   HepLorentzVector::setMetric (savedMetric);
173   return *this;                                   173   return *this;
174                                                   174 
175 } // set ( col1, col2, col3, col4 )               175 } // set ( col1, col2, col3, col4 )
176                                                   176 
177 HepLorentzRotation & HepLorentzRotation::setRo    177 HepLorentzRotation & HepLorentzRotation::setRows
178    (const HepLorentzVector & rrow1,               178    (const HepLorentzVector & rrow1,
179           const HepLorentzVector & rrow2,         179           const HepLorentzVector & rrow2,
180     const HepLorentzVector & rrow3,               180     const HepLorentzVector & rrow3,
181     const HepLorentzVector & rrow4) {             181     const HepLorentzVector & rrow4) {
182   // Set based on using those rows as columns:    182   // Set based on using those rows as columns:
183   set (rrow1, rrow2, rrow3, rrow4);               183   set (rrow1, rrow2, rrow3, rrow4);
184   // Now transpose in place:                      184   // Now transpose in place:
185   double q1, q2, q3;                              185   double q1, q2, q3;
186   q1  = mxy;  q2  = mxz;  q3  = mxt;              186   q1  = mxy;  q2  = mxz;  q3  = mxt;
187   mxy = myx;  mxz = mzx;  mxt = mtx;              187   mxy = myx;  mxz = mzx;  mxt = mtx;
188   myx = q1;   mzx = q2;   mtx = q3;               188   myx = q1;   mzx = q2;   mtx = q3;
189   q1  = myz;  q2  = myt;  q3  = mzt;              189   q1  = myz;  q2  = myt;  q3  = mzt;
190   myz = mzy;  myt = mty;  mzt = mtz;              190   myz = mzy;  myt = mty;  mzt = mtz;
191   mzy = q1;   mty = q2;   mtz = q3;               191   mzy = q1;   mty = q2;   mtz = q3;
192   return *this;                                   192   return *this;
193 } // LorentzTransformation::setRows(row1 ... r    193 } // LorentzTransformation::setRows(row1 ... row4)
194                                                   194 
195 HepLorentzRotation::HepLorentzRotation ( const    195 HepLorentzRotation::HepLorentzRotation ( const HepLorentzVector & ccol1,
196                              const HepLorentzV    196                              const HepLorentzVector & ccol2,
197                              const HepLorentzV    197                              const HepLorentzVector & ccol3,
198                              const HepLorentzV    198                              const HepLorentzVector & ccol4 ) 
199 {                                                 199 {
200   set ( ccol1, ccol2, ccol3, ccol4 );             200   set ( ccol1, ccol2, ccol3, ccol4 );
201 }                                                 201 }
202                                                   202 
203 }  // namespace CLHEP                             203 }  // namespace CLHEP
204                                                   204 
205                                                   205