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