Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/error_propagation/src/G4ErrorMatrix.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 /error_propagation/src/G4ErrorMatrix.cc (Version 11.3.0) and /error_propagation/src/G4ErrorMatrix.cc (Version 9.2.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4ErrorMatrix.cc,v 1.3 2007/06/21 15:04:06 gunter Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-09-02 $
 26 //                                                 28 //
 27 // -------------------------------------------     29 // ------------------------------------------------------------
 28 //      GEANT 4 class implementation file          30 //      GEANT 4 class implementation file
 29 // -------------------------------------------     31 // ------------------------------------------------------------
 30                                                    32 
 31 #include "globals.hh"                              33 #include "globals.hh"
 32                                                    34 
 33 #include <cmath>                                   35 #include <cmath>
 34 #include <iostream>                                36 #include <iostream>
 35                                                    37 
 36 #include "G4ErrorMatrix.hh"                        38 #include "G4ErrorMatrix.hh"
 37 #include "G4ErrorSymMatrix.hh"                     39 #include "G4ErrorSymMatrix.hh"
 38                                                    40 
 39 // Simple operation for all elements               41 // Simple operation for all elements
 40                                                    42 
 41 #define SIMPLE_UOP(OPER)                       <<  43 #define SIMPLE_UOP(OPER)                            \
 42   G4ErrorMatrixIter a = m.begin();             <<  44     G4ErrorMatrixIter a=m.begin();                  \
 43   G4ErrorMatrixIter e = m.end();               <<  45     G4ErrorMatrixIter e=m.end();                    \
 44   for(; a != e; a++)                           <<  46    for(;a!=e; a++) (*a) OPER t;
 45     (*a) OPER t;                               <<  47 
 46                                                <<  48 #define SIMPLE_BOP(OPER)                            \
 47 #define SIMPLE_BOP(OPER)                       <<  49     G4ErrorMatrixIter a=m.begin();                  \
 48   G4ErrorMatrixIter a      = m.begin();        <<  50     G4ErrorMatrixConstIter b=m2.m.begin();          \
 49   G4ErrorMatrixConstIter b = mat2.m.begin();   <<  51     G4ErrorMatrixIter e=m.end();                    \
 50   G4ErrorMatrixIter e      = m.end();          <<  52    for(;a!=e; a++, b++) (*a) OPER (*b);
 51   for(; a != e; a++, b++)                      <<  53 
 52     (*a) OPER(*b);                             <<  54 #define SIMPLE_TOP(OPER)                            \
 53                                                <<  55     G4ErrorMatrixConstIter a=m1.m.begin();          \
 54 #define SIMPLE_TOP(OPER)                       <<  56     G4ErrorMatrixConstIter b=m2.m.begin();          \
 55   G4ErrorMatrixConstIter a = mat1.m.begin();   <<  57     G4ErrorMatrixIter t=mret.m.begin();             \
 56   G4ErrorMatrixConstIter b = mat2.m.begin();   <<  58     G4ErrorMatrixConstIter e=m1.m.end();            \
 57   G4ErrorMatrixIter t      = mret.m.begin();   <<  59    for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
 58   G4ErrorMatrixConstIter e = mat1.m.end();     << 
 59   for(; a != e; a++, b++, t++)                 << 
 60     (*t) = (*a) OPER(*b);                      << 
 61                                                    60 
 62 // Static functions.                               61 // Static functions.
 63                                                    62 
 64 #define CHK_DIM_2(r1, r2, c1, c2, fun)         <<  63 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
 65   if(r1 != r2 || c1 != c2)                     <<  64    if (r1!=r2 || c1!=c2)  { \
 66   {                                            <<  65     G4ErrorMatrix::error("Range error in Matrix function " #fun "(1).");  \
 67     G4ErrorMatrix::error("Range error in Matri <<  66    }
 68   }                                            <<  67 
 69                                                <<  68 #define CHK_DIM_1(c1,r2,fun) \
 70 #define CHK_DIM_1(c1, r2, fun)                 <<  69    if (c1!=r2) { \
 71   if(c1 != r2)                                 <<  70      G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
 72   {                                            <<  71    }
 73     G4ErrorMatrix::error("Range error in Matri << 
 74   }                                            << 
 75                                                    72 
 76 // Constructors. (Default constructors are inl     73 // Constructors. (Default constructors are inlined and in .icc file)
 77                                                    74 
 78 G4ErrorMatrix::G4ErrorMatrix(G4int p, G4int q) <<  75 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q)
 79   : m(p * q)                                   <<  76    : m(p*q), nrow(p), ncol(q)
 80   , nrow(p)                                    << 
 81   , ncol(q)                                    << 
 82 {                                                  77 {
 83   size = nrow * ncol;                              78   size = nrow * ncol;
 84 }                                                  79 }
 85                                                    80 
 86 G4ErrorMatrix::G4ErrorMatrix(G4int p, G4int q, <<  81 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q,G4int init)
 87   : m(p * q)                                   <<  82    : m(p*q), nrow(p), ncol(q)
 88   , nrow(p)                                    << 
 89   , ncol(q)                                    << 
 90 {                                                  83 {
 91   size = nrow * ncol;                          <<  84    size = nrow * ncol;
 92                                                    85 
 93   if(size > 0)                                 <<  86    if (size > 0)
 94   {                                            <<  87    {
 95     switch(init)                               <<  88       switch(init)
 96     {                                          <<  89       {
 97       case 0:                                      90       case 0:
 98         break;                                 <<  91          break;
 99                                                    92 
100       case 1:                                      93       case 1:
101       {                                        <<  94          {
102         if(ncol == nrow)                       <<  95             if ( ncol == nrow )
103         {                                      <<  96             {
104           G4ErrorMatrixIter a = m.begin();     <<  97                G4ErrorMatrixIter a = m.begin();
105           G4ErrorMatrixIter b = m.end();       <<  98                G4ErrorMatrixIter b = m.end();
106           for(; a < b; a += (ncol + 1))        <<  99                for( ; a<b; a+=(ncol+1)) *a = 1.0;
107             *a = 1.0;                          << 100             } else {
108         }                                      << 101                error("Invalid dimension in G4ErrorMatrix(G4int,G4int,1).");
109         else                                   << 102             }
110         {                                      << 103             break;
111           error("Invalid dimension in G4ErrorM << 104          }
112         }                                      << 
113         break;                                 << 
114       }                                        << 
115       default:                                    105       default:
116         error("G4ErrorMatrix: initialization m << 106          error("G4ErrorMatrix: initialization must be either 0 or 1.");
117     }                                          << 107       }
118   }                                            << 108    }
119 }                                                 109 }
120                                                   110 
121 //                                                111 //
122 // Destructor                                     112 // Destructor
123 //                                                113 //
124 G4ErrorMatrix::~G4ErrorMatrix() {}             << 114 G4ErrorMatrix::~G4ErrorMatrix()
                                                   >> 115 {
                                                   >> 116 }
125                                                   117 
126 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatr << 118 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatrix &m1)
127   : m(mat1.size)                               << 119    : m(m1.size), nrow(m1.nrow), ncol(m1.ncol), size(m1.size)
128   , nrow(mat1.nrow)                            << 
129   , ncol(mat1.ncol)                            << 
130   , size(mat1.size)                            << 
131 {                                                 120 {
132   m = mat1.m;                                  << 121    m = m1.m;
133 }                                                 122 }
134                                                   123 
135 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymM << 124 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymMatrix &m1)
136   : m(mat1.nrow * mat1.nrow)                   << 125    : m(m1.nrow*m1.nrow), nrow(m1.nrow), ncol(m1.nrow)
137   , nrow(mat1.nrow)                            << 
138   , ncol(mat1.nrow)                            << 
139 {                                                 126 {
140   size = nrow * ncol;                          << 127    size = nrow * ncol;
141                                                   128 
142   G4int n                    = ncol;           << 129    G4int n = ncol;
143   G4ErrorMatrixConstIter sjk = mat1.m.begin(); << 130    G4ErrorMatrixConstIter sjk = m1.m.begin();
144   G4ErrorMatrixIter m1j      = m.begin();      << 131    G4ErrorMatrixIter m1j = m.begin();
145   G4ErrorMatrixIter mj       = m.begin();      << 132    G4ErrorMatrixIter mj  = m.begin();
146   // j >= k                                    << 133    // j >= k
147   for(G4int j = 1; j <= nrow; j++)             << 134    for(G4int j=1;j<=nrow;j++)
148   {                                            << 135    {
149     G4ErrorMatrixIter mjk = mj;                << 136       G4ErrorMatrixIter mjk = mj;
150     G4ErrorMatrixIter mkj = m1j;               << 137       G4ErrorMatrixIter mkj = m1j;
151     for(G4int k = 1; k <= j; k++)              << 138       for(G4int k=1;k<=j;k++)
152     {                                          << 139       {
153       *(mjk++) = *sjk;                         << 140          *(mjk++) = *sjk;
154       if(j != k)                               << 141          if(j!=k) *mkj = *sjk;
155         *mkj = *sjk;                           << 142          sjk++;
156       sjk++;                                   << 143          mkj += n;
157       mkj += n;                                << 144       }
158     }                                          << 145       mj += n;
159     mj += n;                                   << 146       m1j++;
160     m1j++;                                     << 147    }
161   }                                            << 
162 }                                                 148 }
163                                                   149 
164 //                                                150 //
165 //                                                151 //
166 // Sub matrix                                     152 // Sub matrix
167 //                                                153 //
168 //                                                154 //
169                                                   155 
170 G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row << 156 G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row, G4int max_row,
171                                  G4int max_col << 157                          G4int min_col,G4int max_col) const
172 {                                                 158 {
173   G4ErrorMatrix mret(max_row - min_row + 1, ma << 159   G4ErrorMatrix mret(max_row-min_row+1,max_col-min_col+1);
174   if(max_row > num_row() || max_col > num_col( << 160   if(max_row > num_row() || max_col >num_col())
175   {                                            << 161     { error("G4ErrorMatrix::sub: Index out of range"); }
176     error("G4ErrorMatrix::sub: Index out of ra << 162   G4ErrorMatrixIter a = mret.m.begin();
177   }                                            << 163   G4int nc = num_col();
178   G4ErrorMatrixIter a       = mret.m.begin();  << 
179   G4int nc                  = num_col();       << 
180   G4ErrorMatrixConstIter b1 = m.begin() + (min    164   G4ErrorMatrixConstIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
181                                                << 165   
182   for(G4int irow = 1; irow <= mret.num_row();  << 166   for(G4int irow=1; irow<=mret.num_row(); irow++)
183   {                                               167   {
184     G4ErrorMatrixConstIter brc = b1;              168     G4ErrorMatrixConstIter brc = b1;
185     for(G4int icol = 1; icol <= mret.num_col() << 169     for(G4int icol=1; icol<=mret.num_col(); icol++)
186     {                                             170     {
187       *(a++) = *(brc++);                          171       *(a++) = *(brc++);
188     }                                             172     }
189     b1 += nc;                                     173     b1 += nc;
190   }                                               174   }
191   return mret;                                    175   return mret;
192 }                                                 176 }
193                                                   177 
194 void G4ErrorMatrix::sub(G4int row, G4int col,  << 178 void G4ErrorMatrix::sub(G4int row,G4int col,const G4ErrorMatrix &m1)
195 {                                                 179 {
196   if(row < 1 || row + mat1.num_row() - 1 > num << 180   if(row <1 || row+m1.num_row()-1 > num_row() || 
197      col + mat1.num_col() - 1 > num_col())     << 181      col <1 || col+m1.num_col()-1 > num_col()   )
198   {                                            << 182     { error("G4ErrorMatrix::sub: Index out of range"); }
199     error("G4ErrorMatrix::sub: Index out of ra << 183   G4ErrorMatrixConstIter a = m1.m.begin();
200   }                                            << 184   G4int nc = num_col();
201   G4ErrorMatrixConstIter a = mat1.m.begin();   << 185   G4ErrorMatrixIter b1 = m.begin() + (row - 1) * nc + col - 1;
202   G4int nc                 = num_col();        << 186   
203   G4ErrorMatrixIter b1     = m.begin() + (row  << 187   for(G4int irow=1; irow<=m1.num_row(); irow++)
204                                                << 
205   for(G4int irow = 1; irow <= mat1.num_row();  << 
206   {                                               188   {
207     G4ErrorMatrixIter brc = b1;                   189     G4ErrorMatrixIter brc = b1;
208     for(G4int icol = 1; icol <= mat1.num_col() << 190     for(G4int icol=1; icol<=m1.num_col(); icol++)
209     {                                             191     {
210       *(brc++) = *(a++);                          192       *(brc++) = *(a++);
211     }                                             193     }
212     b1 += nc;                                     194     b1 += nc;
213   }                                               195   }
214 }                                                 196 }
215                                                   197 
216 //                                                198 //
217 // Direct sum of two matricies                    199 // Direct sum of two matricies
218 //                                                200 //
219                                                   201 
220 G4ErrorMatrix dsum(const G4ErrorMatrix& mat1,  << 202 G4ErrorMatrix dsum(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
221 {                                                 203 {
222   G4ErrorMatrix mret(mat1.num_row() + mat2.num << 204   G4ErrorMatrix mret(m1.num_row() + m2.num_row(), m1.num_col() + m2.num_col(),
223                      mat1.num_col() + mat2.num << 205                  0);
224   mret.sub(1, 1, mat1);                        << 206   mret.sub(1,1,m1);
225   mret.sub(mat1.num_row() + 1, mat1.num_col()  << 207   mret.sub(m1.num_row()+1,m1.num_col()+1,m2);
226   return mret;                                    208   return mret;
227 }                                                 209 }
228                                                   210 
229 /* -------------------------------------------    211 /* -----------------------------------------------------------------------
230    This section contains support routines for     212    This section contains support routines for matrix.h. This section contains
231    The two argument functions +,-. They call t    213    The two argument functions +,-. They call the copy constructor and +=,-=.
232    -------------------------------------------    214    ----------------------------------------------------------------------- */
233                                                   215 
234 G4ErrorMatrix G4ErrorMatrix::operator-() const << 216 G4ErrorMatrix G4ErrorMatrix::operator- () const 
235 {                                                 217 {
236   G4ErrorMatrix mat2(nrow, ncol);              << 218    G4ErrorMatrix m2(nrow, ncol);
237   G4ErrorMatrixConstIter a = m.begin();        << 219    G4ErrorMatrixConstIter a=m.begin();
238   G4ErrorMatrixIter b      = mat2.m.begin();   << 220    G4ErrorMatrixIter b=m2.m.begin();
239   G4ErrorMatrixConstIter e = m.end();          << 221    G4ErrorMatrixConstIter e=m.end();
240   for(; a < e; a++, b++)                       << 222    for(;a<e; a++, b++) (*b) = -(*a);
241     (*b) = -(*a);                              << 223    return m2;
242   return mat2;                                 << 
243 }                                                 224 }
244                                                   225 
245 G4ErrorMatrix operator+(const G4ErrorMatrix& m << 226 G4ErrorMatrix operator+(const G4ErrorMatrix &m1,const G4ErrorMatrix &m2)
246 {                                                 227 {
247   G4ErrorMatrix mret(mat1.nrow, mat1.ncol);    << 228   G4ErrorMatrix mret(m1.nrow, m1.ncol);
248   CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 229   CHK_DIM_2(m1.num_row(),m2.num_row(), m1.num_col(),m2.num_col(),+);
249   SIMPLE_TOP(+)                                   230   SIMPLE_TOP(+)
250   return mret;                                    231   return mret;
251 }                                                 232 }
252                                                   233 
253 //                                                234 //
254 // operator -                                     235 // operator -
255 //                                                236 //
256                                                   237 
257 G4ErrorMatrix operator-(const G4ErrorMatrix& m << 238 G4ErrorMatrix operator-(const G4ErrorMatrix &m1,const G4ErrorMatrix &m2)
258 {                                                 239 {
259   G4ErrorMatrix mret(mat1.num_row(), mat1.num_ << 240   G4ErrorMatrix mret(m1.num_row(), m1.num_col());
260   CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 241   CHK_DIM_2(m1.num_row(),m2.num_row(),
                                                   >> 242                          m1.num_col(),m2.num_col(),-);
261   SIMPLE_TOP(-)                                   243   SIMPLE_TOP(-)
262   return mret;                                    244   return mret;
263 }                                                 245 }
264                                                   246 
265 /* -------------------------------------------    247 /* -----------------------------------------------------------------------
266    This section contains support routines for     248    This section contains support routines for matrix.h. This file contains
267    The two argument functions *,/. They call c    249    The two argument functions *,/. They call copy constructor and then /=,*=.
268    -------------------------------------------    250    ----------------------------------------------------------------------- */
269                                                   251 
270 G4ErrorMatrix operator/(const G4ErrorMatrix& m << 252 G4ErrorMatrix operator/(const G4ErrorMatrix &m1,G4double t)
271 {                                                 253 {
272   G4ErrorMatrix mret(mat1);                    << 254   G4ErrorMatrix mret(m1);
273   mret /= t;                                      255   mret /= t;
274   return mret;                                    256   return mret;
275 }                                                 257 }
276                                                   258 
277 G4ErrorMatrix operator*(const G4ErrorMatrix& m << 259 G4ErrorMatrix operator*(const G4ErrorMatrix &m1,G4double t)
278 {                                                 260 {
279   G4ErrorMatrix mret(mat1);                    << 261   G4ErrorMatrix mret(m1);
280   mret *= t;                                      262   mret *= t;
281   return mret;                                    263   return mret;
282 }                                                 264 }
283                                                   265 
284 G4ErrorMatrix operator*(G4double t, const G4Er << 266 G4ErrorMatrix operator*(G4double t,const G4ErrorMatrix &m1)
285 {                                                 267 {
286   G4ErrorMatrix mret(mat1);                    << 268   G4ErrorMatrix mret(m1);
287   mret *= t;                                      269   mret *= t;
288   return mret;                                    270   return mret;
289 }                                                 271 }
290                                                   272 
291 G4ErrorMatrix operator*(const G4ErrorMatrix& m << 273 G4ErrorMatrix operator*(const G4ErrorMatrix &m1,const G4ErrorMatrix &m2)
292 {                                                 274 {
293   // initialize matrix to 0.0                     275   // initialize matrix to 0.0
294   G4ErrorMatrix mret(mat1.nrow, mat2.ncol, 0); << 276   G4ErrorMatrix mret(m1.nrow,m2.ncol,0);
295   CHK_DIM_1(mat1.ncol, mat2.nrow, *);          << 277   CHK_DIM_1(m1.ncol,m2.nrow,*);
296                                                   278 
297   G4int m1cols = mat1.ncol;                    << 279   G4int m1cols = m1.ncol;
298   G4int m2cols = mat2.ncol;                    << 280   G4int m2cols = m2.ncol;
299                                                   281 
300   for(G4int i = 0; i < mat1.nrow; i++)         << 282   for (G4int i=0; i<m1.nrow; i++)
301   {                                               283   {
302     for(G4int j = 0; j < m1cols; j++)          << 284      for (G4int j=0; j<m1cols; j++) 
303     {                                          << 285      {
304       G4double temp        = mat1.m[i * m1cols << 286          G4double temp = m1.m[i*m1cols+j];
305       G4ErrorMatrixIter pt = mret.m.begin() +  << 287          G4ErrorMatrixIter pt = mret.m.begin() + i*m2cols;
306                                                << 288 
307       // Loop over k (the column index in matr << 289         // Loop over k (the column index in matrix m2)
308       G4ErrorMatrixConstIter pb           = ma << 290         G4ErrorMatrixConstIter pb = m2.m.begin() + m2cols*j;
309       const G4ErrorMatrixConstIter pblast = pb << 291         const G4ErrorMatrixConstIter pblast = pb + m2cols;
310       while(pb < pblast)  // Loop checking, 06 << 292         while (pb < pblast)
311       {                                        << 293         {
312         (*pt) += temp * (*pb);                 << 294            (*pt) += temp * (*pb);
313         pb++;                                  << 295            pb++;
314         pt++;                                  << 296            pt++;
315       }                                        << 297         }
316     }                                          << 298      }
317   }                                               299   }
318   return mret;                                    300   return mret;
319 }                                                 301 }
320                                                   302 
321 /* -------------------------------------------    303 /* -----------------------------------------------------------------------
322    This section contains the assignment and in    304    This section contains the assignment and inplace operators =,+=,-=,*=,/=.
323    -------------------------------------------    305    ----------------------------------------------------------------------- */
324                                                   306 
325 G4ErrorMatrix& G4ErrorMatrix::operator+=(const << 307 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorMatrix &m2)
326 {                                                 308 {
327   CHK_DIM_2(num_row(), mat2.num_row(), num_col << 309   CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
328   SIMPLE_BOP(+=)                                  310   SIMPLE_BOP(+=)
329   return (*this);                                 311   return (*this);
330 }                                                 312 }
331                                                   313 
332 G4ErrorMatrix& G4ErrorMatrix::operator-=(const << 314 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorMatrix &m2)
333 {                                                 315 {
334   CHK_DIM_2(num_row(), mat2.num_row(), num_col << 316   CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
335   SIMPLE_BOP(-=)                                  317   SIMPLE_BOP(-=)
336   return (*this);                                 318   return (*this);
337 }                                                 319 }
338                                                   320 
339 G4ErrorMatrix& G4ErrorMatrix::operator/=(G4dou << 321 G4ErrorMatrix & G4ErrorMatrix::operator/=(G4double t)
340 {                                                 322 {
341   SIMPLE_UOP(/=)                                  323   SIMPLE_UOP(/=)
342   return (*this);                                 324   return (*this);
343 }                                                 325 }
344                                                   326 
345 G4ErrorMatrix& G4ErrorMatrix::operator*=(G4dou << 327 G4ErrorMatrix & G4ErrorMatrix::operator*=(G4double t)
346 {                                                 328 {
347   SIMPLE_UOP(*=)                                  329   SIMPLE_UOP(*=)
348   return (*this);                                 330   return (*this);
349 }                                                 331 }
350                                                   332 
351 G4ErrorMatrix& G4ErrorMatrix::operator=(const  << 333 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorMatrix &m1)
352 {                                                 334 {
353   if(&mat1 == this)                            << 335    if(m1.nrow*m1.ncol != size) //??fixme?? m1.size != size
354   {                                            << 336    {
355     return *this;                              << 337       size = m1.nrow * m1.ncol;
356   }                                            << 338       m.resize(size); //??fixme?? if (size < m1.size) m.resize(m1.size);
357                                                << 339    }
358   if(mat1.nrow * mat1.ncol != size)  //??fixme << 340    nrow = m1.nrow;
359   {                                            << 341    ncol = m1.ncol;
360     size = mat1.nrow * mat1.ncol;              << 342    m = m1.m;
361     m.resize(size);  //??fixme?? if (size < ma << 343    return (*this);
362   }                                            << 
363   nrow = mat1.nrow;                            << 
364   ncol = mat1.ncol;                            << 
365   m    = mat1.m;                               << 
366   return (*this);                              << 
367 }                                                 344 }
368                                                   345 
                                                   >> 346 
369 // Print the G4ErrorMatrix.                       347 // Print the G4ErrorMatrix.
370                                                   348 
371 std::ostream& operator<<(std::ostream& os, con << 349 std::ostream& operator<<(std::ostream &s, const G4ErrorMatrix &q)
372 {                                                 350 {
373   os << "\n";                                  << 351   s << "\n";
374                                                   352 
375   // Fixed format needs 3 extra characters for    353   // Fixed format needs 3 extra characters for field,
376   // while scientific needs 7                     354   // while scientific needs 7
377                                                   355 
378   std::size_t width;                           << 356   G4int width;
379   if(os.flags() & std::ios::fixed)             << 357   if(s.flags() & std::ios::fixed)
380   {                                            << 358     { width = s.precision()+3; }
381     width = os.precision() + 3;                << 
382   }                                            << 
383   else                                            359   else
384   {                                            << 360     { width = s.precision()+7; }
385     width = os.precision() + 7;                << 361   for(G4int irow = 1; irow<= q.num_row(); irow++)
386   }                                            << 
387   for(G4int irow = 1; irow <= q.num_row(); ++i << 
388   {                                            << 
389     for(G4int icol = 1; icol <= q.num_col(); + << 
390     {                                             362     {
391       os.width(width);                         << 363       for(G4int icol = 1; icol <= q.num_col(); icol++)
392       os << q(irow, icol) << " ";              << 364         {
                                                   >> 365           s.width(width);
                                                   >> 366           s << q(irow,icol) << " ";
                                                   >> 367         }
                                                   >> 368       s << G4endl;
393     }                                             369     }
394     os << G4endl;                              << 370   return s;
395   }                                            << 
396   return os;                                   << 
397 }                                                 371 }
398                                                   372 
399 G4ErrorMatrix G4ErrorMatrix::T() const            373 G4ErrorMatrix G4ErrorMatrix::T() const
400 {                                                 374 {
401   G4ErrorMatrix mret(ncol, nrow);              << 375    G4ErrorMatrix mret(ncol,nrow);
402   G4ErrorMatrixConstIter pl  = m.end();        << 376     G4ErrorMatrixConstIter pl = m.end();
403   G4ErrorMatrixConstIter pme = m.begin();      << 377     G4ErrorMatrixConstIter pme = m.begin();
404   G4ErrorMatrixIter pt       = mret.m.begin(); << 378     G4ErrorMatrixIter pt = mret.m.begin();
405   G4ErrorMatrixIter ptl      = mret.m.end();   << 379     G4ErrorMatrixIter ptl = mret.m.end();
406   for(; pme < pl; pme++, pt += nrow)           << 380    for (; pme < pl; pme++, pt+=nrow)
407   {                                            << 381    {
408     if(pt >= ptl)                              << 382       if (pt >= ptl) { pt -= (size-1); }
409     {                                          << 383       (*pt) = (*pme);
410       pt -= (size - 1);                        << 384    }
411     }                                          << 385    return mret;
412     (*pt) = (*pme);                            << 
413   }                                            << 
414   return mret;                                 << 
415 }                                                 386 }
416                                                   387 
417 G4ErrorMatrix G4ErrorMatrix::apply(G4double (* << 388 G4ErrorMatrix
                                                   >> 389 G4ErrorMatrix::apply(G4double (*f)(G4double, G4int, G4int)) const
418 {                                                 390 {
419   G4ErrorMatrix mret(num_row(), num_col());    << 391   G4ErrorMatrix mret(num_row(),num_col());
420   G4ErrorMatrixConstIter a = m.cbegin();       << 392   G4ErrorMatrixConstIter a = m.begin();
421   G4ErrorMatrixIter b      = mret.m.begin();   << 393   G4ErrorMatrixIter b = mret.m.begin();
422   for(G4int ir = 1; ir <= num_row(); ++ir)     << 394   for(G4int ir=1;ir<=num_row();ir++)
423   {                                               395   {
424     for(G4int ic = 1; ic <= num_col(); ++ic)   << 396     for(G4int ic=1;ic<=num_col();ic++)
425     {                                             397     {
426       *(b++) = (*f)(*(a++), ir, ic);              398       *(b++) = (*f)(*(a++), ir, ic);
427     }                                             399     }
428   }                                               400   }
429   return mret;                                    401   return mret;
430 }                                                 402 }
431                                                   403 
432 G4int G4ErrorMatrix::dfinv_matrix(G4int* ir)   << 404 G4int G4ErrorMatrix::dfinv_matrix(G4int *ir) {
433 {                                              << 405   if (num_col()!=num_row())
434   if(num_col() != num_row())                   << 406     { error("dfinv_matrix: G4ErrorMatrix is not NxN"); }
435   {                                            << 
436     error("dfinv_matrix: G4ErrorMatrix is not  << 
437   }                                            << 
438   G4int n = num_col();                            407   G4int n = num_col();
439   if(n == 1)                                   << 408   if (n==1) { return 0; }
440   {                                            << 
441     return 0;                                  << 
442   }                                            << 
443                                                   409 
444   G4double s31, s32;                              410   G4double s31, s32;
445   G4double s33, s34;                              411   G4double s33, s34;
446                                                   412 
447   G4ErrorMatrixIter m11 = m.begin();              413   G4ErrorMatrixIter m11 = m.begin();
448   G4ErrorMatrixIter m12 = m11 + 1;                414   G4ErrorMatrixIter m12 = m11 + 1;
449   G4ErrorMatrixIter m21 = m11 + n;                415   G4ErrorMatrixIter m21 = m11 + n;
450   G4ErrorMatrixIter m22 = m12 + n;                416   G4ErrorMatrixIter m22 = m12 + n;
451   *m21                  = -(*m22) * (*m11) * ( << 417   *m21 = -(*m22) * (*m11) * (*m21);
452   *m12                  = -(*m12);             << 418   *m12 = -(*m12);
453   if(n > 2)                                    << 419   if (n>2)
454   {                                               420   {
455     G4ErrorMatrixIter mi    = m.begin() + 2 *  << 421     G4ErrorMatrixIter mi = m.begin() + 2 * n;
456     G4ErrorMatrixIter mii   = m.begin() + 2 *  << 422     G4ErrorMatrixIter mii= m.begin() + 2 * n + 2;
457     G4ErrorMatrixIter mimim = m.begin() + n +     423     G4ErrorMatrixIter mimim = m.begin() + n + 1;
458     for(G4int i = 3; i <= n; i++)              << 424     for (G4int i=3;i<=n;i++)
459     {                                             425     {
460       G4int im2             = i - 2;           << 426       G4int im2 = i - 2;
461       G4ErrorMatrixIter mj  = m.begin();       << 427       G4ErrorMatrixIter mj = m.begin();
462       G4ErrorMatrixIter mji = mj + i - 1;         428       G4ErrorMatrixIter mji = mj + i - 1;
463       G4ErrorMatrixIter mij = mi;                 429       G4ErrorMatrixIter mij = mi;
464       for(G4int j = 1; j <= im2; j++)          << 430       for (G4int j=1;j<=im2;j++)
465       {                                        << 431       { 
466         s31                    = 0.0;          << 432         s31 = 0.0;
467         s32                    = *mji;         << 433         s32 = *mji;
468         G4ErrorMatrixIter mkj  = mj + j - 1;   << 434         G4ErrorMatrixIter mkj = mj + j - 1;
469         G4ErrorMatrixIter mik  = mi + j - 1;   << 435         G4ErrorMatrixIter mik = mi + j - 1;
470         G4ErrorMatrixIter mjkp = mj + j;          436         G4ErrorMatrixIter mjkp = mj + j;
471         G4ErrorMatrixIter mkpi = mj + n + i -     437         G4ErrorMatrixIter mkpi = mj + n + i - 1;
472         for(G4int k = j; k <= im2; k++)        << 438         for (G4int k=j;k<=im2;k++)
473         {                                         439         {
474           s31 += (*mkj) * (*(mik++));             440           s31 += (*mkj) * (*(mik++));
475           s32 += (*(mjkp++)) * (*mkpi);           441           s32 += (*(mjkp++)) * (*mkpi);
476           mkj += n;                               442           mkj += n;
477           mkpi += n;                              443           mkpi += n;
478         }                                         444         }
479         *mij = -(*mii) * (((*(mij - n))) * ((* << 445         *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
480         *mji = -s32;                              446         *mji = -s32;
481         mj += n;                                  447         mj += n;
482         mji += n;                                 448         mji += n;
483         mij++;                                    449         mij++;
484       }                                           450       }
485       *(mii - 1)   = -(*mii) * (*mimim) * (*(m << 451       *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
486       *(mimim + 1) = -(*(mimim + 1));          << 452       *(mimim+1) = -(*(mimim+1));
487       mi += n;                                    453       mi += n;
488       mimim += (n + 1);                        << 454       mimim += (n+1);
489       mii += (n + 1);                          << 455       mii += (n+1);
490     }                                             456     }
491   }                                               457   }
492   G4ErrorMatrixIter mi  = m.begin();           << 458   G4ErrorMatrixIter mi = m.begin();
493   G4ErrorMatrixIter mii = m.begin();              459   G4ErrorMatrixIter mii = m.begin();
494   for(G4int i = 1; i < n; i++)                 << 460   for (G4int i=1;i<n;i++)
495   {                                               461   {
496     G4int ni              = n - i;             << 462     G4int ni = n - i;
497     G4ErrorMatrixIter mij = mi;                   463     G4ErrorMatrixIter mij = mi;
498     G4int j;                                      464     G4int j;
499     for(j = 1; j <= i; j++)                    << 465     for (j=1; j<=i;j++)
500     {                                             466     {
501       s33                       = *mij;        << 467       s33 = *mij;
502       G4ErrorMatrixIter mikj    = mi + n + j - << 468       G4ErrorMatrixIter mikj = mi + n + j - 1;
503       G4ErrorMatrixIter miik    = mii + 1;     << 469       G4ErrorMatrixIter miik = mii + 1;
504       G4ErrorMatrixIter min_end = mi + n;         470       G4ErrorMatrixIter min_end = mi + n;
505       for(; miik < min_end;)                   << 471       for (;miik<min_end;)
506       {                                           472       {
507         s33 += (*mikj) * (*(miik++));             473         s33 += (*mikj) * (*(miik++));
508         mikj += n;                                474         mikj += n;
509       }                                           475       }
510       *(mij++) = s33;                             476       *(mij++) = s33;
511     }                                             477     }
512     for(j = 1; j <= ni; j++)                   << 478     for (j=1;j<=ni;j++)
513     {                                             479     {
514       s34                     = 0.0;           << 480       s34 = 0.0;
515       G4ErrorMatrixIter miik  = mii + j;       << 481       G4ErrorMatrixIter miik = mii + j;
516       G4ErrorMatrixIter mikij = mii + j * n +     482       G4ErrorMatrixIter mikij = mii + j * n + j;
517       for(G4int k = j; k <= ni; k++)           << 483       for (G4int k=j;k<=ni;k++)
518       {                                           484       {
519         s34 += *mikij * (*(miik++));              485         s34 += *mikij * (*(miik++));
520         mikij += n;                               486         mikij += n;
521       }                                           487       }
522       *(mii + j) = s34;                        << 488       *(mii+j) = s34;
523     }                                             489     }
524     mi += n;                                      490     mi += n;
525     mii += (n + 1);                            << 491     mii += (n+1);
526   }                                               492   }
527   G4int nxch = ir[n];                             493   G4int nxch = ir[n];
528   if(nxch == 0)                                << 494   if (nxch==0) return 0;
529     return 0;                                  << 495   for (G4int mm=1;mm<=nxch;mm++)
530   for(G4int mq = 1; mq <= nxch; mq++)          << 
531   {                                               496   {
532     G4int k               = nxch - mq + 1;     << 497     G4int k = nxch - mm + 1;
533     G4int ij              = ir[k];             << 498     G4int ij = ir[k];
534     G4int i               = ij >> 12;          << 499     G4int i = ij >> 12;
535     G4int j               = ij % 4096;         << 500     G4int j = ij%4096;
536     G4ErrorMatrixIter mki = m.begin() + i - 1;    501     G4ErrorMatrixIter mki = m.begin() + i - 1;
537     G4ErrorMatrixIter mkj = m.begin() + j - 1;    502     G4ErrorMatrixIter mkj = m.begin() + j - 1;
538     for(k = 1; k <= n; k++)                    << 503     for (k=1; k<=n;k++)
539     {                                             504     {
540       // 2/24/05 David Sachs fix of improper s    505       // 2/24/05 David Sachs fix of improper swap bug that was present
541       // for many years:                          506       // for many years:
542       G4double ti = *mki;  // 2/24/05          << 507       G4double ti = *mki; // 2/24/05
543       *mki        = *mkj;                      << 508       *mki = *mkj;
544       *mkj        = ti;  // 2/24/05            << 509       *mkj = ti;        // 2/24/05
545       mki += n;                                   510       mki += n;
546       mkj += n;                                   511       mkj += n;
547     }                                             512     }
548   }                                               513   }
549   return 0;                                       514   return 0;
550 }                                                 515 }
551                                                   516 
552 G4int G4ErrorMatrix::dfact_matrix(G4double& de << 517 G4int G4ErrorMatrix::dfact_matrix(G4double &det, G4int *ir)
553 {                                                 518 {
554   if(ncol != nrow)                             << 519   if (ncol!=nrow)
555     error("dfact_matrix: G4ErrorMatrix is not  << 520      error("dfact_matrix: G4ErrorMatrix is not NxN");
556                                                   521 
557   G4int ifail, jfail;                             522   G4int ifail, jfail;
558   G4int n = ncol;                                 523   G4int n = ncol;
559                                                   524 
560   G4double tf;                                    525   G4double tf;
561   G4double g1 = 1.0e-19, g2 = 1.0e19;             526   G4double g1 = 1.0e-19, g2 = 1.0e19;
562                                                   527 
563   G4double p, q, t;                               528   G4double p, q, t;
564   G4double s11, s12;                              529   G4double s11, s12;
565                                                   530 
566   G4double epsilon = 8 * DBL_EPSILON;          << 531   G4double epsilon = 8*DBL_EPSILON;
567   // could be set to zero (like it was before)    532   // could be set to zero (like it was before)
568   // but then the algorithm often doesn't dete    533   // but then the algorithm often doesn't detect
569   // that a matrix is singular                    534   // that a matrix is singular
570                                                   535 
571   G4int normal = 0, imposs = -1;                  536   G4int normal = 0, imposs = -1;
572   G4int jrange = 0, jover = 1, junder = -1;       537   G4int jrange = 0, jover = 1, junder = -1;
573   ifail                 = normal;              << 538   ifail = normal;
574   jfail                 = jrange;              << 539   jfail = jrange;
575   G4int nxch            = 0;                   << 540   G4int nxch = 0;
576   det                   = 1.0;                 << 541   det = 1.0;
577   G4ErrorMatrixIter mj  = m.begin();           << 542   G4ErrorMatrixIter mj = m.begin();
578   G4ErrorMatrixIter mjj = mj;                     543   G4ErrorMatrixIter mjj = mj;
579   for(G4int j = 1; j <= n; j++)                << 544   for (G4int j=1;j<=n;j++)
580   {                                               545   {
581     G4int k = j;                                  546     G4int k = j;
582     p       = (std::fabs(*mjj));               << 547     p = (std::fabs(*mjj));
583     if(j != n)                                 << 548     if (j!=n) {
584     {                                          << 549       G4ErrorMatrixIter mij = mj + n + j - 1; 
585       G4ErrorMatrixIter mij = mj + n + j - 1;  << 550       for (G4int i=j+1;i<=n;i++)
586       for(G4int i = j + 1; i <= n; i++)        << 
587       {                                           551       {
588         q = (std::fabs(*(mij)));                  552         q = (std::fabs(*(mij)));
589         if(q > p)                              << 553         if (q > p)
590         {                                         554         {
591           k = i;                                  555           k = i;
592           p = q;                                  556           p = q;
593         }                                         557         }
594         mij += n;                                 558         mij += n;
595       }                                           559       }
596       if(k == j)                               << 560       if (k==j)
597       {                                           561       {
598         if(p <= epsilon)                       << 562         if (p <= epsilon)
599         {                                         563         {
600           det   = 0;                           << 564           det = 0;
601           ifail = imposs;                         565           ifail = imposs;
602           jfail = jrange;                         566           jfail = jrange;
603           return ifail;                           567           return ifail;
604         }                                         568         }
605         det = -det;  // in this case the sign  << 569         det = -det; // in this case the sign of the determinant
606                      // must not change. So I  << 570                     // must not change. So I change it twice. 
607       }                                           571       }
608       G4ErrorMatrixIter mjl = mj;                 572       G4ErrorMatrixIter mjl = mj;
609       G4ErrorMatrixIter mkl = m.begin() + (k - << 573       G4ErrorMatrixIter mkl = m.begin() + (k-1)*n;
610       for(G4int l = 1; l <= n; l++)            << 574       for (G4int l=1;l<=n;l++)
611       {                                           575       {
612         tf       = *mjl;                       << 576         tf = *mjl;
613         *(mjl++) = *mkl;                          577         *(mjl++) = *mkl;
614         *(mkl++) = tf;                            578         *(mkl++) = tf;
615       }                                           579       }
616       nxch     = nxch + 1;  // this makes the  << 580       nxch = nxch + 1;  // this makes the determinant change its sign
617       ir[nxch] = (((j) << 12) + (k));          << 581       ir[nxch] = (((j)<<12)+(k));
618     }                                             582     }
619     else                                          583     else
620     {                                             584     {
621       if(p <= epsilon)                         << 585       if (p <= epsilon)
622       {                                           586       {
623         det   = 0.0;                           << 587         det = 0.0;
624         ifail = imposs;                           588         ifail = imposs;
625         jfail = jrange;                           589         jfail = jrange;
626         return ifail;                             590         return ifail;
627       }                                           591       }
628     }                                             592     }
629     det *= *mjj;                                  593     det *= *mjj;
630     *mjj = 1.0 / *mjj;                            594     *mjj = 1.0 / *mjj;
631     t    = (std::fabs(det));                   << 595     t = (std::fabs(det));
632     if(t < g1)                                 << 596     if (t < g1)
633     {                                             597     {
634       det = 0.0;                                  598       det = 0.0;
635       if(jfail == jrange)                      << 599       if (jfail == jrange) jfail = junder;
636         jfail = junder;                        << 
637     }                                             600     }
638     else if(t > g2)                            << 601     else if (t > g2)
639     {                                             602     {
640       det = 1.0;                                  603       det = 1.0;
641       if(jfail == jrange)                      << 604       if (jfail==jrange) jfail = jover;
642         jfail = jover;                         << 
643     }                                             605     }
644     if(j != n)                                 << 606     if (j!=n)
645     {                                             607     {
646       G4ErrorMatrixIter mk   = mj + n;         << 608       G4ErrorMatrixIter mk = mj + n;
647       G4ErrorMatrixIter mkjp = mk + j;            609       G4ErrorMatrixIter mkjp = mk + j;
648       G4ErrorMatrixIter mjk  = mj + j;         << 610       G4ErrorMatrixIter mjk = mj + j;
649       for(k = j + 1; k <= n; k++)              << 611       for (k=j+1;k<=n;k++)
650       {                                           612       {
651         s11 = -(*mjk);                         << 613         s11 = - (*mjk);
652         s12 = -(*mkjp);                        << 614         s12 = - (*mkjp);
653         if(j != 1)                             << 615         if (j!=1)
654         {                                         616         {
655           G4ErrorMatrixIter mik  = m.begin() + << 617           G4ErrorMatrixIter mik = m.begin() + k - 1;
656           G4ErrorMatrixIter mijp = m.begin() +    618           G4ErrorMatrixIter mijp = m.begin() + j;
657           G4ErrorMatrixIter mki  = mk;         << 619           G4ErrorMatrixIter mki = mk;
658           G4ErrorMatrixIter mji  = mj;         << 620           G4ErrorMatrixIter mji = mj;
659           for(G4int i = 1; i < j; i++)         << 621           for (G4int i=1;i<j;i++)
660           {                                       622           {
661             s11 += (*mik) * (*(mji++));           623             s11 += (*mik) * (*(mji++));
662             s12 += (*mijp) * (*(mki++));          624             s12 += (*mijp) * (*(mki++));
663             mik += n;                             625             mik += n;
664             mijp += n;                            626             mijp += n;
665           }                                       627           }
666         }                                         628         }
667         *(mjk++) = -s11 * (*mjj);                 629         *(mjk++) = -s11 * (*mjj);
668         *(mkjp)  = -(((*(mjj + 1))) * ((*(mkjp << 630         *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
669         mk += n;                                  631         mk += n;
670         mkjp += n;                                632         mkjp += n;
671       }                                           633       }
672     }                                             634     }
673     mj += n;                                      635     mj += n;
674     mjj += (n + 1);                            << 636     mjj += (n+1);
675   }                                               637   }
676   if(nxch % 2 == 1)                            << 638   if (nxch%2==1) det = -det;
677     det = -det;                                << 639   if (jfail !=jrange) det = 0.0;
678   if(jfail != jrange)                          << 
679     det = 0.0;                                 << 
680   ir[n] = nxch;                                   640   ir[n] = nxch;
681   return 0;                                       641   return 0;
682 }                                                 642 }
683                                                   643 
684 void G4ErrorMatrix::invert(G4int& ierr)        << 644 void G4ErrorMatrix::invert(G4int &ierr)
685 {                                                 645 {
686   if(ncol != nrow)                                646   if(ncol != nrow)
687   {                                            << 647      { error("G4ErrorMatrix::invert: G4ErrorMatrix is not NxN"); }
688     error("G4ErrorMatrix::invert: G4ErrorMatri << 
689   }                                            << 
690                                                   648 
691   static G4ThreadLocal G4int max_array = 20;   << 649   static G4int max_array = 20;
692   static G4ThreadLocal G4int* ir       = 0;    << 650   static G4int *ir = new G4int [max_array+1];
693   if(!ir)                                      << 
694     ir = new G4int[max_array + 1];             << 
695                                                   651 
696   if(ncol > max_array)                         << 652   if (ncol > max_array)
697   {                                               653   {
698     delete[] ir;                               << 654     delete [] ir;
699     max_array = nrow;                             655     max_array = nrow;
700     ir        = new G4int[max_array + 1];      << 656     ir = new G4int [max_array+1];
701   }                                               657   }
702   G4double t1, t2, t3;                            658   G4double t1, t2, t3;
703   G4double det, temp, ss;                      << 659   G4double det, temp, s;
704   G4int ifail;                                    660   G4int ifail;
705   switch(nrow)                                    661   switch(nrow)
706   {                                               662   {
707     case 3:                                    << 663   case 3:
708       G4double c11, c12, c13, c21, c22, c23, c << 664     G4double c11,c12,c13,c21,c22,c23,c31,c32,c33;
709       ifail = 0;                               << 665     ifail = 0;
710       c11   = (*(m.begin() + 4)) * (*(m.begin( << 666     c11 = (*(m.begin()+4)) * (*(m.begin()+8))
711             (*(m.begin() + 5)) * (*(m.begin()  << 667         - (*(m.begin()+5)) * (*(m.begin()+7));
712       c12 = (*(m.begin() + 5)) * (*(m.begin()  << 668     c12 = (*(m.begin()+5)) * (*(m.begin()+6))
713             (*(m.begin() + 3)) * (*(m.begin()  << 669         - (*(m.begin()+3)) * (*(m.begin()+8));
714       c13 = (*(m.begin() + 3)) * (*(m.begin()  << 670     c13 = (*(m.begin()+3)) * (*(m.begin()+7))
715             (*(m.begin() + 4)) * (*(m.begin()  << 671         - (*(m.begin()+4)) * (*(m.begin()+6));
716       c21 = (*(m.begin() + 7)) * (*(m.begin()  << 672     c21 = (*(m.begin()+7)) * (*(m.begin()+2))
717             (*(m.begin() + 8)) * (*(m.begin()  << 673         - (*(m.begin()+8)) * (*(m.begin()+1));
718       c22 = (*(m.begin() + 8)) * (*m.begin())  << 674     c22 = (*(m.begin()+8)) * (*m.begin())
719             (*(m.begin() + 6)) * (*(m.begin()  << 675         - (*(m.begin()+6)) * (*(m.begin()+2));
720       c23 = (*(m.begin() + 6)) * (*(m.begin()  << 676     c23 = (*(m.begin()+6)) * (*(m.begin()+1))
721             (*(m.begin() + 7)) * (*m.begin()); << 677         - (*(m.begin()+7)) * (*m.begin());
722       c31 = (*(m.begin() + 1)) * (*(m.begin()  << 678     c31 = (*(m.begin()+1)) * (*(m.begin()+5))
723             (*(m.begin() + 2)) * (*(m.begin()  << 679         - (*(m.begin()+2)) * (*(m.begin()+4));
724       c32 = (*(m.begin() + 2)) * (*(m.begin()  << 680     c32 = (*(m.begin()+2)) * (*(m.begin()+3))
725             (*m.begin()) * (*(m.begin() + 5)); << 681         - (*m.begin()) * (*(m.begin()+5));
726       c33 = (*m.begin()) * (*(m.begin() + 4))  << 682     c33 = (*m.begin()) * (*(m.begin()+4))
727             (*(m.begin() + 1)) * (*(m.begin()  << 683         - (*(m.begin()+1)) * (*(m.begin()+3));
728       t1 = std::fabs(*m.begin());              << 684     t1 = std::fabs(*m.begin());
729       t2 = std::fabs(*(m.begin() + 3));        << 685     t2 = std::fabs(*(m.begin()+3));
730       t3 = std::fabs(*(m.begin() + 6));        << 686     t3 = std::fabs(*(m.begin()+6));
731       if(t1 >= t2)                             << 687     if (t1 >= t2)
732       {                                        << 688     {
733         if(t3 >= t1)                           << 689       if (t3 >= t1)
734         {                                      << 
735           temp = *(m.begin() + 6);             << 
736           det  = c23 * c12 - c22 * c13;        << 
737         }                                      << 
738         else                                   << 
739         {                                      << 
740           temp = *(m.begin());                 << 
741           det  = c22 * c33 - c23 * c32;        << 
742         }                                      << 
743       }                                        << 
744       else if(t3 >= t2)                        << 
745       {                                           690       {
746         temp = *(m.begin() + 6);               << 691         temp = *(m.begin()+6);
747         det  = c23 * c12 - c22 * c13;          << 692         det = c23*c12-c22*c13;
748       }                                           693       }
749       else                                        694       else
750       {                                           695       {
751         temp = *(m.begin() + 3);               << 696         temp = *(m.begin());
752         det  = c13 * c32 - c12 * c33;          << 697         det = c22*c33-c23*c32;
753       }                                        << 
754       if(det == 0)                             << 
755       {                                        << 
756         ierr = 1;                              << 
757         return;                                << 
758       }                                           698       }
759       {                                        << 699     }
760         ss                   = temp / det;     << 700     else if (t3 >= t2)
761         G4ErrorMatrixIter mq = m.begin();      << 701     {
762         *(mq++)              = ss * c11;       << 702       temp = *(m.begin()+6);
763         *(mq++)              = ss * c21;       << 703       det = c23*c12-c22*c13;
764         *(mq++)              = ss * c31;       << 704     }
765         *(mq++)              = ss * c12;       << 705     else
766         *(mq++)              = ss * c22;       << 706     {
767         *(mq++)              = ss * c32;       << 707       temp = *(m.begin()+3);
768         *(mq++)              = ss * c13;       << 708       det = c13*c32-c12*c33;
769         *(mq++)              = ss * c23;       << 709     }
770         *(mq)                = ss * c33;       << 710     if (det==0)
771       }                                        << 711     {
772       break;                                   << 712       ierr = 1;
773     case 2:                                    << 
774       ifail = 0;                               << 
775       det   = (*m.begin()) * (*(m.begin() + 3) << 
776             (*(m.begin() + 1)) * (*(m.begin()  << 
777       if(det == 0)                             << 
778       {                                        << 
779         ierr = 1;                              << 
780         return;                                << 
781       }                                        << 
782       ss   = 1.0 / det;                        << 
783       temp = ss * (*(m.begin() + 3));          << 
784       *(m.begin() + 1) *= -ss;                 << 
785       *(m.begin() + 2) *= -ss;                 << 
786       *(m.begin() + 3) = ss * (*m.begin());    << 
787       *(m.begin())     = temp;                 << 
788       break;                                   << 
789     case 1:                                    << 
790       ifail = 0;                               << 
791       if((*(m.begin())) == 0)                  << 
792       {                                        << 
793         ierr = 1;                              << 
794         return;                                << 
795       }                                        << 
796       *(m.begin()) = 1.0 / (*(m.begin()));     << 
797       break;                                   << 
798     case 4:                                    << 
799       invertHaywood4(ierr);                    << 
800       return;                                     713       return;
801     case 5:                                    << 714     }
802       invertHaywood5(ierr);                    << 715     {
                                                   >> 716       G4double s = temp/det;
                                                   >> 717       G4ErrorMatrixIter mm = m.begin();
                                                   >> 718       *(mm++) = s*c11;
                                                   >> 719       *(mm++) = s*c21;
                                                   >> 720       *(mm++) = s*c31;
                                                   >> 721       *(mm++) = s*c12;
                                                   >> 722       *(mm++) = s*c22;
                                                   >> 723       *(mm++) = s*c32;
                                                   >> 724       *(mm++) = s*c13;
                                                   >> 725       *(mm++) = s*c23;
                                                   >> 726       *(mm) = s*c33;
                                                   >> 727     }
                                                   >> 728     break;
                                                   >> 729   case 2:
                                                   >> 730     ifail = 0;
                                                   >> 731     det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
                                                   >> 732     if (det==0)
                                                   >> 733     {
                                                   >> 734       ierr = 1;
803       return;                                     735       return;
804     case 6:                                    << 736     }
805       invertHaywood6(ierr);                    << 737     s = 1.0/det;
                                                   >> 738     temp = s*(*(m.begin()+3));
                                                   >> 739     *(m.begin()+1) *= -s;
                                                   >> 740     *(m.begin()+2) *= -s;
                                                   >> 741     *(m.begin()+3) = s*(*m.begin());
                                                   >> 742     *(m.begin()) = temp;
                                                   >> 743     break;
                                                   >> 744   case 1:
                                                   >> 745     ifail = 0;
                                                   >> 746     if ((*(m.begin()))==0)
                                                   >> 747     {
                                                   >> 748       ierr = 1;
806       return;                                     749       return;
807     default:                                   << 750     }
808       ifail = dfact_matrix(det, ir);           << 751     *(m.begin()) = 1.0/(*(m.begin()));
809       if(ifail)                                << 752     break;
810       {                                        << 753   case 4:
811         ierr = 1;                              << 754     invertHaywood4(ierr);
812         return;                                << 755     return;
813       }                                        << 756   case 5:
814       dfinv_matrix(ir);                        << 757     invertHaywood5(ierr);
815       break;                                   << 758     return;
                                                   >> 759   case 6:
                                                   >> 760     invertHaywood6(ierr);
                                                   >> 761     return;
                                                   >> 762   default:
                                                   >> 763     ifail = dfact_matrix(det, ir);
                                                   >> 764     if(ifail) {
                                                   >> 765       ierr = 1;
                                                   >> 766       return;
                                                   >> 767     }
                                                   >> 768     dfinv_matrix(ir);
                                                   >> 769     break;
816   }                                               770   }
817   ierr = 0;                                       771   ierr = 0;
818   return;                                         772   return;
819 }                                                 773 }
820                                                   774 
821 G4double G4ErrorMatrix::determinant() const       775 G4double G4ErrorMatrix::determinant() const
822 {                                                 776 {
823   static G4ThreadLocal G4int max_array = 20;   << 777   static G4int max_array = 20;
824   static G4ThreadLocal G4int* ir       = 0;    << 778   static G4int *ir = new G4int [max_array+1];
825   if(!ir)                                      << 
826     ir = new G4int[max_array + 1];             << 
827   if(ncol != nrow)                                779   if(ncol != nrow)
                                                   >> 780     { error("G4ErrorMatrix::determinant: G4ErrorMatrix is not NxN"); }
                                                   >> 781   if (ncol > max_array)
828   {                                               782   {
829     error("G4ErrorMatrix::determinant: G4Error << 783     delete [] ir;
830   }                                            << 
831   if(ncol > max_array)                         << 
832   {                                            << 
833     delete[] ir;                               << 
834     max_array = nrow;                             784     max_array = nrow;
835     ir        = new G4int[max_array + 1];      << 785     ir = new G4int [max_array+1];
836   }                                               786   }
837   G4double det;                                   787   G4double det;
838   G4ErrorMatrix mt(*this);                        788   G4ErrorMatrix mt(*this);
839   G4int i = mt.dfact_matrix(det, ir);             789   G4int i = mt.dfact_matrix(det, ir);
840   if(i == 0)                                   << 790   if(i==0) return det;
841     return det;                                << 
842   return 0;                                       791   return 0;
843 }                                                 792 }
844                                                   793 
845 G4double G4ErrorMatrix::trace() const             794 G4double G4ErrorMatrix::trace() const
846 {                                                 795 {
847   G4double t = 0.0;                               796   G4double t = 0.0;
848   for(G4ErrorMatrixConstIter d = m.begin(); d  << 797   for (G4ErrorMatrixConstIter d = m.begin(); d < m.end(); d += (ncol+1) )
849   {                                               798   {
850     t += *d;                                      799     t += *d;
851   }                                               800   }
852   return t;                                       801   return t;
853 }                                                 802 }
854                                                   803 
855 void G4ErrorMatrix::error(const char* msg)     << 804 void G4ErrorMatrix::error(const char *s)
856 {                                                 805 {
857   std::ostringstream message;                  << 806   G4cerr << s << G4endl;
858   message << msg;                              << 807   G4cerr << "---Exiting to System." << G4endl;
859   G4Exception("G4ErrorMatrix::error()", "GEANT << 808   abort();
860               message, "Exiting to System.");  << 
861 }                                                 809 }
862                                                   810 
                                                   >> 811 
863 // Aij are indices for a 6x6 matrix.              812 // Aij are indices for a 6x6 matrix.
864 // Mij are indices for a 5x5 matrix.              813 // Mij are indices for a 5x5 matrix.
865 // Fij are indices for a 4x4 matrix.              814 // Fij are indices for a 4x4 matrix.
866                                                   815 
867 #define A00 0                                     816 #define A00 0
868 #define A01 1                                     817 #define A01 1
869 #define A02 2                                     818 #define A02 2
870 #define A03 3                                     819 #define A03 3
871 #define A04 4                                     820 #define A04 4
872 #define A05 5                                     821 #define A05 5
873                                                   822 
874 #define A10 6                                     823 #define A10 6
875 #define A11 7                                     824 #define A11 7
876 #define A12 8                                     825 #define A12 8
877 #define A13 9                                     826 #define A13 9
878 #define A14 10                                    827 #define A14 10
879 #define A15 11                                    828 #define A15 11
880                                                   829 
881 #define A20 12                                    830 #define A20 12
882 #define A21 13                                    831 #define A21 13
883 #define A22 14                                    832 #define A22 14
884 #define A23 15                                    833 #define A23 15
885 #define A24 16                                    834 #define A24 16
886 #define A25 17                                    835 #define A25 17
887                                                   836 
888 #define A30 18                                    837 #define A30 18
889 #define A31 19                                    838 #define A31 19
890 #define A32 20                                    839 #define A32 20
891 #define A33 21                                    840 #define A33 21
892 #define A34 22                                    841 #define A34 22
893 #define A35 23                                    842 #define A35 23
894                                                   843 
895 #define A40 24                                    844 #define A40 24
896 #define A41 25                                    845 #define A41 25
897 #define A42 26                                    846 #define A42 26
898 #define A43 27                                    847 #define A43 27
899 #define A44 28                                    848 #define A44 28
900 #define A45 29                                    849 #define A45 29
901                                                   850 
902 #define A50 30                                    851 #define A50 30
903 #define A51 31                                    852 #define A51 31
904 #define A52 32                                    853 #define A52 32
905 #define A53 33                                    854 #define A53 33
906 #define A54 34                                    855 #define A54 34
907 #define A55 35                                    856 #define A55 35
908                                                   857 
909 #define M00 0                                     858 #define M00 0
910 #define M01 1                                     859 #define M01 1
911 #define M02 2                                     860 #define M02 2
912 #define M03 3                                     861 #define M03 3
913 #define M04 4                                     862 #define M04 4
914                                                   863 
915 #define M10 5                                     864 #define M10 5
916 #define M11 6                                     865 #define M11 6
917 #define M12 7                                     866 #define M12 7
918 #define M13 8                                     867 #define M13 8
919 #define M14 9                                     868 #define M14 9
920                                                   869 
921 #define M20 10                                    870 #define M20 10
922 #define M21 11                                    871 #define M21 11
923 #define M22 12                                    872 #define M22 12
924 #define M23 13                                    873 #define M23 13
925 #define M24 14                                    874 #define M24 14
926                                                   875 
927 #define M30 15                                    876 #define M30 15
928 #define M31 16                                    877 #define M31 16
929 #define M32 17                                    878 #define M32 17
930 #define M33 18                                    879 #define M33 18
931 #define M34 19                                    880 #define M34 19
932                                                   881 
933 #define M40 20                                    882 #define M40 20
934 #define M41 21                                    883 #define M41 21
935 #define M42 22                                    884 #define M42 22
936 #define M43 23                                    885 #define M43 23
937 #define M44 24                                    886 #define M44 24
938                                                   887 
939 #define F00 0                                     888 #define F00 0
940 #define F01 1                                     889 #define F01 1
941 #define F02 2                                     890 #define F02 2
942 #define F03 3                                     891 #define F03 3
943                                                   892 
944 #define F10 4                                     893 #define F10 4
945 #define F11 5                                     894 #define F11 5
946 #define F12 6                                     895 #define F12 6
947 #define F13 7                                     896 #define F13 7
948                                                   897 
949 #define F20 8                                     898 #define F20 8
950 #define F21 9                                     899 #define F21 9
951 #define F22 10                                    900 #define F22 10
952 #define F23 11                                    901 #define F23 11
953                                                   902 
954 #define F30 12                                    903 #define F30 12
955 #define F31 13                                    904 #define F31 13
956 #define F32 14                                    905 #define F32 14
957 #define F33 15                                    906 #define F33 15
958                                                   907 
959 void G4ErrorMatrix::invertHaywood4(G4int& ifai << 908 
                                                   >> 909 void G4ErrorMatrix::invertHaywood4  (G4int & ifail)
960 {                                                 910 {
961   ifail = 0;                                      911   ifail = 0;
962                                                   912 
963   // Find all NECESSARY 2x2 dets:  (18 of them    913   // Find all NECESSARY 2x2 dets:  (18 of them)
964                                                   914 
965   G4double Det2_12_01 = m[F10] * m[F21] - m[F1 << 915   G4double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
966   G4double Det2_12_02 = m[F10] * m[F22] - m[F1 << 916   G4double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
967   G4double Det2_12_03 = m[F10] * m[F23] - m[F1 << 917   G4double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20];
968   G4double Det2_12_13 = m[F11] * m[F23] - m[F1 << 918   G4double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21];
969   G4double Det2_12_23 = m[F12] * m[F23] - m[F1 << 919   G4double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22];
970   G4double Det2_12_12 = m[F11] * m[F22] - m[F1 << 920   G4double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
971   G4double Det2_13_01 = m[F10] * m[F31] - m[F1 << 921   G4double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
972   G4double Det2_13_02 = m[F10] * m[F32] - m[F1 << 922   G4double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
973   G4double Det2_13_03 = m[F10] * m[F33] - m[F1 << 923   G4double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
974   G4double Det2_13_12 = m[F11] * m[F32] - m[F1 << 924   G4double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
975   G4double Det2_13_13 = m[F11] * m[F33] - m[F1 << 925   G4double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
976   G4double Det2_13_23 = m[F12] * m[F33] - m[F1 << 926   G4double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32];
977   G4double Det2_23_01 = m[F20] * m[F31] - m[F2 << 927   G4double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
978   G4double Det2_23_02 = m[F20] * m[F32] - m[F2 << 928   G4double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
979   G4double Det2_23_03 = m[F20] * m[F33] - m[F2 << 929   G4double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
980   G4double Det2_23_12 = m[F21] * m[F32] - m[F2 << 930   G4double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
981   G4double Det2_23_13 = m[F21] * m[F33] - m[F2 << 931   G4double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
982   G4double Det2_23_23 = m[F22] * m[F33] - m[F2 << 932   G4double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
983                                                   933 
984   // Find all NECESSARY 3x3 dets:   (16 of the    934   // Find all NECESSARY 3x3 dets:   (16 of them)
985                                                   935 
986   G4double Det3_012_012 =                      << 936   G4double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02 
987     m[F00] * Det2_12_12 - m[F01] * Det2_12_02  << 937                                 + m[F02]*Det2_12_01;
988   G4double Det3_012_013 =                      << 938   G4double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03 
989     m[F00] * Det2_12_13 - m[F01] * Det2_12_03  << 939                                 + m[F03]*Det2_12_01;
990   G4double Det3_012_023 =                      << 940   G4double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03 
991     m[F00] * Det2_12_23 - m[F02] * Det2_12_03  << 941                                 + m[F03]*Det2_12_02;
992   G4double Det3_012_123 =                      << 942   G4double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13 
993     m[F01] * Det2_12_23 - m[F02] * Det2_12_13  << 943                                 + m[F03]*Det2_12_12;
994   G4double Det3_013_012 =                      << 944   G4double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02 
995     m[F00] * Det2_13_12 - m[F01] * Det2_13_02  << 945                                 + m[F02]*Det2_13_01;
996   G4double Det3_013_013 =                      << 946   G4double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
997     m[F00] * Det2_13_13 - m[F01] * Det2_13_03  << 947                                 + m[F03]*Det2_13_01;
998   G4double Det3_013_023 =                      << 948   G4double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
999     m[F00] * Det2_13_23 - m[F02] * Det2_13_03  << 949                                 + m[F03]*Det2_13_02;
1000   G4double Det3_013_123 =                     << 950   G4double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
1001     m[F01] * Det2_13_23 - m[F02] * Det2_13_13 << 951                                 + m[F03]*Det2_13_12;
1002   G4double Det3_023_012 =                     << 952   G4double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02 
1003     m[F00] * Det2_23_12 - m[F01] * Det2_23_02 << 953                                 + m[F02]*Det2_23_01;
1004   G4double Det3_023_013 =                     << 954   G4double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
1005     m[F00] * Det2_23_13 - m[F01] * Det2_23_03 << 955                                 + m[F03]*Det2_23_01;
1006   G4double Det3_023_023 =                     << 956   G4double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
1007     m[F00] * Det2_23_23 - m[F02] * Det2_23_03 << 957                                 + m[F03]*Det2_23_02;
1008   G4double Det3_023_123 =                     << 958   G4double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
1009     m[F01] * Det2_23_23 - m[F02] * Det2_23_13 << 959                                 + m[F03]*Det2_23_12;
1010   G4double Det3_123_012 =                     << 960   G4double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02 
1011     m[F10] * Det2_23_12 - m[F11] * Det2_23_02 << 961                                 + m[F12]*Det2_23_01;
1012   G4double Det3_123_013 =                     << 962   G4double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03 
1013     m[F10] * Det2_23_13 - m[F11] * Det2_23_03 << 963                                 + m[F13]*Det2_23_01;
1014   G4double Det3_123_023 =                     << 964   G4double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03 
1015     m[F10] * Det2_23_23 - m[F12] * Det2_23_03 << 965                                 + m[F13]*Det2_23_02;
1016   G4double Det3_123_123 =                     << 966   G4double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13 
1017     m[F11] * Det2_23_23 - m[F12] * Det2_23_13 << 967                                 + m[F13]*Det2_23_12;
1018                                                  968 
1019   // Find the 4x4 det:                           969   // Find the 4x4 det:
1020                                                  970 
1021   G4double det = m[F00] * Det3_123_123 - m[F0 << 971   G4double det =    m[F00]*Det3_123_123 
1022                  m[F02] * Det3_123_013 - m[F0 << 972                 - m[F01]*Det3_123_023 
                                                   >> 973                 + m[F02]*Det3_123_013 
                                                   >> 974                 - m[F03]*Det3_123_012;
1023                                                  975 
1024   if(det == 0)                                << 976   if ( det == 0 )
1025   {                                           << 977   {  
1026     ifail = 1;                                   978     ifail = 1;
1027     return;                                      979     return;
1028   }                                           << 980   } 
1029                                                  981 
1030   G4double oneOverDet = 1.0 / det;            << 982   G4double oneOverDet = 1.0/det;
1031   G4double mn1OverDet = -oneOverDet;          << 983   G4double mn1OverDet = - oneOverDet;
1032                                                  984 
1033   m[F00] = Det3_123_123 * oneOverDet;         << 985   m[F00] =  Det3_123_123 * oneOverDet;
1034   m[F01] = Det3_023_123 * mn1OverDet;         << 986   m[F01] =  Det3_023_123 * mn1OverDet;
1035   m[F02] = Det3_013_123 * oneOverDet;         << 987   m[F02] =  Det3_013_123 * oneOverDet;
1036   m[F03] = Det3_012_123 * mn1OverDet;         << 988   m[F03] =  Det3_012_123 * mn1OverDet;
1037                                               << 989 
1038   m[F10] = Det3_123_023 * mn1OverDet;         << 990   m[F10] =  Det3_123_023 * mn1OverDet;
1039   m[F11] = Det3_023_023 * oneOverDet;         << 991   m[F11] =  Det3_023_023 * oneOverDet;
1040   m[F12] = Det3_013_023 * mn1OverDet;         << 992   m[F12] =  Det3_013_023 * mn1OverDet;
1041   m[F13] = Det3_012_023 * oneOverDet;         << 993   m[F13] =  Det3_012_023 * oneOverDet;
1042                                               << 994 
1043   m[F20] = Det3_123_013 * oneOverDet;         << 995   m[F20] =  Det3_123_013 * oneOverDet;
1044   m[F21] = Det3_023_013 * mn1OverDet;         << 996   m[F21] =  Det3_023_013 * mn1OverDet;
1045   m[F22] = Det3_013_013 * oneOverDet;         << 997   m[F22] =  Det3_013_013 * oneOverDet;
1046   m[F23] = Det3_012_013 * mn1OverDet;         << 998   m[F23] =  Det3_012_013 * mn1OverDet;
1047                                               << 999 
1048   m[F30] = Det3_123_012 * mn1OverDet;         << 1000   m[F30] =  Det3_123_012 * mn1OverDet;
1049   m[F31] = Det3_023_012 * oneOverDet;         << 1001   m[F31] =  Det3_023_012 * oneOverDet;
1050   m[F32] = Det3_013_012 * mn1OverDet;         << 1002   m[F32] =  Det3_013_012 * mn1OverDet;
1051   m[F33] = Det3_012_012 * oneOverDet;         << 1003   m[F33] =  Det3_012_012 * oneOverDet;
1052                                                  1004 
1053   return;                                        1005   return;
1054 }                                                1006 }
1055                                                  1007 
1056 void G4ErrorMatrix::invertHaywood5(G4int& ifa << 1008 void G4ErrorMatrix::invertHaywood5  (G4int & ifail)
1057 {                                                1009 {
1058   ifail = 0;                                     1010   ifail = 0;
1059                                                  1011 
1060   // Find all NECESSARY 2x2 dets:  (30 of the    1012   // Find all NECESSARY 2x2 dets:  (30 of them)
1061                                                  1013 
1062   G4double Det2_23_01 = m[M20] * m[M31] - m[M << 1014   G4double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
1063   G4double Det2_23_02 = m[M20] * m[M32] - m[M << 1015   G4double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
1064   G4double Det2_23_03 = m[M20] * m[M33] - m[M << 1016   G4double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
1065   G4double Det2_23_04 = m[M20] * m[M34] - m[M << 1017   G4double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
1066   G4double Det2_23_12 = m[M21] * m[M32] - m[M << 1018   G4double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31];
1067   G4double Det2_23_13 = m[M21] * m[M33] - m[M << 1019   G4double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
1068   G4double Det2_23_14 = m[M21] * m[M34] - m[M << 1020   G4double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31];
1069   G4double Det2_23_23 = m[M22] * m[M33] - m[M << 1021   G4double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
1070   G4double Det2_23_24 = m[M22] * m[M34] - m[M << 1022   G4double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32];
1071   G4double Det2_23_34 = m[M23] * m[M34] - m[M << 1023   G4double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33];
1072   G4double Det2_24_01 = m[M20] * m[M41] - m[M << 1024   G4double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
1073   G4double Det2_24_02 = m[M20] * m[M42] - m[M << 1025   G4double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
1074   G4double Det2_24_03 = m[M20] * m[M43] - m[M << 1026   G4double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
1075   G4double Det2_24_04 = m[M20] * m[M44] - m[M << 1027   G4double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
1076   G4double Det2_24_12 = m[M21] * m[M42] - m[M << 1028   G4double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
1077   G4double Det2_24_13 = m[M21] * m[M43] - m[M << 1029   G4double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
1078   G4double Det2_24_14 = m[M21] * m[M44] - m[M << 1030   G4double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
1079   G4double Det2_24_23 = m[M22] * m[M43] - m[M << 1031   G4double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
1080   G4double Det2_24_24 = m[M22] * m[M44] - m[M << 1032   G4double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
1081   G4double Det2_24_34 = m[M23] * m[M44] - m[M << 1033   G4double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43];
1082   G4double Det2_34_01 = m[M30] * m[M41] - m[M << 1034   G4double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
1083   G4double Det2_34_02 = m[M30] * m[M42] - m[M << 1035   G4double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
1084   G4double Det2_34_03 = m[M30] * m[M43] - m[M << 1036   G4double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
1085   G4double Det2_34_04 = m[M30] * m[M44] - m[M << 1037   G4double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
1086   G4double Det2_34_12 = m[M31] * m[M42] - m[M << 1038   G4double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
1087   G4double Det2_34_13 = m[M31] * m[M43] - m[M << 1039   G4double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
1088   G4double Det2_34_14 = m[M31] * m[M44] - m[M << 1040   G4double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
1089   G4double Det2_34_23 = m[M32] * m[M43] - m[M << 1041   G4double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
1090   G4double Det2_34_24 = m[M32] * m[M44] - m[M << 1042   G4double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
1091   G4double Det2_34_34 = m[M33] * m[M44] - m[M << 1043   G4double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
1092                                                  1044 
1093   // Find all NECESSARY 3x3 dets:   (40 of th    1045   // Find all NECESSARY 3x3 dets:   (40 of them)
1094                                                  1046 
1095   G4double Det3_123_012 =                     << 1047   G4double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02 
1096     m[M10] * Det2_23_12 - m[M11] * Det2_23_02 << 1048                                 + m[M12]*Det2_23_01;
1097   G4double Det3_123_013 =                     << 1049   G4double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03 
1098     m[M10] * Det2_23_13 - m[M11] * Det2_23_03 << 1050                                 + m[M13]*Det2_23_01;
1099   G4double Det3_123_014 =                     << 1051   G4double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04 
1100     m[M10] * Det2_23_14 - m[M11] * Det2_23_04 << 1052                                 + m[M14]*Det2_23_01;
1101   G4double Det3_123_023 =                     << 1053   G4double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03 
1102     m[M10] * Det2_23_23 - m[M12] * Det2_23_03 << 1054                                 + m[M13]*Det2_23_02;
1103   G4double Det3_123_024 =                     << 1055   G4double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04 
1104     m[M10] * Det2_23_24 - m[M12] * Det2_23_04 << 1056                                 + m[M14]*Det2_23_02;
1105   G4double Det3_123_034 =                     << 1057   G4double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04 
1106     m[M10] * Det2_23_34 - m[M13] * Det2_23_04 << 1058                                 + m[M14]*Det2_23_03;
1107   G4double Det3_123_123 =                     << 1059   G4double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13 
1108     m[M11] * Det2_23_23 - m[M12] * Det2_23_13 << 1060                                 + m[M13]*Det2_23_12;
1109   G4double Det3_123_124 =                     << 1061   G4double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14 
1110     m[M11] * Det2_23_24 - m[M12] * Det2_23_14 << 1062                                 + m[M14]*Det2_23_12;
1111   G4double Det3_123_134 =                     << 1063   G4double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14 
1112     m[M11] * Det2_23_34 - m[M13] * Det2_23_14 << 1064                                 + m[M14]*Det2_23_13;
1113   G4double Det3_123_234 =                     << 1065   G4double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24 
1114     m[M12] * Det2_23_34 - m[M13] * Det2_23_24 << 1066                                 + m[M14]*Det2_23_23;
1115   G4double Det3_124_012 =                     << 1067   G4double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02 
1116     m[M10] * Det2_24_12 - m[M11] * Det2_24_02 << 1068                                 + m[M12]*Det2_24_01;
1117   G4double Det3_124_013 =                     << 1069   G4double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03 
1118     m[M10] * Det2_24_13 - m[M11] * Det2_24_03 << 1070                                 + m[M13]*Det2_24_01;
1119   G4double Det3_124_014 =                     << 1071   G4double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04 
1120     m[M10] * Det2_24_14 - m[M11] * Det2_24_04 << 1072                                 + m[M14]*Det2_24_01;
1121   G4double Det3_124_023 =                     << 1073   G4double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03 
1122     m[M10] * Det2_24_23 - m[M12] * Det2_24_03 << 1074                                 + m[M13]*Det2_24_02;
1123   G4double Det3_124_024 =                     << 1075   G4double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04 
1124     m[M10] * Det2_24_24 - m[M12] * Det2_24_04 << 1076                                 + m[M14]*Det2_24_02;
1125   G4double Det3_124_034 =                     << 1077   G4double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04 
1126     m[M10] * Det2_24_34 - m[M13] * Det2_24_04 << 1078                                 + m[M14]*Det2_24_03;
1127   G4double Det3_124_123 =                     << 1079   G4double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13 
1128     m[M11] * Det2_24_23 - m[M12] * Det2_24_13 << 1080                                 + m[M13]*Det2_24_12;
1129   G4double Det3_124_124 =                     << 1081   G4double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14 
1130     m[M11] * Det2_24_24 - m[M12] * Det2_24_14 << 1082                                 + m[M14]*Det2_24_12;
1131   G4double Det3_124_134 =                     << 1083   G4double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14 
1132     m[M11] * Det2_24_34 - m[M13] * Det2_24_14 << 1084                                 + m[M14]*Det2_24_13;
1133   G4double Det3_124_234 =                     << 1085   G4double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24 
1134     m[M12] * Det2_24_34 - m[M13] * Det2_24_24 << 1086                                 + m[M14]*Det2_24_23;
1135   G4double Det3_134_012 =                     << 1087   G4double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02 
1136     m[M10] * Det2_34_12 - m[M11] * Det2_34_02 << 1088                                 + m[M12]*Det2_34_01;
1137   G4double Det3_134_013 =                     << 1089   G4double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03 
1138     m[M10] * Det2_34_13 - m[M11] * Det2_34_03 << 1090                                 + m[M13]*Det2_34_01;
1139   G4double Det3_134_014 =                     << 1091   G4double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04 
1140     m[M10] * Det2_34_14 - m[M11] * Det2_34_04 << 1092                                 + m[M14]*Det2_34_01;
1141   G4double Det3_134_023 =                     << 1093   G4double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03 
1142     m[M10] * Det2_34_23 - m[M12] * Det2_34_03 << 1094                                 + m[M13]*Det2_34_02;
1143   G4double Det3_134_024 =                     << 1095   G4double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04 
1144     m[M10] * Det2_34_24 - m[M12] * Det2_34_04 << 1096                                 + m[M14]*Det2_34_02;
1145   G4double Det3_134_034 =                     << 1097   G4double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04 
1146     m[M10] * Det2_34_34 - m[M13] * Det2_34_04 << 1098                                 + m[M14]*Det2_34_03;
1147   G4double Det3_134_123 =                     << 1099   G4double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13 
1148     m[M11] * Det2_34_23 - m[M12] * Det2_34_13 << 1100                                 + m[M13]*Det2_34_12;
1149   G4double Det3_134_124 =                     << 1101   G4double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14 
1150     m[M11] * Det2_34_24 - m[M12] * Det2_34_14 << 1102                                 + m[M14]*Det2_34_12;
1151   G4double Det3_134_134 =                     << 1103   G4double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14 
1152     m[M11] * Det2_34_34 - m[M13] * Det2_34_14 << 1104                                 + m[M14]*Det2_34_13;
1153   G4double Det3_134_234 =                     << 1105   G4double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24 
1154     m[M12] * Det2_34_34 - m[M13] * Det2_34_24 << 1106                                 + m[M14]*Det2_34_23;
1155   G4double Det3_234_012 =                     << 1107   G4double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02 
1156     m[M20] * Det2_34_12 - m[M21] * Det2_34_02 << 1108                                 + m[M22]*Det2_34_01;
1157   G4double Det3_234_013 =                     << 1109   G4double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03 
1158     m[M20] * Det2_34_13 - m[M21] * Det2_34_03 << 1110                                 + m[M23]*Det2_34_01;
1159   G4double Det3_234_014 =                     << 1111   G4double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04 
1160     m[M20] * Det2_34_14 - m[M21] * Det2_34_04 << 1112                                 + m[M24]*Det2_34_01;
1161   G4double Det3_234_023 =                     << 1113   G4double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03 
1162     m[M20] * Det2_34_23 - m[M22] * Det2_34_03 << 1114                                 + m[M23]*Det2_34_02;
1163   G4double Det3_234_024 =                     << 1115   G4double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04 
1164     m[M20] * Det2_34_24 - m[M22] * Det2_34_04 << 1116                                 + m[M24]*Det2_34_02;
1165   G4double Det3_234_034 =                     << 1117   G4double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04 
1166     m[M20] * Det2_34_34 - m[M23] * Det2_34_04 << 1118                                 + m[M24]*Det2_34_03;
1167   G4double Det3_234_123 =                     << 1119   G4double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13 
1168     m[M21] * Det2_34_23 - m[M22] * Det2_34_13 << 1120                                 + m[M23]*Det2_34_12;
1169   G4double Det3_234_124 =                     << 1121   G4double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14 
1170     m[M21] * Det2_34_24 - m[M22] * Det2_34_14 << 1122                                 + m[M24]*Det2_34_12;
1171   G4double Det3_234_134 =                     << 1123   G4double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14 
1172     m[M21] * Det2_34_34 - m[M23] * Det2_34_14 << 1124                                 + m[M24]*Det2_34_13;
1173   G4double Det3_234_234 =                     << 1125   G4double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24 
1174     m[M22] * Det2_34_34 - m[M23] * Det2_34_24 << 1126                                 + m[M24]*Det2_34_23;
1175                                                  1127 
1176   // Find all NECESSARY 4x4 dets:   (25 of th    1128   // Find all NECESSARY 4x4 dets:   (25 of them)
1177                                                  1129 
1178   G4double Det4_0123_0123 = m[M00] * Det3_123 << 1130   G4double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023 
1179                             m[M02] * Det3_123 << 1131                                 + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
1180   G4double Det4_0123_0124 = m[M00] * Det3_123 << 1132   G4double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024 
1181                             m[M02] * Det3_123 << 1133                                 + m[M02]*Det3_123_014 - m[M04]*Det3_123_012;
1182   G4double Det4_0123_0134 = m[M00] * Det3_123 << 1134   G4double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034 
1183                             m[M03] * Det3_123 << 1135                                 + m[M03]*Det3_123_014 - m[M04]*Det3_123_013;
1184   G4double Det4_0123_0234 = m[M00] * Det3_123 << 1136   G4double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034 
1185                             m[M03] * Det3_123 << 1137                                 + m[M03]*Det3_123_024 - m[M04]*Det3_123_023;
1186   G4double Det4_0123_1234 = m[M01] * Det3_123 << 1138   G4double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134 
1187                             m[M03] * Det3_123 << 1139                                 + m[M03]*Det3_123_124 - m[M04]*Det3_123_123;
1188   G4double Det4_0124_0123 = m[M00] * Det3_124 << 1140   G4double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023 
1189                             m[M02] * Det3_124 << 1141                                 + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
1190   G4double Det4_0124_0124 = m[M00] * Det3_124 << 1142   G4double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024 
1191                             m[M02] * Det3_124 << 1143                                 + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
1192   G4double Det4_0124_0134 = m[M00] * Det3_124 << 1144   G4double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034 
1193                             m[M03] * Det3_124 << 1145                                 + m[M03]*Det3_124_014 - m[M04]*Det3_124_013;
1194   G4double Det4_0124_0234 = m[M00] * Det3_124 << 1146   G4double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034 
1195                             m[M03] * Det3_124 << 1147                                 + m[M03]*Det3_124_024 - m[M04]*Det3_124_023;
1196   G4double Det4_0124_1234 = m[M01] * Det3_124 << 1148   G4double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134 
1197                             m[M03] * Det3_124 << 1149                                 + m[M03]*Det3_124_124 - m[M04]*Det3_124_123;
1198   G4double Det4_0134_0123 = m[M00] * Det3_134 << 1150   G4double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023 
1199                             m[M02] * Det3_134 << 1151                                 + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
1200   G4double Det4_0134_0124 = m[M00] * Det3_134 << 1152   G4double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024 
1201                             m[M02] * Det3_134 << 1153                                 + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
1202   G4double Det4_0134_0134 = m[M00] * Det3_134 << 1154   G4double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034 
1203                             m[M03] * Det3_134 << 1155                                 + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
1204   G4double Det4_0134_0234 = m[M00] * Det3_134 << 1156   G4double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034 
1205                             m[M03] * Det3_134 << 1157                                 + m[M03]*Det3_134_024 - m[M04]*Det3_134_023;
1206   G4double Det4_0134_1234 = m[M01] * Det3_134 << 1158   G4double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134 
1207                             m[M03] * Det3_134 << 1159                                 + m[M03]*Det3_134_124 - m[M04]*Det3_134_123;
1208   G4double Det4_0234_0123 = m[M00] * Det3_234 << 1160   G4double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023 
1209                             m[M02] * Det3_234 << 1161                                 + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
1210   G4double Det4_0234_0124 = m[M00] * Det3_234 << 1162   G4double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024 
1211                             m[M02] * Det3_234 << 1163                                 + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
1212   G4double Det4_0234_0134 = m[M00] * Det3_234 << 1164   G4double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034 
1213                             m[M03] * Det3_234 << 1165                                 + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
1214   G4double Det4_0234_0234 = m[M00] * Det3_234 << 1166   G4double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034 
1215                             m[M03] * Det3_234 << 1167                                 + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
1216   G4double Det4_0234_1234 = m[M01] * Det3_234 << 1168   G4double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134 
1217                             m[M03] * Det3_234 << 1169                                 + m[M03]*Det3_234_124 - m[M04]*Det3_234_123;
1218   G4double Det4_1234_0123 = m[M10] * Det3_234 << 1170   G4double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023 
1219                             m[M12] * Det3_234 << 1171                                 + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
1220   G4double Det4_1234_0124 = m[M10] * Det3_234 << 1172   G4double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024 
1221                             m[M12] * Det3_234 << 1173                                 + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
1222   G4double Det4_1234_0134 = m[M10] * Det3_234 << 1174   G4double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034 
1223                             m[M13] * Det3_234 << 1175                                 + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
1224   G4double Det4_1234_0234 = m[M10] * Det3_234 << 1176   G4double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034 
1225                             m[M13] * Det3_234 << 1177                                 + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
1226   G4double Det4_1234_1234 = m[M11] * Det3_234 << 1178   G4double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134 
1227                             m[M13] * Det3_234 << 1179                                 + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
1228                                                  1180 
1229   // Find the 5x5 det:                           1181   // Find the 5x5 det:
1230                                                  1182 
1231   G4double det = m[M00] * Det4_1234_1234 - m[ << 1183   G4double det =  m[M00]*Det4_1234_1234 
1232                  m[M02] * Det4_1234_0134 - m[ << 1184                 - m[M01]*Det4_1234_0234 
1233                  m[M04] * Det4_1234_0123;     << 1185                 + m[M02]*Det4_1234_0134 
                                                   >> 1186                 - m[M03]*Det4_1234_0124 
                                                   >> 1187                 + m[M04]*Det4_1234_0123;
1234                                                  1188 
1235   if(det == 0)                                << 1189   if ( det == 0 )
1236   {                                           << 1190   {  
1237     ifail = 1;                                   1191     ifail = 1;
1238     return;                                      1192     return;
1239   }                                           << 1193   } 
1240                                                  1194 
1241   G4double oneOverDet = 1.0 / det;            << 1195   G4double oneOverDet = 1.0/det;
1242   G4double mn1OverDet = -oneOverDet;          << 1196   G4double mn1OverDet = - oneOverDet;
1243                                                  1197 
1244   m[M00] = Det4_1234_1234 * oneOverDet;       << 1198   m[M00] =  Det4_1234_1234 * oneOverDet;
1245   m[M01] = Det4_0234_1234 * mn1OverDet;       << 1199   m[M01] =  Det4_0234_1234 * mn1OverDet;
1246   m[M02] = Det4_0134_1234 * oneOverDet;       << 1200   m[M02] =  Det4_0134_1234 * oneOverDet;
1247   m[M03] = Det4_0124_1234 * mn1OverDet;       << 1201   m[M03] =  Det4_0124_1234 * mn1OverDet;
1248   m[M04] = Det4_0123_1234 * oneOverDet;       << 1202   m[M04] =  Det4_0123_1234 * oneOverDet;
1249                                               << 1203 
1250   m[M10] = Det4_1234_0234 * mn1OverDet;       << 1204   m[M10] =  Det4_1234_0234 * mn1OverDet;
1251   m[M11] = Det4_0234_0234 * oneOverDet;       << 1205   m[M11] =  Det4_0234_0234 * oneOverDet;
1252   m[M12] = Det4_0134_0234 * mn1OverDet;       << 1206   m[M12] =  Det4_0134_0234 * mn1OverDet;
1253   m[M13] = Det4_0124_0234 * oneOverDet;       << 1207   m[M13] =  Det4_0124_0234 * oneOverDet;
1254   m[M14] = Det4_0123_0234 * mn1OverDet;       << 1208   m[M14] =  Det4_0123_0234 * mn1OverDet;
1255                                               << 1209 
1256   m[M20] = Det4_1234_0134 * oneOverDet;       << 1210   m[M20] =  Det4_1234_0134 * oneOverDet;
1257   m[M21] = Det4_0234_0134 * mn1OverDet;       << 1211   m[M21] =  Det4_0234_0134 * mn1OverDet;
1258   m[M22] = Det4_0134_0134 * oneOverDet;       << 1212   m[M22] =  Det4_0134_0134 * oneOverDet;
1259   m[M23] = Det4_0124_0134 * mn1OverDet;       << 1213   m[M23] =  Det4_0124_0134 * mn1OverDet;
1260   m[M24] = Det4_0123_0134 * oneOverDet;       << 1214   m[M24] =  Det4_0123_0134 * oneOverDet;
1261                                               << 1215 
1262   m[M30] = Det4_1234_0124 * mn1OverDet;       << 1216   m[M30] =  Det4_1234_0124 * mn1OverDet;
1263   m[M31] = Det4_0234_0124 * oneOverDet;       << 1217   m[M31] =  Det4_0234_0124 * oneOverDet;
1264   m[M32] = Det4_0134_0124 * mn1OverDet;       << 1218   m[M32] =  Det4_0134_0124 * mn1OverDet;
1265   m[M33] = Det4_0124_0124 * oneOverDet;       << 1219   m[M33] =  Det4_0124_0124 * oneOverDet;
1266   m[M34] = Det4_0123_0124 * mn1OverDet;       << 1220   m[M34] =  Det4_0123_0124 * mn1OverDet;
1267                                               << 1221 
1268   m[M40] = Det4_1234_0123 * oneOverDet;       << 1222   m[M40] =  Det4_1234_0123 * oneOverDet;
1269   m[M41] = Det4_0234_0123 * mn1OverDet;       << 1223   m[M41] =  Det4_0234_0123 * mn1OverDet;
1270   m[M42] = Det4_0134_0123 * oneOverDet;       << 1224   m[M42] =  Det4_0134_0123 * oneOverDet;
1271   m[M43] = Det4_0124_0123 * mn1OverDet;       << 1225   m[M43] =  Det4_0124_0123 * mn1OverDet;
1272   m[M44] = Det4_0123_0123 * oneOverDet;       << 1226   m[M44] =  Det4_0123_0123 * oneOverDet;
1273                                                  1227 
1274   return;                                        1228   return;
1275 }                                                1229 }
1276                                                  1230 
1277 void G4ErrorMatrix::invertHaywood6(G4int& ifa << 1231 void G4ErrorMatrix::invertHaywood6  (G4int & ifail)
1278 {                                                1232 {
1279   ifail = 0;                                     1233   ifail = 0;
1280                                                  1234 
1281   // Find all NECESSARY 2x2 dets:  (45 of the    1235   // Find all NECESSARY 2x2 dets:  (45 of them)
1282                                                  1236 
1283   G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1237   G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1284   G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1238   G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1285   G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1239   G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1286   G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1240   G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1287   G4double Det2_34_05 = m[A30] * m[A45] - m[A << 1241   G4double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40];
1288   G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1242   G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1289   G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1243   G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1290   G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1244   G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1291   G4double Det2_34_15 = m[A31] * m[A45] - m[A << 1245   G4double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41];
1292   G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1246   G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1293   G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1247   G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1294   G4double Det2_34_25 = m[A32] * m[A45] - m[A << 1248   G4double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42];
1295   G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1249   G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1296   G4double Det2_34_35 = m[A33] * m[A45] - m[A << 1250   G4double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43];
1297   G4double Det2_34_45 = m[A34] * m[A45] - m[A << 1251   G4double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44];
1298   G4double Det2_35_01 = m[A30] * m[A51] - m[A << 1252   G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1299   G4double Det2_35_02 = m[A30] * m[A52] - m[A << 1253   G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1300   G4double Det2_35_03 = m[A30] * m[A53] - m[A << 1254   G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1301   G4double Det2_35_04 = m[A30] * m[A54] - m[A << 1255   G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1302   G4double Det2_35_05 = m[A30] * m[A55] - m[A << 1256   G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1303   G4double Det2_35_12 = m[A31] * m[A52] - m[A << 1257   G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1304   G4double Det2_35_13 = m[A31] * m[A53] - m[A << 1258   G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1305   G4double Det2_35_14 = m[A31] * m[A54] - m[A << 1259   G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1306   G4double Det2_35_15 = m[A31] * m[A55] - m[A << 1260   G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1307   G4double Det2_35_23 = m[A32] * m[A53] - m[A << 1261   G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1308   G4double Det2_35_24 = m[A32] * m[A54] - m[A << 1262   G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1309   G4double Det2_35_25 = m[A32] * m[A55] - m[A << 1263   G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1310   G4double Det2_35_34 = m[A33] * m[A54] - m[A << 1264   G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1311   G4double Det2_35_35 = m[A33] * m[A55] - m[A << 1265   G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1312   G4double Det2_35_45 = m[A34] * m[A55] - m[A << 1266   G4double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54];
1313   G4double Det2_45_01 = m[A40] * m[A51] - m[A << 1267   G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1314   G4double Det2_45_02 = m[A40] * m[A52] - m[A << 1268   G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1315   G4double Det2_45_03 = m[A40] * m[A53] - m[A << 1269   G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1316   G4double Det2_45_04 = m[A40] * m[A54] - m[A << 1270   G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1317   G4double Det2_45_05 = m[A40] * m[A55] - m[A << 1271   G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1318   G4double Det2_45_12 = m[A41] * m[A52] - m[A << 1272   G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1319   G4double Det2_45_13 = m[A41] * m[A53] - m[A << 1273   G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1320   G4double Det2_45_14 = m[A41] * m[A54] - m[A << 1274   G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1321   G4double Det2_45_15 = m[A41] * m[A55] - m[A << 1275   G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1322   G4double Det2_45_23 = m[A42] * m[A53] - m[A << 1276   G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1323   G4double Det2_45_24 = m[A42] * m[A54] - m[A << 1277   G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1324   G4double Det2_45_25 = m[A42] * m[A55] - m[A << 1278   G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1325   G4double Det2_45_34 = m[A43] * m[A54] - m[A << 1279   G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1326   G4double Det2_45_35 = m[A43] * m[A55] - m[A << 1280   G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1327   G4double Det2_45_45 = m[A44] * m[A55] - m[A << 1281   G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1328                                                  1282 
1329   // Find all NECESSARY 3x3 dets:  (80 of the    1283   // Find all NECESSARY 3x3 dets:  (80 of them)
1330                                                  1284 
1331   G4double Det3_234_012 =                     << 1285   G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
1332     m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1286                                                 + m[A22]*Det2_34_01;
1333   G4double Det3_234_013 =                     << 1287   G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
1334     m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1288                                                 + m[A23]*Det2_34_01;
1335   G4double Det3_234_014 =                     << 1289   G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
1336     m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1290                                                 + m[A24]*Det2_34_01;
1337   G4double Det3_234_015 =                     << 1291   G4double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
1338     m[A20] * Det2_34_15 - m[A21] * Det2_34_05 << 1292                                                 + m[A25]*Det2_34_01;
1339   G4double Det3_234_023 =                     << 1293   G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
1340     m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1294                                                 + m[A23]*Det2_34_02;
1341   G4double Det3_234_024 =                     << 1295   G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
1342     m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1296                                                 + m[A24]*Det2_34_02;
1343   G4double Det3_234_025 =                     << 1297   G4double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05 
1344     m[A20] * Det2_34_25 - m[A22] * Det2_34_05 << 1298                                                 + m[A25]*Det2_34_02;
1345   G4double Det3_234_034 =                     << 1299   G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
1346     m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1300                                                 + m[A24]*Det2_34_03;
1347   G4double Det3_234_035 =                     << 1301   G4double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05                 
1348     m[A20] * Det2_34_35 - m[A23] * Det2_34_05 << 1302                                                 + m[A25]*Det2_34_03;
1349   G4double Det3_234_045 =                     << 1303   G4double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05 
1350     m[A20] * Det2_34_45 - m[A24] * Det2_34_05 << 1304                                                 + m[A25]*Det2_34_04;
1351   G4double Det3_234_123 =                     << 1305   G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
1352     m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1306                                                 + m[A23]*Det2_34_12;
1353   G4double Det3_234_124 =                     << 1307   G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
1354     m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1308                                                 + m[A24]*Det2_34_12;
1355   G4double Det3_234_125 =                     << 1309   G4double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15 
1356     m[A21] * Det2_34_25 - m[A22] * Det2_34_15 << 1310                                                 + m[A25]*Det2_34_12;
1357   G4double Det3_234_134 =                     << 1311   G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
1358     m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1312                                                 + m[A24]*Det2_34_13;
1359   G4double Det3_234_135 =                     << 1313   G4double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15 
1360     m[A21] * Det2_34_35 - m[A23] * Det2_34_15 << 1314                                                 + m[A25]*Det2_34_13;
1361   G4double Det3_234_145 =                     << 1315   G4double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15 
1362     m[A21] * Det2_34_45 - m[A24] * Det2_34_15 << 1316                                                 + m[A25]*Det2_34_14;
1363   G4double Det3_234_234 =                     << 1317   G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
1364     m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1318                                                 + m[A24]*Det2_34_23;
1365   G4double Det3_234_235 =                     << 1319   G4double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
1366     m[A22] * Det2_34_35 - m[A23] * Det2_34_25 << 1320                                                 + m[A25]*Det2_34_23;
1367   G4double Det3_234_245 =                     << 1321   G4double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
1368     m[A22] * Det2_34_45 - m[A24] * Det2_34_25 << 1322                                                 + m[A25]*Det2_34_24;
1369   G4double Det3_234_345 =                     << 1323   G4double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
1370     m[A23] * Det2_34_45 - m[A24] * Det2_34_35 << 1324                                                 + m[A25]*Det2_34_34;
1371   G4double Det3_235_012 =                     << 1325   G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02 
1372     m[A20] * Det2_35_12 - m[A21] * Det2_35_02 << 1326                                                 + m[A22]*Det2_35_01;
1373   G4double Det3_235_013 =                     << 1327   G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03 
1374     m[A20] * Det2_35_13 - m[A21] * Det2_35_03 << 1328                                                 + m[A23]*Det2_35_01;
1375   G4double Det3_235_014 =                     << 1329   G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04 
1376     m[A20] * Det2_35_14 - m[A21] * Det2_35_04 << 1330                                                 + m[A24]*Det2_35_01;
1377   G4double Det3_235_015 =                     << 1331   G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05 
1378     m[A20] * Det2_35_15 - m[A21] * Det2_35_05 << 1332                                                 + m[A25]*Det2_35_01;
1379   G4double Det3_235_023 =                     << 1333   G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03 
1380     m[A20] * Det2_35_23 - m[A22] * Det2_35_03 << 1334                                                 + m[A23]*Det2_35_02;
1381   G4double Det3_235_024 =                     << 1335   G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04 
1382     m[A20] * Det2_35_24 - m[A22] * Det2_35_04 << 1336                                                 + m[A24]*Det2_35_02;
1383   G4double Det3_235_025 =                     << 1337   G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05 
1384     m[A20] * Det2_35_25 - m[A22] * Det2_35_05 << 1338                                                 + m[A25]*Det2_35_02;
1385   G4double Det3_235_034 =                     << 1339   G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04 
1386     m[A20] * Det2_35_34 - m[A23] * Det2_35_04 << 1340                                                 + m[A24]*Det2_35_03;
1387   G4double Det3_235_035 =                     << 1341   G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05 
1388     m[A20] * Det2_35_35 - m[A23] * Det2_35_05 << 1342                                                 + m[A25]*Det2_35_03;
1389   G4double Det3_235_045 =                     << 1343   G4double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05 
1390     m[A20] * Det2_35_45 - m[A24] * Det2_35_05 << 1344                                                 + m[A25]*Det2_35_04;
1391   G4double Det3_235_123 =                     << 1345   G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13 
1392     m[A21] * Det2_35_23 - m[A22] * Det2_35_13 << 1346                                                 + m[A23]*Det2_35_12;
1393   G4double Det3_235_124 =                     << 1347   G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14 
1394     m[A21] * Det2_35_24 - m[A22] * Det2_35_14 << 1348                                                 + m[A24]*Det2_35_12;
1395   G4double Det3_235_125 =                     << 1349   G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15 
1396     m[A21] * Det2_35_25 - m[A22] * Det2_35_15 << 1350                                                 + m[A25]*Det2_35_12;
1397   G4double Det3_235_134 =                     << 1351   G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14 
1398     m[A21] * Det2_35_34 - m[A23] * Det2_35_14 << 1352                                                 + m[A24]*Det2_35_13;
1399   G4double Det3_235_135 =                     << 1353   G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15 
1400     m[A21] * Det2_35_35 - m[A23] * Det2_35_15 << 1354                                                 + m[A25]*Det2_35_13;
1401   G4double Det3_235_145 =                     << 1355   G4double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15 
1402     m[A21] * Det2_35_45 - m[A24] * Det2_35_15 << 1356                                                 + m[A25]*Det2_35_14;
1403   G4double Det3_235_234 =                     << 1357   G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24 
1404     m[A22] * Det2_35_34 - m[A23] * Det2_35_24 << 1358                                                 + m[A24]*Det2_35_23;
1405   G4double Det3_235_235 =                     << 1359   G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25 
1406     m[A22] * Det2_35_35 - m[A23] * Det2_35_25 << 1360                                                 + m[A25]*Det2_35_23;
1407   G4double Det3_235_245 =                     << 1361   G4double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25 
1408     m[A22] * Det2_35_45 - m[A24] * Det2_35_25 << 1362                                                 + m[A25]*Det2_35_24;
1409   G4double Det3_235_345 =                     << 1363   G4double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35 
1410     m[A23] * Det2_35_45 - m[A24] * Det2_35_35 << 1364                                                 + m[A25]*Det2_35_34;
1411   G4double Det3_245_012 =                     << 1365   G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02 
1412     m[A20] * Det2_45_12 - m[A21] * Det2_45_02 << 1366                                                 + m[A22]*Det2_45_01;
1413   G4double Det3_245_013 =                     << 1367   G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03 
1414     m[A20] * Det2_45_13 - m[A21] * Det2_45_03 << 1368                                                 + m[A23]*Det2_45_01;
1415   G4double Det3_245_014 =                     << 1369   G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04 
1416     m[A20] * Det2_45_14 - m[A21] * Det2_45_04 << 1370                                                 + m[A24]*Det2_45_01;
1417   G4double Det3_245_015 =                     << 1371   G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05 
1418     m[A20] * Det2_45_15 - m[A21] * Det2_45_05 << 1372                                                 + m[A25]*Det2_45_01;
1419   G4double Det3_245_023 =                     << 1373   G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03 
1420     m[A20] * Det2_45_23 - m[A22] * Det2_45_03 << 1374                                                 + m[A23]*Det2_45_02;
1421   G4double Det3_245_024 =                     << 1375   G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04 
1422     m[A20] * Det2_45_24 - m[A22] * Det2_45_04 << 1376                                                 + m[A24]*Det2_45_02;
1423   G4double Det3_245_025 =                     << 1377   G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05 
1424     m[A20] * Det2_45_25 - m[A22] * Det2_45_05 << 1378                                                 + m[A25]*Det2_45_02;
1425   G4double Det3_245_034 =                     << 1379   G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04 
1426     m[A20] * Det2_45_34 - m[A23] * Det2_45_04 << 1380                                                 + m[A24]*Det2_45_03;
1427   G4double Det3_245_035 =                     << 1381   G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05 
1428     m[A20] * Det2_45_35 - m[A23] * Det2_45_05 << 1382                                                 + m[A25]*Det2_45_03;
1429   G4double Det3_245_045 =                     << 1383   G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05 
1430     m[A20] * Det2_45_45 - m[A24] * Det2_45_05 << 1384                                                 + m[A25]*Det2_45_04;
1431   G4double Det3_245_123 =                     << 1385   G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13 
1432     m[A21] * Det2_45_23 - m[A22] * Det2_45_13 << 1386                                                 + m[A23]*Det2_45_12;
1433   G4double Det3_245_124 =                     << 1387   G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14 
1434     m[A21] * Det2_45_24 - m[A22] * Det2_45_14 << 1388                                                 + m[A24]*Det2_45_12;
1435   G4double Det3_245_125 =                     << 1389   G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15 
1436     m[A21] * Det2_45_25 - m[A22] * Det2_45_15 << 1390                                                 + m[A25]*Det2_45_12;
1437   G4double Det3_245_134 =                     << 1391   G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14 
1438     m[A21] * Det2_45_34 - m[A23] * Det2_45_14 << 1392                                                 + m[A24]*Det2_45_13;
1439   G4double Det3_245_135 =                     << 1393   G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15 
1440     m[A21] * Det2_45_35 - m[A23] * Det2_45_15 << 1394                                                 + m[A25]*Det2_45_13;
1441   G4double Det3_245_145 =                     << 1395   G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15 
1442     m[A21] * Det2_45_45 - m[A24] * Det2_45_15 << 1396                                                 + m[A25]*Det2_45_14;
1443   G4double Det3_245_234 =                     << 1397   G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24 
1444     m[A22] * Det2_45_34 - m[A23] * Det2_45_24 << 1398                                                 + m[A24]*Det2_45_23;
1445   G4double Det3_245_235 =                     << 1399   G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25 
1446     m[A22] * Det2_45_35 - m[A23] * Det2_45_25 << 1400                                                 + m[A25]*Det2_45_23;
1447   G4double Det3_245_245 =                     << 1401   G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25 
1448     m[A22] * Det2_45_45 - m[A24] * Det2_45_25 << 1402                                                 + m[A25]*Det2_45_24;
1449   G4double Det3_245_345 =                     << 1403   G4double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35 
1450     m[A23] * Det2_45_45 - m[A24] * Det2_45_35 << 1404                                                 + m[A25]*Det2_45_34;
1451   G4double Det3_345_012 =                     << 1405   G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02 
1452     m[A30] * Det2_45_12 - m[A31] * Det2_45_02 << 1406                                                 + m[A32]*Det2_45_01;
1453   G4double Det3_345_013 =                     << 1407   G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03 
1454     m[A30] * Det2_45_13 - m[A31] * Det2_45_03 << 1408                                                 + m[A33]*Det2_45_01;
1455   G4double Det3_345_014 =                     << 1409   G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04 
1456     m[A30] * Det2_45_14 - m[A31] * Det2_45_04 << 1410                                                 + m[A34]*Det2_45_01;
1457   G4double Det3_345_015 =                     << 1411   G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05 
1458     m[A30] * Det2_45_15 - m[A31] * Det2_45_05 << 1412                                                 + m[A35]*Det2_45_01;
1459   G4double Det3_345_023 =                     << 1413   G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03 
1460     m[A30] * Det2_45_23 - m[A32] * Det2_45_03 << 1414                                                 + m[A33]*Det2_45_02;
1461   G4double Det3_345_024 =                     << 1415   G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04 
1462     m[A30] * Det2_45_24 - m[A32] * Det2_45_04 << 1416                                                 + m[A34]*Det2_45_02;
1463   G4double Det3_345_025 =                     << 1417   G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05 
1464     m[A30] * Det2_45_25 - m[A32] * Det2_45_05 << 1418                                                 + m[A35]*Det2_45_02;
1465   G4double Det3_345_034 =                     << 1419   G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04 
1466     m[A30] * Det2_45_34 - m[A33] * Det2_45_04 << 1420                                                 + m[A34]*Det2_45_03;
1467   G4double Det3_345_035 =                     << 1421   G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05 
1468     m[A30] * Det2_45_35 - m[A33] * Det2_45_05 << 1422                                                 + m[A35]*Det2_45_03;
1469   G4double Det3_345_045 =                     << 1423   G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05 
1470     m[A30] * Det2_45_45 - m[A34] * Det2_45_05 << 1424                                                 + m[A35]*Det2_45_04;
1471   G4double Det3_345_123 =                     << 1425   G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13 
1472     m[A31] * Det2_45_23 - m[A32] * Det2_45_13 << 1426                                                 + m[A33]*Det2_45_12;
1473   G4double Det3_345_124 =                     << 1427   G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14 
1474     m[A31] * Det2_45_24 - m[A32] * Det2_45_14 << 1428                                                 + m[A34]*Det2_45_12;
1475   G4double Det3_345_125 =                     << 1429   G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15 
1476     m[A31] * Det2_45_25 - m[A32] * Det2_45_15 << 1430                                                 + m[A35]*Det2_45_12;
1477   G4double Det3_345_134 =                     << 1431   G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14 
1478     m[A31] * Det2_45_34 - m[A33] * Det2_45_14 << 1432                                                 + m[A34]*Det2_45_13;
1479   G4double Det3_345_135 =                     << 1433   G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15 
1480     m[A31] * Det2_45_35 - m[A33] * Det2_45_15 << 1434                                                 + m[A35]*Det2_45_13;
1481   G4double Det3_345_145 =                     << 1435   G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15 
1482     m[A31] * Det2_45_45 - m[A34] * Det2_45_15 << 1436                                                 + m[A35]*Det2_45_14;
1483   G4double Det3_345_234 =                     << 1437   G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24 
1484     m[A32] * Det2_45_34 - m[A33] * Det2_45_24 << 1438                                                 + m[A34]*Det2_45_23;
1485   G4double Det3_345_235 =                     << 1439   G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25 
1486     m[A32] * Det2_45_35 - m[A33] * Det2_45_25 << 1440                                                 + m[A35]*Det2_45_23;
1487   G4double Det3_345_245 =                     << 1441   G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25 
1488     m[A32] * Det2_45_45 - m[A34] * Det2_45_25 << 1442                                                 + m[A35]*Det2_45_24;
1489   G4double Det3_345_345 =                     << 1443   G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35 
1490     m[A33] * Det2_45_45 - m[A34] * Det2_45_35 << 1444                                                 + m[A35]*Det2_45_34;
1491                                                  1445 
1492   // Find all NECESSARY 4x4 dets:  (75 of the    1446   // Find all NECESSARY 4x4 dets:  (75 of them)
1493                                                  1447 
1494   G4double Det4_1234_0123 = m[A10] * Det3_234 << 1448   G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
1495                             m[A12] * Det3_234 << 1449                         + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1496   G4double Det4_1234_0124 = m[A10] * Det3_234 << 1450   G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
1497                             m[A12] * Det3_234 << 1451                         + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1498   G4double Det4_1234_0125 = m[A10] * Det3_234 << 1452   G4double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025 
1499                             m[A12] * Det3_234 << 1453                         + m[A12]*Det3_234_015 - m[A15]*Det3_234_012;
1500   G4double Det4_1234_0134 = m[A10] * Det3_234 << 1454   G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
1501                             m[A13] * Det3_234 << 1455                         + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1502   G4double Det4_1234_0135 = m[A10] * Det3_234 << 1456   G4double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
1503                             m[A13] * Det3_234 << 1457                         + m[A13]*Det3_234_015 - m[A15]*Det3_234_013;
1504   G4double Det4_1234_0145 = m[A10] * Det3_234 << 1458   G4double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
1505                             m[A14] * Det3_234 << 1459                         + m[A14]*Det3_234_015 - m[A15]*Det3_234_014;
1506   G4double Det4_1234_0234 = m[A10] * Det3_234 << 1460   G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
1507                             m[A13] * Det3_234 << 1461                         + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1508   G4double Det4_1234_0235 = m[A10] * Det3_234 << 1462   G4double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035 
1509                             m[A13] * Det3_234 << 1463                         + m[A13]*Det3_234_025 - m[A15]*Det3_234_023;
1510   G4double Det4_1234_0245 = m[A10] * Det3_234 << 1464   G4double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045 
1511                             m[A14] * Det3_234 << 1465                         + m[A14]*Det3_234_025 - m[A15]*Det3_234_024;
1512   G4double Det4_1234_0345 = m[A10] * Det3_234 << 1466   G4double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045 
1513                             m[A14] * Det3_234 << 1467                         + m[A14]*Det3_234_035 - m[A15]*Det3_234_034;
1514   G4double Det4_1234_1234 = m[A11] * Det3_234 << 1468   G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
1515                             m[A13] * Det3_234 << 1469                         + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1516   G4double Det4_1234_1235 = m[A11] * Det3_234 << 1470   G4double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135 
1517                             m[A13] * Det3_234 << 1471                         + m[A13]*Det3_234_125 - m[A15]*Det3_234_123;
1518   G4double Det4_1234_1245 = m[A11] * Det3_234 << 1472   G4double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145 
1519                             m[A14] * Det3_234 << 1473                         + m[A14]*Det3_234_125 - m[A15]*Det3_234_124;
1520   G4double Det4_1234_1345 = m[A11] * Det3_234 << 1474   G4double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145 
1521                             m[A14] * Det3_234 << 1475                         + m[A14]*Det3_234_135 - m[A15]*Det3_234_134;
1522   G4double Det4_1234_2345 = m[A12] * Det3_234 << 1476   G4double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245 
1523                             m[A14] * Det3_234 << 1477                         + m[A14]*Det3_234_235 - m[A15]*Det3_234_234;
1524   G4double Det4_1235_0123 = m[A10] * Det3_235 << 1478   G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023 
1525                             m[A12] * Det3_235 << 1479                         + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1526   G4double Det4_1235_0124 = m[A10] * Det3_235 << 1480   G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024 
1527                             m[A12] * Det3_235 << 1481                         + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1528   G4double Det4_1235_0125 = m[A10] * Det3_235 << 1482   G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025 
1529                             m[A12] * Det3_235 << 1483                         + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1530   G4double Det4_1235_0134 = m[A10] * Det3_235 << 1484   G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034 
1531                             m[A13] * Det3_235 << 1485                         + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1532   G4double Det4_1235_0135 = m[A10] * Det3_235 << 1486   G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035 
1533                             m[A13] * Det3_235 << 1487                         + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1534   G4double Det4_1235_0145 = m[A10] * Det3_235 << 1488   G4double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045 
1535                             m[A14] * Det3_235 << 1489                         + m[A14]*Det3_235_015 - m[A15]*Det3_235_014;
1536   G4double Det4_1235_0234 = m[A10] * Det3_235 << 1490   G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034 
1537                             m[A13] * Det3_235 << 1491                         + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1538   G4double Det4_1235_0235 = m[A10] * Det3_235 << 1492   G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035 
1539                             m[A13] * Det3_235 << 1493                         + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1540   G4double Det4_1235_0245 = m[A10] * Det3_235 << 1494   G4double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045 
1541                             m[A14] * Det3_235 << 1495                         + m[A14]*Det3_235_025 - m[A15]*Det3_235_024;
1542   G4double Det4_1235_0345 = m[A10] * Det3_235 << 1496   G4double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045 
1543                             m[A14] * Det3_235 << 1497                         + m[A14]*Det3_235_035 - m[A15]*Det3_235_034;
1544   G4double Det4_1235_1234 = m[A11] * Det3_235 << 1498   G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134 
1545                             m[A13] * Det3_235 << 1499                         + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1546   G4double Det4_1235_1235 = m[A11] * Det3_235 << 1500   G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135 
1547                             m[A13] * Det3_235 << 1501                         + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1548   G4double Det4_1235_1245 = m[A11] * Det3_235 << 1502   G4double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145 
1549                             m[A14] * Det3_235 << 1503                         + m[A14]*Det3_235_125 - m[A15]*Det3_235_124;
1550   G4double Det4_1235_1345 = m[A11] * Det3_235 << 1504   G4double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145 
1551                             m[A14] * Det3_235 << 1505                         + m[A14]*Det3_235_135 - m[A15]*Det3_235_134;
1552   G4double Det4_1235_2345 = m[A12] * Det3_235 << 1506   G4double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245 
1553                             m[A14] * Det3_235 << 1507                         + m[A14]*Det3_235_235 - m[A15]*Det3_235_234;
1554   G4double Det4_1245_0123 = m[A10] * Det3_245 << 1508   G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023 
1555                             m[A12] * Det3_245 << 1509                         + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1556   G4double Det4_1245_0124 = m[A10] * Det3_245 << 1510   G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024 
1557                             m[A12] * Det3_245 << 1511                         + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1558   G4double Det4_1245_0125 = m[A10] * Det3_245 << 1512   G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025 
1559                             m[A12] * Det3_245 << 1513                         + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1560   G4double Det4_1245_0134 = m[A10] * Det3_245 << 1514   G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034 
1561                             m[A13] * Det3_245 << 1515                         + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1562   G4double Det4_1245_0135 = m[A10] * Det3_245 << 1516   G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035 
1563                             m[A13] * Det3_245 << 1517                         + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1564   G4double Det4_1245_0145 = m[A10] * Det3_245 << 1518   G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045 
1565                             m[A14] * Det3_245 << 1519                         + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1566   G4double Det4_1245_0234 = m[A10] * Det3_245 << 1520   G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034 
1567                             m[A13] * Det3_245 << 1521                         + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1568   G4double Det4_1245_0235 = m[A10] * Det3_245 << 1522   G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035 
1569                             m[A13] * Det3_245 << 1523                         + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1570   G4double Det4_1245_0245 = m[A10] * Det3_245 << 1524   G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045 
1571                             m[A14] * Det3_245 << 1525                         + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1572   G4double Det4_1245_0345 = m[A10] * Det3_245 << 1526   G4double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045 
1573                             m[A14] * Det3_245 << 1527                         + m[A14]*Det3_245_035 - m[A15]*Det3_245_034;
1574   G4double Det4_1245_1234 = m[A11] * Det3_245 << 1528   G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134 
1575                             m[A13] * Det3_245 << 1529                         + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1576   G4double Det4_1245_1235 = m[A11] * Det3_245 << 1530   G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135 
1577                             m[A13] * Det3_245 << 1531                         + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1578   G4double Det4_1245_1245 = m[A11] * Det3_245 << 1532   G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145 
1579                             m[A14] * Det3_245 << 1533                         + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1580   G4double Det4_1245_1345 = m[A11] * Det3_245 << 1534   G4double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145 
1581                             m[A14] * Det3_245 << 1535                         + m[A14]*Det3_245_135 - m[A15]*Det3_245_134;
1582   G4double Det4_1245_2345 = m[A12] * Det3_245 << 1536   G4double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245 
1583                             m[A14] * Det3_245 << 1537                         + m[A14]*Det3_245_235 - m[A15]*Det3_245_234;
1584   G4double Det4_1345_0123 = m[A10] * Det3_345 << 1538   G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023 
1585                             m[A12] * Det3_345 << 1539                         + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1586   G4double Det4_1345_0124 = m[A10] * Det3_345 << 1540   G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024 
1587                             m[A12] * Det3_345 << 1541                         + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1588   G4double Det4_1345_0125 = m[A10] * Det3_345 << 1542   G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025 
1589                             m[A12] * Det3_345 << 1543                         + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1590   G4double Det4_1345_0134 = m[A10] * Det3_345 << 1544   G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034 
1591                             m[A13] * Det3_345 << 1545                         + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1592   G4double Det4_1345_0135 = m[A10] * Det3_345 << 1546   G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035 
1593                             m[A13] * Det3_345 << 1547                         + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1594   G4double Det4_1345_0145 = m[A10] * Det3_345 << 1548   G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045 
1595                             m[A14] * Det3_345 << 1549                         + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1596   G4double Det4_1345_0234 = m[A10] * Det3_345 << 1550   G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034 
1597                             m[A13] * Det3_345 << 1551                         + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1598   G4double Det4_1345_0235 = m[A10] * Det3_345 << 1552   G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035 
1599                             m[A13] * Det3_345 << 1553                         + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1600   G4double Det4_1345_0245 = m[A10] * Det3_345 << 1554   G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045 
1601                             m[A14] * Det3_345 << 1555                         + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1602   G4double Det4_1345_0345 = m[A10] * Det3_345 << 1556   G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045 
1603                             m[A14] * Det3_345 << 1557                         + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1604   G4double Det4_1345_1234 = m[A11] * Det3_345 << 1558   G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134 
1605                             m[A13] * Det3_345 << 1559                         + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1606   G4double Det4_1345_1235 = m[A11] * Det3_345 << 1560   G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135 
1607                             m[A13] * Det3_345 << 1561                         + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1608   G4double Det4_1345_1245 = m[A11] * Det3_345 << 1562   G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145 
1609                             m[A14] * Det3_345 << 1563                         + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1610   G4double Det4_1345_1345 = m[A11] * Det3_345 << 1564   G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145 
1611                             m[A14] * Det3_345 << 1565                         + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1612   G4double Det4_1345_2345 = m[A12] * Det3_345 << 1566   G4double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245 
1613                             m[A14] * Det3_345 << 1567                         + m[A14]*Det3_345_235 - m[A15]*Det3_345_234;
1614   G4double Det4_2345_0123 = m[A20] * Det3_345 << 1568   G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023 
1615                             m[A22] * Det3_345 << 1569                         + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1616   G4double Det4_2345_0124 = m[A20] * Det3_345 << 1570   G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024 
1617                             m[A22] * Det3_345 << 1571                         + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1618   G4double Det4_2345_0125 = m[A20] * Det3_345 << 1572   G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025 
1619                             m[A22] * Det3_345 << 1573                         + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1620   G4double Det4_2345_0134 = m[A20] * Det3_345 << 1574   G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034 
1621                             m[A23] * Det3_345 << 1575                         + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1622   G4double Det4_2345_0135 = m[A20] * Det3_345 << 1576   G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035 
1623                             m[A23] * Det3_345 << 1577                         + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1624   G4double Det4_2345_0145 = m[A20] * Det3_345 << 1578   G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045 
1625                             m[A24] * Det3_345 << 1579                         + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1626   G4double Det4_2345_0234 = m[A20] * Det3_345 << 1580   G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034 
1627                             m[A23] * Det3_345 << 1581                         + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1628   G4double Det4_2345_0235 = m[A20] * Det3_345 << 1582   G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035 
1629                             m[A23] * Det3_345 << 1583                         + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1630   G4double Det4_2345_0245 = m[A20] * Det3_345 << 1584   G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045 
1631                             m[A24] * Det3_345 << 1585                         + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1632   G4double Det4_2345_0345 = m[A20] * Det3_345 << 1586   G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045 
1633                             m[A24] * Det3_345 << 1587                         + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1634   G4double Det4_2345_1234 = m[A21] * Det3_345 << 1588   G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134 
1635                             m[A23] * Det3_345 << 1589                         + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1636   G4double Det4_2345_1235 = m[A21] * Det3_345 << 1590   G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135 
1637                             m[A23] * Det3_345 << 1591                         + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1638   G4double Det4_2345_1245 = m[A21] * Det3_345 << 1592   G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145 
1639                             m[A24] * Det3_345 << 1593                         + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1640   G4double Det4_2345_1345 = m[A21] * Det3_345 << 1594   G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145 
1641                             m[A24] * Det3_345 << 1595                         + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1642   G4double Det4_2345_2345 = m[A22] * Det3_345 << 1596   G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245 
1643                             m[A24] * Det3_345 << 1597                         + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1644                                                  1598 
1645   // Find all NECESSARY 5x5 dets:  (36 of the    1599   // Find all NECESSARY 5x5 dets:  (36 of them)
1646                                                  1600 
1647   G4double Det5_01234_01234 =                 << 1601   G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234 
1648     m[A00] * Det4_1234_1234 - m[A01] * Det4_1 << 1602     + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1649     m[A02] * Det4_1234_0134 - m[A03] * Det4_1 << 1603   G4double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
1650   G4double Det5_01234_01235 =                 << 1604     + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123;
1651     m[A00] * Det4_1234_1235 - m[A01] * Det4_1 << 1605   G4double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
1652     m[A02] * Det4_1234_0135 - m[A03] * Det4_1 << 1606     + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124;
1653   G4double Det5_01234_01245 =                 << 1607   G4double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
1654     m[A00] * Det4_1234_1245 - m[A01] * Det4_1 << 1608     + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134;
1655     m[A02] * Det4_1234_0145 - m[A04] * Det4_1 << 1609   G4double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
1656   G4double Det5_01234_01345 =                 << 1610     + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234;
1657     m[A00] * Det4_1234_1345 - m[A01] * Det4_1 << 1611   G4double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
1658     m[A03] * Det4_1234_0145 - m[A04] * Det4_1 << 1612     + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234;
1659   G4double Det5_01234_02345 =                 << 1613   G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234 
1660     m[A00] * Det4_1234_2345 - m[A02] * Det4_1 << 1614     + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1661     m[A03] * Det4_1234_0245 - m[A04] * Det4_1 << 1615   G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235 
1662   G4double Det5_01234_12345 =                 << 1616     + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1663     m[A01] * Det4_1234_2345 - m[A02] * Det4_1 << 1617   G4double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245 
1664     m[A03] * Det4_1234_1245 - m[A04] * Det4_1 << 1618     + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124;
1665   G4double Det5_01235_01234 =                 << 1619   G4double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345 
1666     m[A00] * Det4_1235_1234 - m[A01] * Det4_1 << 1620     + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134;
1667     m[A02] * Det4_1235_0134 - m[A03] * Det4_1 << 1621   G4double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345 
1668   G4double Det5_01235_01235 =                 << 1622     + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234;
1669     m[A00] * Det4_1235_1235 - m[A01] * Det4_1 << 1623   G4double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345 
1670     m[A02] * Det4_1235_0135 - m[A03] * Det4_1 << 1624     + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234;
1671   G4double Det5_01235_01245 =                 << 1625   G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234 
1672     m[A00] * Det4_1235_1245 - m[A01] * Det4_1 << 1626     + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1673     m[A02] * Det4_1235_0145 - m[A04] * Det4_1 << 1627   G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235 
1674   G4double Det5_01235_01345 =                 << 1628     + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1675     m[A00] * Det4_1235_1345 - m[A01] * Det4_1 << 1629   G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245 
1676     m[A03] * Det4_1235_0145 - m[A04] * Det4_1 << 1630     + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1677   G4double Det5_01235_02345 =                 << 1631   G4double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345 
1678     m[A00] * Det4_1235_2345 - m[A02] * Det4_1 << 1632     + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134;
1679     m[A03] * Det4_1235_0245 - m[A04] * Det4_1 << 1633   G4double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345 
1680   G4double Det5_01235_12345 =                 << 1634     + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234;
1681     m[A01] * Det4_1235_2345 - m[A02] * Det4_1 << 1635   G4double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345 
1682     m[A03] * Det4_1235_1245 - m[A04] * Det4_1 << 1636     + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234;
1683   G4double Det5_01245_01234 =                 << 1637   G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234 
1684     m[A00] * Det4_1245_1234 - m[A01] * Det4_1 << 1638     + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1685     m[A02] * Det4_1245_0134 - m[A03] * Det4_1 << 1639   G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235 
1686   G4double Det5_01245_01235 =                 << 1640     + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1687     m[A00] * Det4_1245_1235 - m[A01] * Det4_1 << 1641   G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245 
1688     m[A02] * Det4_1245_0135 - m[A03] * Det4_1 << 1642     + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1689   G4double Det5_01245_01245 =                 << 1643   G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345 
1690     m[A00] * Det4_1245_1245 - m[A01] * Det4_1 << 1644     + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1691     m[A02] * Det4_1245_0145 - m[A04] * Det4_1 << 1645   G4double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345 
1692   G4double Det5_01245_01345 =                 << 1646     + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234;
1693     m[A00] * Det4_1245_1345 - m[A01] * Det4_1 << 1647   G4double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345 
1694     m[A03] * Det4_1245_0145 - m[A04] * Det4_1 << 1648     + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234; //
1695   G4double Det5_01245_02345 =                 << 1649   G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234 
1696     m[A00] * Det4_1245_2345 - m[A02] * Det4_1 << 1650     + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1697     m[A03] * Det4_1245_0245 - m[A04] * Det4_1 << 1651   G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235 
1698   G4double Det5_01245_12345 =                 << 1652     + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1699     m[A01] * Det4_1245_2345 - m[A02] * Det4_1 << 1653   G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245 
1700     m[A03] * Det4_1245_1245 - m[A04] * Det4_1 << 1654     + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1701   G4double Det5_01345_01234 =                 << 1655   G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345 
1702     m[A00] * Det4_1345_1234 - m[A01] * Det4_1 << 1656     + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1703     m[A02] * Det4_1345_0134 - m[A03] * Det4_1 << 1657   G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345 
1704   G4double Det5_01345_01235 =                 << 1658     + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1705     m[A00] * Det4_1345_1235 - m[A01] * Det4_1 << 1659   G4double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345 
1706     m[A02] * Det4_1345_0135 - m[A03] * Det4_1 << 1660     + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234;
1707   G4double Det5_01345_01245 =                 << 1661   G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234 
1708     m[A00] * Det4_1345_1245 - m[A01] * Det4_1 << 1662     + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1709     m[A02] * Det4_1345_0145 - m[A04] * Det4_1 << 1663   G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235 
1710   G4double Det5_01345_01345 =                 << 1664     + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1711     m[A00] * Det4_1345_1345 - m[A01] * Det4_1 << 1665   G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245 
1712     m[A03] * Det4_1345_0145 - m[A04] * Det4_1 << 1666     + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1713   G4double Det5_01345_02345 =                 << 1667   G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345 
1714     m[A00] * Det4_1345_2345 - m[A02] * Det4_1 << 1668     + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1715     m[A03] * Det4_1345_0245 - m[A04] * Det4_1 << 1669   G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345 
1716   G4double Det5_01345_12345 =                 << 1670     + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1717     m[A01] * Det4_1345_2345 - m[A02] * Det4_1 << 1671   G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345 
1718     m[A03] * Det4_1345_1245 - m[A04] * Det4_1 << 1672     + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1719     m[A05] * Det4_1345_1234;  //              << 1673 
1720   G4double Det5_02345_01234 =                 << 1674   // Find the determinant 
1721     m[A00] * Det4_2345_1234 - m[A01] * Det4_2 << 1675 
1722     m[A02] * Det4_2345_0134 - m[A03] * Det4_2 << 1676   G4double det =  m[A00]*Det5_12345_12345 
1723   G4double Det5_02345_01235 =                 << 1677                 - m[A01]*Det5_12345_02345 
1724     m[A00] * Det4_2345_1235 - m[A01] * Det4_2 << 1678                 + m[A02]*Det5_12345_01345 
1725     m[A02] * Det4_2345_0135 - m[A03] * Det4_2 << 1679                 - m[A03]*Det5_12345_01245 
1726   G4double Det5_02345_01245 =                 << 1680                 + m[A04]*Det5_12345_01235 
1727     m[A00] * Det4_2345_1245 - m[A01] * Det4_2 << 1681                 - m[A05]*Det5_12345_01234;
1728     m[A02] * Det4_2345_0145 - m[A04] * Det4_2 << 
1729   G4double Det5_02345_01345 =                 << 
1730     m[A00] * Det4_2345_1345 - m[A01] * Det4_2 << 
1731     m[A03] * Det4_2345_0145 - m[A04] * Det4_2 << 
1732   G4double Det5_02345_02345 =                 << 
1733     m[A00] * Det4_2345_2345 - m[A02] * Det4_2 << 
1734     m[A03] * Det4_2345_0245 - m[A04] * Det4_2 << 
1735   G4double Det5_02345_12345 =                 << 
1736     m[A01] * Det4_2345_2345 - m[A02] * Det4_2 << 
1737     m[A03] * Det4_2345_1245 - m[A04] * Det4_2 << 
1738   G4double Det5_12345_01234 =                 << 
1739     m[A10] * Det4_2345_1234 - m[A11] * Det4_2 << 
1740     m[A12] * Det4_2345_0134 - m[A13] * Det4_2 << 
1741   G4double Det5_12345_01235 =                 << 
1742     m[A10] * Det4_2345_1235 - m[A11] * Det4_2 << 
1743     m[A12] * Det4_2345_0135 - m[A13] * Det4_2 << 
1744   G4double Det5_12345_01245 =                 << 
1745     m[A10] * Det4_2345_1245 - m[A11] * Det4_2 << 
1746     m[A12] * Det4_2345_0145 - m[A14] * Det4_2 << 
1747   G4double Det5_12345_01345 =                 << 
1748     m[A10] * Det4_2345_1345 - m[A11] * Det4_2 << 
1749     m[A13] * Det4_2345_0145 - m[A14] * Det4_2 << 
1750   G4double Det5_12345_02345 =                 << 
1751     m[A10] * Det4_2345_2345 - m[A12] * Det4_2 << 
1752     m[A13] * Det4_2345_0245 - m[A14] * Det4_2 << 
1753   G4double Det5_12345_12345 =                 << 
1754     m[A11] * Det4_2345_2345 - m[A12] * Det4_2 << 
1755     m[A13] * Det4_2345_1245 - m[A14] * Det4_2 << 
1756                                               << 
1757   // Find the determinant                     << 
1758                                               << 
1759   G4double det = m[A00] * Det5_12345_12345 -  << 
1760                  m[A02] * Det5_12345_01345 -  << 
1761                  m[A04] * Det5_12345_01235 -  << 
1762                                                  1682 
1763   if(det == 0)                                << 1683   if ( det == 0 )
1764   {                                           << 1684   {  
1765     ifail = 1;                                   1685     ifail = 1;
1766     return;                                      1686     return;
1767   }                                           << 1687   } 
1768                                                  1688 
1769   G4double oneOverDet = 1.0 / det;            << 1689   G4double oneOverDet = 1.0/det;
1770   G4double mn1OverDet = -oneOverDet;          << 1690   G4double mn1OverDet = - oneOverDet;
1771                                                  1691 
1772   m[A00] = Det5_12345_12345 * oneOverDet;     << 1692   m[A00] =  Det5_12345_12345*oneOverDet;
1773   m[A01] = Det5_02345_12345 * mn1OverDet;     << 1693   m[A01] =  Det5_02345_12345*mn1OverDet;
1774   m[A02] = Det5_01345_12345 * oneOverDet;     << 1694   m[A02] =  Det5_01345_12345*oneOverDet;
1775   m[A03] = Det5_01245_12345 * mn1OverDet;     << 1695   m[A03] =  Det5_01245_12345*mn1OverDet;
1776   m[A04] = Det5_01235_12345 * oneOverDet;     << 1696   m[A04] =  Det5_01235_12345*oneOverDet;
1777   m[A05] = Det5_01234_12345 * mn1OverDet;     << 1697   m[A05] =  Det5_01234_12345*mn1OverDet;
1778                                               << 1698 
1779   m[A10] = Det5_12345_02345 * mn1OverDet;     << 1699   m[A10] =  Det5_12345_02345*mn1OverDet;
1780   m[A11] = Det5_02345_02345 * oneOverDet;     << 1700   m[A11] =  Det5_02345_02345*oneOverDet;
1781   m[A12] = Det5_01345_02345 * mn1OverDet;     << 1701   m[A12] =  Det5_01345_02345*mn1OverDet;
1782   m[A13] = Det5_01245_02345 * oneOverDet;     << 1702   m[A13] =  Det5_01245_02345*oneOverDet;
1783   m[A14] = Det5_01235_02345 * mn1OverDet;     << 1703   m[A14] =  Det5_01235_02345*mn1OverDet;
1784   m[A15] = Det5_01234_02345 * oneOverDet;     << 1704   m[A15] =  Det5_01234_02345*oneOverDet;
1785                                               << 1705 
1786   m[A20] = Det5_12345_01345 * oneOverDet;     << 1706   m[A20] =  Det5_12345_01345*oneOverDet;
1787   m[A21] = Det5_02345_01345 * mn1OverDet;     << 1707   m[A21] =  Det5_02345_01345*mn1OverDet;
1788   m[A22] = Det5_01345_01345 * oneOverDet;     << 1708   m[A22] =  Det5_01345_01345*oneOverDet;
1789   m[A23] = Det5_01245_01345 * mn1OverDet;     << 1709   m[A23] =  Det5_01245_01345*mn1OverDet;
1790   m[A24] = Det5_01235_01345 * oneOverDet;     << 1710   m[A24] =  Det5_01235_01345*oneOverDet;
1791   m[A25] = Det5_01234_01345 * mn1OverDet;     << 1711   m[A25] =  Det5_01234_01345*mn1OverDet;
1792                                               << 1712 
1793   m[A30] = Det5_12345_01245 * mn1OverDet;     << 1713   m[A30] =  Det5_12345_01245*mn1OverDet;
1794   m[A31] = Det5_02345_01245 * oneOverDet;     << 1714   m[A31] =  Det5_02345_01245*oneOverDet;
1795   m[A32] = Det5_01345_01245 * mn1OverDet;     << 1715   m[A32] =  Det5_01345_01245*mn1OverDet;
1796   m[A33] = Det5_01245_01245 * oneOverDet;     << 1716   m[A33] =  Det5_01245_01245*oneOverDet;
1797   m[A34] = Det5_01235_01245 * mn1OverDet;     << 1717   m[A34] =  Det5_01235_01245*mn1OverDet;
1798   m[A35] = Det5_01234_01245 * oneOverDet;     << 1718   m[A35] =  Det5_01234_01245*oneOverDet;
1799                                               << 1719 
1800   m[A40] = Det5_12345_01235 * oneOverDet;     << 1720   m[A40] =  Det5_12345_01235*oneOverDet;
1801   m[A41] = Det5_02345_01235 * mn1OverDet;     << 1721   m[A41] =  Det5_02345_01235*mn1OverDet;
1802   m[A42] = Det5_01345_01235 * oneOverDet;     << 1722   m[A42] =  Det5_01345_01235*oneOverDet;
1803   m[A43] = Det5_01245_01235 * mn1OverDet;     << 1723   m[A43] =  Det5_01245_01235*mn1OverDet;
1804   m[A44] = Det5_01235_01235 * oneOverDet;     << 1724   m[A44] =  Det5_01235_01235*oneOverDet;
1805   m[A45] = Det5_01234_01235 * mn1OverDet;     << 1725   m[A45] =  Det5_01234_01235*mn1OverDet;
1806                                               << 1726 
1807   m[A50] = Det5_12345_01234 * mn1OverDet;     << 1727   m[A50] =  Det5_12345_01234*mn1OverDet;
1808   m[A51] = Det5_02345_01234 * oneOverDet;     << 1728   m[A51] =  Det5_02345_01234*oneOverDet;
1809   m[A52] = Det5_01345_01234 * mn1OverDet;     << 1729   m[A52] =  Det5_01345_01234*mn1OverDet;
1810   m[A53] = Det5_01245_01234 * oneOverDet;     << 1730   m[A53] =  Det5_01245_01234*oneOverDet;
1811   m[A54] = Det5_01235_01234 * mn1OverDet;     << 1731   m[A54] =  Det5_01235_01234*mn1OverDet;
1812   m[A55] = Det5_01234_01234 * oneOverDet;     << 1732   m[A55] =  Det5_01234_01234*oneOverDet;
1813                                                  1733 
1814   return;                                        1734   return;
1815 }                                                1735 }
1816                                                  1736