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