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 7.1)


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