G4ErrorMatrix.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4ErrorMatrix.cc 69766 2013-05-14 14:33:55Z gcosmo $
00027 //
00028 // ------------------------------------------------------------
00029 //      GEANT 4 class implementation file
00030 // ------------------------------------------------------------
00031 
00032 #include "globals.hh"
00033 
00034 #include <cmath>
00035 #include <iostream>
00036 
00037 #include "G4ErrorMatrix.hh"
00038 #include "G4ErrorSymMatrix.hh"
00039 
00040 // Simple operation for all elements
00041 
00042 #define SIMPLE_UOP(OPER)                            \
00043     G4ErrorMatrixIter a=m.begin();                  \
00044     G4ErrorMatrixIter e=m.end();                    \
00045    for(;a!=e; a++) (*a) OPER t;
00046 
00047 #define SIMPLE_BOP(OPER)                            \
00048     G4ErrorMatrixIter a=m.begin();                  \
00049     G4ErrorMatrixConstIter b=mat2.m.begin();          \
00050     G4ErrorMatrixIter e=m.end();                    \
00051    for(;a!=e; a++, b++) (*a) OPER (*b);
00052 
00053 #define SIMPLE_TOP(OPER)                            \
00054     G4ErrorMatrixConstIter a=mat1.m.begin();          \
00055     G4ErrorMatrixConstIter b=mat2.m.begin();          \
00056     G4ErrorMatrixIter t=mret.m.begin();             \
00057     G4ErrorMatrixConstIter e=mat1.m.end();            \
00058    for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
00059 
00060 // Static functions.
00061 
00062 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
00063    if (r1!=r2 || c1!=c2)  { \
00064     G4ErrorMatrix::error("Range error in Matrix function " #fun "(1).");  \
00065    }
00066 
00067 #define CHK_DIM_1(c1,r2,fun) \
00068    if (c1!=r2) { \
00069      G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
00070    }
00071 
00072 // Constructors. (Default constructors are inlined and in .icc file)
00073 
00074 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q)
00075    : m(p*q), nrow(p), ncol(q)
00076 {
00077   size = nrow * ncol;
00078 }
00079 
00080 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q,G4int init)
00081    : m(p*q), nrow(p), ncol(q)
00082 {
00083    size = nrow * ncol;
00084 
00085    if (size > 0)
00086    {
00087       switch(init)
00088       {
00089       case 0:
00090          break;
00091 
00092       case 1:
00093          {
00094             if ( ncol == nrow )
00095             {
00096                G4ErrorMatrixIter a = m.begin();
00097                G4ErrorMatrixIter b = m.end();
00098                for( ; a<b; a+=(ncol+1)) *a = 1.0;
00099             } else {
00100                error("Invalid dimension in G4ErrorMatrix(G4int,G4int,1).");
00101             }
00102             break;
00103          }
00104       default:
00105          error("G4ErrorMatrix: initialization must be either 0 or 1.");
00106       }
00107    }
00108 }
00109 
00110 //
00111 // Destructor
00112 //
00113 G4ErrorMatrix::~G4ErrorMatrix()
00114 {
00115 }
00116 
00117 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatrix &mat1)
00118    : m(mat1.size), nrow(mat1.nrow), ncol(mat1.ncol), size(mat1.size)
00119 {
00120    m = mat1.m;
00121 }
00122 
00123 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymMatrix &mat1)
00124    : m(mat1.nrow*mat1.nrow), nrow(mat1.nrow), ncol(mat1.nrow)
00125 {
00126    size = nrow * ncol;
00127 
00128    G4int n = ncol;
00129    G4ErrorMatrixConstIter sjk = mat1.m.begin();
00130    G4ErrorMatrixIter m1j = m.begin();
00131    G4ErrorMatrixIter mj  = m.begin();
00132    // j >= k
00133    for(G4int j=1;j<=nrow;j++)
00134    {
00135       G4ErrorMatrixIter mjk = mj;
00136       G4ErrorMatrixIter mkj = m1j;
00137       for(G4int k=1;k<=j;k++)
00138       {
00139          *(mjk++) = *sjk;
00140          if(j!=k) *mkj = *sjk;
00141          sjk++;
00142          mkj += n;
00143       }
00144       mj += n;
00145       m1j++;
00146    }
00147 }
00148 
00149 //
00150 //
00151 // Sub matrix
00152 //
00153 //
00154 
00155 G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row, G4int max_row,
00156                          G4int min_col,G4int max_col) const
00157 {
00158   G4ErrorMatrix mret(max_row-min_row+1,max_col-min_col+1);
00159   if(max_row > num_row() || max_col >num_col())
00160     { error("G4ErrorMatrix::sub: Index out of range"); }
00161   G4ErrorMatrixIter a = mret.m.begin();
00162   G4int nc = num_col();
00163   G4ErrorMatrixConstIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
00164   
00165   for(G4int irow=1; irow<=mret.num_row(); irow++)
00166   {
00167     G4ErrorMatrixConstIter brc = b1;
00168     for(G4int icol=1; icol<=mret.num_col(); icol++)
00169     {
00170       *(a++) = *(brc++);
00171     }
00172     b1 += nc;
00173   }
00174   return mret;
00175 }
00176 
00177 void G4ErrorMatrix::sub(G4int row,G4int col,const G4ErrorMatrix &mat1)
00178 {
00179   if(row <1 || row+mat1.num_row()-1 > num_row() || 
00180      col <1 || col+mat1.num_col()-1 > num_col()   )
00181     { error("G4ErrorMatrix::sub: Index out of range"); }
00182   G4ErrorMatrixConstIter a = mat1.m.begin();
00183   G4int nc = num_col();
00184   G4ErrorMatrixIter b1 = m.begin() + (row - 1) * nc + col - 1;
00185   
00186   for(G4int irow=1; irow<=mat1.num_row(); irow++)
00187   {
00188     G4ErrorMatrixIter brc = b1;
00189     for(G4int icol=1; icol<=mat1.num_col(); icol++)
00190     {
00191       *(brc++) = *(a++);
00192     }
00193     b1 += nc;
00194   }
00195 }
00196 
00197 //
00198 // Direct sum of two matricies
00199 //
00200 
00201 G4ErrorMatrix dsum(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
00202 {
00203   G4ErrorMatrix mret(mat1.num_row() + mat2.num_row(),
00204                      mat1.num_col() + mat2.num_col(), 0);
00205   mret.sub(1,1,mat1);
00206   mret.sub(mat1.num_row()+1,mat1.num_col()+1,mat2);
00207   return mret;
00208 }
00209 
00210 /* -----------------------------------------------------------------------
00211    This section contains support routines for matrix.h. This section contains
00212    The two argument functions +,-. They call the copy constructor and +=,-=.
00213    ----------------------------------------------------------------------- */
00214 
00215 G4ErrorMatrix G4ErrorMatrix::operator- () const 
00216 {
00217    G4ErrorMatrix mat2(nrow, ncol);
00218    G4ErrorMatrixConstIter a=m.begin();
00219    G4ErrorMatrixIter b=mat2.m.begin();
00220    G4ErrorMatrixConstIter e=m.end();
00221    for(;a<e; a++, b++) (*b) = -(*a);
00222    return mat2;
00223 }
00224 
00225 G4ErrorMatrix operator+(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
00226 {
00227   G4ErrorMatrix mret(mat1.nrow, mat1.ncol);
00228   CHK_DIM_2(mat1.num_row(),mat2.num_row(), mat1.num_col(),mat2.num_col(),+);
00229   SIMPLE_TOP(+)
00230   return mret;
00231 }
00232 
00233 //
00234 // operator -
00235 //
00236 
00237 G4ErrorMatrix operator-(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
00238 {
00239   G4ErrorMatrix mret(mat1.num_row(), mat1.num_col());
00240   CHK_DIM_2(mat1.num_row(),mat2.num_row(),
00241                          mat1.num_col(),mat2.num_col(),-);
00242   SIMPLE_TOP(-)
00243   return mret;
00244 }
00245 
00246 /* -----------------------------------------------------------------------
00247    This section contains support routines for matrix.h. This file contains
00248    The two argument functions *,/. They call copy constructor and then /=,*=.
00249    ----------------------------------------------------------------------- */
00250 
00251 G4ErrorMatrix operator/(const G4ErrorMatrix &mat1,G4double t)
00252 {
00253   G4ErrorMatrix mret(mat1);
00254   mret /= t;
00255   return mret;
00256 }
00257 
00258 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,G4double t)
00259 {
00260   G4ErrorMatrix mret(mat1);
00261   mret *= t;
00262   return mret;
00263 }
00264 
00265 G4ErrorMatrix operator*(G4double t,const G4ErrorMatrix &mat1)
00266 {
00267   G4ErrorMatrix mret(mat1);
00268   mret *= t;
00269   return mret;
00270 }
00271 
00272 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
00273 {
00274   // initialize matrix to 0.0
00275   G4ErrorMatrix mret(mat1.nrow,mat2.ncol,0);
00276   CHK_DIM_1(mat1.ncol,mat2.nrow,*);
00277 
00278   G4int m1cols = mat1.ncol;
00279   G4int m2cols = mat2.ncol;
00280 
00281   for (G4int i=0; i<mat1.nrow; i++)
00282   {
00283      for (G4int j=0; j<m1cols; j++) 
00284      {
00285          G4double temp = mat1.m[i*m1cols+j];
00286          G4ErrorMatrixIter pt = mret.m.begin() + i*m2cols;
00287 
00288         // Loop over k (the column index in matrix mat2)
00289         G4ErrorMatrixConstIter pb = mat2.m.begin() + m2cols*j;
00290         const G4ErrorMatrixConstIter pblast = pb + m2cols;
00291         while (pb < pblast)
00292         {
00293            (*pt) += temp * (*pb);
00294            pb++;
00295            pt++;
00296         }
00297      }
00298   }
00299   return mret;
00300 }
00301 
00302 /* -----------------------------------------------------------------------
00303    This section contains the assignment and inplace operators =,+=,-=,*=,/=.
00304    ----------------------------------------------------------------------- */
00305 
00306 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorMatrix &mat2)
00307 {
00308   CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
00309   SIMPLE_BOP(+=)
00310   return (*this);
00311 }
00312 
00313 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorMatrix &mat2)
00314 {
00315   CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
00316   SIMPLE_BOP(-=)
00317   return (*this);
00318 }
00319 
00320 G4ErrorMatrix & G4ErrorMatrix::operator/=(G4double t)
00321 {
00322   SIMPLE_UOP(/=)
00323   return (*this);
00324 }
00325 
00326 G4ErrorMatrix & G4ErrorMatrix::operator*=(G4double t)
00327 {
00328   SIMPLE_UOP(*=)
00329   return (*this);
00330 }
00331 
00332 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorMatrix &mat1)
00333 {
00334    if (&mat1 == this) { return *this; }
00335 
00336    if(mat1.nrow*mat1.ncol != size) //??fixme?? mat1.size != size
00337    {
00338       size = mat1.nrow * mat1.ncol;
00339       m.resize(size); //??fixme?? if (size < mat1.size) m.resize(mat1.size);
00340    }
00341    nrow = mat1.nrow;
00342    ncol = mat1.ncol;
00343    m = mat1.m;
00344    return (*this);
00345 }
00346 
00347 
00348 // Print the G4ErrorMatrix.
00349 
00350 std::ostream& operator<<(std::ostream &os, const G4ErrorMatrix &q)
00351 {
00352   os << "\n";
00353 
00354   // Fixed format needs 3 extra characters for field,
00355   // while scientific needs 7
00356 
00357   G4int width;
00358   if(os.flags() & std::ios::fixed)
00359     { width = os.precision()+3; }
00360   else
00361     { width = os.precision()+7; }
00362   for(G4int irow = 1; irow<= q.num_row(); irow++)
00363     {
00364       for(G4int icol = 1; icol <= q.num_col(); icol++)
00365         {
00366           os.width(width);
00367           os << q(irow,icol) << " ";
00368         }
00369       os << G4endl;
00370     }
00371   return os;
00372 }
00373 
00374 G4ErrorMatrix G4ErrorMatrix::T() const
00375 {
00376    G4ErrorMatrix mret(ncol,nrow);
00377     G4ErrorMatrixConstIter pl = m.end();
00378     G4ErrorMatrixConstIter pme = m.begin();
00379     G4ErrorMatrixIter pt = mret.m.begin();
00380     G4ErrorMatrixIter ptl = mret.m.end();
00381    for (; pme < pl; pme++, pt+=nrow)
00382    {
00383       if (pt >= ptl) { pt -= (size-1); }
00384       (*pt) = (*pme);
00385    }
00386    return mret;
00387 }
00388 
00389 G4ErrorMatrix
00390 G4ErrorMatrix::apply(G4double (*f)(G4double, G4int, G4int)) const
00391 {
00392   G4ErrorMatrix mret(num_row(),num_col());
00393   G4ErrorMatrixConstIter a = m.begin();
00394   G4ErrorMatrixIter b = mret.m.begin();
00395   for(G4int ir=1;ir<=num_row();ir++)
00396   {
00397     for(G4int ic=1;ic<=num_col();ic++)
00398     {
00399       *(b++) = (*f)(*(a++), ir, ic);
00400     }
00401   }
00402   return mret;
00403 }
00404 
00405 G4int G4ErrorMatrix::dfinv_matrix(G4int *ir) {
00406   if (num_col()!=num_row())
00407     { error("dfinv_matrix: G4ErrorMatrix is not NxN"); }
00408   G4int n = num_col();
00409   if (n==1) { return 0; }
00410 
00411   G4double s31, s32;
00412   G4double s33, s34;
00413 
00414   G4ErrorMatrixIter m11 = m.begin();
00415   G4ErrorMatrixIter m12 = m11 + 1;
00416   G4ErrorMatrixIter m21 = m11 + n;
00417   G4ErrorMatrixIter m22 = m12 + n;
00418   *m21 = -(*m22) * (*m11) * (*m21);
00419   *m12 = -(*m12);
00420   if (n>2)
00421   {
00422     G4ErrorMatrixIter mi = m.begin() + 2 * n;
00423     G4ErrorMatrixIter mii= m.begin() + 2 * n + 2;
00424     G4ErrorMatrixIter mimim = m.begin() + n + 1;
00425     for (G4int i=3;i<=n;i++)
00426     {
00427       G4int im2 = i - 2;
00428       G4ErrorMatrixIter mj = m.begin();
00429       G4ErrorMatrixIter mji = mj + i - 1;
00430       G4ErrorMatrixIter mij = mi;
00431       for (G4int j=1;j<=im2;j++)
00432       { 
00433         s31 = 0.0;
00434         s32 = *mji;
00435         G4ErrorMatrixIter mkj = mj + j - 1;
00436         G4ErrorMatrixIter mik = mi + j - 1;
00437         G4ErrorMatrixIter mjkp = mj + j;
00438         G4ErrorMatrixIter mkpi = mj + n + i - 1;
00439         for (G4int k=j;k<=im2;k++)
00440         {
00441           s31 += (*mkj) * (*(mik++));
00442           s32 += (*(mjkp++)) * (*mkpi);
00443           mkj += n;
00444           mkpi += n;
00445         }
00446         *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
00447         *mji = -s32;
00448         mj += n;
00449         mji += n;
00450         mij++;
00451       }
00452       *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
00453       *(mimim+1) = -(*(mimim+1));
00454       mi += n;
00455       mimim += (n+1);
00456       mii += (n+1);
00457     }
00458   }
00459   G4ErrorMatrixIter mi = m.begin();
00460   G4ErrorMatrixIter mii = m.begin();
00461   for (G4int i=1;i<n;i++)
00462   {
00463     G4int ni = n - i;
00464     G4ErrorMatrixIter mij = mi;
00465     G4int j;
00466     for (j=1; j<=i;j++)
00467     {
00468       s33 = *mij;
00469       G4ErrorMatrixIter mikj = mi + n + j - 1;
00470       G4ErrorMatrixIter miik = mii + 1;
00471       G4ErrorMatrixIter min_end = mi + n;
00472       for (;miik<min_end;)
00473       {
00474         s33 += (*mikj) * (*(miik++));
00475         mikj += n;
00476       }
00477       *(mij++) = s33;
00478     }
00479     for (j=1;j<=ni;j++)
00480     {
00481       s34 = 0.0;
00482       G4ErrorMatrixIter miik = mii + j;
00483       G4ErrorMatrixIter mikij = mii + j * n + j;
00484       for (G4int k=j;k<=ni;k++)
00485       {
00486         s34 += *mikij * (*(miik++));
00487         mikij += n;
00488       }
00489       *(mii+j) = s34;
00490     }
00491     mi += n;
00492     mii += (n+1);
00493   }
00494   G4int nxch = ir[n];
00495   if (nxch==0) return 0;
00496   for (G4int mq=1;mq<=nxch;mq++)
00497   {
00498     G4int k = nxch - mq + 1;
00499     G4int ij = ir[k];
00500     G4int i = ij >> 12;
00501     G4int j = ij%4096;
00502     G4ErrorMatrixIter mki = m.begin() + i - 1;
00503     G4ErrorMatrixIter mkj = m.begin() + j - 1;
00504     for (k=1; k<=n;k++)
00505     {
00506       // 2/24/05 David Sachs fix of improper swap bug that was present
00507       // for many years:
00508       G4double ti = *mki; // 2/24/05
00509       *mki = *mkj;
00510       *mkj = ti;        // 2/24/05
00511       mki += n;
00512       mkj += n;
00513     }
00514   }
00515   return 0;
00516 }
00517 
00518 G4int G4ErrorMatrix::dfact_matrix(G4double &det, G4int *ir)
00519 {
00520   if (ncol!=nrow)
00521      error("dfact_matrix: G4ErrorMatrix is not NxN");
00522 
00523   G4int ifail, jfail;
00524   G4int n = ncol;
00525 
00526   G4double tf;
00527   G4double g1 = 1.0e-19, g2 = 1.0e19;
00528 
00529   G4double p, q, t;
00530   G4double s11, s12;
00531 
00532   G4double epsilon = 8*DBL_EPSILON;
00533   // could be set to zero (like it was before)
00534   // but then the algorithm often doesn't detect
00535   // that a matrix is singular
00536 
00537   G4int normal = 0, imposs = -1;
00538   G4int jrange = 0, jover = 1, junder = -1;
00539   ifail = normal;
00540   jfail = jrange;
00541   G4int nxch = 0;
00542   det = 1.0;
00543   G4ErrorMatrixIter mj = m.begin();
00544   G4ErrorMatrixIter mjj = mj;
00545   for (G4int j=1;j<=n;j++)
00546   {
00547     G4int k = j;
00548     p = (std::fabs(*mjj));
00549     if (j!=n) {
00550       G4ErrorMatrixIter mij = mj + n + j - 1; 
00551       for (G4int i=j+1;i<=n;i++)
00552       {
00553         q = (std::fabs(*(mij)));
00554         if (q > p)
00555         {
00556           k = i;
00557           p = q;
00558         }
00559         mij += n;
00560       }
00561       if (k==j)
00562       {
00563         if (p <= epsilon)
00564         {
00565           det = 0;
00566           ifail = imposs;
00567           jfail = jrange;
00568           return ifail;
00569         }
00570         det = -det; // in this case the sign of the determinant
00571                     // must not change. So I change it twice. 
00572       }
00573       G4ErrorMatrixIter mjl = mj;
00574       G4ErrorMatrixIter mkl = m.begin() + (k-1)*n;
00575       for (G4int l=1;l<=n;l++)
00576       {
00577         tf = *mjl;
00578         *(mjl++) = *mkl;
00579         *(mkl++) = tf;
00580       }
00581       nxch = nxch + 1;  // this makes the determinant change its sign
00582       ir[nxch] = (((j)<<12)+(k));
00583     }
00584     else
00585     {
00586       if (p <= epsilon)
00587       {
00588         det = 0.0;
00589         ifail = imposs;
00590         jfail = jrange;
00591         return ifail;
00592       }
00593     }
00594     det *= *mjj;
00595     *mjj = 1.0 / *mjj;
00596     t = (std::fabs(det));
00597     if (t < g1)
00598     {
00599       det = 0.0;
00600       if (jfail == jrange) jfail = junder;
00601     }
00602     else if (t > g2)
00603     {
00604       det = 1.0;
00605       if (jfail==jrange) jfail = jover;
00606     }
00607     if (j!=n)
00608     {
00609       G4ErrorMatrixIter mk = mj + n;
00610       G4ErrorMatrixIter mkjp = mk + j;
00611       G4ErrorMatrixIter mjk = mj + j;
00612       for (k=j+1;k<=n;k++)
00613       {
00614         s11 = - (*mjk);
00615         s12 = - (*mkjp);
00616         if (j!=1)
00617         {
00618           G4ErrorMatrixIter mik = m.begin() + k - 1;
00619           G4ErrorMatrixIter mijp = m.begin() + j;
00620           G4ErrorMatrixIter mki = mk;
00621           G4ErrorMatrixIter mji = mj;
00622           for (G4int i=1;i<j;i++)
00623           {
00624             s11 += (*mik) * (*(mji++));
00625             s12 += (*mijp) * (*(mki++));
00626             mik += n;
00627             mijp += n;
00628           }
00629         }
00630         *(mjk++) = -s11 * (*mjj);
00631         *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
00632         mk += n;
00633         mkjp += n;
00634       }
00635     }
00636     mj += n;
00637     mjj += (n+1);
00638   }
00639   if (nxch%2==1) det = -det;
00640   if (jfail !=jrange) det = 0.0;
00641   ir[n] = nxch;
00642   return 0;
00643 }
00644 
00645 void G4ErrorMatrix::invert(G4int &ierr)
00646 {
00647   if(ncol != nrow)
00648      { error("G4ErrorMatrix::invert: G4ErrorMatrix is not NxN"); }
00649 
00650   static G4int max_array = 20;
00651   static G4int *ir = new G4int [max_array+1];
00652 
00653   if (ncol > max_array)
00654   {
00655     delete [] ir;
00656     max_array = nrow;
00657     ir = new G4int [max_array+1];
00658   }
00659   G4double t1, t2, t3;
00660   G4double det, temp, ss;
00661   G4int ifail;
00662   switch(nrow)
00663   {
00664   case 3:
00665     G4double c11,c12,c13,c21,c22,c23,c31,c32,c33;
00666     ifail = 0;
00667     c11 = (*(m.begin()+4)) * (*(m.begin()+8))
00668         - (*(m.begin()+5)) * (*(m.begin()+7));
00669     c12 = (*(m.begin()+5)) * (*(m.begin()+6))
00670         - (*(m.begin()+3)) * (*(m.begin()+8));
00671     c13 = (*(m.begin()+3)) * (*(m.begin()+7))
00672         - (*(m.begin()+4)) * (*(m.begin()+6));
00673     c21 = (*(m.begin()+7)) * (*(m.begin()+2))
00674         - (*(m.begin()+8)) * (*(m.begin()+1));
00675     c22 = (*(m.begin()+8)) * (*m.begin())
00676         - (*(m.begin()+6)) * (*(m.begin()+2));
00677     c23 = (*(m.begin()+6)) * (*(m.begin()+1))
00678         - (*(m.begin()+7)) * (*m.begin());
00679     c31 = (*(m.begin()+1)) * (*(m.begin()+5))
00680         - (*(m.begin()+2)) * (*(m.begin()+4));
00681     c32 = (*(m.begin()+2)) * (*(m.begin()+3))
00682         - (*m.begin()) * (*(m.begin()+5));
00683     c33 = (*m.begin()) * (*(m.begin()+4))
00684         - (*(m.begin()+1)) * (*(m.begin()+3));
00685     t1 = std::fabs(*m.begin());
00686     t2 = std::fabs(*(m.begin()+3));
00687     t3 = std::fabs(*(m.begin()+6));
00688     if (t1 >= t2)
00689     {
00690       if (t3 >= t1)
00691       {
00692         temp = *(m.begin()+6);
00693         det = c23*c12-c22*c13;
00694       }
00695       else
00696       {
00697         temp = *(m.begin());
00698         det = c22*c33-c23*c32;
00699       }
00700     }
00701     else if (t3 >= t2)
00702     {
00703       temp = *(m.begin()+6);
00704       det = c23*c12-c22*c13;
00705     }
00706     else
00707     {
00708       temp = *(m.begin()+3);
00709       det = c13*c32-c12*c33;
00710     }
00711     if (det==0)
00712     {
00713       ierr = 1;
00714       return;
00715     }
00716     {
00717       ss = temp/det;
00718       G4ErrorMatrixIter mq = m.begin();
00719       *(mq++) = ss*c11;
00720       *(mq++) = ss*c21;
00721       *(mq++) = ss*c31;
00722       *(mq++) = ss*c12;
00723       *(mq++) = ss*c22;
00724       *(mq++) = ss*c32;
00725       *(mq++) = ss*c13;
00726       *(mq++) = ss*c23;
00727       *(mq) = ss*c33;
00728     }
00729     break;
00730   case 2:
00731     ifail = 0;
00732     det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
00733     if (det==0)
00734     {
00735       ierr = 1;
00736       return;
00737     }
00738     ss = 1.0/det;
00739     temp = ss*(*(m.begin()+3));
00740     *(m.begin()+1) *= -ss;
00741     *(m.begin()+2) *= -ss;
00742     *(m.begin()+3) = ss*(*m.begin());
00743     *(m.begin()) = temp;
00744     break;
00745   case 1:
00746     ifail = 0;
00747     if ((*(m.begin()))==0)
00748     {
00749       ierr = 1;
00750       return;
00751     }
00752     *(m.begin()) = 1.0/(*(m.begin()));
00753     break;
00754   case 4:
00755     invertHaywood4(ierr);
00756     return;
00757   case 5:
00758     invertHaywood5(ierr);
00759     return;
00760   case 6:
00761     invertHaywood6(ierr);
00762     return;
00763   default:
00764     ifail = dfact_matrix(det, ir);
00765     if(ifail) {
00766       ierr = 1;
00767       return;
00768     }
00769     dfinv_matrix(ir);
00770     break;
00771   }
00772   ierr = 0;
00773   return;
00774 }
00775 
00776 G4double G4ErrorMatrix::determinant() const
00777 {
00778   static G4int max_array = 20;
00779   static G4int *ir = new G4int [max_array+1];
00780   if(ncol != nrow)
00781     { error("G4ErrorMatrix::determinant: G4ErrorMatrix is not NxN"); }
00782   if (ncol > max_array)
00783   {
00784     delete [] ir;
00785     max_array = nrow;
00786     ir = new G4int [max_array+1];
00787   }
00788   G4double det;
00789   G4ErrorMatrix mt(*this);
00790   G4int i = mt.dfact_matrix(det, ir);
00791   if(i==0) return det;
00792   return 0;
00793 }
00794 
00795 G4double G4ErrorMatrix::trace() const
00796 {
00797   G4double t = 0.0;
00798   for (G4ErrorMatrixConstIter d = m.begin(); d < m.end(); d += (ncol+1) )
00799   {
00800     t += *d;
00801   }
00802   return t;
00803 }
00804 
00805 void G4ErrorMatrix::error(const char *msg)
00806 {
00807   G4cerr << msg << G4endl;
00808   G4cerr << "---Exiting to System." << G4endl;
00809   abort();
00810 }
00811 
00812 
00813 // Aij are indices for a 6x6 matrix.
00814 // Mij are indices for a 5x5 matrix.
00815 // Fij are indices for a 4x4 matrix.
00816 
00817 #define A00 0
00818 #define A01 1
00819 #define A02 2
00820 #define A03 3
00821 #define A04 4
00822 #define A05 5
00823 
00824 #define A10 6
00825 #define A11 7
00826 #define A12 8
00827 #define A13 9
00828 #define A14 10
00829 #define A15 11
00830 
00831 #define A20 12
00832 #define A21 13
00833 #define A22 14
00834 #define A23 15
00835 #define A24 16
00836 #define A25 17
00837 
00838 #define A30 18
00839 #define A31 19
00840 #define A32 20
00841 #define A33 21
00842 #define A34 22
00843 #define A35 23
00844 
00845 #define A40 24
00846 #define A41 25
00847 #define A42 26
00848 #define A43 27
00849 #define A44 28
00850 #define A45 29
00851 
00852 #define A50 30
00853 #define A51 31
00854 #define A52 32
00855 #define A53 33
00856 #define A54 34
00857 #define A55 35
00858 
00859 #define M00 0
00860 #define M01 1
00861 #define M02 2
00862 #define M03 3
00863 #define M04 4
00864 
00865 #define M10 5
00866 #define M11 6
00867 #define M12 7
00868 #define M13 8
00869 #define M14 9
00870 
00871 #define M20 10
00872 #define M21 11
00873 #define M22 12
00874 #define M23 13
00875 #define M24 14
00876 
00877 #define M30 15
00878 #define M31 16
00879 #define M32 17
00880 #define M33 18
00881 #define M34 19
00882 
00883 #define M40 20
00884 #define M41 21
00885 #define M42 22
00886 #define M43 23
00887 #define M44 24
00888 
00889 #define F00 0
00890 #define F01 1
00891 #define F02 2
00892 #define F03 3
00893 
00894 #define F10 4
00895 #define F11 5
00896 #define F12 6
00897 #define F13 7
00898 
00899 #define F20 8
00900 #define F21 9
00901 #define F22 10
00902 #define F23 11
00903 
00904 #define F30 12
00905 #define F31 13
00906 #define F32 14
00907 #define F33 15
00908 
00909 
00910 void G4ErrorMatrix::invertHaywood4  (G4int & ifail)
00911 {
00912   ifail = 0;
00913 
00914   // Find all NECESSARY 2x2 dets:  (18 of them)
00915 
00916   G4double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
00917   G4double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
00918   G4double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20];
00919   G4double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21];
00920   G4double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22];
00921   G4double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
00922   G4double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
00923   G4double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
00924   G4double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
00925   G4double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
00926   G4double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
00927   G4double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32];
00928   G4double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
00929   G4double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
00930   G4double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
00931   G4double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
00932   G4double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
00933   G4double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
00934 
00935   // Find all NECESSARY 3x3 dets:   (16 of them)
00936 
00937   G4double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02 
00938                                 + m[F02]*Det2_12_01;
00939   G4double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03 
00940                                 + m[F03]*Det2_12_01;
00941   G4double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03 
00942                                 + m[F03]*Det2_12_02;
00943   G4double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13 
00944                                 + m[F03]*Det2_12_12;
00945   G4double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02 
00946                                 + m[F02]*Det2_13_01;
00947   G4double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
00948                                 + m[F03]*Det2_13_01;
00949   G4double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
00950                                 + m[F03]*Det2_13_02;
00951   G4double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
00952                                 + m[F03]*Det2_13_12;
00953   G4double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02 
00954                                 + m[F02]*Det2_23_01;
00955   G4double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
00956                                 + m[F03]*Det2_23_01;
00957   G4double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
00958                                 + m[F03]*Det2_23_02;
00959   G4double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
00960                                 + m[F03]*Det2_23_12;
00961   G4double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02 
00962                                 + m[F12]*Det2_23_01;
00963   G4double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03 
00964                                 + m[F13]*Det2_23_01;
00965   G4double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03 
00966                                 + m[F13]*Det2_23_02;
00967   G4double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13 
00968                                 + m[F13]*Det2_23_12;
00969 
00970   // Find the 4x4 det:
00971 
00972   G4double det =    m[F00]*Det3_123_123 
00973                 - m[F01]*Det3_123_023 
00974                 + m[F02]*Det3_123_013 
00975                 - m[F03]*Det3_123_012;
00976 
00977   if ( det == 0 )
00978   {  
00979     ifail = 1;
00980     return;
00981   } 
00982 
00983   G4double oneOverDet = 1.0/det;
00984   G4double mn1OverDet = - oneOverDet;
00985 
00986   m[F00] =  Det3_123_123 * oneOverDet;
00987   m[F01] =  Det3_023_123 * mn1OverDet;
00988   m[F02] =  Det3_013_123 * oneOverDet;
00989   m[F03] =  Det3_012_123 * mn1OverDet;
00990 
00991   m[F10] =  Det3_123_023 * mn1OverDet;
00992   m[F11] =  Det3_023_023 * oneOverDet;
00993   m[F12] =  Det3_013_023 * mn1OverDet;
00994   m[F13] =  Det3_012_023 * oneOverDet;
00995 
00996   m[F20] =  Det3_123_013 * oneOverDet;
00997   m[F21] =  Det3_023_013 * mn1OverDet;
00998   m[F22] =  Det3_013_013 * oneOverDet;
00999   m[F23] =  Det3_012_013 * mn1OverDet;
01000 
01001   m[F30] =  Det3_123_012 * mn1OverDet;
01002   m[F31] =  Det3_023_012 * oneOverDet;
01003   m[F32] =  Det3_013_012 * mn1OverDet;
01004   m[F33] =  Det3_012_012 * oneOverDet;
01005 
01006   return;
01007 }
01008 
01009 void G4ErrorMatrix::invertHaywood5  (G4int & ifail)
01010 {
01011   ifail = 0;
01012 
01013   // Find all NECESSARY 2x2 dets:  (30 of them)
01014 
01015   G4double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
01016   G4double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
01017   G4double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
01018   G4double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
01019   G4double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31];
01020   G4double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
01021   G4double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31];
01022   G4double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
01023   G4double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32];
01024   G4double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33];
01025   G4double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
01026   G4double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
01027   G4double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
01028   G4double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
01029   G4double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
01030   G4double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
01031   G4double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
01032   G4double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
01033   G4double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
01034   G4double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43];
01035   G4double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
01036   G4double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
01037   G4double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
01038   G4double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
01039   G4double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
01040   G4double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
01041   G4double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
01042   G4double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
01043   G4double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
01044   G4double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
01045 
01046   // Find all NECESSARY 3x3 dets:   (40 of them)
01047 
01048   G4double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02 
01049                                 + m[M12]*Det2_23_01;
01050   G4double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03 
01051                                 + m[M13]*Det2_23_01;
01052   G4double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04 
01053                                 + m[M14]*Det2_23_01;
01054   G4double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03 
01055                                 + m[M13]*Det2_23_02;
01056   G4double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04 
01057                                 + m[M14]*Det2_23_02;
01058   G4double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04 
01059                                 + m[M14]*Det2_23_03;
01060   G4double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13 
01061                                 + m[M13]*Det2_23_12;
01062   G4double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14 
01063                                 + m[M14]*Det2_23_12;
01064   G4double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14 
01065                                 + m[M14]*Det2_23_13;
01066   G4double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24 
01067                                 + m[M14]*Det2_23_23;
01068   G4double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02 
01069                                 + m[M12]*Det2_24_01;
01070   G4double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03 
01071                                 + m[M13]*Det2_24_01;
01072   G4double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04 
01073                                 + m[M14]*Det2_24_01;
01074   G4double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03 
01075                                 + m[M13]*Det2_24_02;
01076   G4double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04 
01077                                 + m[M14]*Det2_24_02;
01078   G4double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04 
01079                                 + m[M14]*Det2_24_03;
01080   G4double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13 
01081                                 + m[M13]*Det2_24_12;
01082   G4double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14 
01083                                 + m[M14]*Det2_24_12;
01084   G4double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14 
01085                                 + m[M14]*Det2_24_13;
01086   G4double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24 
01087                                 + m[M14]*Det2_24_23;
01088   G4double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02 
01089                                 + m[M12]*Det2_34_01;
01090   G4double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03 
01091                                 + m[M13]*Det2_34_01;
01092   G4double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04 
01093                                 + m[M14]*Det2_34_01;
01094   G4double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03 
01095                                 + m[M13]*Det2_34_02;
01096   G4double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04 
01097                                 + m[M14]*Det2_34_02;
01098   G4double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04 
01099                                 + m[M14]*Det2_34_03;
01100   G4double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13 
01101                                 + m[M13]*Det2_34_12;
01102   G4double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14 
01103                                 + m[M14]*Det2_34_12;
01104   G4double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14 
01105                                 + m[M14]*Det2_34_13;
01106   G4double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24 
01107                                 + m[M14]*Det2_34_23;
01108   G4double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02 
01109                                 + m[M22]*Det2_34_01;
01110   G4double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03 
01111                                 + m[M23]*Det2_34_01;
01112   G4double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04 
01113                                 + m[M24]*Det2_34_01;
01114   G4double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03 
01115                                 + m[M23]*Det2_34_02;
01116   G4double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04 
01117                                 + m[M24]*Det2_34_02;
01118   G4double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04 
01119                                 + m[M24]*Det2_34_03;
01120   G4double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13 
01121                                 + m[M23]*Det2_34_12;
01122   G4double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14 
01123                                 + m[M24]*Det2_34_12;
01124   G4double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14 
01125                                 + m[M24]*Det2_34_13;
01126   G4double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24 
01127                                 + m[M24]*Det2_34_23;
01128 
01129   // Find all NECESSARY 4x4 dets:   (25 of them)
01130 
01131   G4double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023 
01132                                 + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
01133   G4double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024 
01134                                 + m[M02]*Det3_123_014 - m[M04]*Det3_123_012;
01135   G4double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034 
01136                                 + m[M03]*Det3_123_014 - m[M04]*Det3_123_013;
01137   G4double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034 
01138                                 + m[M03]*Det3_123_024 - m[M04]*Det3_123_023;
01139   G4double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134 
01140                                 + m[M03]*Det3_123_124 - m[M04]*Det3_123_123;
01141   G4double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023 
01142                                 + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
01143   G4double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024 
01144                                 + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
01145   G4double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034 
01146                                 + m[M03]*Det3_124_014 - m[M04]*Det3_124_013;
01147   G4double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034 
01148                                 + m[M03]*Det3_124_024 - m[M04]*Det3_124_023;
01149   G4double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134 
01150                                 + m[M03]*Det3_124_124 - m[M04]*Det3_124_123;
01151   G4double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023 
01152                                 + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
01153   G4double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024 
01154                                 + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
01155   G4double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034 
01156                                 + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
01157   G4double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034 
01158                                 + m[M03]*Det3_134_024 - m[M04]*Det3_134_023;
01159   G4double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134 
01160                                 + m[M03]*Det3_134_124 - m[M04]*Det3_134_123;
01161   G4double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023 
01162                                 + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
01163   G4double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024 
01164                                 + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
01165   G4double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034 
01166                                 + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
01167   G4double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034 
01168                                 + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
01169   G4double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134 
01170                                 + m[M03]*Det3_234_124 - m[M04]*Det3_234_123;
01171   G4double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023 
01172                                 + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
01173   G4double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024 
01174                                 + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
01175   G4double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034 
01176                                 + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
01177   G4double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034 
01178                                 + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
01179   G4double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134 
01180                                 + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
01181 
01182   // Find the 5x5 det:
01183 
01184   G4double det =  m[M00]*Det4_1234_1234 
01185                 - m[M01]*Det4_1234_0234 
01186                 + m[M02]*Det4_1234_0134 
01187                 - m[M03]*Det4_1234_0124 
01188                 + m[M04]*Det4_1234_0123;
01189 
01190   if ( det == 0 )
01191   {  
01192     ifail = 1;
01193     return;
01194   } 
01195 
01196   G4double oneOverDet = 1.0/det;
01197   G4double mn1OverDet = - oneOverDet;
01198 
01199   m[M00] =  Det4_1234_1234 * oneOverDet;
01200   m[M01] =  Det4_0234_1234 * mn1OverDet;
01201   m[M02] =  Det4_0134_1234 * oneOverDet;
01202   m[M03] =  Det4_0124_1234 * mn1OverDet;
01203   m[M04] =  Det4_0123_1234 * oneOverDet;
01204 
01205   m[M10] =  Det4_1234_0234 * mn1OverDet;
01206   m[M11] =  Det4_0234_0234 * oneOverDet;
01207   m[M12] =  Det4_0134_0234 * mn1OverDet;
01208   m[M13] =  Det4_0124_0234 * oneOverDet;
01209   m[M14] =  Det4_0123_0234 * mn1OverDet;
01210 
01211   m[M20] =  Det4_1234_0134 * oneOverDet;
01212   m[M21] =  Det4_0234_0134 * mn1OverDet;
01213   m[M22] =  Det4_0134_0134 * oneOverDet;
01214   m[M23] =  Det4_0124_0134 * mn1OverDet;
01215   m[M24] =  Det4_0123_0134 * oneOverDet;
01216 
01217   m[M30] =  Det4_1234_0124 * mn1OverDet;
01218   m[M31] =  Det4_0234_0124 * oneOverDet;
01219   m[M32] =  Det4_0134_0124 * mn1OverDet;
01220   m[M33] =  Det4_0124_0124 * oneOverDet;
01221   m[M34] =  Det4_0123_0124 * mn1OverDet;
01222 
01223   m[M40] =  Det4_1234_0123 * oneOverDet;
01224   m[M41] =  Det4_0234_0123 * mn1OverDet;
01225   m[M42] =  Det4_0134_0123 * oneOverDet;
01226   m[M43] =  Det4_0124_0123 * mn1OverDet;
01227   m[M44] =  Det4_0123_0123 * oneOverDet;
01228 
01229   return;
01230 }
01231 
01232 void G4ErrorMatrix::invertHaywood6  (G4int & ifail)
01233 {
01234   ifail = 0;
01235 
01236   // Find all NECESSARY 2x2 dets:  (45 of them)
01237 
01238   G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
01239   G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
01240   G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
01241   G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
01242   G4double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40];
01243   G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
01244   G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
01245   G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
01246   G4double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41];
01247   G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
01248   G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
01249   G4double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42];
01250   G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
01251   G4double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43];
01252   G4double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44];
01253   G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
01254   G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
01255   G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
01256   G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
01257   G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
01258   G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
01259   G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
01260   G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
01261   G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
01262   G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
01263   G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
01264   G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
01265   G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
01266   G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
01267   G4double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54];
01268   G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
01269   G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
01270   G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
01271   G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
01272   G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
01273   G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
01274   G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
01275   G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
01276   G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
01277   G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
01278   G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
01279   G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
01280   G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
01281   G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
01282   G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
01283 
01284   // Find all NECESSARY 3x3 dets:  (80 of them)
01285 
01286   G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
01287                                                 + m[A22]*Det2_34_01;
01288   G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
01289                                                 + m[A23]*Det2_34_01;
01290   G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
01291                                                 + m[A24]*Det2_34_01;
01292   G4double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
01293                                                 + m[A25]*Det2_34_01;
01294   G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
01295                                                 + m[A23]*Det2_34_02;
01296   G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
01297                                                 + m[A24]*Det2_34_02;
01298   G4double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05 
01299                                                 + m[A25]*Det2_34_02;
01300   G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
01301                                                 + m[A24]*Det2_34_03;
01302   G4double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05                 
01303                                                 + m[A25]*Det2_34_03;
01304   G4double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05 
01305                                                 + m[A25]*Det2_34_04;
01306   G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
01307                                                 + m[A23]*Det2_34_12;
01308   G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
01309                                                 + m[A24]*Det2_34_12;
01310   G4double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15 
01311                                                 + m[A25]*Det2_34_12;
01312   G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
01313                                                 + m[A24]*Det2_34_13;
01314   G4double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15 
01315                                                 + m[A25]*Det2_34_13;
01316   G4double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15 
01317                                                 + m[A25]*Det2_34_14;
01318   G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
01319                                                 + m[A24]*Det2_34_23;
01320   G4double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
01321                                                 + m[A25]*Det2_34_23;
01322   G4double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
01323                                                 + m[A25]*Det2_34_24;
01324   G4double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
01325                                                 + m[A25]*Det2_34_34;
01326   G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02 
01327                                                 + m[A22]*Det2_35_01;
01328   G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03 
01329                                                 + m[A23]*Det2_35_01;
01330   G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04 
01331                                                 + m[A24]*Det2_35_01;
01332   G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05 
01333                                                 + m[A25]*Det2_35_01;
01334   G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03 
01335                                                 + m[A23]*Det2_35_02;
01336   G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04 
01337                                                 + m[A24]*Det2_35_02;
01338   G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05 
01339                                                 + m[A25]*Det2_35_02;
01340   G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04 
01341                                                 + m[A24]*Det2_35_03;
01342   G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05 
01343                                                 + m[A25]*Det2_35_03;
01344   G4double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05 
01345                                                 + m[A25]*Det2_35_04;
01346   G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13 
01347                                                 + m[A23]*Det2_35_12;
01348   G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14 
01349                                                 + m[A24]*Det2_35_12;
01350   G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15 
01351                                                 + m[A25]*Det2_35_12;
01352   G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14 
01353                                                 + m[A24]*Det2_35_13;
01354   G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15 
01355                                                 + m[A25]*Det2_35_13;
01356   G4double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15 
01357                                                 + m[A25]*Det2_35_14;
01358   G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24 
01359                                                 + m[A24]*Det2_35_23;
01360   G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25 
01361                                                 + m[A25]*Det2_35_23;
01362   G4double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25 
01363                                                 + m[A25]*Det2_35_24;
01364   G4double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35 
01365                                                 + m[A25]*Det2_35_34;
01366   G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02 
01367                                                 + m[A22]*Det2_45_01;
01368   G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03 
01369                                                 + m[A23]*Det2_45_01;
01370   G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04 
01371                                                 + m[A24]*Det2_45_01;
01372   G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05 
01373                                                 + m[A25]*Det2_45_01;
01374   G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03 
01375                                                 + m[A23]*Det2_45_02;
01376   G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04 
01377                                                 + m[A24]*Det2_45_02;
01378   G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05 
01379                                                 + m[A25]*Det2_45_02;
01380   G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04 
01381                                                 + m[A24]*Det2_45_03;
01382   G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05 
01383                                                 + m[A25]*Det2_45_03;
01384   G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05 
01385                                                 + m[A25]*Det2_45_04;
01386   G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13 
01387                                                 + m[A23]*Det2_45_12;
01388   G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14 
01389                                                 + m[A24]*Det2_45_12;
01390   G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15 
01391                                                 + m[A25]*Det2_45_12;
01392   G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14 
01393                                                 + m[A24]*Det2_45_13;
01394   G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15 
01395                                                 + m[A25]*Det2_45_13;
01396   G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15 
01397                                                 + m[A25]*Det2_45_14;
01398   G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24 
01399                                                 + m[A24]*Det2_45_23;
01400   G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25 
01401                                                 + m[A25]*Det2_45_23;
01402   G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25 
01403                                                 + m[A25]*Det2_45_24;
01404   G4double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35 
01405                                                 + m[A25]*Det2_45_34;
01406   G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02 
01407                                                 + m[A32]*Det2_45_01;
01408   G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03 
01409                                                 + m[A33]*Det2_45_01;
01410   G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04 
01411                                                 + m[A34]*Det2_45_01;
01412   G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05 
01413                                                 + m[A35]*Det2_45_01;
01414   G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03 
01415                                                 + m[A33]*Det2_45_02;
01416   G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04 
01417                                                 + m[A34]*Det2_45_02;
01418   G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05 
01419                                                 + m[A35]*Det2_45_02;
01420   G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04 
01421                                                 + m[A34]*Det2_45_03;
01422   G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05 
01423                                                 + m[A35]*Det2_45_03;
01424   G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05 
01425                                                 + m[A35]*Det2_45_04;
01426   G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13 
01427                                                 + m[A33]*Det2_45_12;
01428   G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14 
01429                                                 + m[A34]*Det2_45_12;
01430   G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15 
01431                                                 + m[A35]*Det2_45_12;
01432   G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14 
01433                                                 + m[A34]*Det2_45_13;
01434   G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15 
01435                                                 + m[A35]*Det2_45_13;
01436   G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15 
01437                                                 + m[A35]*Det2_45_14;
01438   G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24 
01439                                                 + m[A34]*Det2_45_23;
01440   G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25 
01441                                                 + m[A35]*Det2_45_23;
01442   G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25 
01443                                                 + m[A35]*Det2_45_24;
01444   G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35 
01445                                                 + m[A35]*Det2_45_34;
01446 
01447   // Find all NECESSARY 4x4 dets:  (75 of them)
01448 
01449   G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
01450                         + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
01451   G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
01452                         + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
01453   G4double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025 
01454                         + m[A12]*Det3_234_015 - m[A15]*Det3_234_012;
01455   G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
01456                         + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
01457   G4double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
01458                         + m[A13]*Det3_234_015 - m[A15]*Det3_234_013;
01459   G4double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
01460                         + m[A14]*Det3_234_015 - m[A15]*Det3_234_014;
01461   G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
01462                         + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
01463   G4double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035 
01464                         + m[A13]*Det3_234_025 - m[A15]*Det3_234_023;
01465   G4double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045 
01466                         + m[A14]*Det3_234_025 - m[A15]*Det3_234_024;
01467   G4double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045 
01468                         + m[A14]*Det3_234_035 - m[A15]*Det3_234_034;
01469   G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
01470                         + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
01471   G4double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135 
01472                         + m[A13]*Det3_234_125 - m[A15]*Det3_234_123;
01473   G4double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145 
01474                         + m[A14]*Det3_234_125 - m[A15]*Det3_234_124;
01475   G4double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145 
01476                         + m[A14]*Det3_234_135 - m[A15]*Det3_234_134;
01477   G4double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245 
01478                         + m[A14]*Det3_234_235 - m[A15]*Det3_234_234;
01479   G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023 
01480                         + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
01481   G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024 
01482                         + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
01483   G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025 
01484                         + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
01485   G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034 
01486                         + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
01487   G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035 
01488                         + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
01489   G4double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045 
01490                         + m[A14]*Det3_235_015 - m[A15]*Det3_235_014;
01491   G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034 
01492                         + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
01493   G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035 
01494                         + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
01495   G4double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045 
01496                         + m[A14]*Det3_235_025 - m[A15]*Det3_235_024;
01497   G4double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045 
01498                         + m[A14]*Det3_235_035 - m[A15]*Det3_235_034;
01499   G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134 
01500                         + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
01501   G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135 
01502                         + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
01503   G4double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145 
01504                         + m[A14]*Det3_235_125 - m[A15]*Det3_235_124;
01505   G4double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145 
01506                         + m[A14]*Det3_235_135 - m[A15]*Det3_235_134;
01507   G4double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245 
01508                         + m[A14]*Det3_235_235 - m[A15]*Det3_235_234;
01509   G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023 
01510                         + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
01511   G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024 
01512                         + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
01513   G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025 
01514                         + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
01515   G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034 
01516                         + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
01517   G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035 
01518                         + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
01519   G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045 
01520                         + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
01521   G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034 
01522                         + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
01523   G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035 
01524                         + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
01525   G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045 
01526                         + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
01527   G4double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045 
01528                         + m[A14]*Det3_245_035 - m[A15]*Det3_245_034;
01529   G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134 
01530                         + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
01531   G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135 
01532                         + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
01533   G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145 
01534                         + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
01535   G4double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145 
01536                         + m[A14]*Det3_245_135 - m[A15]*Det3_245_134;
01537   G4double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245 
01538                         + m[A14]*Det3_245_235 - m[A15]*Det3_245_234;
01539   G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023 
01540                         + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
01541   G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024 
01542                         + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
01543   G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025 
01544                         + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
01545   G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034 
01546                         + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
01547   G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035 
01548                         + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
01549   G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045 
01550                         + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
01551   G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034 
01552                         + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
01553   G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035 
01554                         + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
01555   G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045 
01556                         + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
01557   G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045 
01558                         + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
01559   G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134 
01560                         + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
01561   G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135 
01562                         + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
01563   G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145 
01564                         + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
01565   G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145 
01566                         + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
01567   G4double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245 
01568                         + m[A14]*Det3_345_235 - m[A15]*Det3_345_234;
01569   G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023 
01570                         + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
01571   G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024 
01572                         + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
01573   G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025 
01574                         + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
01575   G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034 
01576                         + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
01577   G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035 
01578                         + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
01579   G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045 
01580                         + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
01581   G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034 
01582                         + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
01583   G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035 
01584                         + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
01585   G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045 
01586                         + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
01587   G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045 
01588                         + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
01589   G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134 
01590                         + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
01591   G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135 
01592                         + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
01593   G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145 
01594                         + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
01595   G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145 
01596                         + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
01597   G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245 
01598                         + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
01599 
01600   // Find all NECESSARY 5x5 dets:  (36 of them)
01601 
01602   G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234 
01603     + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
01604   G4double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
01605     + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123;
01606   G4double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
01607     + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124;
01608   G4double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
01609     + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134;
01610   G4double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
01611     + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234;
01612   G4double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
01613     + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234;
01614   G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234 
01615     + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
01616   G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235 
01617     + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
01618   G4double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245 
01619     + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124;
01620   G4double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345 
01621     + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134;
01622   G4double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345 
01623     + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234;
01624   G4double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345 
01625     + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234;
01626   G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234 
01627     + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
01628   G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235 
01629     + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
01630   G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245 
01631     + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
01632   G4double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345 
01633     + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134;
01634   G4double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345 
01635     + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234;
01636   G4double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345 
01637     + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234;
01638   G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234 
01639     + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
01640   G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235 
01641     + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
01642   G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245 
01643     + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
01644   G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345 
01645     + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
01646   G4double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345 
01647     + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234;
01648   G4double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345 
01649     + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234; //
01650   G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234 
01651     + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
01652   G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235 
01653     + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
01654   G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245 
01655     + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
01656   G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345 
01657     + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
01658   G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345 
01659     + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
01660   G4double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345 
01661     + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234;
01662   G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234 
01663     + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
01664   G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235 
01665     + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
01666   G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245 
01667     + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
01668   G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345 
01669     + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
01670   G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345 
01671     + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
01672   G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345 
01673     + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
01674 
01675   // Find the determinant 
01676 
01677   G4double det =  m[A00]*Det5_12345_12345 
01678                 - m[A01]*Det5_12345_02345 
01679                 + m[A02]*Det5_12345_01345 
01680                 - m[A03]*Det5_12345_01245 
01681                 + m[A04]*Det5_12345_01235 
01682                 - m[A05]*Det5_12345_01234;
01683 
01684   if ( det == 0 )
01685   {  
01686     ifail = 1;
01687     return;
01688   } 
01689 
01690   G4double oneOverDet = 1.0/det;
01691   G4double mn1OverDet = - oneOverDet;
01692 
01693   m[A00] =  Det5_12345_12345*oneOverDet;
01694   m[A01] =  Det5_02345_12345*mn1OverDet;
01695   m[A02] =  Det5_01345_12345*oneOverDet;
01696   m[A03] =  Det5_01245_12345*mn1OverDet;
01697   m[A04] =  Det5_01235_12345*oneOverDet;
01698   m[A05] =  Det5_01234_12345*mn1OverDet;
01699 
01700   m[A10] =  Det5_12345_02345*mn1OverDet;
01701   m[A11] =  Det5_02345_02345*oneOverDet;
01702   m[A12] =  Det5_01345_02345*mn1OverDet;
01703   m[A13] =  Det5_01245_02345*oneOverDet;
01704   m[A14] =  Det5_01235_02345*mn1OverDet;
01705   m[A15] =  Det5_01234_02345*oneOverDet;
01706 
01707   m[A20] =  Det5_12345_01345*oneOverDet;
01708   m[A21] =  Det5_02345_01345*mn1OverDet;
01709   m[A22] =  Det5_01345_01345*oneOverDet;
01710   m[A23] =  Det5_01245_01345*mn1OverDet;
01711   m[A24] =  Det5_01235_01345*oneOverDet;
01712   m[A25] =  Det5_01234_01345*mn1OverDet;
01713 
01714   m[A30] =  Det5_12345_01245*mn1OverDet;
01715   m[A31] =  Det5_02345_01245*oneOverDet;
01716   m[A32] =  Det5_01345_01245*mn1OverDet;
01717   m[A33] =  Det5_01245_01245*oneOverDet;
01718   m[A34] =  Det5_01235_01245*mn1OverDet;
01719   m[A35] =  Det5_01234_01245*oneOverDet;
01720 
01721   m[A40] =  Det5_12345_01235*oneOverDet;
01722   m[A41] =  Det5_02345_01235*mn1OverDet;
01723   m[A42] =  Det5_01345_01235*oneOverDet;
01724   m[A43] =  Det5_01245_01235*mn1OverDet;
01725   m[A44] =  Det5_01235_01235*oneOverDet;
01726   m[A45] =  Det5_01234_01235*mn1OverDet;
01727 
01728   m[A50] =  Det5_12345_01234*mn1OverDet;
01729   m[A51] =  Det5_02345_01234*oneOverDet;
01730   m[A52] =  Det5_01345_01234*mn1OverDet;
01731   m[A53] =  Det5_01245_01234*oneOverDet;
01732   m[A54] =  Det5_01235_01234*mn1OverDet;
01733   m[A55] =  Det5_01234_01234*oneOverDet;
01734 
01735   return;
01736 }

Generated on Mon May 27 17:48:11 2013 for Geant4 by  doxygen 1.4.7