Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4ErrorSymMatrix.cc,v 1.3 2007-06-21 15:04:10 gunter Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
 26 //                                                 28 //
 27 // -------------------------------------------     29 // ------------------------------------------------------------
 28 //      GEANT 4 class implementation file          30 //      GEANT 4 class implementation file
 29 // -------------------------------------------     31 // ------------------------------------------------------------
 30                                                    32 
 31 #include "globals.hh"                              33 #include "globals.hh"
 32 #include <iostream>                                34 #include <iostream>
 33 #include <cmath>                                   35 #include <cmath>
 34                                                    36 
 35 #include "G4ErrorSymMatrix.hh"                     37 #include "G4ErrorSymMatrix.hh"
 36 #include "G4ErrorMatrix.hh"                        38 #include "G4ErrorMatrix.hh"
 37                                                    39 
 38 // Simple operation for all elements               40 // Simple operation for all elements
 39                                                    41 
 40 #define SIMPLE_UOP(OPER)                       <<  42 #define SIMPLE_UOP(OPER)                      \
 41   G4ErrorMatrixIter a = m.begin();             <<  43     G4ErrorMatrixIter a=m.begin();            \
 42   G4ErrorMatrixIter e = m.begin() + num_size() <<  44     G4ErrorMatrixIter e=m.begin()+num_size(); \
 43   for(; a < e; a++)                            <<  45    for(;a<e; a++) (*a) OPER t;
 44     (*a) OPER t;                               <<  46 
 45                                                <<  47 #define SIMPLE_BOP(OPER)                           \
 46 #define SIMPLE_BOP(OPER)                       <<  48     G4ErrorMatrixIter a=m.begin();                 \
 47   G4ErrorMatrixIter a      = m.begin();        <<  49     G4ErrorMatrixConstIter b=m2.m.begin();         \
 48   G4ErrorMatrixConstIter b = mat2.m.begin();   <<  50     G4ErrorMatrixConstIter e=m.begin()+num_size(); \
 49   G4ErrorMatrixConstIter e = m.begin() + num_s <<  51    for(;a<e; a++, b++) (*a) OPER (*b);
 50   for(; a < e; a++, b++)                       <<  52 
 51     (*a) OPER(*b);                             <<  53 #define SIMPLE_TOP(OPER)                                 \
 52                                                <<  54     G4ErrorMatrixConstIter a=m1.m.begin();               \
 53 #define SIMPLE_TOP(OPER)                       <<  55     G4ErrorMatrixConstIter b=m2.m.begin();               \
 54   G4ErrorMatrixConstIter a = mat1.m.begin();   <<  56     G4ErrorMatrixIter t=mret.m.begin();                  \
 55   G4ErrorMatrixConstIter b = mat2.m.begin();   <<  57     G4ErrorMatrixConstIter e=m1.m.begin()+m1.num_size(); \
 56   G4ErrorMatrixIter t      = mret.m.begin();   <<  58    for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
 57   G4ErrorMatrixConstIter e = mat1.m.begin() +  <<  59 
 58   for(; a < e; a++, b++, t++)                  <<  60 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
 59     (*t) = (*a) OPER(*b);                      <<  61    if (r1!=r2 || c1!=c2)  { \
 60                                                <<  62     G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
 61 #define CHK_DIM_2(r1, r2, c1, c2, fun)         <<  63    }
 62   if(r1 != r2 || c1 != c2)                     <<  64 
 63   {                                            <<  65 #define CHK_DIM_1(c1,r2,fun) \
 64     G4ErrorMatrix::error("Range error in Matri <<  66    if (c1!=r2) { \
 65   }                                            <<  67      G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
 66                                                <<  68    }
 67 #define CHK_DIM_1(c1, r2, fun)                 << 
 68   if(c1 != r2)                                 << 
 69   {                                            << 
 70     G4ErrorMatrix::error("Range error in Matri << 
 71   }                                            << 
 72                                                    69 
 73 // Constructors. (Default constructors are inl     70 // Constructors. (Default constructors are inlined and in .icc file)
 74                                                    71 
 75 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)        72 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)
 76   : m(p * (p + 1) / 2)                         <<  73    : m(p*(p+1)/2), nrow(p)
 77   , nrow(p)                                    << 
 78 {                                                  74 {
 79   size = nrow * (nrow + 1) / 2;                <<  75    size = nrow * (nrow+1) / 2;
 80   m.assign(size, 0);                           <<  76    m.assign(size,0);
 81 }                                                  77 }
 82                                                    78 
 83 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4     79 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4int init)
 84   : m(p * (p + 1) / 2)                         <<  80    : m(p*(p+1)/2), nrow(p)
 85   , nrow(p)                                    << 
 86 {                                                  81 {
 87   size = nrow * (nrow + 1) / 2;                <<  82    size = nrow * (nrow+1) / 2;
 88                                                    83 
 89   m.assign(size, 0);                           <<  84    m.assign(size,0);
 90   switch(init)                                 <<  85    switch(init)
 91   {                                            <<  86    {
 92     case 0:                                    <<  87    case 0:
 93       break;                                       88       break;
 94                                                <<  89       
 95     case 1:                                    <<  90    case 1:
 96     {                                          << 
 97       G4ErrorMatrixIter a = m.begin();         << 
 98       for(G4int i = 1; i <= nrow; i++)         << 
 99       {                                            91       {
100         *a = 1.0;                              <<  92          G4ErrorMatrixIter a = m.begin();
101         a += (i + 1);                          <<  93          for(G4int i=1;i<=nrow;i++)
                                                   >>  94          {
                                                   >>  95             *a = 1.0;
                                                   >>  96             a += (i+1);
                                                   >>  97          }
                                                   >>  98          break;
102       }                                            99       }
103       break;                                   << 100    default:
104     }                                          << 
105     default:                                   << 
106       G4ErrorMatrix::error("G4ErrorSymMatrix:     101       G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1.");
107   }                                            << 102    }
108 }                                                 103 }
109                                                   104 
110 //                                                105 //
111 // Destructor                                     106 // Destructor
112 //                                                107 //
113                                                   108 
114 G4ErrorSymMatrix::~G4ErrorSymMatrix() {}       << 109 G4ErrorSymMatrix::~G4ErrorSymMatrix()
                                                   >> 110 {
                                                   >> 111 }
115                                                   112 
116 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4Err << 113 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4ErrorSymMatrix &m1)
117   : m(mat1.size)                               << 114    : m(m1.size), nrow(m1.nrow), size(m1.size)
118   , nrow(mat1.nrow)                            << 
119   , size(mat1.size)                            << 
120 {                                                 115 {
121   m = mat1.m;                                  << 116    m = m1.m;
122 }                                                 117 }
123                                                   118 
124 //                                                119 //
125 //                                                120 //
126 // Sub matrix                                     121 // Sub matrix
127 //                                                122 //
128 //                                                123 //
129                                                   124 
130 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m    125 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) const
131 {                                                 126 {
132   G4ErrorSymMatrix mret(max_row - min_row + 1) << 127   G4ErrorSymMatrix mret(max_row-min_row+1);
133   if(max_row > num_row())                         128   if(max_row > num_row())
134   {                                            << 129     { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
135     G4ErrorMatrix::error("G4ErrorSymMatrix::su << 130   G4ErrorMatrixIter a = mret.m.begin();
136   }                                            << 131   G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
137   G4ErrorMatrixIter a       = mret.m.begin();  << 132   for(G4int irow=1; irow<=mret.num_row(); irow++)
138   G4ErrorMatrixConstIter b1 = m.begin() + (min << 
139   for(G4int irow = 1; irow <= mret.num_row();  << 
140   {                                               133   {
141     G4ErrorMatrixConstIter b = b1;                134     G4ErrorMatrixConstIter b = b1;
142     for(G4int icol = 1; icol <= irow; icol++)  << 135     for(G4int icol=1; icol<=irow; icol++)
143     {                                             136     {
144       *(a++) = *(b++);                            137       *(a++) = *(b++);
145     }                                             138     }
146     b1 += irow + min_row - 1;                  << 139     b1 += irow+min_row-1;
147   }                                               140   }
148   return mret;                                    141   return mret;
149 }                                                 142 }
150                                                   143 
151 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m << 144 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) 
152 {                                                 145 {
153   G4ErrorSymMatrix mret(max_row - min_row + 1) << 146   G4ErrorSymMatrix mret(max_row-min_row+1);
154   if(max_row > num_row())                         147   if(max_row > num_row())
155   {                                            << 148     { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
156     G4ErrorMatrix::error("G4ErrorSymMatrix::su << 149   G4ErrorMatrixIter a = mret.m.begin();
157   }                                            << 150   G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
158   G4ErrorMatrixIter a  = mret.m.begin();       << 151   for(G4int irow=1; irow<=mret.num_row(); irow++)
159   G4ErrorMatrixIter b1 = m.begin() + (min_row  << 
160   for(G4int irow = 1; irow <= mret.num_row();  << 
161   {                                               152   {
162     G4ErrorMatrixIter b = b1;                     153     G4ErrorMatrixIter b = b1;
163     for(G4int icol = 1; icol <= irow; icol++)  << 154     for(G4int icol=1; icol<=irow; icol++)
164     {                                             155     {
165       *(a++) = *(b++);                            156       *(a++) = *(b++);
166     }                                             157     }
167     b1 += irow + min_row - 1;                  << 158     b1 += irow+min_row-1;
168   }                                               159   }
169   return mret;                                    160   return mret;
170 }                                                 161 }
171                                                   162 
172 void G4ErrorSymMatrix::sub(G4int row, const G4 << 163 void G4ErrorSymMatrix::sub(G4int row,const G4ErrorSymMatrix &m1)
173 {                                                 164 {
174   if(row < 1 || row + mat1.num_row() - 1 > num << 165   if(row <1 || row+m1.num_row()-1 > num_row() )
175   {                                            << 166     { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
176     G4ErrorMatrix::error("G4ErrorSymMatrix::su << 167   G4ErrorMatrixConstIter a = m1.m.begin();
177   }                                            << 168   G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
178   G4ErrorMatrixConstIter a = mat1.m.begin();   << 169   for(G4int irow=1; irow<=m1.num_row(); irow++)
179   G4ErrorMatrixIter b1     = m.begin() + (row  << 
180   for(G4int irow = 1; irow <= mat1.num_row();  << 
181   {                                               170   {
182     G4ErrorMatrixIter b = b1;                     171     G4ErrorMatrixIter b = b1;
183     for(G4int icol = 1; icol <= irow; icol++)  << 172     for(G4int icol=1; icol<=irow; icol++)
184     {                                             173     {
185       *(b++) = *(a++);                            174       *(b++) = *(a++);
186     }                                             175     }
187     b1 += irow + row - 1;                      << 176     b1 += irow+row-1;
188   }                                               177   }
189 }                                                 178 }
190                                                   179 
191 //                                                180 //
192 // Direct sum of two matricies                    181 // Direct sum of two matricies
193 //                                                182 //
194                                                   183 
195 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix&  << 184 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &m1,
196                       const G4ErrorSymMatrix&  << 185                       const G4ErrorSymMatrix &m2)
197 {                                                 186 {
198   G4ErrorSymMatrix mret(mat1.num_row() + mat2. << 187   G4ErrorSymMatrix mret(m1.num_row() + m2.num_row(), 0);
199   mret.sub(1, mat1);                           << 188   mret.sub(1,m1);
200   mret.sub(mat1.num_row() + 1, mat2);          << 189   mret.sub(m1.num_row()+1,m2);
201   return mret;                                    190   return mret;
202 }                                                 191 }
203                                                   192 
204 /* -------------------------------------------    193 /* -----------------------------------------------------------------------
205    This section contains support routines for     194    This section contains support routines for matrix.h. This section contains
206    The two argument functions +,-. They call t    195    The two argument functions +,-. They call the copy constructor and +=,-=.
207    -------------------------------------------    196    ----------------------------------------------------------------------- */
208                                                   197 
209 G4ErrorSymMatrix G4ErrorSymMatrix::operator-() << 198 G4ErrorSymMatrix G4ErrorSymMatrix::operator- () const 
210 {                                                 199 {
211   G4ErrorSymMatrix mat2(nrow);                 << 200   G4ErrorSymMatrix m2(nrow);
212   G4ErrorMatrixConstIter a = m.begin();        << 201   G4ErrorMatrixConstIter a=m.begin();
213   G4ErrorMatrixIter b      = mat2.m.begin();   << 202   G4ErrorMatrixIter b=m2.m.begin();
214   G4ErrorMatrixConstIter e = m.begin() + num_s << 203   G4ErrorMatrixConstIter e=m.begin()+num_size();
215   for(; a < e; a++, b++)                       << 204   for(;a<e; a++, b++) { (*b) = -(*a); }
216   {                                            << 205   return m2;
217     (*b) = -(*a);                              << 
218   }                                            << 
219   return mat2;                                 << 
220 }                                                 206 }
221                                                   207 
222 G4ErrorMatrix operator+(const G4ErrorMatrix& m << 208 G4ErrorMatrix operator+(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
223 {                                                 209 {
224   G4ErrorMatrix mret(mat1);                    << 210   G4ErrorMatrix mret(m1);
225   CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 211   CHK_DIM_2(m1.num_row(),m2.num_row(), m1.num_col(),m2.num_col(),+);
226   mret += mat2;                                << 212   mret += m2;
227   return mret;                                    213   return mret;
228 }                                                 214 }
229                                                   215 
230 G4ErrorMatrix operator+(const G4ErrorSymMatrix << 216 G4ErrorMatrix operator+(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
231 {                                                 217 {
232   G4ErrorMatrix mret(mat2);                    << 218   G4ErrorMatrix mret(m2);
233   CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 219   CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),+);
234   mret += mat1;                                << 220   mret += m1;
235   return mret;                                    221   return mret;
236 }                                                 222 }
237                                                   223 
238 G4ErrorSymMatrix operator+(const G4ErrorSymMat << 224 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &m1,
239                            const G4ErrorSymMat << 225                            const G4ErrorSymMatrix &m2)
240 {                                                 226 {
241   G4ErrorSymMatrix mret(mat1.nrow);            << 227   G4ErrorSymMatrix mret(m1.nrow);
242   CHK_DIM_1(mat1.nrow, mat2.nrow, +);          << 228   CHK_DIM_1(m1.nrow, m2.nrow,+);
243   SIMPLE_TOP(+)                                   229   SIMPLE_TOP(+)
244   return mret;                                    230   return mret;
245 }                                                 231 }
246                                                   232 
247 //                                                233 //
248 // operator -                                     234 // operator -
249 //                                                235 //
250                                                   236 
251 G4ErrorMatrix operator-(const G4ErrorMatrix& m << 237 G4ErrorMatrix operator-(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
252 {                                                 238 {
253   G4ErrorMatrix mret(mat1);                    << 239   G4ErrorMatrix mret(m1);
254   CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 240   CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),-);
255   mret -= mat2;                                << 241   mret -= m2;
256   return mret;                                    242   return mret;
257 }                                                 243 }
258                                                   244 
259 G4ErrorMatrix operator-(const G4ErrorSymMatrix << 245 G4ErrorMatrix operator-(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
260 {                                                 246 {
261   G4ErrorMatrix mret(mat1);                    << 247   G4ErrorMatrix mret(m1);
262   CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma << 248   CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),m2.num_col(),-);
263   mret -= mat2;                                << 249   mret -= m2;
264   return mret;                                    250   return mret;
265 }                                                 251 }
266                                                   252 
267 G4ErrorSymMatrix operator-(const G4ErrorSymMat << 253 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &m1,
268                            const G4ErrorSymMat << 254                            const G4ErrorSymMatrix &m2)
269 {                                                 255 {
270   G4ErrorSymMatrix mret(mat1.num_row());       << 256   G4ErrorSymMatrix mret(m1.num_row());
271   CHK_DIM_1(mat1.num_row(), mat2.num_row(), -) << 257   CHK_DIM_1(m1.num_row(),m2.num_row(),-);
272   SIMPLE_TOP(-)                                   258   SIMPLE_TOP(-)
273   return mret;                                    259   return mret;
274 }                                                 260 }
275                                                   261 
276 /* -------------------------------------------    262 /* -----------------------------------------------------------------------
277    This section contains support routines for     263    This section contains support routines for matrix.h. This file contains
278    The two argument functions *,/. They call c    264    The two argument functions *,/. They call copy constructor and then /=,*=.
279    -------------------------------------------    265    ----------------------------------------------------------------------- */
280                                                   266 
281 G4ErrorSymMatrix operator/(const G4ErrorSymMat << 267 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &m1,G4double t)
282 {                                                 268 {
283   G4ErrorSymMatrix mret(mat1);                 << 269   G4ErrorSymMatrix mret(m1);
284   mret /= t;                                      270   mret /= t;
285   return mret;                                    271   return mret;
286 }                                                 272 }
287                                                   273 
288 G4ErrorSymMatrix operator*(const G4ErrorSymMat << 274 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &m1,G4double t)
289 {                                                 275 {
290   G4ErrorSymMatrix mret(mat1);                 << 276   G4ErrorSymMatrix mret(m1);
291   mret *= t;                                      277   mret *= t;
292   return mret;                                    278   return mret;
293 }                                                 279 }
294                                                   280 
295 G4ErrorSymMatrix operator*(G4double t, const G << 281 G4ErrorSymMatrix operator*(G4double t,const G4ErrorSymMatrix &m1)
296 {                                                 282 {
297   G4ErrorSymMatrix mret(mat1);                 << 283   G4ErrorSymMatrix mret(m1);
298   mret *= t;                                      284   mret *= t;
299   return mret;                                    285   return mret;
300 }                                                 286 }
301                                                   287 
302 G4ErrorMatrix operator*(const G4ErrorMatrix& m << 288 G4ErrorMatrix operator*(const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
303 {                                                 289 {
304   G4ErrorMatrix mret(mat1.num_row(), mat2.num_ << 290   G4ErrorMatrix mret(m1.num_row(),m2.num_col());
305   CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 291   CHK_DIM_1(m1.num_col(),m2.num_row(),*);
306   G4ErrorMatrixConstIter mit1, mit2, sp, snp;  << 292   G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
307   G4double temp;                                  293   G4double temp;
308   G4ErrorMatrixIter mir = mret.m.begin();      << 294   G4ErrorMatrixIter mir=mret.m.begin();
309   for(mit1 = mat1.m.begin();                   << 295   for(mit1=m1.m.begin();
310       mit1 < mat1.m.begin() + mat1.num_row() * << 296       mit1<m1.m.begin()+m1.num_row()*m1.num_col();
311   {                                            << 297       mit1 = mit2)
312     snp = mat2.m.begin();                      << 298     {
313     for(int step = 1; step <= mat2.num_row();  << 299       snp=m2.m.begin();
314     {                                          << 300       for(int step=1;step<=m2.num_row();++step)
315       mit2 = mit1;                             << 301         {
316       sp   = snp;                              << 302             mit2=mit1;
317       snp += step;                             << 303             sp=snp;
318       temp = 0;                                << 304             snp+=step;
319       while(sp < snp)  // Loop checking, 06.08 << 305             temp=0;
320       {                                        << 306             while(sp<snp)
321         temp += *(sp++) * (*(mit2++));         << 307             temp+=*(sp++)*(*(mit2++));
322       }                                        << 308             if( step<m2.num_row() ) { // only if we aren't on the last row
323       if(step < mat2.num_row())                << 309                 sp+=step-1;
324       {  // only if we aren't on the last row  << 310                 for(int stept=step+1;stept<=m2.num_row();stept++)
325         sp += step - 1;                        << 311                 {
326         for(int stept = step + 1; stept <= mat << 312                    temp+=*sp*(*(mit2++));
327         {                                      << 313                    if(stept<m2.num_row()) sp+=stept;
328           temp += *sp * (*(mit2++));           << 314                 }
329           if(stept < mat2.num_row())           << 315             } // if(step
330             sp += stept;                       << 316             *(mir++)=temp;
331         }                                      << 317         } // for(step
332       }  // if(step                            << 318     } // for(mit1
333       *(mir++) = temp;                         << 
334     }  // for(step                             << 
335   }    // for(mit1                             << 
336   return mret;                                    319   return mret;
337 }                                                 320 }
338                                                   321 
339 G4ErrorMatrix operator*(const G4ErrorSymMatrix << 322 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
340 {                                                 323 {
341   G4ErrorMatrix mret(mat1.num_row(), mat2.num_ << 324   G4ErrorMatrix mret(m1.num_row(),m2.num_col());
342   CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 325   CHK_DIM_1(m1.num_col(),m2.num_row(),*);
343   G4int step, stept;                           << 326   G4int step,stept;
344   G4ErrorMatrixConstIter mit1, mit2, sp, snp;  << 327   G4ErrorMatrixConstIter mit1,mit2,sp,snp;
345   G4double temp;                                  328   G4double temp;
346   G4ErrorMatrixIter mir = mret.m.begin();      << 329   G4ErrorMatrixIter mir=mret.m.begin();
347   for(step = 1, snp = mat1.m.begin(); step <=  << 330   for(step=1,snp=m1.m.begin();step<=m1.num_row();snp+=step++)
348   {                                               331   {
349     for(mit1 = mat2.m.begin(); mit1 < mat2.m.b << 332     for(mit1=m2.m.begin();mit1<m2.m.begin()+m2.num_col();mit1++)
350     {                                             333     {
351       mit2 = mit1;                             << 334       mit2=mit1;
352       sp   = snp;                              << 335       sp=snp;
353       temp = 0;                                << 336       temp=0;
354       while(sp < snp + step)  // Loop checking << 337       while(sp<snp+step)
355       {                                           338       {
356         temp += *mit2 * (*(sp++));             << 339         temp+=*mit2*(*(sp++));
357         mit2 += mat2.num_col();                << 340         mit2+=m2.num_col();
358       }                                           341       }
359       sp += step - 1;                          << 342       sp+=step-1;
360       for(stept = step + 1; stept <= mat1.num_ << 343       for(stept=step+1;stept<=m1.num_row();stept++)
361       {                                           344       {
362         temp += *mit2 * (*sp);                 << 345         temp+=*mit2*(*sp);
363         mit2 += mat2.num_col();                << 346         mit2+=m2.num_col();
364         sp += stept;                           << 347         sp+=stept;
365       }                                           348       }
366       *(mir++) = temp;                         << 349       *(mir++)=temp;
367     }                                             350     }
368   }                                               351   }
369   return mret;                                    352   return mret;
370 }                                                 353 }
371                                                   354 
372 G4ErrorMatrix operator*(const G4ErrorSymMatrix << 355 G4ErrorMatrix operator*(const G4ErrorSymMatrix &m1, const G4ErrorSymMatrix &m2)
373                         const G4ErrorSymMatrix << 
374 {                                                 356 {
375   G4ErrorMatrix mret(mat1.num_row(), mat1.num_ << 357   G4ErrorMatrix mret(m1.num_row(),m1.num_row());
376   CHK_DIM_1(mat1.num_col(), mat2.num_row(), *) << 358   CHK_DIM_1(m1.num_col(),m2.num_row(),*);
377   G4int step1, stept1, step2, stept2;          << 359   G4int step1,stept1,step2,stept2;
378   G4ErrorMatrixConstIter snp1, sp1, snp2, sp2; << 360   G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
379   G4double temp;                                  361   G4double temp;
380   G4ErrorMatrixIter mr = mret.m.begin();          362   G4ErrorMatrixIter mr = mret.m.begin();
381   for(step1 = 1, snp1 = mat1.m.begin(); step1  << 363   for(step1=1,snp1=m1.m.begin();step1<=m1.num_row();snp1+=step1++)
382       snp1 += step1++)                         << 
383   {                                               364   {
384     for(step2 = 1, snp2 = mat2.m.begin(); step << 365     for(step2=1,snp2=m2.m.begin();step2<=m2.num_row();)
385     {                                             366     {
386       sp1 = snp1;                              << 367       sp1=snp1;
387       sp2 = snp2;                              << 368       sp2=snp2;
388       snp2 += step2;                           << 369       snp2+=step2;
389       temp = 0;                                << 370       temp=0;
390       if(step1 < step2)                        << 371       if(step1<step2)
391       {                                        << 372       {
392         while(sp1 < snp1 + step1)  // Loop che << 373         while(sp1<snp1+step1)
393         {                                      << 374           { temp+=(*(sp1++))*(*(sp2++)); }
394           temp += (*(sp1++)) * (*(sp2++));     << 375         sp1+=step1-1;
395         }                                      << 376         for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
396         sp1 += step1 - 1;                      << 377           { temp+=(*sp1)*(*(sp2++)); }
397         for(stept1 = step1 + 1; stept1 != step << 378         sp2+=step2-1;
398         {                                      << 379         for(stept2=++step2;stept2<=m2.num_row();sp1+=stept1++,sp2+=stept2++)
399           temp += (*sp1) * (*(sp2++));         << 380           { temp+=(*sp1)*(*sp2); }
400         }                                      << 
401         sp2 += step2 - 1;                      << 
402         for(stept2 = ++step2; stept2 <= mat2.n << 
403             sp1 += stept1++, sp2 += stept2++)  << 
404         {                                      << 
405           temp += (*sp1) * (*sp2);             << 
406         }                                      << 
407       }                                           381       }
408       else                                        382       else
409       {                                           383       {
410         while(sp2 < snp2)  // Loop checking, 0 << 384         while(sp2<snp2)
411         {                                      << 385           { temp+=(*(sp1++))*(*(sp2++)); }
412           temp += (*(sp1++)) * (*(sp2++));     << 386         sp2+=step2-1;
413         }                                      << 387         for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
414         sp2 += step2 - 1;                      << 388           { temp+=(*(sp1++))*(*sp2); }
415         for(stept2 = ++step2; stept2 != step1  << 389         sp1+=step1-1;
416         {                                      << 390         for(stept1=step1+1;stept1<=m1.num_row();sp1+=stept1++,sp2+=stept2++)
417           temp += (*(sp1++)) * (*sp2);         << 391           { temp+=(*sp1)*(*sp2); }
418         }                                      << 
419         sp1 += step1 - 1;                      << 
420         for(stept1 = step1 + 1; stept1 <= mat1 << 
421             sp1 += stept1++, sp2 += stept2++)  << 
422         {                                      << 
423           temp += (*sp1) * (*sp2);             << 
424         }                                      << 
425       }                                           392       }
426       *(mr++) = temp;                          << 393       *(mr++)=temp;
427     }                                             394     }
428   }                                               395   }
429   return mret;                                    396   return mret;
430 }                                                 397 }
431                                                   398 
432 /* -------------------------------------------    399 /* -----------------------------------------------------------------------
433    This section contains the assignment and in    400    This section contains the assignment and inplace operators =,+=,-=,*=,/=.
434    -------------------------------------------    401    ----------------------------------------------------------------------- */
435                                                   402 
436 G4ErrorMatrix& G4ErrorMatrix::operator+=(const << 403 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorSymMatrix &m2)
437 {                                                 404 {
438   CHK_DIM_2(num_row(), mat2.num_row(), num_col << 405   CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
439   G4int n                    = num_col();      << 406   G4int n = num_col();
440   G4ErrorMatrixConstIter sjk = mat2.m.begin(); << 407   G4ErrorMatrixConstIter sjk = m2.m.begin();
441   G4ErrorMatrixIter m1j      = m.begin();      << 408   G4ErrorMatrixIter m1j = m.begin();
442   G4ErrorMatrixIter mj       = m.begin();      << 409   G4ErrorMatrixIter mj = m.begin();
443   // j >= k                                       410   // j >= k
444   for(G4int j = 1; j <= num_row(); j++)        << 411   for(G4int j=1;j<=num_row();j++)
445   {                                               412   {
446     G4ErrorMatrixIter mjk = mj;                   413     G4ErrorMatrixIter mjk = mj;
447     G4ErrorMatrixIter mkj = m1j;                  414     G4ErrorMatrixIter mkj = m1j;
448     for(G4int k = 1; k <= j; k++)              << 415     for(G4int k=1;k<=j;k++)
449     {                                             416     {
450       *(mjk++) += *sjk;                           417       *(mjk++) += *sjk;
451       if(j != k)                               << 418       if(j!=k) *mkj += *sjk;
452         *mkj += *sjk;                          << 
453       sjk++;                                      419       sjk++;
454       mkj += n;                                   420       mkj += n;
455     }                                             421     }
456     mj += n;                                      422     mj += n;
457     m1j++;                                        423     m1j++;
458   }                                               424   }
459   return (*this);                                 425   return (*this);
460 }                                                 426 }
461                                                   427 
462 G4ErrorSymMatrix& G4ErrorSymMatrix::operator+= << 428 G4ErrorSymMatrix & G4ErrorSymMatrix::operator+=(const G4ErrorSymMatrix &m2)
463 {                                                 429 {
464   CHK_DIM_2(num_row(), mat2.num_row(), num_col << 430   CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
465   SIMPLE_BOP(+=)                                  431   SIMPLE_BOP(+=)
466   return (*this);                                 432   return (*this);
467 }                                                 433 }
468                                                   434 
469 G4ErrorMatrix& G4ErrorMatrix::operator-=(const << 435 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorSymMatrix &m2)
470 {                                                 436 {
471   CHK_DIM_2(num_row(), mat2.num_row(), num_col << 437   CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
472   G4int n                    = num_col();      << 438   G4int n = num_col();
473   G4ErrorMatrixConstIter sjk = mat2.m.begin(); << 439   G4ErrorMatrixConstIter sjk = m2.m.begin();
474   G4ErrorMatrixIter m1j      = m.begin();      << 440   G4ErrorMatrixIter m1j = m.begin();
475   G4ErrorMatrixIter mj       = m.begin();      << 441   G4ErrorMatrixIter mj = m.begin();
476   // j >= k                                       442   // j >= k
477   for(G4int j = 1; j <= num_row(); j++)        << 443   for(G4int j=1;j<=num_row();j++)
478   {                                               444   {
479     G4ErrorMatrixIter mjk = mj;                   445     G4ErrorMatrixIter mjk = mj;
480     G4ErrorMatrixIter mkj = m1j;                  446     G4ErrorMatrixIter mkj = m1j;
481     for(G4int k = 1; k <= j; k++)              << 447     for(G4int k=1;k<=j;k++)
482     {                                             448     {
483       *(mjk++) -= *sjk;                           449       *(mjk++) -= *sjk;
484       if(j != k)                               << 450       if(j!=k) *mkj -= *sjk;
485         *mkj -= *sjk;                          << 
486       sjk++;                                      451       sjk++;
487       mkj += n;                                   452       mkj += n;
488     }                                             453     }
489     mj += n;                                      454     mj += n;
490     m1j++;                                        455     m1j++;
491   }                                               456   }
492   return (*this);                                 457   return (*this);
493 }                                                 458 }
494                                                   459 
495 G4ErrorSymMatrix& G4ErrorSymMatrix::operator-= << 460 G4ErrorSymMatrix & G4ErrorSymMatrix::operator-=(const G4ErrorSymMatrix &m2)
496 {                                                 461 {
497   CHK_DIM_2(num_row(), mat2.num_row(), num_col << 462   CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
498   SIMPLE_BOP(-=)                                  463   SIMPLE_BOP(-=)
499   return (*this);                                 464   return (*this);
500 }                                                 465 }
501                                                   466 
502 G4ErrorSymMatrix& G4ErrorSymMatrix::operator/= << 467 G4ErrorSymMatrix & G4ErrorSymMatrix::operator/=(G4double t)
503 {                                                 468 {
504   SIMPLE_UOP(/=)                                  469   SIMPLE_UOP(/=)
505   return (*this);                                 470   return (*this);
506 }                                                 471 }
507                                                   472 
508 G4ErrorSymMatrix& G4ErrorSymMatrix::operator*= << 473 G4ErrorSymMatrix & G4ErrorSymMatrix::operator*=(G4double t)
509 {                                                 474 {
510   SIMPLE_UOP(*=)                                  475   SIMPLE_UOP(*=)
511   return (*this);                                 476   return (*this);
512 }                                                 477 }
513                                                   478 
514 G4ErrorMatrix& G4ErrorMatrix::operator=(const  << 479 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorSymMatrix &m1)
515 {                                                 480 {
516   if(mat1.nrow * mat1.nrow != size)            << 481    if(m1.nrow*m1.nrow != size)
517   {                                            << 482    {
518     size = mat1.nrow * mat1.nrow;              << 483       size = m1.nrow * m1.nrow;
519     m.resize(size);                            << 484       m.resize(size);
520   }                                            << 485    }
521   nrow                       = mat1.nrow;      << 486    nrow = m1.nrow;
522   ncol                       = mat1.nrow;      << 487    ncol = m1.nrow;
523   G4int n                    = ncol;           << 488    G4int n = ncol;
524   G4ErrorMatrixConstIter sjk = mat1.m.begin(); << 489    G4ErrorMatrixConstIter sjk = m1.m.begin();
525   G4ErrorMatrixIter m1j      = m.begin();      << 490    G4ErrorMatrixIter m1j = m.begin();
526   G4ErrorMatrixIter mj       = m.begin();      << 491    G4ErrorMatrixIter mj = m.begin();
527   // j >= k                                    << 492    // j >= k
528   for(G4int j = 1; j <= num_row(); j++)        << 493    for(G4int j=1;j<=num_row();j++)
529   {                                            << 494    {
530     G4ErrorMatrixIter mjk = mj;                << 495       G4ErrorMatrixIter mjk = mj;
531     G4ErrorMatrixIter mkj = m1j;               << 496       G4ErrorMatrixIter mkj = m1j;
532     for(G4int k = 1; k <= j; k++)              << 497       for(G4int k=1;k<=j;k++)
533     {                                          << 498       {
534       *(mjk++) = *sjk;                         << 499          *(mjk++) = *sjk;
535       if(j != k)                               << 500          if(j!=k) *mkj = *sjk;
536         *mkj = *sjk;                           << 501          sjk++;
537       sjk++;                                   << 502          mkj += n;
538       mkj += n;                                << 503       }
539     }                                          << 504       mj += n;
540     mj += n;                                   << 505       m1j++;
541     m1j++;                                     << 506    }
542   }                                            << 507    return (*this);
543   return (*this);                              << 508 }
544 }                                              << 509 
545                                                << 510 G4ErrorSymMatrix & G4ErrorSymMatrix::operator=(const G4ErrorSymMatrix &m1)
546 G4ErrorSymMatrix& G4ErrorSymMatrix::operator=( << 511 {
547 {                                              << 512    if(m1.nrow != nrow)
548   if(&mat1 == this)                            << 513    {
549   {                                            << 514       nrow = m1.nrow;
550     return *this;                              << 515       size = m1.size;
551   }                                            << 516       m.resize(size);
552   if(mat1.nrow != nrow)                        << 517    }
553   {                                            << 518    m = m1.m;
554     nrow = mat1.nrow;                          << 519    return (*this);
555     size = mat1.size;                          << 
556     m.resize(size);                            << 
557   }                                            << 
558   m = mat1.m;                                  << 
559   return (*this);                              << 
560 }                                                 520 }
561                                                   521 
562 // Print the Matrix.                              522 // Print the Matrix.
563                                                   523 
564 std::ostream& operator<<(std::ostream& os, con << 524 std::ostream& operator<<(std::ostream &s, const G4ErrorSymMatrix &q)
565 {                                                 525 {
566   os << G4endl;                                << 526   s << G4endl;
567                                                   527 
568   // Fixed format needs 3 extra characters for    528   // Fixed format needs 3 extra characters for field,
569   // while scientific needs 7                     529   // while scientific needs 7
570                                                   530 
571   std::size_t width;                           << 531   G4int width;
572   if(os.flags() & std::ios::fixed)             << 532   if(s.flags() & std::ios::fixed)
573   {                                               533   {
574     width = os.precision() + 3;                << 534     width = s.precision()+3;
575   }                                               535   }
576   else                                            536   else
577   {                                               537   {
578     width = os.precision() + 7;                << 538     width = s.precision()+7;
579   }                                               539   }
580   for(G4int irow = 1; irow <= q.num_row(); ++i << 540   for(G4int irow = 1; irow<= q.num_row(); irow++)
581   {                                               541   {
582     for(G4int icol = 1; icol <= q.num_col(); + << 542     for(G4int icol = 1; icol <= q.num_col(); icol++)
583     {                                             543     {
584       os.width(width);                         << 544       s.width(width);
585       os << q(irow, icol) << " ";              << 545       s << q(irow,icol) << " ";
586     }                                             546     }
587     os << G4endl;                              << 547     s << G4endl;
588   }                                               548   }
589   return os;                                   << 549   return s;
590 }                                                 550 }
591                                                   551 
592 G4ErrorSymMatrix G4ErrorSymMatrix::apply(G4dou << 552 G4ErrorSymMatrix G4ErrorSymMatrix::
593                                                << 553 apply(G4double (*f)(G4double, G4int, G4int)) const
594 {                                                 554 {
595   G4ErrorSymMatrix mret(num_row());               555   G4ErrorSymMatrix mret(num_row());
596   G4ErrorMatrixConstIter a = m.cbegin();       << 556   G4ErrorMatrixConstIter a = m.begin();
597   G4ErrorMatrixIter b      = mret.m.begin();   << 557   G4ErrorMatrixIter b = mret.m.begin();
598   for(G4int ir = 1; ir <= num_row(); ++ir)     << 558   for(G4int ir=1;ir<=num_row();ir++)
599   {                                               559   {
600     for(G4int ic = 1; ic <= ir; ++ic)          << 560     for(G4int ic=1;ic<=ir;ic++)
601     {                                             561     {
602       *(b++) = (*f)(*(a++), ir, ic);              562       *(b++) = (*f)(*(a++), ir, ic);
603     }                                             563     }
604   }                                               564   }
605   return mret;                                    565   return mret;
606 }                                                 566 }
607                                                   567 
608 void G4ErrorSymMatrix::assign(const G4ErrorMat << 568 void G4ErrorSymMatrix::assign (const G4ErrorMatrix &m1)
609 {                                                 569 {
610   if(mat1.nrow != nrow)                        << 570    if(m1.nrow != nrow)
611   {                                            << 571    {
612     nrow = mat1.nrow;                          << 572       nrow = m1.nrow;
613     size = nrow * (nrow + 1) / 2;              << 573       size = nrow * (nrow+1) / 2;
614     m.resize(size);                            << 574       m.resize(size);
615   }                                            << 575    }
616   G4ErrorMatrixConstIter a = mat1.m.begin();   << 576    G4ErrorMatrixConstIter a = m1.m.begin();
617   G4ErrorMatrixIter b      = m.begin();        << 577    G4ErrorMatrixIter b = m.begin();
618   for(G4int r = 1; r <= nrow; r++)             << 578    for(G4int r=1;r<=nrow;r++)
619   {                                            << 579    {
620     G4ErrorMatrixConstIter d = a;              << 580       G4ErrorMatrixConstIter d = a;
621     for(G4int c = 1; c <= r; c++)              << 581       for(G4int c=1;c<=r;c++)
622     {                                          << 582       {
623       *(b++) = *(d++);                         << 583          *(b++) = *(d++);
624     }                                          << 584       }
625     a += nrow;                                 << 585       a += nrow;
626   }                                            << 586    }
627 }                                                 587 }
628                                                   588 
629 G4ErrorSymMatrix G4ErrorSymMatrix::similarity( << 589 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorMatrix &m1) const
630 {                                                 590 {
631   G4ErrorSymMatrix mret(mat1.num_row());       << 591   G4ErrorSymMatrix mret(m1.num_row());
632   G4ErrorMatrix temp = mat1 * (*this);         << 592   G4ErrorMatrix temp = m1*(*this);
633                                                   593 
634   // If mat1*(*this) has correct dimensions, t << 594 // If m1*(*this) has correct dimensions, then so will the m1.T multiplication.
635   // multiplication. So there is no need to ch << 595 // So there is no need to check dimensions again.
636                                                   596 
637   G4int n                  = mat1.num_col();   << 597   G4int n = m1.num_col();
638   G4ErrorMatrixIter mr     = mret.m.begin();   << 598   G4ErrorMatrixIter mr = mret.m.begin();
639   G4ErrorMatrixIter tempr1 = temp.m.begin();      599   G4ErrorMatrixIter tempr1 = temp.m.begin();
640   for(G4int r = 1; r <= mret.num_row(); r++)   << 600   for(G4int r=1;r<=mret.num_row();r++)
641   {                                               601   {
642     G4ErrorMatrixConstIter m1c1 = mat1.m.begin << 602     G4ErrorMatrixConstIter m1c1 = m1.m.begin();
643     for(G4int c = 1; c <= r; c++)              << 603     for(G4int c=1;c<=r;c++)
644     {                                             604     {
645       G4double tmp                = 0.0;       << 605       G4double tmp = 0.0;
646       G4ErrorMatrixIter tempri    = tempr1;    << 606       G4ErrorMatrixIter tempri = tempr1;
647       G4ErrorMatrixConstIter m1ci = m1c1;         607       G4ErrorMatrixConstIter m1ci = m1c1;
648       for(G4int i = 1; i <= mat1.num_col(); i+ << 608       for(G4int i=1;i<=m1.num_col();i++)
649       {                                           609       {
650         tmp += (*(tempri++)) * (*(m1ci++));    << 610         tmp+=(*(tempri++))*(*(m1ci++));
651       }                                           611       }
652       *(mr++) = tmp;                              612       *(mr++) = tmp;
653       m1c1 += n;                                  613       m1c1 += n;
654     }                                             614     }
655     tempr1 += n;                                  615     tempr1 += n;
656   }                                               616   }
657   return mret;                                    617   return mret;
658 }                                                 618 }
659                                                   619 
660 G4ErrorSymMatrix G4ErrorSymMatrix::similarity( << 620 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorSymMatrix &m1) const
661   const G4ErrorSymMatrix& mat1) const          << 
662 {                                                 621 {
663   G4ErrorSymMatrix mret(mat1.num_row());       << 622   G4ErrorSymMatrix mret(m1.num_row());
664   G4ErrorMatrix temp       = mat1 * (*this);   << 623   G4ErrorMatrix temp = m1*(*this);
665   G4int n                  = mat1.num_col();   << 624   G4int n = m1.num_col();
666   G4ErrorMatrixIter mr     = mret.m.begin();   << 625   G4ErrorMatrixIter mr = mret.m.begin();
667   G4ErrorMatrixIter tempr1 = temp.m.begin();      626   G4ErrorMatrixIter tempr1 = temp.m.begin();
668   for(G4int r = 1; r <= mret.num_row(); r++)   << 627   for(G4int r=1;r<=mret.num_row();r++)
669   {                                               628   {
670     G4ErrorMatrixConstIter m1c1 = mat1.m.begin << 629     G4ErrorMatrixConstIter m1c1 = m1.m.begin();
671     G4int c;                                      630     G4int c;
672     for(c = 1; c <= r; c++)                    << 631     for(c=1;c<=r;c++)
673     {                                             632     {
674       G4double tmp                = 0.0;       << 633       G4double tmp = 0.0;
675       G4ErrorMatrixIter tempri    = tempr1;    << 634       G4ErrorMatrixIter tempri = tempr1;
676       G4ErrorMatrixConstIter m1ci = m1c1;         635       G4ErrorMatrixConstIter m1ci = m1c1;
677       G4int i;                                    636       G4int i;
678       for(i = 1; i < c; i++)                   << 637       for(i=1;i<c;i++)
679       {                                           638       {
680         tmp += (*(tempri++)) * (*(m1ci++));    << 639         tmp+=(*(tempri++))*(*(m1ci++));
681       }                                           640       }
682       for(i = c; i <= mat1.num_col(); i++)     << 641       for(i=c;i<=m1.num_col();i++)
683       {                                           642       {
684         tmp += (*(tempri++)) * (*(m1ci));      << 643         tmp+=(*(tempri++))*(*(m1ci));
685         m1ci += i;                                644         m1ci += i;
686       }                                           645       }
687       *(mr++) = tmp;                              646       *(mr++) = tmp;
688       m1c1 += c;                                  647       m1c1 += c;
689     }                                             648     }
690     tempr1 += n;                                  649     tempr1 += n;
691   }                                               650   }
692   return mret;                                    651   return mret;
693 }                                                 652 }
694                                                   653 
695 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT << 654 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT(const G4ErrorMatrix &m1) const
696 {                                                 655 {
697   G4ErrorSymMatrix mret(mat1.num_col());       << 656   G4ErrorSymMatrix mret(m1.num_col());
698   G4ErrorMatrix temp       = (*this) * mat1;   << 657   G4ErrorMatrix temp = (*this)*m1;
699   G4int n                  = mat1.num_col();   << 658   G4int n = m1.num_col();
700   G4ErrorMatrixIter mrc    = mret.m.begin();   << 659   G4ErrorMatrixIter mrc = mret.m.begin();
701   G4ErrorMatrixIter temp1r = temp.m.begin();      660   G4ErrorMatrixIter temp1r = temp.m.begin();
702   for(G4int r = 1; r <= mret.num_row(); r++)   << 661   for(G4int r=1;r<=mret.num_row();r++)
703   {                                               662   {
704     G4ErrorMatrixConstIter m11c = mat1.m.begin << 663     G4ErrorMatrixConstIter m11c = m1.m.begin();
705     for(G4int c = 1; c <= r; c++)              << 664     for(G4int c=1;c<=r;c++)
706     {                                             665     {
707       G4double tmp                = 0.0;       << 666       G4double tmp = 0.0;
708       G4ErrorMatrixIter tempir    = temp1r;    << 667       G4ErrorMatrixIter tempir = temp1r;
709       G4ErrorMatrixConstIter m1ic = m11c;         668       G4ErrorMatrixConstIter m1ic = m11c;
710       for(G4int i = 1; i <= mat1.num_row(); i+ << 669       for(G4int i=1;i<=m1.num_row();i++)
711       {                                           670       {
712         tmp += (*(tempir)) * (*(m1ic));        << 671         tmp+=(*(tempir))*(*(m1ic));
713         tempir += n;                              672         tempir += n;
714         m1ic += n;                                673         m1ic += n;
715       }                                           674       }
716       *(mrc++) = tmp;                             675       *(mrc++) = tmp;
717       m11c++;                                     676       m11c++;
718     }                                             677     }
719     temp1r++;                                     678     temp1r++;
720   }                                               679   }
721   return mret;                                    680   return mret;
722 }                                                 681 }
723                                                   682 
724 void G4ErrorSymMatrix::invert(G4int& ifail)    << 683 void G4ErrorSymMatrix::invert(G4int &ifail)
725 {                                              << 684 {  
726   ifail = 0;                                      685   ifail = 0;
727                                                   686 
728   switch(nrow)                                    687   switch(nrow)
729   {                                               688   {
730     case 3:                                    << 689   case 3:
731     {                                             690     {
732       G4double det, temp;                         691       G4double det, temp;
733       G4double t1, t2, t3;                        692       G4double t1, t2, t3;
734       G4double c11, c12, c13, c22, c23, c33;   << 693       G4double c11,c12,c13,c22,c23,c33;
735       c11 = (*(m.begin() + 2)) * (*(m.begin()  << 694       c11 = (*(m.begin()+2)) * (*(m.begin()+5))
736             (*(m.begin() + 4)) * (*(m.begin()  << 695           - (*(m.begin()+4)) * (*(m.begin()+4));
737       c12 = (*(m.begin() + 4)) * (*(m.begin()  << 696       c12 = (*(m.begin()+4)) * (*(m.begin()+3))
738             (*(m.begin() + 1)) * (*(m.begin()  << 697           - (*(m.begin()+1)) * (*(m.begin()+5));
739       c13 = (*(m.begin() + 1)) * (*(m.begin()  << 698       c13 = (*(m.begin()+1)) * (*(m.begin()+4))
740             (*(m.begin() + 2)) * (*(m.begin()  << 699           - (*(m.begin()+2)) * (*(m.begin()+3));
741       c22 = (*(m.begin() + 5)) * (*m.begin())  << 700       c22 = (*(m.begin()+5)) * (*m.begin())
742             (*(m.begin() + 3)) * (*(m.begin()  << 701           - (*(m.begin()+3)) * (*(m.begin()+3));
743       c23 = (*(m.begin() + 3)) * (*(m.begin()  << 702       c23 = (*(m.begin()+3)) * (*(m.begin()+1))
744             (*(m.begin() + 4)) * (*m.begin()); << 703           - (*(m.begin()+4)) * (*m.begin());
745       c33 = (*m.begin()) * (*(m.begin() + 2))  << 704       c33 = (*m.begin()) * (*(m.begin()+2))
746             (*(m.begin() + 1)) * (*(m.begin()  << 705           - (*(m.begin()+1)) * (*(m.begin()+1));
747       t1 = std::fabs(*m.begin());                 706       t1 = std::fabs(*m.begin());
748       t2 = std::fabs(*(m.begin() + 1));        << 707       t2 = std::fabs(*(m.begin()+1));
749       t3 = std::fabs(*(m.begin() + 3));        << 708       t3 = std::fabs(*(m.begin()+3));
750       if(t1 >= t2)                             << 709       if (t1 >= t2)
751       {                                           710       {
752         if(t3 >= t1)                           << 711         if (t3 >= t1)
753         {                                         712         {
754           temp = *(m.begin() + 3);             << 713           temp = *(m.begin()+3);
755           det  = c23 * c12 - c22 * c13;        << 714           det = c23*c12-c22*c13;
756         }                                         715         }
757         else                                      716         else
758         {                                         717         {
759           temp = *m.begin();                      718           temp = *m.begin();
760           det  = c22 * c33 - c23 * c23;        << 719           det = c22*c33-c23*c23;
761         }                                         720         }
762       }                                           721       }
763       else if(t3 >= t2)                        << 722       else if (t3 >= t2)
764       {                                           723       {
765         temp = *(m.begin() + 3);               << 724         temp = *(m.begin()+3);
766         det  = c23 * c12 - c22 * c13;          << 725         det = c23*c12-c22*c13;
767       }                                           726       }
768       else                                        727       else
769       {                                           728       {
770         temp = *(m.begin() + 1);               << 729         temp = *(m.begin()+1);
771         det  = c13 * c23 - c12 * c33;          << 730         det = c13*c23-c12*c33;
772       }                                           731       }
773       if(det == 0)                             << 732       if (det==0)
774       {                                           733       {
775         ifail = 1;                                734         ifail = 1;
776         return;                                   735         return;
777       }                                           736       }
778       {                                           737       {
779         G4double ss          = temp / det;     << 738         G4double s = temp/det;
780         G4ErrorMatrixIter mq = m.begin();      << 739         G4ErrorMatrixIter mm = m.begin();
781         *(mq++)              = ss * c11;       << 740         *(mm++) = s*c11;
782         *(mq++)              = ss * c12;       << 741         *(mm++) = s*c12;
783         *(mq++)              = ss * c22;       << 742         *(mm++) = s*c22;
784         *(mq++)              = ss * c13;       << 743         *(mm++) = s*c13;
785         *(mq++)              = ss * c23;       << 744         *(mm++) = s*c23;
786         *(mq)                = ss * c33;       << 745         *(mm) = s*c33;
787       }                                           746       }
788     }                                             747     }
789     break;                                        748     break;
790     case 2:                                    << 749  case 2:
791     {                                             750     {
792       G4double det, temp, ss;                  << 751       G4double det, temp, s;
793       det = (*m.begin()) * (*(m.begin() + 2))  << 752       det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
794             (*(m.begin() + 1)) * (*(m.begin()  << 753       if (det==0)
795       if(det == 0)                             << 
796       {                                           754       {
797         ifail = 1;                                755         ifail = 1;
798         return;                                   756         return;
799       }                                           757       }
800       ss = 1.0 / det;                          << 758       s = 1.0/det;
801       *(m.begin() + 1) *= -ss;                 << 759       *(m.begin()+1) *= -s;
802       temp             = ss * (*(m.begin() + 2 << 760       temp = s*(*(m.begin()+2));
803       *(m.begin() + 2) = ss * (*m.begin());    << 761       *(m.begin()+2) = s*(*m.begin());
804       *m.begin()       = temp;                 << 762       *m.begin() = temp;
805       break;                                      763       break;
806     }                                             764     }
807     case 1:                                    << 765  case 1:
808     {                                             766     {
809       if((*m.begin()) == 0)                    << 767       if ((*m.begin())==0)
810       {                                           768       {
811         ifail = 1;                                769         ifail = 1;
812         return;                                   770         return;
813       }                                           771       }
814       *m.begin() = 1.0 / (*m.begin());         << 772       *m.begin() = 1.0/(*m.begin());
815       break;                                      773       break;
816     }                                             774     }
817     case 5:                                    << 775  case 5:
818     {                                             776     {
819       invert5(ifail);                             777       invert5(ifail);
820       return;                                     778       return;
821     }                                             779     }
822     case 6:                                    << 780  case 6:
823     {                                             781     {
824       invert6(ifail);                             782       invert6(ifail);
825       return;                                     783       return;
826     }                                             784     }
827     case 4:                                    << 785  case 4:
828     {                                             786     {
829       invert4(ifail);                             787       invert4(ifail);
830       return;                                     788       return;
831     }                                             789     }
832     default:                                   << 790  default:
833     {                                             791     {
834       invertBunchKaufman(ifail);                  792       invertBunchKaufman(ifail);
835       return;                                     793       return;
836     }                                             794     }
837   }                                               795   }
838   return;  // inversion successful             << 796   return; // inversion successful
839 }                                                 797 }
840                                                   798 
841 G4double G4ErrorSymMatrix::determinant() const    799 G4double G4ErrorSymMatrix::determinant() const
842 {                                                 800 {
843   static const G4int max_array = 20;              801   static const G4int max_array = 20;
844                                                   802 
845   // ir must point to an array which is ***1 l    803   // ir must point to an array which is ***1 longer than*** nrow
846                                                   804 
847   static std::vector<G4int> ir_vec(max_array + << 805   static std::vector<G4int> ir_vec (max_array+1); 
848   if(ir_vec.size() <= static_cast<unsigned int << 806   if (ir_vec.size() <= static_cast<unsigned int>(nrow))
849   {                                               807   {
850     ir_vec.resize(nrow + 1);                   << 808     ir_vec.resize(nrow+1);
851   }                                               809   }
852   G4int* ir = &ir_vec[0];                      << 810   G4int * ir = &ir_vec[0];   
853                                                   811 
854   G4double det;                                   812   G4double det;
855   G4ErrorMatrix mt(*this);                        813   G4ErrorMatrix mt(*this);
856   G4int i = mt.dfact_matrix(det, ir);             814   G4int i = mt.dfact_matrix(det, ir);
857   if(i == 0)                                   << 815   if(i==0) { return det; }
858   {                                            << 
859     return det;                                << 
860   }                                            << 
861   return 0.0;                                     816   return 0.0;
862 }                                                 817 }
863                                                   818 
864 G4double G4ErrorSymMatrix::trace() const          819 G4double G4ErrorSymMatrix::trace() const
865 {                                                 820 {
866   G4double t = 0.0;                            << 821    G4double t = 0.0;
867   for(G4int i = 0; i < nrow; i++)              << 822    for (G4int i=0; i<nrow; i++) 
868   {                                            << 823      { t += *(m.begin() + (i+3)*i/2); }
869     t += *(m.begin() + (i + 3) * i / 2);       << 824    return t;
870   }                                            << 
871   return t;                                    << 
872 }                                                 825 }
873                                                   826 
874 void G4ErrorSymMatrix::invertBunchKaufman(G4in << 827 void G4ErrorSymMatrix::invertBunchKaufman(G4int &ifail)
875 {                                                 828 {
876   // Bunch-Kaufman diagonal pivoting method       829   // Bunch-Kaufman diagonal pivoting method
877   // It is decribed in J.R. Bunch, L. Kaufman  << 830   // It is decribed in J.R. Bunch, L. Kaufman (1977). 
878   // "Some Stable Methods for Calculating Iner << 831   // "Some Stable Methods for Calculating Inertia and Solving Symmetric 
879   // Linear Systems", Math. Comp. 31, p. 162-1 << 832   // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub, 
880   // Charles F. van Loan, "Matrix Computations << 833   // Charles F. van Loan, "Matrix Computations" (the second edition 
881   // has a bug.) and implemented in "lapack"      834   // has a bug.) and implemented in "lapack"
882   // Mario Stanke, 09/97                          835   // Mario Stanke, 09/97
883                                                   836 
884   G4int i, j, k, ss;                           << 837   G4int i, j, k, s;
885   G4int pivrow;                                   838   G4int pivrow;
886                                                   839 
887   // Establish the two working-space arrays ne    840   // Establish the two working-space arrays needed:  x and piv are
888   // used as pointers to arrays of doubles and    841   // used as pointers to arrays of doubles and ints respectively, each
889   // of length nrow.  We do not want to reallo    842   // of length nrow.  We do not want to reallocate each time through
890   // unless the size needs to grow.  We do not    843   // unless the size needs to grow.  We do not want to leak memory, even
891   // by having a new without a delete that is     844   // by having a new without a delete that is only done once.
892                                                << 845   
893   static const G4int max_array                 << 846   static const G4int max_array = 25;
894   static G4ThreadLocal std::vector<G4double>*  << 847   static std::vector<G4double> xvec (max_array);
895   if(!xvec)                                    << 848   static std::vector<G4int>    pivv (max_array);
896     xvec = new std::vector<G4double>(max_array << 849   typedef std::vector<G4int>::iterator pivIter; 
897   static G4ThreadLocal std::vector<G4int>* piv << 850   if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow);
898   if(!pivv)                                    << 851   if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow);
899     pivv = new std::vector<G4int>(max_array);  << 852     // Note - resize should do  nothing if the size is already larger than nrow,
900   typedef std::vector<G4int>::iterator pivIter << 853     //        but on VC++ there are indications that it does so we check.
901   if(xvec->size() < static_cast<unsigned int>( << 854     // Note - the data elements in a vector are guaranteed to be contiguous,
902     xvec->resize(nrow);                        << 855     //        so x[i] and piv[i] are optimally fast.
903   if(pivv->size() < static_cast<unsigned int>( << 856   G4ErrorMatrixIter   x   = xvec.begin();
904     pivv->resize(nrow);                        << 857     // x[i] is used as helper storage, needs to have at least size nrow.
905   // Note - resize should do  nothing if the s << 858   pivIter piv = pivv.begin();
906   //        but on VC++ there are indications  << 859     // piv[i] is used to store details of exchanges
907   // Note - the data elements in a vector are  << 860       
908   //        so x[i] and piv[i] are optimally f << 
909   G4ErrorMatrixIter x = xvec->begin();         << 
910   // x[i] is used as helper storage, needs to  << 
911   pivIter piv = pivv->begin();                 << 
912   // piv[i] is used to store details of exchan << 
913                                                << 
914   G4double temp1, temp2;                          861   G4double temp1, temp2;
915   G4ErrorMatrixIter ip, mjj, iq;                  862   G4ErrorMatrixIter ip, mjj, iq;
916   G4double lambda, sigma;                         863   G4double lambda, sigma;
917   const G4double alpha   = .6404;  // = (1+sqr << 864   const G4double alpha = .6404; // = (1+sqrt(17))/8
918   const G4double epsilon = 32 * DBL_EPSILON;   << 865   const G4double epsilon = 32*DBL_EPSILON;
919   // whenever a sum of two doubles is below or << 866     // whenever a sum of two doubles is below or equal to epsilon
920   // it is set to zero.                        << 867     // it is set to zero.
921   // this constant could be set to zero but th << 868     // this constant could be set to zero but then the algorithm
922   // doesn't neccessarily detect that a matrix << 869     // doesn't neccessarily detect that a matrix is singular
923                                                << 870   
924   for(i = 0; i < nrow; i++)                    << 871   for (i = 0; i < nrow; i++)
925   {                                               872   {
926     piv[i] = i + 1;                            << 873     piv[i] = i+1;
927   }                                               874   }
928                                                << 875       
929   ifail = 0;                                      876   ifail = 0;
930                                                << 877       
931   // compute the factorization P*A*P^T = L * D << 878   // compute the factorization P*A*P^T = L * D * L^T 
932   // L is unit lower triangular, D is direct s    879   // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
933   // L and D^-1 are stored in A = *this, P is     880   // L and D^-1 are stored in A = *this, P is stored in piv[]
934                                                << 881         
935   for(j = 1; j < nrow; j += ss)  // main loop  << 882   for (j=1; j < nrow; j+=s)  // main loop over columns
936   {                                               883   {
937     mjj    = m.begin() + j * (j - 1) / 2 + j - << 884           mjj = m.begin() + j*(j-1)/2 + j-1;
938     lambda = 0;  // compute lambda = max of A( << 885           lambda = 0;           // compute lambda = max of A(j+1:n,j)
939     pivrow = j + 1;                            << 886           pivrow = j+1;
940     ip     = m.begin() + (j + 1) * j / 2 + j - << 887           ip = m.begin() + (j+1)*j/2 + j-1;
941     for(i = j + 1; i <= nrow; ip += i++)       << 888           for (i=j+1; i <= nrow ; ip += i++)
942     {                                          << 
943       if(std::fabs(*ip) > lambda)              << 
944       {                                        << 
945         lambda = std::fabs(*ip);               << 
946         pivrow = i;                            << 
947       }                                        << 
948     }                                          << 
949     if(lambda == 0)                            << 
950     {                                          << 
951       if(*mjj == 0)                            << 
952       {                                        << 
953         ifail = 1;                             << 
954         return;                                << 
955       }                                        << 
956       ss   = 1;                                << 
957       *mjj = 1. / *mjj;                        << 
958     }                                          << 
959     else                                       << 
960     {                                          << 
961       if(std::fabs(*mjj) >= lambda * alpha)    << 
962       {                                        << 
963         ss     = 1;                            << 
964         pivrow = j;                            << 
965       }                                        << 
966       else                                     << 
967       {                                        << 
968         sigma = 0;  // compute sigma = max A(p << 
969         ip    = m.begin() + pivrow * (pivrow - << 
970         for(k = j; k < pivrow; k++)            << 
971         {                                      << 
972           if(std::fabs(*ip) > sigma)           << 
973             sigma = std::fabs(*ip);            << 
974           ip++;                                << 
975         }                                      << 
976         if(sigma * std::fabs(*mjj) >= alpha *  << 
977         {                                      << 
978           ss     = 1;                          << 
979           pivrow = j;                          << 
980         }                                      << 
981         else if(std::fabs(*(m.begin() + pivrow << 
982                             1)) >= alpha * sig << 
983         {                                      << 
984           ss = 1;                              << 
985         }                                      << 
986         else                                   << 
987         {                                      << 
988           ss = 2;                              << 
989         }                                      << 
990       }                                        << 
991       if(pivrow == j)  // no permutation necce << 
992       {                                        << 
993         piv[j - 1] = pivrow;                   << 
994         if(*mjj == 0)                          << 
995         {                                      << 
996           ifail = 1;                           << 
997           return;                              << 
998         }                                      << 
999         temp2 = *mjj = 1. / *mjj;  // invert D << 
1000                                               << 
1001         // update A(j+1:n, j+1,n)             << 
1002         for(i = j + 1; i <= nrow; i++)        << 
1003         {                                     << 
1004           temp1 = *(m.begin() + i * (i - 1) / << 
1005           ip    = m.begin() + i * (i - 1) / 2 << 
1006           for(k = j + 1; k <= i; k++)         << 
1007           {                                      889           {
1008             *ip -= temp1 * *(m.begin() + k *  << 890             if (std::fabs(*ip) > lambda)
1009             if(std::fabs(*ip) <= epsilon)     << 
1010             {                                 << 
1011               *ip = 0;                        << 
1012             }                                 << 
1013             ip++;                             << 
1014           }                                   << 
1015         }                                     << 
1016         // update L                           << 
1017         ip = m.begin() + (j + 1) * j / 2 + j  << 
1018         for(i = j + 1; i <= nrow; ip += i++)  << 
1019         {                                     << 
1020           *ip *= temp2;                       << 
1021         }                                     << 
1022       }                                       << 
1023       else if(ss == 1)  // 1x1 pivot          << 
1024       {                                       << 
1025         piv[j - 1] = pivrow;                  << 
1026                                               << 
1027         // interchange rows and columns j and << 
1028         // submatrix (j:n,j:n)                << 
1029         ip = m.begin() + pivrow * (pivrow - 1 << 
1030         for(i = j + 1; i < pivrow; i++, ip++) << 
1031         {                                     << 
1032           temp1 = *(m.begin() + i * (i - 1) / << 
1033           *(m.begin() + i * (i - 1) / 2 + j - << 
1034           *ip                                 << 
1035         }                                     << 
1036         temp1 = *mjj;                         << 
1037         *mjj  = *(m.begin() + pivrow * (pivro << 
1038         *(m.begin() + pivrow * (pivrow - 1) / << 
1039         ip = m.begin() + (pivrow + 1) * pivro << 
1040         iq = ip + pivrow - j;                 << 
1041         for(i = pivrow + 1; i <= nrow; ip +=  << 
1042         {                                     << 
1043           temp1 = *iq;                        << 
1044           *iq   = *ip;                        << 
1045           *ip   = temp1;                      << 
1046         }                                     << 
1047                                               << 
1048         if(*mjj == 0)                         << 
1049         {                                     << 
1050           ifail = 1;                          << 
1051           return;                             << 
1052         }                                     << 
1053         temp2 = *mjj = 1. / *mjj;  // invert  << 
1054                                               << 
1055         // update A(j+1:n, j+1:n)             << 
1056         for(i = j + 1; i <= nrow; i++)        << 
1057         {                                     << 
1058           temp1 = *(m.begin() + i * (i - 1) / << 
1059           ip    = m.begin() + i * (i - 1) / 2 << 
1060           for(k = j + 1; k <= i; k++)         << 
1061           {                                   << 
1062             *ip -= temp1 * *(m.begin() + k *  << 
1063             if(std::fabs(*ip) <= epsilon)     << 
1064             {                                 << 
1065               *ip = 0;                        << 
1066             }                                 << 
1067             ip++;                             << 
1068           }                                   << 
1069         }                                     << 
1070         // update L                           << 
1071         ip = m.begin() + (j + 1) * j / 2 + j  << 
1072         for(i = j + 1; i <= nrow; ip += i++)  << 
1073         {                                     << 
1074           *ip *= temp2;                       << 
1075         }                                     << 
1076       }                                       << 
1077       else  // ss=2, ie use a 2x2 pivot       << 
1078       {                                       << 
1079         piv[j - 1] = -pivrow;                 << 
1080         piv[j]     = 0;  // that means this i << 
1081                                               << 
1082         if(j + 1 != pivrow)                   << 
1083         {                                     << 
1084           // interchange rows and columns j+1 << 
1085           // submatrix (j:n,j:n)              << 
1086           ip = m.begin() + pivrow * (pivrow - << 
1087           for(i = j + 2; i < pivrow; i++, ip+ << 
1088           {                                   << 
1089             temp1 = *(m.begin() + i * (i - 1) << 
1090             *(m.begin() + i * (i - 1) / 2 + j << 
1091             *ip                               << 
1092           }                                   << 
1093           temp1 = *(mjj + j + 1);             << 
1094           *(mjj + j + 1) =                    << 
1095             *(m.begin() + pivrow * (pivrow -  << 
1096           *(m.begin() + pivrow * (pivrow - 1) << 
1097           temp1                               << 
1098           *(mjj + j) = *(m.begin() + pivrow * << 
1099           *(m.begin() + pivrow * (pivrow - 1) << 
1100           ip = m.begin() + (pivrow + 1) * piv << 
1101           iq = ip + pivrow - (j + 1);         << 
1102           for(i = pivrow + 1; i <= nrow; ip + << 
1103           {                                   << 
1104             temp1 = *iq;                      << 
1105             *iq   = *ip;                      << 
1106             *ip   = temp1;                    << 
1107           }                                   << 
1108         }                                     << 
1109         // invert D(j:j+1,j:j+1)              << 
1110         temp2 = *mjj * *(mjj + j + 1) - *(mjj << 
1111         if(temp2 == 0)                        << 
1112         {                                     << 
1113           G4Exception("G4ErrorSymMatrix::bunc << 
1114                       "GEANT4e-Notification", << 
1115                       "Error in pivot choice! << 
1116         }                                     << 
1117         temp2 = 1. / temp2;                   << 
1118                                               << 
1119         // this quotient is guaranteed to exi << 
1120         // of the pivot                       << 
1121                                               << 
1122         temp1          = *mjj;                << 
1123         *mjj           = *(mjj + j + 1) * tem << 
1124         *(mjj + j + 1) = temp1 * temp2;       << 
1125         *(mjj + j)     = -*(mjj + j) * temp2; << 
1126                                               << 
1127         if(j < nrow - 1)  // otherwise do not << 
1128         {                                     << 
1129           // update A(j+2:n, j+2:n)           << 
1130           for(i = j + 2; i <= nrow; i++)      << 
1131           {                                   << 
1132             ip    = m.begin() + i * (i - 1) / << 
1133             temp1 = *ip * *mjj + *(ip + 1) *  << 
1134             if(std::fabs(temp1) <= epsilon)   << 
1135             {                                 << 
1136               temp1 = 0;                      << 
1137             }                                 << 
1138             temp2 = *ip * *(mjj + j) + *(ip + << 
1139             if(std::fabs(temp2) <= epsilon)   << 
1140             {                                 << 
1141               temp2 = 0;                      << 
1142             }                                 << 
1143             for(k = j + 2; k <= i; k++)       << 
1144             {                                 << 
1145               ip = m.begin() + i * (i - 1) /  << 
1146               iq = m.begin() + k * (k - 1) /  << 
1147               *ip -= temp1 * *iq + temp2 * *( << 
1148               if(std::fabs(*ip) <= epsilon)   << 
1149               {                                  891               {
1150                 *ip = 0;                      << 892                 lambda = std::fabs(*ip);
                                                   >> 893                 pivrow = i;
1151               }                                  894               }
1152             }                                 << 
1153           }                                      895           }
1154           // update L                         << 896           if (lambda == 0 )
1155           for(i = j + 2; i <= nrow; i++)      << 
1156           {                                   << 
1157             ip    = m.begin() + i * (i - 1) / << 
1158             temp1 = *ip * *mjj + *(ip + 1) *  << 
1159             if(std::fabs(temp1) <= epsilon)   << 
1160             {                                    897             {
1161               temp1 = 0;                      << 898               if (*mjj == 0)
                                                   >> 899                 {
                                                   >> 900                   ifail = 1;
                                                   >> 901                   return;
                                                   >> 902                 }
                                                   >> 903               s=1;
                                                   >> 904               *mjj = 1./ *mjj;
1162             }                                    905             }
1163             *(ip + 1) = *ip * *(mjj + j) + *( << 906           else
1164             if(std::fabs(*(ip + 1)) <= epsilo << 
1165             {                                    907             {
1166               *(ip + 1) = 0;                  << 908               if (std::fabs(*mjj) >= lambda*alpha)
                                                   >> 909                 {
                                                   >> 910                   s=1;
                                                   >> 911                   pivrow=j;
                                                   >> 912                 }
                                                   >> 913               else
                                                   >> 914                 {
                                                   >> 915                   sigma = 0;  // compute sigma = max A(pivrow, j:pivrow-1)
                                                   >> 916                   ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
                                                   >> 917                   for (k=j; k < pivrow; k++)
                                                   >> 918                     {
                                                   >> 919                       if (std::fabs(*ip) > sigma)
                                                   >> 920                         sigma = std::fabs(*ip);
                                                   >> 921                       ip++;
                                                   >> 922                     }
                                                   >> 923                   if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
                                                   >> 924                     {
                                                   >> 925                       s=1;
                                                   >> 926                       pivrow = j;
                                                   >> 927                     }
                                                   >> 928                   else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1)) 
                                                   >> 929                                 >= alpha * sigma)
                                                   >> 930                     { s=1; }
                                                   >> 931                   else
                                                   >> 932                     { s=2; }
                                                   >> 933                 }
                                                   >> 934               if (pivrow == j)  // no permutation neccessary
                                                   >> 935                 {
                                                   >> 936                   piv[j-1] = pivrow;
                                                   >> 937                   if (*mjj == 0)
                                                   >> 938                     {
                                                   >> 939                       ifail=1;
                                                   >> 940                       return;
                                                   >> 941                     }
                                                   >> 942                   temp2 = *mjj = 1./ *mjj; // invert D(j,j)
                                                   >> 943                   
                                                   >> 944                   // update A(j+1:n, j+1,n)
                                                   >> 945                   for (i=j+1; i <= nrow; i++)
                                                   >> 946                     {
                                                   >> 947                       temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
                                                   >> 948                       ip = m.begin()+i*(i-1)/2+j;
                                                   >> 949                       for (k=j+1; k<=i; k++)
                                                   >> 950                         {
                                                   >> 951                           *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
                                                   >> 952                           if (std::fabs(*ip) <= epsilon)
                                                   >> 953                             { *ip=0; }
                                                   >> 954                           ip++;
                                                   >> 955                         }
                                                   >> 956                     }
                                                   >> 957                   // update L 
                                                   >> 958                   ip = m.begin() + (j+1)*j/2 + j-1; 
                                                   >> 959                   for (i=j+1; i <= nrow; ip += i++)
                                                   >> 960                   {
                                                   >> 961                     *ip *= temp2;
                                                   >> 962                   }
                                                   >> 963                 }
                                                   >> 964               else if (s==1) // 1x1 pivot 
                                                   >> 965                 {
                                                   >> 966                   piv[j-1] = pivrow;
                                                   >> 967                   
                                                   >> 968                   // interchange rows and columns j and pivrow in
                                                   >> 969                   // submatrix (j:n,j:n)
                                                   >> 970                   ip = m.begin() + pivrow*(pivrow-1)/2 + j;
                                                   >> 971                   for (i=j+1; i < pivrow; i++, ip++)
                                                   >> 972                     {
                                                   >> 973                       temp1 = *(m.begin() + i*(i-1)/2 + j-1);
                                                   >> 974                       *(m.begin() + i*(i-1)/2 + j-1)= *ip;
                                                   >> 975                       *ip = temp1;
                                                   >> 976                     }
                                                   >> 977                   temp1 = *mjj;
                                                   >> 978                   *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
                                                   >> 979                   *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
                                                   >> 980                   ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
                                                   >> 981                   iq = ip + pivrow-j;
                                                   >> 982                   for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
                                                   >> 983                     {
                                                   >> 984                       temp1 = *iq;
                                                   >> 985                       *iq = *ip;
                                                   >> 986                       *ip = temp1;
                                                   >> 987                     }
                                                   >> 988                   
                                                   >> 989                   if (*mjj == 0)
                                                   >> 990                     {
                                                   >> 991                       ifail = 1;
                                                   >> 992                       return;
                                                   >> 993                     }
                                                   >> 994                   temp2 = *mjj = 1./ *mjj; // invert D(j,j)
                                                   >> 995                   
                                                   >> 996                   // update A(j+1:n, j+1:n)
                                                   >> 997                   for (i = j+1; i <= nrow; i++)
                                                   >> 998                     {
                                                   >> 999                       temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
                                                   >> 1000                       ip = m.begin()+i*(i-1)/2+j;
                                                   >> 1001                       for (k=j+1; k<=i; k++)
                                                   >> 1002                         {
                                                   >> 1003                           *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
                                                   >> 1004                           if (std::fabs(*ip) <= epsilon)
                                                   >> 1005                             { *ip=0; }
                                                   >> 1006                           ip++;
                                                   >> 1007                         }
                                                   >> 1008                     }
                                                   >> 1009                   // update L
                                                   >> 1010                   ip = m.begin() + (j+1)*j/2 + j-1;
                                                   >> 1011                   for (i=j+1; i<=nrow; ip += i++)
                                                   >> 1012                   {
                                                   >> 1013                     *ip *= temp2;
                                                   >> 1014                   }
                                                   >> 1015                 }
                                                   >> 1016               else // s=2, ie use a 2x2 pivot
                                                   >> 1017                 {
                                                   >> 1018                   piv[j-1] = -pivrow;
                                                   >> 1019                   piv[j] = 0; // that means this is the second row of a 2x2 pivot
                                                   >> 1020                   
                                                   >> 1021                   if (j+1 != pivrow) 
                                                   >> 1022                     {
                                                   >> 1023                       // interchange rows and columns j+1 and pivrow in
                                                   >> 1024                       // submatrix (j:n,j:n) 
                                                   >> 1025                       ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
                                                   >> 1026                       for (i=j+2; i < pivrow; i++, ip++)
                                                   >> 1027                         {
                                                   >> 1028                           temp1 = *(m.begin() + i*(i-1)/2 + j);
                                                   >> 1029                           *(m.begin() + i*(i-1)/2 + j) = *ip;
                                                   >> 1030                           *ip = temp1;
                                                   >> 1031                         }
                                                   >> 1032                       temp1 = *(mjj + j + 1);
                                                   >> 1033                       *(mjj + j + 1) = 
                                                   >> 1034                         *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
                                                   >> 1035                       *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
                                                   >> 1036                       temp1 = *(mjj + j);
                                                   >> 1037                       *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
                                                   >> 1038                       *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
                                                   >> 1039                       ip = m.begin() + (pivrow+1)*pivrow/2 + j;
                                                   >> 1040                       iq = ip + pivrow-(j+1);
                                                   >> 1041                       for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
                                                   >> 1042                         {
                                                   >> 1043                           temp1 = *iq;
                                                   >> 1044                           *iq = *ip;
                                                   >> 1045                           *ip = temp1;
                                                   >> 1046                         }
                                                   >> 1047                     } 
                                                   >> 1048                   // invert D(j:j+1,j:j+1)
                                                   >> 1049                   temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j); 
                                                   >> 1050                   if (temp2 == 0)
                                                   >> 1051                   {
                                                   >> 1052                     G4cerr
                                                   >> 1053                       << "G4ErrorSymMatrix::bunch_invert: error in pivot choice"
                                                   >> 1054                       << G4endl;
                                                   >> 1055                   }
                                                   >> 1056                   temp2 = 1. / temp2;
                                                   >> 1057 
                                                   >> 1058                   // this quotient is guaranteed to exist by the choice 
                                                   >> 1059                   // of the pivot
                                                   >> 1060 
                                                   >> 1061                   temp1 = *mjj;
                                                   >> 1062                   *mjj = *(mjj + j + 1) * temp2;
                                                   >> 1063                   *(mjj + j + 1) = temp1 * temp2;
                                                   >> 1064                   *(mjj + j) = - *(mjj + j) * temp2;
                                                   >> 1065                   
                                                   >> 1066                   if (j < nrow-1) // otherwise do nothing
                                                   >> 1067                     {
                                                   >> 1068                       // update A(j+2:n, j+2:n)
                                                   >> 1069                       for (i=j+2; i <= nrow ; i++)
                                                   >> 1070                         {
                                                   >> 1071                           ip = m.begin() + i*(i-1)/2 + j-1;
                                                   >> 1072                           temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
                                                   >> 1073                           if (std::fabs(temp1 ) <= epsilon)
                                                   >> 1074                             { temp1 = 0; }
                                                   >> 1075                           temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
                                                   >> 1076                           if (std::fabs(temp2 ) <= epsilon)
                                                   >> 1077                             { temp2 = 0; }
                                                   >> 1078                           for (k = j+2; k <= i ; k++)
                                                   >> 1079                             {
                                                   >> 1080                               ip = m.begin() + i*(i-1)/2 + k-1;
                                                   >> 1081                               iq = m.begin() + k*(k-1)/2 + j-1;
                                                   >> 1082                               *ip -= temp1 * *iq + temp2 * *(iq+1);
                                                   >> 1083                               if (std::fabs(*ip) <= epsilon)
                                                   >> 1084                                 { *ip = 0; }
                                                   >> 1085                             }
                                                   >> 1086                         }
                                                   >> 1087                       // update L
                                                   >> 1088                       for (i=j+2; i <= nrow ; i++)
                                                   >> 1089                         {
                                                   >> 1090                           ip = m.begin() + i*(i-1)/2 + j-1;
                                                   >> 1091                           temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
                                                   >> 1092                           if (std::fabs(temp1) <= epsilon)
                                                   >> 1093                             { temp1 = 0; }
                                                   >> 1094                           *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
                                                   >> 1095                           if (std::fabs(*(ip+1)) <= epsilon)
                                                   >> 1096                             { *(ip+1) = 0; }
                                                   >> 1097                           *ip = temp1;
                                                   >> 1098                         }
                                                   >> 1099                     }
                                                   >> 1100                 }
1167             }                                    1101             }
1168             *ip = temp1;                      << 1102   } // end of main loop over columns
1169           }                                   << 
1170         }                                     << 
1171       }                                       << 
1172     }                                         << 
1173   }  // end of main loop over columns         << 
1174                                                  1103 
1175   if(j == nrow)  // the the last pivot is 1x1 << 1104   if (j == nrow) // the the last pivot is 1x1
1176   {                                              1105   {
1177     mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1106     mjj = m.begin() + j*(j-1)/2 + j-1;
1178     if(*mjj == 0)                             << 1107     if (*mjj == 0)
1179     {                                            1108     {
1180       ifail = 1;                                 1109       ifail = 1;
1181       return;                                    1110       return;
1182     }                                            1111     }
1183     else                                         1112     else
1184     {                                            1113     {
1185       *mjj = 1. / *mjj;                          1114       *mjj = 1. / *mjj;
1186     }                                            1115     }
1187   }  // end of last pivot code                << 1116   } // end of last pivot code
1188                                                  1117 
1189   // computing the inverse from the factoriza    1118   // computing the inverse from the factorization
1190                                               << 1119          
1191   for(j = nrow; j >= 1; j -= ss)  // loop ove << 1120   for (j = nrow ; j >= 1 ; j -= s) // loop over columns
1192   {                                              1121   {
1193     mjj = m.begin() + j * (j - 1) / 2 + j - 1 << 1122           mjj = m.begin() + j*(j-1)/2 + j-1;
1194     if(piv[j - 1] > 0)  // 1x1 pivot, compute << 1123           if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
1195     {                                         << 1124             {
1196       ss = 1;                                 << 1125               s = 1; 
1197       if(j < nrow)                            << 1126               if (j < nrow)
1198       {                                       << 1127                 {
1199         ip = m.begin() + (j + 1) * j / 2 + j  << 1128                   ip = m.begin() + (j+1)*j/2 + j-1;
1200         for(i = 0; i < nrow - j; ip += 1 + j  << 1129                   for (i=0; i < nrow-j; ip += 1+j+i++)
1201         {                                     << 1130                   {
1202           x[i] = *ip;                         << 1131                     x[i] = *ip;
1203         }                                     << 1132                   }
1204         for(i = j + 1; i <= nrow; i++)        << 1133                   for (i=j+1; i<=nrow ; i++)
1205         {                                     << 1134                   {
1206           temp2 = 0;                          << 1135                     temp2=0;
1207           ip    = m.begin() + i * (i - 1) / 2 << 1136                     ip = m.begin() + i*(i-1)/2 + j;
1208           for(k = 0; k <= i - j - 1; k++)     << 1137                     for (k=0; k <= i-j-1; k++)
1209           {                                   << 1138                       { temp2 += *ip++ * x[k]; }
1210             temp2 += *ip++ * x[k];            << 1139                     for (ip += i-1; k < nrow-j; ip += 1+j+k++) 
1211           }                                   << 1140                       { temp2 += *ip * x[k]; }
1212           for(ip += i - 1; k < nrow - j; ip + << 1141                     *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1213           {                                   << 1142                   }
1214             temp2 += *ip * x[k];              << 1143                   temp2 = 0;
1215           }                                   << 1144                   ip = m.begin() + (j+1)*j/2 + j-1;
1216           *(m.begin() + i * (i - 1) / 2 + j - << 1145                   for (k=0; k < nrow-j; ip += 1+j+k++)
1217         }                                     << 1146                     { temp2 += x[k] * *ip; }
1218         temp2 = 0;                            << 1147                   *mjj -= temp2;
1219         ip    = m.begin() + (j + 1) * j / 2 + << 1148                 }
1220         for(k = 0; k < nrow - j; ip += 1 + j  << 1149             }
1221         {                                     << 1150           else //2x2 pivot, compute columns j and j-1 of the inverse
1222           temp2 += x[k] * *ip;                << 1151             {
1223         }                                     << 1152               if (piv[j-1] != 0)
1224         *mjj -= temp2;                        << 1153                 { G4cerr << "error in piv" << piv[j-1] << G4endl; }
1225       }                                       << 1154               s=2; 
1226     }                                         << 1155               if (j < nrow)
1227     else  // 2x2 pivot, compute columns j and << 1156                 {
1228     {                                         << 1157                   ip = m.begin() + (j+1)*j/2 + j-1;
1229       if(piv[j - 1] != 0)                     << 1158                   for (i=0; i < nrow-j; ip += 1+j+i++)
1230       {                                       << 1159                   {
1231         std::ostringstream message;           << 1160                     x[i] = *ip;
1232         message << "Error in pivot: " << piv[ << 1161                   }
1233         G4Exception("G4ErrorSymMatrix::invert << 1162                   for (i=j+1; i<=nrow ; i++)
1234                     "GEANT4e-Notification", J << 1163                   {
1235       }                                       << 1164                     temp2 = 0;
1236       ss = 2;                                 << 1165                     ip = m.begin() + i*(i-1)/2 + j;
1237       if(j < nrow)                            << 1166                     for (k=0; k <= i-j-1; k++)
1238       {                                       << 1167                       { temp2 += *ip++ * x[k]; }
1239         ip = m.begin() + (j + 1) * j / 2 + j  << 1168                     for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1240         for(i = 0; i < nrow - j; ip += 1 + j  << 1169                       { temp2 += *ip * x[k]; }
1241         {                                     << 1170                     *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1242           x[i] = *ip;                         << 1171                   }    
1243         }                                     << 1172                   temp2 = 0;
1244         for(i = j + 1; i <= nrow; i++)        << 1173                   ip = m.begin() + (j+1)*j/2 + j-1;
1245         {                                     << 1174                   for (k=0; k < nrow-j; ip += 1+j+k++)
1246           temp2 = 0;                          << 1175                     { temp2 += x[k] * *ip; }
1247           ip    = m.begin() + i * (i - 1) / 2 << 1176                   *mjj -= temp2;
1248           for(k = 0; k <= i - j - 1; k++)     << 1177                   temp2 = 0;
1249           {                                   << 1178                   ip = m.begin() + (j+1)*j/2 + j-2;
1250             temp2 += *ip++ * x[k];            << 1179                   for (i=j+1; i <= nrow; ip += i++)
1251           }                                   << 1180                     { temp2 += *ip * *(ip+1); }
1252           for(ip += i - 1; k < nrow - j; ip + << 1181                   *(mjj-1) -= temp2;
1253           {                                   << 1182                   ip = m.begin() + (j+1)*j/2 + j-2;
1254             temp2 += *ip * x[k];              << 1183                   for (i=0; i < nrow-j; ip += 1+j+i++)
1255           }                                   << 1184                   {
1256           *(m.begin() + i * (i - 1) / 2 + j - << 1185                     x[i] = *ip;
1257         }                                     << 1186                   }
1258         temp2 = 0;                            << 1187                   for (i=j+1; i <= nrow ; i++)
1259         ip    = m.begin() + (j + 1) * j / 2 + << 1188                   {
1260         for(k = 0; k < nrow - j; ip += 1 + j  << 1189                     temp2 = 0;
1261         {                                     << 1190                     ip = m.begin() + i*(i-1)/2 + j;
1262           temp2 += x[k] * *ip;                << 1191                     for (k=0; k <= i-j-1; k++)
1263         }                                     << 1192                       { temp2 += *ip++ * x[k]; }
1264         *mjj -= temp2;                        << 1193                     for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1265         temp2 = 0;                            << 1194                       { temp2 += *ip * x[k]; }
1266         ip    = m.begin() + (j + 1) * j / 2 + << 1195                     *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1267         for(i = j + 1; i <= nrow; ip += i++)  << 1196                   }
1268         {                                     << 1197                   temp2 = 0;
1269           temp2 += *ip * *(ip + 1);           << 1198                   ip = m.begin() + (j+1)*j/2 + j-2;
1270         }                                     << 1199                   for (k=0; k < nrow-j; ip += 1+j+k++)
1271         *(mjj - 1) -= temp2;                  << 1200                     { temp2 += x[k] * *ip; }
1272         ip = m.begin() + (j + 1) * j / 2 + j  << 1201                   *(mjj-j) -= temp2;
1273         for(i = 0; i < nrow - j; ip += 1 + j  << 1202                 }
1274         {                                     << 1203             }  
1275           x[i] = *ip;                         << 1204           
1276         }                                     << 1205           // interchange rows and columns j and piv[j-1] 
1277         for(i = j + 1; i <= nrow; i++)        << 1206           // or rows and columns j and -piv[j-2]
1278         {                                     << 1207           
1279           temp2 = 0;                          << 1208           pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1280           ip    = m.begin() + i * (i - 1) / 2 << 1209           ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1281           for(k = 0; k <= i - j - 1; k++)     << 1210           for (i=j+1;i < pivrow; i++, ip++)
1282           {                                   << 1211             {
1283             temp2 += *ip++ * x[k];            << 1212               temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1284           }                                   << 1213               *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1285           for(ip += i - 1; k < nrow - j; ip + << 1214               *ip = temp1;
1286           {                                   << 1215             }
1287             temp2 += *ip * x[k];              << 1216           temp1 = *mjj;
1288           }                                   << 1217           *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1289           *(m.begin() + i * (i - 1) / 2 + j - << 1218           *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1290         }                                     << 1219           if (s==2)
1291         temp2 = 0;                            << 1220             {
1292         ip    = m.begin() + (j + 1) * j / 2 + << 1221               temp1 = *(mjj-1);
1293         for(k = 0; k < nrow - j; ip += 1 + j  << 1222               *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1294         {                                     << 1223               *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1295           temp2 += x[k] * *ip;                << 1224             }
1296         }                                     << 1225           
1297         *(mjj - j) -= temp2;                  << 1226           ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;  // &A(i,j)
1298       }                                       << 1227           iq = ip + pivrow-j;
1299     }                                         << 1228           for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1300                                               << 1229             {
1301     // interchange rows and columns j and piv << 1230               temp1 = *iq;
1302     // or rows and columns j and -piv[j-2]    << 1231               *iq = *ip;
1303                                               << 1232               *ip = temp1;
1304     pivrow = (piv[j - 1] == 0) ? -piv[j - 2]  << 1233             } 
1305     ip     = m.begin() + pivrow * (pivrow - 1 << 1234   } // end of loop over columns (in computing inverse from factorization)
1306     for(i = j + 1; i < pivrow; i++, ip++)     << 1235 
1307     {                                         << 1236   return; // inversion successful
1308       temp1 = *(m.begin() + i * (i - 1) / 2 + << 1237 }
1309       *(m.begin() + i * (i - 1) / 2 + j - 1)  << 1238 
1310       *ip                                     << 1239 G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0;
1311     }                                         << 1240 G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0;
1312     temp1 = *mjj;                             << 1241 G4double G4ErrorSymMatrix::adjustment5x5 = 0.0;
1313     *mjj  = *(m.begin() + pivrow * (pivrow -  << 1242 G4double G4ErrorSymMatrix::adjustment6x6 = 0.0;
1314     *(m.begin() + pivrow * (pivrow - 1) / 2 + << 1243 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
1315     if(ss == 2)                               << 1244 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
1316     {                                         << 1245 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005;
1317       temp1      = *(mjj - 1);                << 1246 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002;
1318       *(mjj - 1) = *(m.begin() + pivrow * (pi << 
1319       *(m.begin() + pivrow * (pivrow - 1) / 2 << 
1320     }                                         << 
1321                                               << 
1322     ip = m.begin() + (pivrow + 1) * pivrow /  << 
1323     iq = ip + pivrow - j;                     << 
1324     for(i = pivrow + 1; i <= nrow; ip += i, i << 
1325     {                                         << 
1326       temp1 = *iq;                            << 
1327       *iq   = *ip;                            << 
1328       *ip   = temp1;                          << 
1329     }                                         << 
1330   }  // end of loop over columns (in computin << 
1331                                               << 
1332   return;  // inversion successful            << 
1333 }                                             << 
1334                                               << 
1335 G4ThreadLocal G4double G4ErrorSymMatrix::posD << 
1336 G4ThreadLocal G4double G4ErrorSymMatrix::posD << 
1337 G4ThreadLocal G4double G4ErrorSymMatrix::adju << 
1338 G4ThreadLocal G4double G4ErrorSymMatrix::adju << 
1339 const G4double G4ErrorSymMatrix::CHOLESKY_THR << 
1340 const G4double G4ErrorSymMatrix::CHOLESKY_THR << 
1341 const G4double G4ErrorSymMatrix::CHOLESKY_CRE << 
1342 const G4double G4ErrorSymMatrix::CHOLESKY_CRE << 
1343                                                  1247 
1344 // Aij are indices for a 6x6 symmetric matrix    1248 // Aij are indices for a 6x6 symmetric matrix.
1345 //     The indices for 5x5 or 4x4 symmetric m << 1249 //     The indices for 5x5 or 4x4 symmetric matrices are the same, 
1346 //     ignoring all combinations with an inde    1250 //     ignoring all combinations with an index which is inapplicable.
1347                                                  1251 
1348 #define A00 0                                    1252 #define A00 0
1349 #define A01 1                                    1253 #define A01 1
1350 #define A02 3                                    1254 #define A02 3
1351 #define A03 6                                    1255 #define A03 6
1352 #define A04 10                                   1256 #define A04 10
1353 #define A05 15                                   1257 #define A05 15
1354                                                  1258 
1355 #define A10 1                                    1259 #define A10 1
1356 #define A11 2                                    1260 #define A11 2
1357 #define A12 4                                    1261 #define A12 4
1358 #define A13 7                                    1262 #define A13 7
1359 #define A14 11                                   1263 #define A14 11
1360 #define A15 16                                   1264 #define A15 16
1361                                                  1265 
1362 #define A20 3                                    1266 #define A20 3
1363 #define A21 4                                    1267 #define A21 4
1364 #define A22 5                                    1268 #define A22 5
1365 #define A23 8                                    1269 #define A23 8
1366 #define A24 12                                   1270 #define A24 12
1367 #define A25 17                                   1271 #define A25 17
1368                                                  1272 
1369 #define A30 6                                    1273 #define A30 6
1370 #define A31 7                                    1274 #define A31 7
1371 #define A32 8                                    1275 #define A32 8
1372 #define A33 9                                    1276 #define A33 9
1373 #define A34 13                                   1277 #define A34 13
1374 #define A35 18                                   1278 #define A35 18
1375                                                  1279 
1376 #define A40 10                                   1280 #define A40 10
1377 #define A41 11                                   1281 #define A41 11
1378 #define A42 12                                   1282 #define A42 12
1379 #define A43 13                                   1283 #define A43 13
1380 #define A44 14                                   1284 #define A44 14
1381 #define A45 19                                   1285 #define A45 19
1382                                                  1286 
1383 #define A50 15                                   1287 #define A50 15
1384 #define A51 16                                   1288 #define A51 16
1385 #define A52 17                                   1289 #define A52 17
1386 #define A53 18                                   1290 #define A53 18
1387 #define A54 19                                   1291 #define A54 19
1388 #define A55 20                                   1292 #define A55 20
1389                                                  1293 
1390 void G4ErrorSymMatrix::invert5(G4int& ifail)  << 1294 void G4ErrorSymMatrix::invert5(G4int & ifail)
1391 {                                             << 1295 { 
1392   if(posDefFraction5x5 >= CHOLESKY_THRESHOLD_ << 1296       if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
1393   {                                           << 1297       {
1394     invertCholesky5(ifail);                   << 1298         invertCholesky5(ifail);
1395     posDefFraction5x5 = .9 * posDefFraction5x << 1299         posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
1396     if(ifail != 0)  // Cholesky failed -- inv << 1300         if (ifail!=0)    // Cholesky failed -- invert using Haywood
1397     {                                         << 1301         {
1398       invertHaywood5(ifail);                  << 1302           invertHaywood5(ifail);      
1399     }                                         << 1303         }
1400   }                                           << 1304       }
1401   else                                        << 1305       else
1402   {                                           << 
1403     if(posDefFraction5x5 + adjustment5x5 >= C << 
1404     {                                         << 
1405       invertCholesky5(ifail);                 << 
1406       posDefFraction5x5 = .9 * posDefFraction << 
1407       if(ifail != 0)  // Cholesky failed -- i << 
1408       {                                          1306       {
1409         invertHaywood5(ifail);                << 1307         if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
1410         adjustment5x5 = 0;                    << 1308         {
                                                   >> 1309           invertCholesky5(ifail);
                                                   >> 1310           posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
                                                   >> 1311           if (ifail!=0)    // Cholesky failed -- invert using Haywood
                                                   >> 1312           {
                                                   >> 1313             invertHaywood5(ifail);      
                                                   >> 1314             adjustment5x5 = 0;
                                                   >> 1315           }
                                                   >> 1316         }
                                                   >> 1317         else
                                                   >> 1318         {
                                                   >> 1319           invertHaywood5(ifail);      
                                                   >> 1320           adjustment5x5 += CHOLESKY_CREEP_5x5;
                                                   >> 1321         }
1411       }                                          1322       }
1412     }                                         << 1323       return;
1413     else                                      << 
1414     {                                         << 
1415       invertHaywood5(ifail);                  << 
1416       adjustment5x5 += CHOLESKY_CREEP_5x5;    << 
1417     }                                         << 
1418   }                                           << 
1419   return;                                     << 
1420 }                                                1324 }
1421                                                  1325 
1422 void G4ErrorSymMatrix::invert6(G4int& ifail)  << 1326 void G4ErrorSymMatrix::invert6(G4int & ifail)
1423 {                                                1327 {
1424   if(posDefFraction6x6 >= CHOLESKY_THRESHOLD_ << 1328       if (posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6)
1425   {                                           << 
1426     invertCholesky6(ifail);                   << 
1427     posDefFraction6x6 = .9 * posDefFraction6x << 
1428     if(ifail != 0)  // Cholesky failed -- inv << 
1429     {                                         << 
1430       invertHaywood6(ifail);                  << 
1431     }                                         << 
1432   }                                           << 
1433   else                                        << 
1434   {                                           << 
1435     if(posDefFraction6x6 + adjustment6x6 >= C << 
1436     {                                         << 
1437       invertCholesky6(ifail);                 << 
1438       posDefFraction6x6 = .9 * posDefFraction << 
1439       if(ifail != 0)  // Cholesky failed -- i << 
1440       {                                          1329       {
1441         invertHaywood6(ifail);                << 1330         invertCholesky6(ifail);
1442         adjustment6x6 = 0;                    << 1331         posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
                                                   >> 1332         if (ifail!=0)    // Cholesky failed -- invert using Haywood
                                                   >> 1333         {
                                                   >> 1334           invertHaywood6(ifail);      
                                                   >> 1335         }
1443       }                                          1336       }
1444     }                                         << 1337       else
1445     else                                      << 1338       {
1446     {                                         << 1339         if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
1447       invertHaywood6(ifail);                  << 1340         {
1448       adjustment6x6 += CHOLESKY_CREEP_6x6;    << 1341           invertCholesky6(ifail);
1449     }                                         << 1342           posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
1450   }                                           << 1343           if (ifail!=0)    // Cholesky failed -- invert using Haywood
1451   return;                                     << 1344           {
                                                   >> 1345             invertHaywood6(ifail);      
                                                   >> 1346             adjustment6x6 = 0;
                                                   >> 1347           }
                                                   >> 1348         }
                                                   >> 1349         else
                                                   >> 1350         {
                                                   >> 1351           invertHaywood6(ifail);      
                                                   >> 1352           adjustment6x6 += CHOLESKY_CREEP_6x6;
                                                   >> 1353         }
                                                   >> 1354       }
                                                   >> 1355       return;
1452 }                                                1356 }
1453                                                  1357 
1454 void G4ErrorSymMatrix::invertHaywood5(G4int&  << 1358 void G4ErrorSymMatrix::invertHaywood5  (G4int & ifail)
1455 {                                                1359 {
1456   ifail = 0;                                     1360   ifail = 0;
1457                                                  1361 
1458   // Find all NECESSARY 2x2 dets:  (25 of the    1362   // Find all NECESSARY 2x2 dets:  (25 of them)
1459                                                  1363 
1460   G4double Det2_23_01 = m[A20] * m[A31] - m[A << 1364   G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
1461   G4double Det2_23_02 = m[A20] * m[A32] - m[A << 1365   G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
1462   G4double Det2_23_03 = m[A20] * m[A33] - m[A << 1366   G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
1463   G4double Det2_23_12 = m[A21] * m[A32] - m[A << 1367   G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
1464   G4double Det2_23_13 = m[A21] * m[A33] - m[A << 1368   G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
1465   G4double Det2_23_23 = m[A22] * m[A33] - m[A << 1369   G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
1466   G4double Det2_24_01 = m[A20] * m[A41] - m[A << 1370   G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
1467   G4double Det2_24_02 = m[A20] * m[A42] - m[A << 1371   G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
1468   G4double Det2_24_03 = m[A20] * m[A43] - m[A << 1372   G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
1469   G4double Det2_24_04 = m[A20] * m[A44] - m[A << 1373   G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
1470   G4double Det2_24_12 = m[A21] * m[A42] - m[A << 1374   G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
1471   G4double Det2_24_13 = m[A21] * m[A43] - m[A << 1375   G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
1472   G4double Det2_24_14 = m[A21] * m[A44] - m[A << 1376   G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
1473   G4double Det2_24_23 = m[A22] * m[A43] - m[A << 1377   G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
1474   G4double Det2_24_24 = m[A22] * m[A44] - m[A << 1378   G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
1475   G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1379   G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1476   G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1380   G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1477   G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1381   G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1478   G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1382   G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1479   G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1383   G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1480   G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1384   G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1481   G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1385   G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1482   G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1386   G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1483   G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1387   G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1484   G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1388   G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1485                                                  1389 
1486   // Find all NECESSARY 3x3 dets:   (30 of th    1390   // Find all NECESSARY 3x3 dets:   (30 of them)
1487                                                  1391 
1488   G4double Det3_123_012 =                     << 1392   G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02 
1489     m[A10] * Det2_23_12 - m[A11] * Det2_23_02 << 1393                                 + m[A12]*Det2_23_01;
1490   G4double Det3_123_013 =                     << 1394   G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03 
1491     m[A10] * Det2_23_13 - m[A11] * Det2_23_03 << 1395                                 + m[A13]*Det2_23_01;
1492   G4double Det3_123_023 =                     << 1396   G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03 
1493     m[A10] * Det2_23_23 - m[A12] * Det2_23_03 << 1397                                 + m[A13]*Det2_23_02;
1494   G4double Det3_123_123 =                     << 1398   G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13 
1495     m[A11] * Det2_23_23 - m[A12] * Det2_23_13 << 1399                                 + m[A13]*Det2_23_12;
1496   G4double Det3_124_012 =                     << 1400   G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02 
1497     m[A10] * Det2_24_12 - m[A11] * Det2_24_02 << 1401                                 + m[A12]*Det2_24_01;
1498   G4double Det3_124_013 =                     << 1402   G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03 
1499     m[A10] * Det2_24_13 - m[A11] * Det2_24_03 << 1403                                 + m[A13]*Det2_24_01;
1500   G4double Det3_124_014 =                     << 1404   G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04 
1501     m[A10] * Det2_24_14 - m[A11] * Det2_24_04 << 1405                                 + m[A14]*Det2_24_01;
1502   G4double Det3_124_023 =                     << 1406   G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03 
1503     m[A10] * Det2_24_23 - m[A12] * Det2_24_03 << 1407                                 + m[A13]*Det2_24_02;
1504   G4double Det3_124_024 =                     << 1408   G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04 
1505     m[A10] * Det2_24_24 - m[A12] * Det2_24_04 << 1409                                 + m[A14]*Det2_24_02;
1506   G4double Det3_124_123 =                     << 1410   G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13 
1507     m[A11] * Det2_24_23 - m[A12] * Det2_24_13 << 1411                                 + m[A13]*Det2_24_12;
1508   G4double Det3_124_124 =                     << 1412   G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14 
1509     m[A11] * Det2_24_24 - m[A12] * Det2_24_14 << 1413                                 + m[A14]*Det2_24_12;
1510   G4double Det3_134_012 =                     << 1414   G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02 
1511     m[A10] * Det2_34_12 - m[A11] * Det2_34_02 << 1415                                 + m[A12]*Det2_34_01;
1512   G4double Det3_134_013 =                     << 1416   G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03 
1513     m[A10] * Det2_34_13 - m[A11] * Det2_34_03 << 1417                                 + m[A13]*Det2_34_01;
1514   G4double Det3_134_014 =                     << 1418   G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04 
1515     m[A10] * Det2_34_14 - m[A11] * Det2_34_04 << 1419                                 + m[A14]*Det2_34_01;
1516   G4double Det3_134_023 =                     << 1420   G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03 
1517     m[A10] * Det2_34_23 - m[A12] * Det2_34_03 << 1421                                 + m[A13]*Det2_34_02;
1518   G4double Det3_134_024 =                     << 1422   G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04 
1519     m[A10] * Det2_34_24 - m[A12] * Det2_34_04 << 1423                                 + m[A14]*Det2_34_02;
1520   G4double Det3_134_034 =                     << 1424   G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04 
1521     m[A10] * Det2_34_34 - m[A13] * Det2_34_04 << 1425                                 + m[A14]*Det2_34_03;
1522   G4double Det3_134_123 =                     << 1426   G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13 
1523     m[A11] * Det2_34_23 - m[A12] * Det2_34_13 << 1427                                 + m[A13]*Det2_34_12;
1524   G4double Det3_134_124 =                     << 1428   G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14 
1525     m[A11] * Det2_34_24 - m[A12] * Det2_34_14 << 1429                                 + m[A14]*Det2_34_12;
1526   G4double Det3_134_134 =                     << 1430   G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14 
1527     m[A11] * Det2_34_34 - m[A13] * Det2_34_14 << 1431                                 + m[A14]*Det2_34_13;
1528   G4double Det3_234_012 =                     << 1432   G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
1529     m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1433                                 + m[A22]*Det2_34_01;
1530   G4double Det3_234_013 =                     << 1434   G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
1531     m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1435                                 + m[A23]*Det2_34_01;
1532   G4double Det3_234_014 =                     << 1436   G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
1533     m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1437                                 + m[A24]*Det2_34_01;
1534   G4double Det3_234_023 =                     << 1438   G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
1535     m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1439                                 + m[A23]*Det2_34_02;
1536   G4double Det3_234_024 =                     << 1440   G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
1537     m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1441                                 + m[A24]*Det2_34_02;
1538   G4double Det3_234_034 =                     << 1442   G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
1539     m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1443                                 + m[A24]*Det2_34_03;
1540   G4double Det3_234_123 =                     << 1444   G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
1541     m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1445                                 + m[A23]*Det2_34_12;
1542   G4double Det3_234_124 =                     << 1446   G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
1543     m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1447                                 + m[A24]*Det2_34_12;
1544   G4double Det3_234_134 =                     << 1448   G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
1545     m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1449                                 + m[A24]*Det2_34_13;
1546   G4double Det3_234_234 =                     << 1450   G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
1547     m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1451                                 + m[A24]*Det2_34_23;
1548                                                  1452 
1549   // Find all NECESSARY 4x4 dets:   (15 of th    1453   // Find all NECESSARY 4x4 dets:   (15 of them)
1550                                                  1454 
1551   G4double Det4_0123_0123 = m[A00] * Det3_123 << 1455   G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023 
1552                             m[A02] * Det3_123 << 1456                                 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
1553   G4double Det4_0124_0123 = m[A00] * Det3_124 << 1457   G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023 
1554                             m[A02] * Det3_124 << 1458                                 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
1555   G4double Det4_0124_0124 = m[A00] * Det3_124 << 1459   G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024 
1556                             m[A02] * Det3_124 << 1460                                 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
1557   G4double Det4_0134_0123 = m[A00] * Det3_134 << 1461   G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023 
1558                             m[A02] * Det3_134 << 1462                                 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
1559   G4double Det4_0134_0124 = m[A00] * Det3_134 << 1463   G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024 
1560                             m[A02] * Det3_134 << 1464                                 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
1561   G4double Det4_0134_0134 = m[A00] * Det3_134 << 1465   G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034 
1562                             m[A03] * Det3_134 << 1466                                 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
1563   G4double Det4_0234_0123 = m[A00] * Det3_234 << 1467   G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023 
1564                             m[A02] * Det3_234 << 1468                                 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
1565   G4double Det4_0234_0124 = m[A00] * Det3_234 << 1469   G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024 
1566                             m[A02] * Det3_234 << 1470                                 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
1567   G4double Det4_0234_0134 = m[A00] * Det3_234 << 1471   G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034 
1568                             m[A03] * Det3_234 << 1472                                 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
1569   G4double Det4_0234_0234 = m[A00] * Det3_234 << 1473   G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034 
1570                             m[A03] * Det3_234 << 1474                                 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
1571   G4double Det4_1234_0123 = m[A10] * Det3_234 << 1475   G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
1572                             m[A12] * Det3_234 << 1476                                 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1573   G4double Det4_1234_0124 = m[A10] * Det3_234 << 1477   G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
1574                             m[A12] * Det3_234 << 1478                                 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1575   G4double Det4_1234_0134 = m[A10] * Det3_234 << 1479   G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
1576                             m[A13] * Det3_234 << 1480                                 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1577   G4double Det4_1234_0234 = m[A10] * Det3_234 << 1481   G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
1578                             m[A13] * Det3_234 << 1482                                 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1579   G4double Det4_1234_1234 = m[A11] * Det3_234 << 1483   G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
1580                             m[A13] * Det3_234 << 1484                                 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1581                                                  1485 
1582   // Find the 5x5 det:                           1486   // Find the 5x5 det:
1583                                                  1487 
1584   G4double det = m[A00] * Det4_1234_1234 - m[ << 1488   G4double det =  m[A00]*Det4_1234_1234 
1585                  m[A02] * Det4_1234_0134 - m[ << 1489                 - m[A01]*Det4_1234_0234 
1586                  m[A04] * Det4_1234_0123;     << 1490                 + m[A02]*Det4_1234_0134 
                                                   >> 1491                 - m[A03]*Det4_1234_0124 
                                                   >> 1492                 + m[A04]*Det4_1234_0123;
1587                                                  1493 
1588   if(det == 0)                                << 1494   if ( det == 0 )
1589   {                                           << 1495   {  
1590     ifail = 1;                                   1496     ifail = 1;
1591     return;                                      1497     return;
1592   }                                           << 1498   } 
1593                                                  1499 
1594   G4double oneOverDet = 1.0 / det;            << 1500   G4double oneOverDet = 1.0/det;
1595   G4double mn1OverDet = -oneOverDet;          << 1501   G4double mn1OverDet = - oneOverDet;
1596                                                  1502 
1597   m[A00] = Det4_1234_1234 * oneOverDet;       << 1503   m[A00] =  Det4_1234_1234 * oneOverDet;
1598   m[A01] = Det4_1234_0234 * mn1OverDet;       << 1504   m[A01] =  Det4_1234_0234 * mn1OverDet;
1599   m[A02] = Det4_1234_0134 * oneOverDet;       << 1505   m[A02] =  Det4_1234_0134 * oneOverDet;
1600   m[A03] = Det4_1234_0124 * mn1OverDet;       << 1506   m[A03] =  Det4_1234_0124 * mn1OverDet;
1601   m[A04] = Det4_1234_0123 * oneOverDet;       << 1507   m[A04] =  Det4_1234_0123 * oneOverDet;
1602                                               << 1508 
1603   m[A11] = Det4_0234_0234 * oneOverDet;       << 1509   m[A11] =  Det4_0234_0234 * oneOverDet;
1604   m[A12] = Det4_0234_0134 * mn1OverDet;       << 1510   m[A12] =  Det4_0234_0134 * mn1OverDet;
1605   m[A13] = Det4_0234_0124 * oneOverDet;       << 1511   m[A13] =  Det4_0234_0124 * oneOverDet;
1606   m[A14] = Det4_0234_0123 * mn1OverDet;       << 1512   m[A14] =  Det4_0234_0123 * mn1OverDet;
1607                                               << 1513 
1608   m[A22] = Det4_0134_0134 * oneOverDet;       << 1514   m[A22] =  Det4_0134_0134 * oneOverDet;
1609   m[A23] = Det4_0134_0124 * mn1OverDet;       << 1515   m[A23] =  Det4_0134_0124 * mn1OverDet;
1610   m[A24] = Det4_0134_0123 * oneOverDet;       << 1516   m[A24] =  Det4_0134_0123 * oneOverDet;
1611                                                  1517 
1612   m[A33] = Det4_0124_0124 * oneOverDet;       << 1518   m[A33] =  Det4_0124_0124 * oneOverDet;
1613   m[A34] = Det4_0124_0123 * mn1OverDet;       << 1519   m[A34] =  Det4_0124_0123 * mn1OverDet;
1614                                                  1520 
1615   m[A44] = Det4_0123_0123 * oneOverDet;       << 1521   m[A44] =  Det4_0123_0123 * oneOverDet;
1616                                                  1522 
1617   return;                                        1523   return;
1618 }                                                1524 }
1619                                                  1525 
1620 void G4ErrorSymMatrix::invertHaywood6(G4int&  << 1526 void G4ErrorSymMatrix::invertHaywood6  (G4int & ifail)
1621 {                                                1527 {
1622   ifail = 0;                                     1528   ifail = 0;
1623                                                  1529 
1624   // Find all NECESSARY 2x2 dets:  (39 of the    1530   // Find all NECESSARY 2x2 dets:  (39 of them)
1625                                                  1531 
1626   G4double Det2_34_01 = m[A30] * m[A41] - m[A << 1532   G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1627   G4double Det2_34_02 = m[A30] * m[A42] - m[A << 1533   G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1628   G4double Det2_34_03 = m[A30] * m[A43] - m[A << 1534   G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1629   G4double Det2_34_04 = m[A30] * m[A44] - m[A << 1535   G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1630   G4double Det2_34_12 = m[A31] * m[A42] - m[A << 1536   G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1631   G4double Det2_34_13 = m[A31] * m[A43] - m[A << 1537   G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1632   G4double Det2_34_14 = m[A31] * m[A44] - m[A << 1538   G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1633   G4double Det2_34_23 = m[A32] * m[A43] - m[A << 1539   G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1634   G4double Det2_34_24 = m[A32] * m[A44] - m[A << 1540   G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1635   G4double Det2_34_34 = m[A33] * m[A44] - m[A << 1541   G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1636   G4double Det2_35_01 = m[A30] * m[A51] - m[A << 1542   G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1637   G4double Det2_35_02 = m[A30] * m[A52] - m[A << 1543   G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1638   G4double Det2_35_03 = m[A30] * m[A53] - m[A << 1544   G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1639   G4double Det2_35_04 = m[A30] * m[A54] - m[A << 1545   G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1640   G4double Det2_35_05 = m[A30] * m[A55] - m[A << 1546   G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1641   G4double Det2_35_12 = m[A31] * m[A52] - m[A << 1547   G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1642   G4double Det2_35_13 = m[A31] * m[A53] - m[A << 1548   G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1643   G4double Det2_35_14 = m[A31] * m[A54] - m[A << 1549   G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1644   G4double Det2_35_15 = m[A31] * m[A55] - m[A << 1550   G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1645   G4double Det2_35_23 = m[A32] * m[A53] - m[A << 1551   G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1646   G4double Det2_35_24 = m[A32] * m[A54] - m[A << 1552   G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1647   G4double Det2_35_25 = m[A32] * m[A55] - m[A << 1553   G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1648   G4double Det2_35_34 = m[A33] * m[A54] - m[A << 1554   G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1649   G4double Det2_35_35 = m[A33] * m[A55] - m[A << 1555   G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1650   G4double Det2_45_01 = m[A40] * m[A51] - m[A << 1556   G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1651   G4double Det2_45_02 = m[A40] * m[A52] - m[A << 1557   G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1652   G4double Det2_45_03 = m[A40] * m[A53] - m[A << 1558   G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1653   G4double Det2_45_04 = m[A40] * m[A54] - m[A << 1559   G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1654   G4double Det2_45_05 = m[A40] * m[A55] - m[A << 1560   G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1655   G4double Det2_45_12 = m[A41] * m[A52] - m[A << 1561   G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1656   G4double Det2_45_13 = m[A41] * m[A53] - m[A << 1562   G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1657   G4double Det2_45_14 = m[A41] * m[A54] - m[A << 1563   G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1658   G4double Det2_45_15 = m[A41] * m[A55] - m[A << 1564   G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1659   G4double Det2_45_23 = m[A42] * m[A53] - m[A << 1565   G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1660   G4double Det2_45_24 = m[A42] * m[A54] - m[A << 1566   G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1661   G4double Det2_45_25 = m[A42] * m[A55] - m[A << 1567   G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1662   G4double Det2_45_34 = m[A43] * m[A54] - m[A << 1568   G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1663   G4double Det2_45_35 = m[A43] * m[A55] - m[A << 1569   G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1664   G4double Det2_45_45 = m[A44] * m[A55] - m[A << 1570   G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1665                                                  1571 
1666   // Find all NECESSARY 3x3 dets:  (65 of the    1572   // Find all NECESSARY 3x3 dets:  (65 of them)
1667                                                  1573 
1668   G4double Det3_234_012 =                     << 1574   G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
1669     m[A20] * Det2_34_12 - m[A21] * Det2_34_02 << 1575                                                 + m[A22]*Det2_34_01;
1670   G4double Det3_234_013 =                     << 1576   G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
1671     m[A20] * Det2_34_13 - m[A21] * Det2_34_03 << 1577                                                 + m[A23]*Det2_34_01;
1672   G4double Det3_234_014 =                     << 1578   G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
1673     m[A20] * Det2_34_14 - m[A21] * Det2_34_04 << 1579                                                 + m[A24]*Det2_34_01;
1674   G4double Det3_234_023 =                     << 1580   G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
1675     m[A20] * Det2_34_23 - m[A22] * Det2_34_03 << 1581                                                 + m[A23]*Det2_34_02;
1676   G4double Det3_234_024 =                     << 1582   G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
1677     m[A20] * Det2_34_24 - m[A22] * Det2_34_04 << 1583                                                 + m[A24]*Det2_34_02;
1678   G4double Det3_234_034 =                     << 1584   G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
1679     m[A20] * Det2_34_34 - m[A23] * Det2_34_04 << 1585                                                 + m[A24]*Det2_34_03;
1680   G4double Det3_234_123 =                     << 1586   G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
1681     m[A21] * Det2_34_23 - m[A22] * Det2_34_13 << 1587                                                 + m[A23]*Det2_34_12;
1682   G4double Det3_234_124 =                     << 1588   G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
1683     m[A21] * Det2_34_24 - m[A22] * Det2_34_14 << 1589                                                 + m[A24]*Det2_34_12;
1684   G4double Det3_234_134 =                     << 1590   G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
1685     m[A21] * Det2_34_34 - m[A23] * Det2_34_14 << 1591                                                 + m[A24]*Det2_34_13;
1686   G4double Det3_234_234 =                     << 1592   G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
1687     m[A22] * Det2_34_34 - m[A23] * Det2_34_24 << 1593                                                 + m[A24]*Det2_34_23;
1688   G4double Det3_235_012 =                     << 1594   G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02 
1689     m[A20] * Det2_35_12 - m[A21] * Det2_35_02 << 1595                                                 + m[A22]*Det2_35_01;
1690   G4double Det3_235_013 =                     << 1596   G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03 
1691     m[A20] * Det2_35_13 - m[A21] * Det2_35_03 << 1597                                                 + m[A23]*Det2_35_01;
1692   G4double Det3_235_014 =                     << 1598   G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04 
1693     m[A20] * Det2_35_14 - m[A21] * Det2_35_04 << 1599                                                 + m[A24]*Det2_35_01;
1694   G4double Det3_235_015 =                     << 1600   G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05 
1695     m[A20] * Det2_35_15 - m[A21] * Det2_35_05 << 1601                                                 + m[A25]*Det2_35_01;
1696   G4double Det3_235_023 =                     << 1602   G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03 
1697     m[A20] * Det2_35_23 - m[A22] * Det2_35_03 << 1603                                                 + m[A23]*Det2_35_02;
1698   G4double Det3_235_024 =                     << 1604   G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04 
1699     m[A20] * Det2_35_24 - m[A22] * Det2_35_04 << 1605                                                 + m[A24]*Det2_35_02;
1700   G4double Det3_235_025 =                     << 1606   G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05 
1701     m[A20] * Det2_35_25 - m[A22] * Det2_35_05 << 1607                                                 + m[A25]*Det2_35_02;
1702   G4double Det3_235_034 =                     << 1608   G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04 
1703     m[A20] * Det2_35_34 - m[A23] * Det2_35_04 << 1609                                                 + m[A24]*Det2_35_03;
1704   G4double Det3_235_035 =                     << 1610   G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05 
1705     m[A20] * Det2_35_35 - m[A23] * Det2_35_05 << 1611                                                 + m[A25]*Det2_35_03;
1706   G4double Det3_235_123 =                     << 1612   G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13 
1707     m[A21] * Det2_35_23 - m[A22] * Det2_35_13 << 1613                                                 + m[A23]*Det2_35_12;
1708   G4double Det3_235_124 =                     << 1614   G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14 
1709     m[A21] * Det2_35_24 - m[A22] * Det2_35_14 << 1615                                                 + m[A24]*Det2_35_12;
1710   G4double Det3_235_125 =                     << 1616   G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15 
1711     m[A21] * Det2_35_25 - m[A22] * Det2_35_15 << 1617                                                 + m[A25]*Det2_35_12;
1712   G4double Det3_235_134 =                     << 1618   G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14 
1713     m[A21] * Det2_35_34 - m[A23] * Det2_35_14 << 1619                                                 + m[A24]*Det2_35_13;
1714   G4double Det3_235_135 =                     << 1620   G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15 
1715     m[A21] * Det2_35_35 - m[A23] * Det2_35_15 << 1621                                                 + m[A25]*Det2_35_13;
1716   G4double Det3_235_234 =                     << 1622   G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24 
1717     m[A22] * Det2_35_34 - m[A23] * Det2_35_24 << 1623                                                 + m[A24]*Det2_35_23;
1718   G4double Det3_235_235 =                     << 1624   G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25 
1719     m[A22] * Det2_35_35 - m[A23] * Det2_35_25 << 1625                                                 + m[A25]*Det2_35_23;
1720   G4double Det3_245_012 =                     << 1626   G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02 
1721     m[A20] * Det2_45_12 - m[A21] * Det2_45_02 << 1627                                                 + m[A22]*Det2_45_01;
1722   G4double Det3_245_013 =                     << 1628   G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03 
1723     m[A20] * Det2_45_13 - m[A21] * Det2_45_03 << 1629                                                 + m[A23]*Det2_45_01;
1724   G4double Det3_245_014 =                     << 1630   G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04 
1725     m[A20] * Det2_45_14 - m[A21] * Det2_45_04 << 1631                                                 + m[A24]*Det2_45_01;
1726   G4double Det3_245_015 =                     << 1632   G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05 
1727     m[A20] * Det2_45_15 - m[A21] * Det2_45_05 << 1633                                                 + m[A25]*Det2_45_01;
1728   G4double Det3_245_023 =                     << 1634   G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03 
1729     m[A20] * Det2_45_23 - m[A22] * Det2_45_03 << 1635                                                 + m[A23]*Det2_45_02;
1730   G4double Det3_245_024 =                     << 1636   G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04 
1731     m[A20] * Det2_45_24 - m[A22] * Det2_45_04 << 1637                                                 + m[A24]*Det2_45_02;
1732   G4double Det3_245_025 =                     << 1638   G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05 
1733     m[A20] * Det2_45_25 - m[A22] * Det2_45_05 << 1639                                                 + m[A25]*Det2_45_02;
1734   G4double Det3_245_034 =                     << 1640   G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04 
1735     m[A20] * Det2_45_34 - m[A23] * Det2_45_04 << 1641                                                 + m[A24]*Det2_45_03;
1736   G4double Det3_245_035 =                     << 1642   G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05 
1737     m[A20] * Det2_45_35 - m[A23] * Det2_45_05 << 1643                                                 + m[A25]*Det2_45_03;
1738   G4double Det3_245_045 =                     << 1644   G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05 
1739     m[A20] * Det2_45_45 - m[A24] * Det2_45_05 << 1645                                                 + m[A25]*Det2_45_04;
1740   G4double Det3_245_123 =                     << 1646   G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13 
1741     m[A21] * Det2_45_23 - m[A22] * Det2_45_13 << 1647                                                 + m[A23]*Det2_45_12;
1742   G4double Det3_245_124 =                     << 1648   G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14 
1743     m[A21] * Det2_45_24 - m[A22] * Det2_45_14 << 1649                                                 + m[A24]*Det2_45_12;
1744   G4double Det3_245_125 =                     << 1650   G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15 
1745     m[A21] * Det2_45_25 - m[A22] * Det2_45_15 << 1651                                                 + m[A25]*Det2_45_12;
1746   G4double Det3_245_134 =                     << 1652   G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14 
1747     m[A21] * Det2_45_34 - m[A23] * Det2_45_14 << 1653                                                 + m[A24]*Det2_45_13;
1748   G4double Det3_245_135 =                     << 1654   G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15 
1749     m[A21] * Det2_45_35 - m[A23] * Det2_45_15 << 1655                                                 + m[A25]*Det2_45_13;
1750   G4double Det3_245_145 =                     << 1656   G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15 
1751     m[A21] * Det2_45_45 - m[A24] * Det2_45_15 << 1657                                                 + m[A25]*Det2_45_14;
1752   G4double Det3_245_234 =                     << 1658   G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24 
1753     m[A22] * Det2_45_34 - m[A23] * Det2_45_24 << 1659                                                 + m[A24]*Det2_45_23;
1754   G4double Det3_245_235 =                     << 1660   G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25 
1755     m[A22] * Det2_45_35 - m[A23] * Det2_45_25 << 1661                                                 + m[A25]*Det2_45_23;
1756   G4double Det3_245_245 =                     << 1662   G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25 
1757     m[A22] * Det2_45_45 - m[A24] * Det2_45_25 << 1663                                                 + m[A25]*Det2_45_24;
1758   G4double Det3_345_012 =                     << 1664   G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02 
1759     m[A30] * Det2_45_12 - m[A31] * Det2_45_02 << 1665                                                 + m[A32]*Det2_45_01;
1760   G4double Det3_345_013 =                     << 1666   G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03 
1761     m[A30] * Det2_45_13 - m[A31] * Det2_45_03 << 1667                                                 + m[A33]*Det2_45_01;
1762   G4double Det3_345_014 =                     << 1668   G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04 
1763     m[A30] * Det2_45_14 - m[A31] * Det2_45_04 << 1669                                                 + m[A34]*Det2_45_01;
1764   G4double Det3_345_015 =                     << 1670   G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05 
1765     m[A30] * Det2_45_15 - m[A31] * Det2_45_05 << 1671                                                 + m[A35]*Det2_45_01;
1766   G4double Det3_345_023 =                     << 1672   G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03 
1767     m[A30] * Det2_45_23 - m[A32] * Det2_45_03 << 1673                                                 + m[A33]*Det2_45_02;
1768   G4double Det3_345_024 =                     << 1674   G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04 
1769     m[A30] * Det2_45_24 - m[A32] * Det2_45_04 << 1675                                                 + m[A34]*Det2_45_02;
1770   G4double Det3_345_025 =                     << 1676   G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05 
1771     m[A30] * Det2_45_25 - m[A32] * Det2_45_05 << 1677                                                 + m[A35]*Det2_45_02;
1772   G4double Det3_345_034 =                     << 1678   G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04 
1773     m[A30] * Det2_45_34 - m[A33] * Det2_45_04 << 1679                                                 + m[A34]*Det2_45_03;
1774   G4double Det3_345_035 =                     << 1680   G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05 
1775     m[A30] * Det2_45_35 - m[A33] * Det2_45_05 << 1681                                                 + m[A35]*Det2_45_03;
1776   G4double Det3_345_045 =                     << 1682   G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05 
1777     m[A30] * Det2_45_45 - m[A34] * Det2_45_05 << 1683                                                 + m[A35]*Det2_45_04;
1778   G4double Det3_345_123 =                     << 1684   G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13 
1779     m[A31] * Det2_45_23 - m[A32] * Det2_45_13 << 1685                                                 + m[A33]*Det2_45_12;
1780   G4double Det3_345_124 =                     << 1686   G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14 
1781     m[A31] * Det2_45_24 - m[A32] * Det2_45_14 << 1687                                                 + m[A34]*Det2_45_12;
1782   G4double Det3_345_125 =                     << 1688   G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15 
1783     m[A31] * Det2_45_25 - m[A32] * Det2_45_15 << 1689                                                 + m[A35]*Det2_45_12;
1784   G4double Det3_345_134 =                     << 1690   G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14 
1785     m[A31] * Det2_45_34 - m[A33] * Det2_45_14 << 1691                                                 + m[A34]*Det2_45_13;
1786   G4double Det3_345_135 =                     << 1692   G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15 
1787     m[A31] * Det2_45_35 - m[A33] * Det2_45_15 << 1693                                                 + m[A35]*Det2_45_13;
1788   G4double Det3_345_145 =                     << 1694   G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15 
1789     m[A31] * Det2_45_45 - m[A34] * Det2_45_15 << 1695                                                 + m[A35]*Det2_45_14;
1790   G4double Det3_345_234 =                     << 1696   G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24 
1791     m[A32] * Det2_45_34 - m[A33] * Det2_45_24 << 1697                                                 + m[A34]*Det2_45_23;
1792   G4double Det3_345_235 =                     << 1698   G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25 
1793     m[A32] * Det2_45_35 - m[A33] * Det2_45_25 << 1699                                                 + m[A35]*Det2_45_23;
1794   G4double Det3_345_245 =                     << 1700   G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25 
1795     m[A32] * Det2_45_45 - m[A34] * Det2_45_25 << 1701                                                 + m[A35]*Det2_45_24;
1796   G4double Det3_345_345 =                     << 1702   G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35 
1797     m[A33] * Det2_45_45 - m[A34] * Det2_45_35 << 1703                                                 + m[A35]*Det2_45_34;
1798                                                  1704 
1799   // Find all NECESSARY 4x4 dets:  (55 of the    1705   // Find all NECESSARY 4x4 dets:  (55 of them)
1800                                                  1706 
1801   G4double Det4_1234_0123 = m[A10] * Det3_234 << 1707   G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
1802                             m[A12] * Det3_234 << 1708                         + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1803   G4double Det4_1234_0124 = m[A10] * Det3_234 << 1709   G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
1804                             m[A12] * Det3_234 << 1710                         + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1805   G4double Det4_1234_0134 = m[A10] * Det3_234 << 1711   G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
1806                             m[A13] * Det3_234 << 1712                         + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1807   G4double Det4_1234_0234 = m[A10] * Det3_234 << 1713   G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
1808                             m[A13] * Det3_234 << 1714                         + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1809   G4double Det4_1234_1234 = m[A11] * Det3_234 << 1715   G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
1810                             m[A13] * Det3_234 << 1716                         + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1811   G4double Det4_1235_0123 = m[A10] * Det3_235 << 1717   G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023 
1812                             m[A12] * Det3_235 << 1718                         + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1813   G4double Det4_1235_0124 = m[A10] * Det3_235 << 1719   G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024 
1814                             m[A12] * Det3_235 << 1720                         + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1815   G4double Det4_1235_0125 = m[A10] * Det3_235 << 1721   G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025 
1816                             m[A12] * Det3_235 << 1722                         + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1817   G4double Det4_1235_0134 = m[A10] * Det3_235 << 1723   G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034 
1818                             m[A13] * Det3_235 << 1724                         + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1819   G4double Det4_1235_0135 = m[A10] * Det3_235 << 1725   G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035 
1820                             m[A13] * Det3_235 << 1726                         + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1821   G4double Det4_1235_0234 = m[A10] * Det3_235 << 1727   G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034 
1822                             m[A13] * Det3_235 << 1728                         + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1823   G4double Det4_1235_0235 = m[A10] * Det3_235 << 1729   G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035 
1824                             m[A13] * Det3_235 << 1730                         + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1825   G4double Det4_1235_1234 = m[A11] * Det3_235 << 1731   G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134 
1826                             m[A13] * Det3_235 << 1732                         + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1827   G4double Det4_1235_1235 = m[A11] * Det3_235 << 1733   G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135 
1828                             m[A13] * Det3_235 << 1734                         + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1829   G4double Det4_1245_0123 = m[A10] * Det3_245 << 1735   G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023 
1830                             m[A12] * Det3_245 << 1736                         + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1831   G4double Det4_1245_0124 = m[A10] * Det3_245 << 1737   G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024 
1832                             m[A12] * Det3_245 << 1738                         + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1833   G4double Det4_1245_0125 = m[A10] * Det3_245 << 1739   G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025 
1834                             m[A12] * Det3_245 << 1740                         + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1835   G4double Det4_1245_0134 = m[A10] * Det3_245 << 1741   G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034 
1836                             m[A13] * Det3_245 << 1742                         + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1837   G4double Det4_1245_0135 = m[A10] * Det3_245 << 1743   G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035 
1838                             m[A13] * Det3_245 << 1744                         + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1839   G4double Det4_1245_0145 = m[A10] * Det3_245 << 1745   G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045 
1840                             m[A14] * Det3_245 << 1746                         + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1841   G4double Det4_1245_0234 = m[A10] * Det3_245 << 1747   G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034 
1842                             m[A13] * Det3_245 << 1748                         + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1843   G4double Det4_1245_0235 = m[A10] * Det3_245 << 1749   G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035 
1844                             m[A13] * Det3_245 << 1750                         + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1845   G4double Det4_1245_0245 = m[A10] * Det3_245 << 1751   G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045 
1846                             m[A14] * Det3_245 << 1752                         + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1847   G4double Det4_1245_1234 = m[A11] * Det3_245 << 1753   G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134 
1848                             m[A13] * Det3_245 << 1754                         + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1849   G4double Det4_1245_1235 = m[A11] * Det3_245 << 1755   G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135 
1850                             m[A13] * Det3_245 << 1756                         + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1851   G4double Det4_1245_1245 = m[A11] * Det3_245 << 1757   G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145 
1852                             m[A14] * Det3_245 << 1758                         + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1853   G4double Det4_1345_0123 = m[A10] * Det3_345 << 1759   G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023 
1854                             m[A12] * Det3_345 << 1760                         + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1855   G4double Det4_1345_0124 = m[A10] * Det3_345 << 1761   G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024 
1856                             m[A12] * Det3_345 << 1762                         + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1857   G4double Det4_1345_0125 = m[A10] * Det3_345 << 1763   G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025 
1858                             m[A12] * Det3_345 << 1764                         + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1859   G4double Det4_1345_0134 = m[A10] * Det3_345 << 1765   G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034 
1860                             m[A13] * Det3_345 << 1766                         + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1861   G4double Det4_1345_0135 = m[A10] * Det3_345 << 1767   G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035 
1862                             m[A13] * Det3_345 << 1768                         + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1863   G4double Det4_1345_0145 = m[A10] * Det3_345 << 1769   G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045 
1864                             m[A14] * Det3_345 << 1770                         + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1865   G4double Det4_1345_0234 = m[A10] * Det3_345 << 1771   G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034 
1866                             m[A13] * Det3_345 << 1772                         + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1867   G4double Det4_1345_0235 = m[A10] * Det3_345 << 1773   G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035 
1868                             m[A13] * Det3_345 << 1774                         + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1869   G4double Det4_1345_0245 = m[A10] * Det3_345 << 1775   G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045 
1870                             m[A14] * Det3_345 << 1776                         + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1871   G4double Det4_1345_0345 = m[A10] * Det3_345 << 1777   G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045 
1872                             m[A14] * Det3_345 << 1778                         + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1873   G4double Det4_1345_1234 = m[A11] * Det3_345 << 1779   G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134 
1874                             m[A13] * Det3_345 << 1780                         + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1875   G4double Det4_1345_1235 = m[A11] * Det3_345 << 1781   G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135 
1876                             m[A13] * Det3_345 << 1782                         + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1877   G4double Det4_1345_1245 = m[A11] * Det3_345 << 1783   G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145 
1878                             m[A14] * Det3_345 << 1784                         + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1879   G4double Det4_1345_1345 = m[A11] * Det3_345 << 1785   G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145 
1880                             m[A14] * Det3_345 << 1786                         + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1881   G4double Det4_2345_0123 = m[A20] * Det3_345 << 1787   G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023 
1882                             m[A22] * Det3_345 << 1788                         + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1883   G4double Det4_2345_0124 = m[A20] * Det3_345 << 1789   G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024 
1884                             m[A22] * Det3_345 << 1790                         + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1885   G4double Det4_2345_0125 = m[A20] * Det3_345 << 1791   G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025 
1886                             m[A22] * Det3_345 << 1792                         + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1887   G4double Det4_2345_0134 = m[A20] * Det3_345 << 1793   G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034 
1888                             m[A23] * Det3_345 << 1794                         + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1889   G4double Det4_2345_0135 = m[A20] * Det3_345 << 1795   G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035 
1890                             m[A23] * Det3_345 << 1796                         + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1891   G4double Det4_2345_0145 = m[A20] * Det3_345 << 1797   G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045 
1892                             m[A24] * Det3_345 << 1798                         + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1893   G4double Det4_2345_0234 = m[A20] * Det3_345 << 1799   G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034 
1894                             m[A23] * Det3_345 << 1800                         + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1895   G4double Det4_2345_0235 = m[A20] * Det3_345 << 1801   G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035 
1896                             m[A23] * Det3_345 << 1802                         + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1897   G4double Det4_2345_0245 = m[A20] * Det3_345 << 1803   G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045 
1898                             m[A24] * Det3_345 << 1804                         + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1899   G4double Det4_2345_0345 = m[A20] * Det3_345 << 1805   G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045 
1900                             m[A24] * Det3_345 << 1806                         + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1901   G4double Det4_2345_1234 = m[A21] * Det3_345 << 1807   G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134 
1902                             m[A23] * Det3_345 << 1808                         + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1903   G4double Det4_2345_1235 = m[A21] * Det3_345 << 1809   G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135 
1904                             m[A23] * Det3_345 << 1810                         + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1905   G4double Det4_2345_1245 = m[A21] * Det3_345 << 1811   G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145 
1906                             m[A24] * Det3_345 << 1812                         + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1907   G4double Det4_2345_1345 = m[A21] * Det3_345 << 1813   G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145 
1908                             m[A24] * Det3_345 << 1814                         + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1909   G4double Det4_2345_2345 = m[A22] * Det3_345 << 1815   G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245 
1910                             m[A24] * Det3_345 << 1816                         + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1911                                                  1817 
1912   // Find all NECESSARY 5x5 dets:  (19 of the    1818   // Find all NECESSARY 5x5 dets:  (19 of them)
1913                                                  1819 
1914   G4double Det5_01234_01234 =                 << 1820   G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234 
1915     m[A00] * Det4_1234_1234 - m[A01] * Det4_1 << 1821     + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1916     m[A02] * Det4_1234_0134 - m[A03] * Det4_1 << 1822   G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234 
1917   G4double Det5_01235_01234 =                 << 1823     + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1918     m[A00] * Det4_1235_1234 - m[A01] * Det4_1 << 1824   G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235 
1919     m[A02] * Det4_1235_0134 - m[A03] * Det4_1 << 1825     + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1920   G4double Det5_01235_01235 =                 << 1826   G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234 
1921     m[A00] * Det4_1235_1235 - m[A01] * Det4_1 << 1827     + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1922     m[A02] * Det4_1235_0135 - m[A03] * Det4_1 << 1828   G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235 
1923   G4double Det5_01245_01234 =                 << 1829     + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1924     m[A00] * Det4_1245_1234 - m[A01] * Det4_1 << 1830   G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245 
1925     m[A02] * Det4_1245_0134 - m[A03] * Det4_1 << 1831     + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1926   G4double Det5_01245_01235 =                 << 1832   G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234 
1927     m[A00] * Det4_1245_1235 - m[A01] * Det4_1 << 1833     + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1928     m[A02] * Det4_1245_0135 - m[A03] * Det4_1 << 1834   G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235 
1929   G4double Det5_01245_01245 =                 << 1835     + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1930     m[A00] * Det4_1245_1245 - m[A01] * Det4_1 << 1836   G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245 
1931     m[A02] * Det4_1245_0145 - m[A04] * Det4_1 << 1837     + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1932   G4double Det5_01345_01234 =                 << 1838   G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345 
1933     m[A00] * Det4_1345_1234 - m[A01] * Det4_1 << 1839     + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1934     m[A02] * Det4_1345_0134 - m[A03] * Det4_1 << 1840   G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234 
1935   G4double Det5_01345_01235 =                 << 1841     + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1936     m[A00] * Det4_1345_1235 - m[A01] * Det4_1 << 1842   G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235 
1937     m[A02] * Det4_1345_0135 - m[A03] * Det4_1 << 1843     + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1938   G4double Det5_01345_01245 =                 << 1844   G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245 
1939     m[A00] * Det4_1345_1245 - m[A01] * Det4_1 << 1845     + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1940     m[A02] * Det4_1345_0145 - m[A04] * Det4_1 << 1846   G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345 
1941   G4double Det5_01345_01345 =                 << 1847     + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1942     m[A00] * Det4_1345_1345 - m[A01] * Det4_1 << 1848   G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345 
1943     m[A03] * Det4_1345_0145 - m[A04] * Det4_1 << 1849     + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1944   G4double Det5_02345_01234 =                 << 1850   G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234 
1945     m[A00] * Det4_2345_1234 - m[A01] * Det4_2 << 1851     + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1946     m[A02] * Det4_2345_0134 - m[A03] * Det4_2 << 1852   G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235 
1947   G4double Det5_02345_01235 =                 << 1853     + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1948     m[A00] * Det4_2345_1235 - m[A01] * Det4_2 << 1854   G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245 
1949     m[A02] * Det4_2345_0135 - m[A03] * Det4_2 << 1855     + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1950   G4double Det5_02345_01245 =                 << 1856   G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345 
1951     m[A00] * Det4_2345_1245 - m[A01] * Det4_2 << 1857     + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1952     m[A02] * Det4_2345_0145 - m[A04] * Det4_2 << 1858   G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345 
1953   G4double Det5_02345_01345 =                 << 1859     + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1954     m[A00] * Det4_2345_1345 - m[A01] * Det4_2 << 1860   G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345 
1955     m[A03] * Det4_2345_0145 - m[A04] * Det4_2 << 1861     + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1956   G4double Det5_02345_02345 =                 << 1862 
1957     m[A00] * Det4_2345_2345 - m[A02] * Det4_2 << 1863   // Find the determinant 
1958     m[A03] * Det4_2345_0245 - m[A04] * Det4_2 << 1864 
1959   G4double Det5_12345_01234 =                 << 1865   G4double det =  m[A00]*Det5_12345_12345 
1960     m[A10] * Det4_2345_1234 - m[A11] * Det4_2 << 1866                 - m[A01]*Det5_12345_02345 
1961     m[A12] * Det4_2345_0134 - m[A13] * Det4_2 << 1867                 + m[A02]*Det5_12345_01345 
1962   G4double Det5_12345_01235 =                 << 1868                 - m[A03]*Det5_12345_01245 
1963     m[A10] * Det4_2345_1235 - m[A11] * Det4_2 << 1869                 + m[A04]*Det5_12345_01235 
1964     m[A12] * Det4_2345_0135 - m[A13] * Det4_2 << 1870                 - m[A05]*Det5_12345_01234;
1965   G4double Det5_12345_01245 =                 << 
1966     m[A10] * Det4_2345_1245 - m[A11] * Det4_2 << 
1967     m[A12] * Det4_2345_0145 - m[A14] * Det4_2 << 
1968   G4double Det5_12345_01345 =                 << 
1969     m[A10] * Det4_2345_1345 - m[A11] * Det4_2 << 
1970     m[A13] * Det4_2345_0145 - m[A14] * Det4_2 << 
1971   G4double Det5_12345_02345 =                 << 
1972     m[A10] * Det4_2345_2345 - m[A12] * Det4_2 << 
1973     m[A13] * Det4_2345_0245 - m[A14] * Det4_2 << 
1974   G4double Det5_12345_12345 =                 << 
1975     m[A11] * Det4_2345_2345 - m[A12] * Det4_2 << 
1976     m[A13] * Det4_2345_1245 - m[A14] * Det4_2 << 
1977                                               << 
1978   // Find the determinant                     << 
1979                                               << 
1980   G4double det = m[A00] * Det5_12345_12345 -  << 
1981                  m[A02] * Det5_12345_01345 -  << 
1982                  m[A04] * Det5_12345_01235 -  << 
1983                                                  1871 
1984   if(det == 0)                                << 1872   if ( det == 0 )
1985   {                                           << 1873   {  
1986     ifail = 1;                                   1874     ifail = 1;
1987     return;                                      1875     return;
1988   }                                           << 1876   } 
1989                                                  1877 
1990   G4double oneOverDet = 1.0 / det;            << 1878   G4double oneOverDet = 1.0/det;
1991   G4double mn1OverDet = -oneOverDet;          << 1879   G4double mn1OverDet = - oneOverDet;
1992                                                  1880 
1993   m[A00] = Det5_12345_12345 * oneOverDet;     << 1881   m[A00] =  Det5_12345_12345*oneOverDet;
1994   m[A01] = Det5_12345_02345 * mn1OverDet;     << 1882   m[A01] =  Det5_12345_02345*mn1OverDet;
1995   m[A02] = Det5_12345_01345 * oneOverDet;     << 1883   m[A02] =  Det5_12345_01345*oneOverDet;
1996   m[A03] = Det5_12345_01245 * mn1OverDet;     << 1884   m[A03] =  Det5_12345_01245*mn1OverDet;
1997   m[A04] = Det5_12345_01235 * oneOverDet;     << 1885   m[A04] =  Det5_12345_01235*oneOverDet;
1998   m[A05] = Det5_12345_01234 * mn1OverDet;     << 1886   m[A05] =  Det5_12345_01234*mn1OverDet;
1999                                               << 1887 
2000   m[A11] = Det5_02345_02345 * oneOverDet;     << 1888   m[A11] =  Det5_02345_02345*oneOverDet;
2001   m[A12] = Det5_02345_01345 * mn1OverDet;     << 1889   m[A12] =  Det5_02345_01345*mn1OverDet;
2002   m[A13] = Det5_02345_01245 * oneOverDet;     << 1890   m[A13] =  Det5_02345_01245*oneOverDet;
2003   m[A14] = Det5_02345_01235 * mn1OverDet;     << 1891   m[A14] =  Det5_02345_01235*mn1OverDet;
2004   m[A15] = Det5_02345_01234 * oneOverDet;     << 1892   m[A15] =  Det5_02345_01234*oneOverDet;
2005                                               << 1893 
2006   m[A22] = Det5_01345_01345 * oneOverDet;     << 1894   m[A22] =  Det5_01345_01345*oneOverDet;
2007   m[A23] = Det5_01345_01245 * mn1OverDet;     << 1895   m[A23] =  Det5_01345_01245*mn1OverDet;
2008   m[A24] = Det5_01345_01235 * oneOverDet;     << 1896   m[A24] =  Det5_01345_01235*oneOverDet;
2009   m[A25] = Det5_01345_01234 * mn1OverDet;     << 1897   m[A25] =  Det5_01345_01234*mn1OverDet;
2010                                               << 1898 
2011   m[A33] = Det5_01245_01245 * oneOverDet;     << 1899   m[A33] =  Det5_01245_01245*oneOverDet;
2012   m[A34] = Det5_01245_01235 * mn1OverDet;     << 1900   m[A34] =  Det5_01245_01235*mn1OverDet;
2013   m[A35] = Det5_01245_01234 * oneOverDet;     << 1901   m[A35] =  Det5_01245_01234*oneOverDet;
2014                                                  1902 
2015   m[A44] = Det5_01235_01235 * oneOverDet;     << 1903   m[A44] =  Det5_01235_01235*oneOverDet;
2016   m[A45] = Det5_01235_01234 * mn1OverDet;     << 1904   m[A45] =  Det5_01235_01234*mn1OverDet;
2017                                                  1905 
2018   m[A55] = Det5_01234_01234 * oneOverDet;     << 1906   m[A55] =  Det5_01234_01234*oneOverDet;
2019                                                  1907 
2020   return;                                        1908   return;
2021 }                                                1909 }
2022                                                  1910 
2023 void G4ErrorSymMatrix::invertCholesky5(G4int& << 1911 void G4ErrorSymMatrix::invertCholesky5 (G4int & ifail)
2024 {                                                1912 {
2025   // Invert by                                << 1913 
                                                   >> 1914   // Invert by 
2026   //                                             1915   //
2027   // a) decomposing M = G*G^T with G lower tr    1916   // a) decomposing M = G*G^T with G lower triangular
2028   //    (if M is not positive definite this w    1917   //    (if M is not positive definite this will fail, leaving this unchanged)
2029   // b) inverting G to form H                    1918   // b) inverting G to form H
2030   // c) multiplying H^T * H to get M^-1.         1919   // c) multiplying H^T * H to get M^-1.
2031   //                                             1920   //
2032   // If the matrix is pos. def. it is inverte    1921   // If the matrix is pos. def. it is inverted and 1 is returned.
2033   // If the matrix is not pos. def. it remain    1922   // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2034                                                  1923 
2035   G4double h10;  // below-diagonal elements o << 1924   G4double h10;                           // below-diagonal elements of H
2036   G4double h20, h21;                             1925   G4double h20, h21;
2037   G4double h30, h31, h32;                        1926   G4double h30, h31, h32;
2038   G4double h40, h41, h42, h43;                   1927   G4double h40, h41, h42, h43;
2039                                                  1928 
2040   G4double h00, h11, h22, h33, h44;  // 1/dia << 1929   G4double h00, h11, h22, h33, h44;       // 1/diagonal elements of G = 
2041                                      // diago << 1930                                           // diagonal elements of H
2042                                                  1931 
2043   G4double g10;  // below-diagonal elements o << 1932   G4double g10;                           // below-diagonal elements of G
2044   G4double g20, g21;                             1933   G4double g20, g21;
2045   G4double g30, g31, g32;                        1934   G4double g30, g31, g32;
2046   G4double g40, g41, g42, g43;                   1935   G4double g40, g41, g42, g43;
2047                                                  1936 
2048   ifail = 1;  // We start by assuing we won't    1937   ifail = 1;  // We start by assuing we won't succeed...
2049                                                  1938 
2050   // Form G -- compute diagonal members of H     1939   // Form G -- compute diagonal members of H directly rather than of G
2051   //-------                                      1940   //-------
2052                                                  1941 
2053   // Scale first column by 1/sqrt(A00)           1942   // Scale first column by 1/sqrt(A00)
2054                                                  1943 
2055   h00 = m[A00];                               << 1944   h00 = m[A00]; 
2056   if(h00 <= 0)                                << 1945   if (h00 <= 0) { return; }
2057   {                                           << 
2058     return;                                   << 
2059   }                                           << 
2060   h00 = 1.0 / std::sqrt(h00);                    1946   h00 = 1.0 / std::sqrt(h00);
2061                                                  1947 
2062   g10 = m[A10] * h00;                            1948   g10 = m[A10] * h00;
2063   g20 = m[A20] * h00;                            1949   g20 = m[A20] * h00;
2064   g30 = m[A30] * h00;                            1950   g30 = m[A30] * h00;
2065   g40 = m[A40] * h00;                            1951   g40 = m[A40] * h00;
2066                                                  1952 
2067   // Form G11 (actually, just h11)               1953   // Form G11 (actually, just h11)
2068                                                  1954 
2069   h11 = m[A11] - (g10 * g10);                    1955   h11 = m[A11] - (g10 * g10);
2070   if(h11 <= 0)                                << 1956   if (h11 <= 0) { return; }
2071   {                                           << 
2072     return;                                   << 
2073   }                                           << 
2074   h11 = 1.0 / std::sqrt(h11);                    1957   h11 = 1.0 / std::sqrt(h11);
2075                                                  1958 
2076   // Subtract inter-column column dot product    1959   // Subtract inter-column column dot products from rest of column 1 and
2077   // scale to get column 1 of G               << 1960   // scale to get column 1 of G 
2078                                                  1961 
2079   g21 = (m[A21] - (g10 * g20)) * h11;            1962   g21 = (m[A21] - (g10 * g20)) * h11;
2080   g31 = (m[A31] - (g10 * g30)) * h11;            1963   g31 = (m[A31] - (g10 * g30)) * h11;
2081   g41 = (m[A41] - (g10 * g40)) * h11;            1964   g41 = (m[A41] - (g10 * g40)) * h11;
2082                                                  1965 
2083   // Form G22 (actually, just h22)               1966   // Form G22 (actually, just h22)
2084                                                  1967 
2085   h22 = m[A22] - (g20 * g20) - (g21 * g21);      1968   h22 = m[A22] - (g20 * g20) - (g21 * g21);
2086   if(h22 <= 0)                                << 1969   if (h22 <= 0) { return; }
2087   {                                           << 
2088     return;                                   << 
2089   }                                           << 
2090   h22 = 1.0 / std::sqrt(h22);                    1970   h22 = 1.0 / std::sqrt(h22);
2091                                                  1971 
2092   // Subtract inter-column column dot product    1972   // Subtract inter-column column dot products from rest of column 2 and
2093   // scale to get column 2 of G               << 1973   // scale to get column 2 of G 
2094                                                  1974 
2095   g32 = (m[A32] - (g20 * g30) - (g21 * g31))     1975   g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2096   g42 = (m[A42] - (g20 * g40) - (g21 * g41))     1976   g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2097                                                  1977 
2098   // Form G33 (actually, just h33)               1978   // Form G33 (actually, just h33)
2099                                                  1979 
2100   h33 = m[A33] - (g30 * g30) - (g31 * g31) -     1980   h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2101   if(h33 <= 0)                                << 1981   if (h33 <= 0) { return; }
2102   {                                           << 
2103     return;                                   << 
2104   }                                           << 
2105   h33 = 1.0 / std::sqrt(h33);                    1982   h33 = 1.0 / std::sqrt(h33);
2106                                                  1983 
2107   // Subtract inter-column column dot product    1984   // Subtract inter-column column dot product from A43 and scale to get G43
2108                                                  1985 
2109   g43 = (m[A43] - (g30 * g40) - (g31 * g41) -    1986   g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2110                                                  1987 
2111   // Finally form h44 - if this is possible i    1988   // Finally form h44 - if this is possible inversion succeeds
2112                                                  1989 
2113   h44 = m[A44] - (g40 * g40) - (g41 * g41) -     1990   h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2114   if(h44 <= 0)                                << 1991   if (h44 <= 0) { return; }
2115   {                                           << 
2116     return;                                   << 
2117   }                                           << 
2118   h44 = 1.0 / std::sqrt(h44);                    1992   h44 = 1.0 / std::sqrt(h44);
2119                                                  1993 
2120   // Form H = 1/G -- diagonal members of H ar    1994   // Form H = 1/G -- diagonal members of H are already correct
2121   //-------------                                1995   //-------------
2122                                                  1996 
2123   // The order here is dictated by speed cons    1997   // The order here is dictated by speed considerations
2124                                                  1998 
2125   h43 = -h33 * g43 * h44;                     << 1999   h43 = -h33 *  g43 * h44;
2126   h32 = -h22 * g32 * h33;                     << 2000   h32 = -h22 *  g32 * h33;
2127   h42 = -h22 * (g32 * h43 + g42 * h44);          2001   h42 = -h22 * (g32 * h43 + g42 * h44);
2128   h21 = -h11 * g21 * h22;                     << 2002   h21 = -h11 *  g21 * h22;
2129   h31 = -h11 * (g21 * h32 + g31 * h33);          2003   h31 = -h11 * (g21 * h32 + g31 * h33);
2130   h41 = -h11 * (g21 * h42 + g31 * h43 + g41 *    2004   h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2131   h10 = -h00 * g10 * h11;                     << 2005   h10 = -h00 *  g10 * h11;
2132   h20 = -h00 * (g10 * h21 + g20 * h22);          2006   h20 = -h00 * (g10 * h21 + g20 * h22);
2133   h30 = -h00 * (g10 * h31 + g20 * h32 + g30 *    2007   h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2134   h40 = -h00 * (g10 * h41 + g20 * h42 + g30 *    2008   h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2135                                                  2009 
2136   // Change this to its inverse = H^T*H          2010   // Change this to its inverse = H^T*H
2137   //------------------------------------         2011   //------------------------------------
2138                                                  2012 
2139   m[A00] = h00 * h00 + h10 * h10 + h20 * h20     2013   m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
2140   m[A01] = h10 * h11 + h20 * h21 + h30 * h31     2014   m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
2141   m[A11] = h11 * h11 + h21 * h21 + h31 * h31     2015   m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
2142   m[A02] = h20 * h22 + h30 * h32 + h40 * h42;    2016   m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
2143   m[A12] = h21 * h22 + h31 * h32 + h41 * h42;    2017   m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
2144   m[A22] = h22 * h22 + h32 * h32 + h42 * h42;    2018   m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
2145   m[A03] = h30 * h33 + h40 * h43;                2019   m[A03] = h30 * h33 + h40 * h43;
2146   m[A13] = h31 * h33 + h41 * h43;                2020   m[A13] = h31 * h33 + h41 * h43;
2147   m[A23] = h32 * h33 + h42 * h43;                2021   m[A23] = h32 * h33 + h42 * h43;
2148   m[A33] = h33 * h33 + h43 * h43;                2022   m[A33] = h33 * h33 + h43 * h43;
2149   m[A04] = h40 * h44;                            2023   m[A04] = h40 * h44;
2150   m[A14] = h41 * h44;                            2024   m[A14] = h41 * h44;
2151   m[A24] = h42 * h44;                            2025   m[A24] = h42 * h44;
2152   m[A34] = h43 * h44;                            2026   m[A34] = h43 * h44;
2153   m[A44] = h44 * h44;                            2027   m[A44] = h44 * h44;
2154                                                  2028 
2155   ifail = 0;                                     2029   ifail = 0;
2156   return;                                        2030   return;
2157 }                                                2031 }
2158                                                  2032 
2159 void G4ErrorSymMatrix::invertCholesky6(G4int& << 2033 void G4ErrorSymMatrix::invertCholesky6 (G4int & ifail)
2160 {                                                2034 {
2161   // Invert by                                << 2035   // Invert by 
2162   //                                             2036   //
2163   // a) decomposing M = G*G^T with G lower tr    2037   // a) decomposing M = G*G^T with G lower triangular
2164   //    (if M is not positive definite this w    2038   //    (if M is not positive definite this will fail, leaving this unchanged)
2165   // b) inverting G to form H                    2039   // b) inverting G to form H
2166   // c) multiplying H^T * H to get M^-1.         2040   // c) multiplying H^T * H to get M^-1.
2167   //                                             2041   //
2168   // If the matrix is pos. def. it is inverte    2042   // If the matrix is pos. def. it is inverted and 1 is returned.
2169   // If the matrix is not pos. def. it remain    2043   // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2170                                                  2044 
2171   G4double h10;  // below-diagonal elements o << 2045   G4double h10;                           // below-diagonal elements of H
2172   G4double h20, h21;                             2046   G4double h20, h21;
2173   G4double h30, h31, h32;                        2047   G4double h30, h31, h32;
2174   G4double h40, h41, h42, h43;                   2048   G4double h40, h41, h42, h43;
2175   G4double h50, h51, h52, h53, h54;              2049   G4double h50, h51, h52, h53, h54;
2176                                                  2050 
2177   G4double h00, h11, h22, h33, h44, h55;  //  << 2051   G4double h00, h11, h22, h33, h44, h55;  // 1/diagonal elements of G = 
2178                                           //     2052                                           // diagonal elements of H
2179                                                  2053 
2180   G4double g10;  // below-diagonal elements o << 2054   G4double g10;                           // below-diagonal elements of G
2181   G4double g20, g21;                             2055   G4double g20, g21;
2182   G4double g30, g31, g32;                        2056   G4double g30, g31, g32;
2183   G4double g40, g41, g42, g43;                   2057   G4double g40, g41, g42, g43;
2184   G4double g50, g51, g52, g53, g54;              2058   G4double g50, g51, g52, g53, g54;
2185                                                  2059 
2186   ifail = 1;  // We start by assuing we won't    2060   ifail = 1;  // We start by assuing we won't succeed...
2187                                                  2061 
2188   // Form G -- compute diagonal members of H     2062   // Form G -- compute diagonal members of H directly rather than of G
2189   //-------                                      2063   //-------
2190                                                  2064 
2191   // Scale first column by 1/sqrt(A00)           2065   // Scale first column by 1/sqrt(A00)
2192                                                  2066 
2193   h00 = m[A00];                               << 2067   h00 = m[A00]; 
2194   if(h00 <= 0)                                << 2068   if (h00 <= 0) { return; }
2195   {                                           << 
2196     return;                                   << 
2197   }                                           << 
2198   h00 = 1.0 / std::sqrt(h00);                    2069   h00 = 1.0 / std::sqrt(h00);
2199                                                  2070 
2200   g10 = m[A10] * h00;                            2071   g10 = m[A10] * h00;
2201   g20 = m[A20] * h00;                            2072   g20 = m[A20] * h00;
2202   g30 = m[A30] * h00;                            2073   g30 = m[A30] * h00;
2203   g40 = m[A40] * h00;                            2074   g40 = m[A40] * h00;
2204   g50 = m[A50] * h00;                            2075   g50 = m[A50] * h00;
2205                                                  2076 
2206   // Form G11 (actually, just h11)               2077   // Form G11 (actually, just h11)
2207                                                  2078 
2208   h11 = m[A11] - (g10 * g10);                    2079   h11 = m[A11] - (g10 * g10);
2209   if(h11 <= 0)                                << 2080   if (h11 <= 0) { return; }
2210   {                                           << 
2211     return;                                   << 
2212   }                                           << 
2213   h11 = 1.0 / std::sqrt(h11);                    2081   h11 = 1.0 / std::sqrt(h11);
2214                                                  2082 
2215   // Subtract inter-column column dot product    2083   // Subtract inter-column column dot products from rest of column 1 and
2216   // scale to get column 1 of G               << 2084   // scale to get column 1 of G 
2217                                                  2085 
2218   g21 = (m[A21] - (g10 * g20)) * h11;            2086   g21 = (m[A21] - (g10 * g20)) * h11;
2219   g31 = (m[A31] - (g10 * g30)) * h11;            2087   g31 = (m[A31] - (g10 * g30)) * h11;
2220   g41 = (m[A41] - (g10 * g40)) * h11;            2088   g41 = (m[A41] - (g10 * g40)) * h11;
2221   g51 = (m[A51] - (g10 * g50)) * h11;            2089   g51 = (m[A51] - (g10 * g50)) * h11;
2222                                                  2090 
2223   // Form G22 (actually, just h22)               2091   // Form G22 (actually, just h22)
2224                                                  2092 
2225   h22 = m[A22] - (g20 * g20) - (g21 * g21);      2093   h22 = m[A22] - (g20 * g20) - (g21 * g21);
2226   if(h22 <= 0)                                << 2094   if (h22 <= 0) { return; }
2227   {                                           << 
2228     return;                                   << 
2229   }                                           << 
2230   h22 = 1.0 / std::sqrt(h22);                    2095   h22 = 1.0 / std::sqrt(h22);
2231                                                  2096 
2232   // Subtract inter-column column dot product    2097   // Subtract inter-column column dot products from rest of column 2 and
2233   // scale to get column 2 of G               << 2098   // scale to get column 2 of G 
2234                                                  2099 
2235   g32 = (m[A32] - (g20 * g30) - (g21 * g31))     2100   g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2236   g42 = (m[A42] - (g20 * g40) - (g21 * g41))     2101   g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2237   g52 = (m[A52] - (g20 * g50) - (g21 * g51))     2102   g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
2238                                                  2103 
2239   // Form G33 (actually, just h33)               2104   // Form G33 (actually, just h33)
2240                                                  2105 
2241   h33 = m[A33] - (g30 * g30) - (g31 * g31) -     2106   h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2242   if(h33 <= 0)                                << 2107   if (h33 <= 0) { return; }
2243   {                                           << 
2244     return;                                   << 
2245   }                                           << 
2246   h33 = 1.0 / std::sqrt(h33);                    2108   h33 = 1.0 / std::sqrt(h33);
2247                                                  2109 
2248   // Subtract inter-column column dot product    2110   // Subtract inter-column column dot products from rest of column 3 and
2249   // scale to get column 3 of G               << 2111   // scale to get column 3 of G 
2250                                                  2112 
2251   g43 = (m[A43] - (g30 * g40) - (g31 * g41) -    2113   g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2252   g53 = (m[A53] - (g30 * g50) - (g31 * g51) -    2114   g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
2253                                                  2115 
2254   // Form G44 (actually, just h44)               2116   // Form G44 (actually, just h44)
2255                                                  2117 
2256   h44 = m[A44] - (g40 * g40) - (g41 * g41) -     2118   h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2257   if(h44 <= 0)                                << 2119   if (h44 <= 0) { return; }
2258   {                                           << 
2259     return;                                   << 
2260   }                                           << 
2261   h44 = 1.0 / std::sqrt(h44);                    2120   h44 = 1.0 / std::sqrt(h44);
2262                                                  2121 
2263   // Subtract inter-column column dot product    2122   // Subtract inter-column column dot product from M54 and scale to get G54
2264                                                  2123 
2265   g54 = (m[A54] - (g40 * g50) - (g41 * g51) -    2124   g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
2266                                                  2125 
2267   // Finally form h55 - if this is possible i    2126   // Finally form h55 - if this is possible inversion succeeds
2268                                                  2127 
2269   h55 = m[A55] - (g50 * g50) - (g51 * g51) -  << 2128   h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
2270         (g54 * g54);                          << 2129   if (h55 <= 0) { return; }
2271   if(h55 <= 0)                                << 
2272   {                                           << 
2273     return;                                   << 
2274   }                                           << 
2275   h55 = 1.0 / std::sqrt(h55);                    2130   h55 = 1.0 / std::sqrt(h55);
2276                                                  2131 
2277   // Form H = 1/G -- diagonal members of H ar    2132   // Form H = 1/G -- diagonal members of H are already correct
2278   //-------------                                2133   //-------------
2279                                                  2134 
2280   // The order here is dictated by speed cons    2135   // The order here is dictated by speed considerations
2281                                                  2136 
2282   h54 = -h44 * g54 * h55;                     << 2137   h54 = -h44 *  g54 * h55;
2283   h43 = -h33 * g43 * h44;                     << 2138   h43 = -h33 *  g43 * h44;
2284   h53 = -h33 * (g43 * h54 + g53 * h55);          2139   h53 = -h33 * (g43 * h54 + g53 * h55);
2285   h32 = -h22 * g32 * h33;                     << 2140   h32 = -h22 *  g32 * h33;
2286   h42 = -h22 * (g32 * h43 + g42 * h44);          2141   h42 = -h22 * (g32 * h43 + g42 * h44);
2287   h52 = -h22 * (g32 * h53 + g42 * h54 + g52 *    2142   h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
2288   h21 = -h11 * g21 * h22;                     << 2143   h21 = -h11 *  g21 * h22;
2289   h31 = -h11 * (g21 * h32 + g31 * h33);          2144   h31 = -h11 * (g21 * h32 + g31 * h33);
2290   h41 = -h11 * (g21 * h42 + g31 * h43 + g41 *    2145   h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2291   h51 = -h11 * (g21 * h52 + g31 * h53 + g41 *    2146   h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
2292   h10 = -h00 * g10 * h11;                     << 2147   h10 = -h00 *  g10 * h11;
2293   h20 = -h00 * (g10 * h21 + g20 * h22);          2148   h20 = -h00 * (g10 * h21 + g20 * h22);
2294   h30 = -h00 * (g10 * h31 + g20 * h32 + g30 *    2149   h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2295   h40 = -h00 * (g10 * h41 + g20 * h42 + g30 *    2150   h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2296   h50 = -h00 * (g10 * h51 + g20 * h52 + g30 *    2151   h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
2297                                                  2152 
2298   // Change this to its inverse = H^T*H          2153   // Change this to its inverse = H^T*H
2299   //------------------------------------         2154   //------------------------------------
2300                                                  2155 
2301   m[A00] =                                    << 2156   m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
2302     h00 * h00 + h10 * h10 + h20 * h20 + h30 * << 
2303   m[A01] = h10 * h11 + h20 * h21 + h30 * h31     2157   m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
2304   m[A11] = h11 * h11 + h21 * h21 + h31 * h31     2158   m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
2305   m[A02] = h20 * h22 + h30 * h32 + h40 * h42     2159   m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
2306   m[A12] = h21 * h22 + h31 * h32 + h41 * h42     2160   m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
2307   m[A22] = h22 * h22 + h32 * h32 + h42 * h42     2161   m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
2308   m[A03] = h30 * h33 + h40 * h43 + h50 * h53;    2162   m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
2309   m[A13] = h31 * h33 + h41 * h43 + h51 * h53;    2163   m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
2310   m[A23] = h32 * h33 + h42 * h43 + h52 * h53;    2164   m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
2311   m[A33] = h33 * h33 + h43 * h43 + h53 * h53;    2165   m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
2312   m[A04] = h40 * h44 + h50 * h54;                2166   m[A04] = h40 * h44 + h50 * h54;
2313   m[A14] = h41 * h44 + h51 * h54;                2167   m[A14] = h41 * h44 + h51 * h54;
2314   m[A24] = h42 * h44 + h52 * h54;                2168   m[A24] = h42 * h44 + h52 * h54;
2315   m[A34] = h43 * h44 + h53 * h54;                2169   m[A34] = h43 * h44 + h53 * h54;
2316   m[A44] = h44 * h44 + h54 * h54;                2170   m[A44] = h44 * h44 + h54 * h54;
2317   m[A05] = h50 * h55;                            2171   m[A05] = h50 * h55;
2318   m[A15] = h51 * h55;                            2172   m[A15] = h51 * h55;
2319   m[A25] = h52 * h55;                            2173   m[A25] = h52 * h55;
2320   m[A35] = h53 * h55;                            2174   m[A35] = h53 * h55;
2321   m[A45] = h54 * h55;                            2175   m[A45] = h54 * h55;
2322   m[A55] = h55 * h55;                            2176   m[A55] = h55 * h55;
2323                                                  2177 
2324   ifail = 0;                                     2178   ifail = 0;
2325   return;                                        2179   return;
2326 }                                                2180 }
2327                                                  2181 
2328 void G4ErrorSymMatrix::invert4(G4int& ifail)  << 2182 void G4ErrorSymMatrix::invert4  (G4int & ifail)
2329 {                                                2183 {
2330   ifail = 0;                                     2184   ifail = 0;
2331                                                  2185 
2332   // Find all NECESSARY 2x2 dets:  (14 of the    2186   // Find all NECESSARY 2x2 dets:  (14 of them)
2333                                                  2187 
2334   G4double Det2_12_01 = m[A10] * m[A21] - m[A << 2188   G4double Det2_12_01 = m[A10]*m[A21] - m[A11]*m[A20];
2335   G4double Det2_12_02 = m[A10] * m[A22] - m[A << 2189   G4double Det2_12_02 = m[A10]*m[A22] - m[A12]*m[A20];
2336   G4double Det2_12_12 = m[A11] * m[A22] - m[A << 2190   G4double Det2_12_12 = m[A11]*m[A22] - m[A12]*m[A21];
2337   G4double Det2_13_01 = m[A10] * m[A31] - m[A << 2191   G4double Det2_13_01 = m[A10]*m[A31] - m[A11]*m[A30];
2338   G4double Det2_13_02 = m[A10] * m[A32] - m[A << 2192   G4double Det2_13_02 = m[A10]*m[A32] - m[A12]*m[A30];
2339   G4double Det2_13_03 = m[A10] * m[A33] - m[A << 2193   G4double Det2_13_03 = m[A10]*m[A33] - m[A13]*m[A30];
2340   G4double Det2_13_12 = m[A11] * m[A32] - m[A << 2194   G4double Det2_13_12 = m[A11]*m[A32] - m[A12]*m[A31];
2341   G4double Det2_13_13 = m[A11] * m[A33] - m[A << 2195   G4double Det2_13_13 = m[A11]*m[A33] - m[A13]*m[A31];
2342   G4double Det2_23_01 = m[A20] * m[A31] - m[A << 2196   G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
2343   G4double Det2_23_02 = m[A20] * m[A32] - m[A << 2197   G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
2344   G4double Det2_23_03 = m[A20] * m[A33] - m[A << 2198   G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
2345   G4double Det2_23_12 = m[A21] * m[A32] - m[A << 2199   G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
2346   G4double Det2_23_13 = m[A21] * m[A33] - m[A << 2200   G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
2347   G4double Det2_23_23 = m[A22] * m[A33] - m[A << 2201   G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
2348                                                  2202 
2349   // Find all NECESSARY 3x3 dets:   (10 of th    2203   // Find all NECESSARY 3x3 dets:   (10 of them)
2350                                                  2204 
2351   G4double Det3_012_012 =                     << 2205   G4double Det3_012_012 = m[A00]*Det2_12_12 - m[A01]*Det2_12_02 
2352     m[A00] * Det2_12_12 - m[A01] * Det2_12_02 << 2206                                 + m[A02]*Det2_12_01;
2353   G4double Det3_013_012 =                     << 2207   G4double Det3_013_012 = m[A00]*Det2_13_12 - m[A01]*Det2_13_02 
2354     m[A00] * Det2_13_12 - m[A01] * Det2_13_02 << 2208                                 + m[A02]*Det2_13_01;
2355   G4double Det3_013_013 =                     << 2209   G4double Det3_013_013 = m[A00]*Det2_13_13 - m[A01]*Det2_13_03
2356     m[A00] * Det2_13_13 - m[A01] * Det2_13_03 << 2210                                 + m[A03]*Det2_13_01;
2357   G4double Det3_023_012 =                     << 2211   G4double Det3_023_012 = m[A00]*Det2_23_12 - m[A01]*Det2_23_02 
2358     m[A00] * Det2_23_12 - m[A01] * Det2_23_02 << 2212                                 + m[A02]*Det2_23_01;
2359   G4double Det3_023_013 =                     << 2213   G4double Det3_023_013 = m[A00]*Det2_23_13 - m[A01]*Det2_23_03
2360     m[A00] * Det2_23_13 - m[A01] * Det2_23_03 << 2214                                 + m[A03]*Det2_23_01;
2361   G4double Det3_023_023 =                     << 2215   G4double Det3_023_023 = m[A00]*Det2_23_23 - m[A02]*Det2_23_03
2362     m[A00] * Det2_23_23 - m[A02] * Det2_23_03 << 2216                                 + m[A03]*Det2_23_02;
2363   G4double Det3_123_012 =                     << 2217   G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02 
2364     m[A10] * Det2_23_12 - m[A11] * Det2_23_02 << 2218                                 + m[A12]*Det2_23_01;
2365   G4double Det3_123_013 =                     << 2219   G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03 
2366     m[A10] * Det2_23_13 - m[A11] * Det2_23_03 << 2220                                 + m[A13]*Det2_23_01;
2367   G4double Det3_123_023 =                     << 2221   G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03 
2368     m[A10] * Det2_23_23 - m[A12] * Det2_23_03 << 2222                                 + m[A13]*Det2_23_02;
2369   G4double Det3_123_123 =                     << 2223   G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13 
2370     m[A11] * Det2_23_23 - m[A12] * Det2_23_13 << 2224                                 + m[A13]*Det2_23_12;
2371                                                  2225 
2372   // Find the 4x4 det:                           2226   // Find the 4x4 det:
2373                                                  2227 
2374   G4double det = m[A00] * Det3_123_123 - m[A0 << 2228   G4double det =  m[A00]*Det3_123_123 
2375                  m[A02] * Det3_123_013 - m[A0 << 2229                 - m[A01]*Det3_123_023 
                                                   >> 2230                 + m[A02]*Det3_123_013 
                                                   >> 2231                 - m[A03]*Det3_123_012;
2376                                                  2232 
2377   if(det == 0)                                << 2233   if ( det == 0 )
2378   {                                           << 2234   {  
2379     ifail = 1;                                   2235     ifail = 1;
2380     return;                                      2236     return;
2381   }                                           << 2237   } 
                                                   >> 2238 
                                                   >> 2239   G4double oneOverDet = 1.0/det;
                                                   >> 2240   G4double mn1OverDet = - oneOverDet;
2382                                                  2241 
2383   G4double oneOverDet = 1.0 / det;            << 2242   m[A00] =  Det3_123_123 * oneOverDet;
2384   G4double mn1OverDet = -oneOverDet;          << 2243   m[A01] =  Det3_123_023 * mn1OverDet;
                                                   >> 2244   m[A02] =  Det3_123_013 * oneOverDet;
                                                   >> 2245   m[A03] =  Det3_123_012 * mn1OverDet;
2385                                                  2246 
2386   m[A00] = Det3_123_123 * oneOverDet;         << 
2387   m[A01] = Det3_123_023 * mn1OverDet;         << 
2388   m[A02] = Det3_123_013 * oneOverDet;         << 
2389   m[A03] = Det3_123_012 * mn1OverDet;         << 
2390                                               << 
2391   m[A11] = Det3_023_023 * oneOverDet;         << 
2392   m[A12] = Det3_023_013 * mn1OverDet;         << 
2393   m[A13] = Det3_023_012 * oneOverDet;         << 
2394                                                  2247 
2395   m[A22] = Det3_013_013 * oneOverDet;         << 2248   m[A11] =  Det3_023_023 * oneOverDet;
2396   m[A23] = Det3_013_012 * mn1OverDet;         << 2249   m[A12] =  Det3_023_013 * mn1OverDet;
                                                   >> 2250   m[A13] =  Det3_023_012 * oneOverDet;
2397                                                  2251 
2398   m[A33] = Det3_012_012 * oneOverDet;         << 2252   m[A22] =  Det3_013_013 * oneOverDet;
                                                   >> 2253   m[A23] =  Det3_013_012 * mn1OverDet;
                                                   >> 2254 
                                                   >> 2255   m[A33] =  Det3_012_012 * oneOverDet;
2399                                                  2256 
2400   return;                                        2257   return;
2401 }                                                2258 }
2402                                                  2259 
2403 void G4ErrorSymMatrix::invertHaywood4(G4int&  << 2260 void G4ErrorSymMatrix::invertHaywood4  (G4int & ifail)
2404 {                                                2261 {
2405   invert4(ifail);  // For the 4x4 case, the m << 2262   invert4(ifail); // For the 4x4 case, the method we use for invert is already
2406                    // the Haywood method.     << 2263                   // the Haywood method.
2407 }                                                2264 }
                                                   >> 2265 
2408                                                  2266