Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ErrorMatrix.cc
Go to the documentation of this file.
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
78G4ErrorMatrix::G4ErrorMatrix(G4int p, G4int q)
79 : m(p * q)
80 , nrow(p)
81 , ncol(q)
82{
83 size = nrow * ncol;
84}
85
86G4ErrorMatrix::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//
124G4ErrorMatrix::~G4ErrorMatrix() {}
125
126G4ErrorMatrix::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
135G4ErrorMatrix::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
170G4ErrorMatrix 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
194void 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
220G4ErrorMatrix 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
234G4ErrorMatrix 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
245G4ErrorMatrix 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
257G4ErrorMatrix 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
270G4ErrorMatrix operator/(const G4ErrorMatrix& mat1, G4double t)
271{
272 G4ErrorMatrix mret(mat1);
273 mret /= t;
274 return mret;
275}
276
277G4ErrorMatrix operator*(const G4ErrorMatrix& mat1, G4double t)
278{
279 G4ErrorMatrix mret(mat1);
280 mret *= t;
281 return mret;
282}
283
284G4ErrorMatrix operator*(G4double t, const G4ErrorMatrix& mat1)
285{
286 G4ErrorMatrix mret(mat1);
287 mret *= t;
288 return mret;
289}
290
291G4ErrorMatrix 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
325G4ErrorMatrix& 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
332G4ErrorMatrix& 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
339G4ErrorMatrix& G4ErrorMatrix::operator/=(G4double t)
340{
341 SIMPLE_UOP(/=)
342 return (*this);
343}
344
345G4ErrorMatrix& G4ErrorMatrix::operator*=(G4double t)
346{
347 SIMPLE_UOP(*=)
348 return (*this);
349}
350
351G4ErrorMatrix& 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
371std::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
399G4ErrorMatrix 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
417G4ErrorMatrix 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
432G4int 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
552G4int 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
684void 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
821G4double 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
845G4double 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
855void 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
959void 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
1056void 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
1277void 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}
epsilon
G4double epsilon(G4double density, G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:95
A41
#define A41
Definition: G4ErrorMatrix.cc:896
F01
#define F01
Definition: G4ErrorMatrix.cc:940
F22
#define F22
Definition: G4ErrorMatrix.cc:951
F00
#define F00
Definition: G4ErrorMatrix.cc:939
A20
#define A20
Definition: G4ErrorMatrix.cc:881
A01
#define A01
Definition: G4ErrorMatrix.cc:868
A53
#define A53
Definition: G4ErrorMatrix.cc:905
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: G4ErrorMatrix.cc:64
M00
#define M00
Definition: G4ErrorMatrix.cc:909
F02
#define F02
Definition: G4ErrorMatrix.cc:941
A23
#define A23
Definition: G4ErrorMatrix.cc:884
M34
#define M34
Definition: G4ErrorMatrix.cc:931
A45
#define A45
Definition: G4ErrorMatrix.cc:900
A13
#define A13
Definition: G4ErrorMatrix.cc:877
A00
#define A00
Definition: G4ErrorMatrix.cc:867
M03
#define M03
Definition: G4ErrorMatrix.cc:912
A54
#define A54
Definition: G4ErrorMatrix.cc:906
dsum
G4ErrorMatrix dsum(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
Definition: G4ErrorMatrix.cc:220
F31
#define F31
Definition: G4ErrorMatrix.cc:955
A40
#define A40
Definition: G4ErrorMatrix.cc:895
F32
#define F32
Definition: G4ErrorMatrix.cc:956
M10
#define M10
Definition: G4ErrorMatrix.cc:915
F21
#define F21
Definition: G4ErrorMatrix.cc:950
F03
#define F03
Definition: G4ErrorMatrix.cc:942
M41
#define M41
Definition: G4ErrorMatrix.cc:934
M20
#define M20
Definition: G4ErrorMatrix.cc:921
M42
#define M42
Definition: G4ErrorMatrix.cc:935
A25
#define A25
Definition: G4ErrorMatrix.cc:886
M13
#define M13
Definition: G4ErrorMatrix.cc:918
operator*
G4ErrorMatrix operator*(const G4ErrorMatrix &mat1, G4double t)
Definition: G4ErrorMatrix.cc:277
A02
#define A02
Definition: G4ErrorMatrix.cc:869
F30
#define F30
Definition: G4ErrorMatrix.cc:954
operator<<
std::ostream & operator<<(std::ostream &os, const G4ErrorMatrix &q)
Definition: G4ErrorMatrix.cc:371
A24
#define A24
Definition: G4ErrorMatrix.cc:885
M33
#define M33
Definition: G4ErrorMatrix.cc:930
A22
#define A22
Definition: G4ErrorMatrix.cc:883
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: G4ErrorMatrix.cc:47
M04
#define M04
Definition: G4ErrorMatrix.cc:913
A04
#define A04
Definition: G4ErrorMatrix.cc:871
operator-
G4ErrorMatrix operator-(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
Definition: G4ErrorMatrix.cc:257
A30
#define A30
Definition: G4ErrorMatrix.cc:888
M23
#define M23
Definition: G4ErrorMatrix.cc:924
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorMatrix.cc:41
A12
#define A12
Definition: G4ErrorMatrix.cc:876
A55
#define A55
Definition: G4ErrorMatrix.cc:907
M11
#define M11
Definition: G4ErrorMatrix.cc:916
A35
#define A35
Definition: G4ErrorMatrix.cc:893
M44
#define M44
Definition: G4ErrorMatrix.cc:937
A50
#define A50
Definition: G4ErrorMatrix.cc:902
A03
#define A03
Definition: G4ErrorMatrix.cc:870
A14
#define A14
Definition: G4ErrorMatrix.cc:878
F11
#define F11
Definition: G4ErrorMatrix.cc:945
M31
#define M31
Definition: G4ErrorMatrix.cc:928
F10
#define F10
Definition: G4ErrorMatrix.cc:944
A51
#define A51
Definition: G4ErrorMatrix.cc:903
M02
#define M02
Definition: G4ErrorMatrix.cc:911
M21
#define M21
Definition: G4ErrorMatrix.cc:922
A21
#define A21
Definition: G4ErrorMatrix.cc:882
M01
#define M01
Definition: G4ErrorMatrix.cc:910
A11
#define A11
Definition: G4ErrorMatrix.cc:875
A10
#define A10
Definition: G4ErrorMatrix.cc:874
A44
#define A44
Definition: G4ErrorMatrix.cc:899
A05
#define A05
Definition: G4ErrorMatrix.cc:872
A32
#define A32
Definition: G4ErrorMatrix.cc:890
A31
#define A31
Definition: G4ErrorMatrix.cc:889
F33
#define F33
Definition: G4ErrorMatrix.cc:957
M24
#define M24
Definition: G4ErrorMatrix.cc:925
F20
#define F20
Definition: G4ErrorMatrix.cc:949
A33
#define A33
Definition: G4ErrorMatrix.cc:891
M12
#define M12
Definition: G4ErrorMatrix.cc:917
operator/
G4ErrorMatrix operator/(const G4ErrorMatrix &mat1, G4double t)
Definition: G4ErrorMatrix.cc:270
M22
#define M22
Definition: G4ErrorMatrix.cc:923
A42
#define A42
Definition: G4ErrorMatrix.cc:897
operator+
G4ErrorMatrix operator+(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
Definition: G4ErrorMatrix.cc:245
F12
#define F12
Definition: G4ErrorMatrix.cc:946
A52
#define A52
Definition: G4ErrorMatrix.cc:904
M43
#define M43
Definition: G4ErrorMatrix.cc:936
F23
#define F23
Definition: G4ErrorMatrix.cc:952
A15
#define A15
Definition: G4ErrorMatrix.cc:879
A34
#define A34
Definition: G4ErrorMatrix.cc:892
M14
#define M14
Definition: G4ErrorMatrix.cc:919
M30
#define M30
Definition: G4ErrorMatrix.cc:927
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: G4ErrorMatrix.cc:54
M40
#define M40
Definition: G4ErrorMatrix.cc:933
A43
#define A43
Definition: G4ErrorMatrix.cc:898
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorMatrix.cc:70
M32
#define M32
Definition: G4ErrorMatrix.cc:929
F13
#define F13
Definition: G4ErrorMatrix.cc:947
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:43
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4ErrorMatrix::apply
G4ErrorMatrix apply(G4double(*f)(G4double, G4int, G4int)) const
Definition: G4ErrorMatrix.cc:417
G4ErrorMatrix::invertHaywood4
virtual void invertHaywood4(G4int &ierr)
Definition: G4ErrorMatrix.cc:959
G4ErrorMatrix::operator-
G4ErrorMatrix operator-() const
Definition: G4ErrorMatrix.cc:234
G4ErrorMatrix::invert
virtual void invert(G4int &ierr)
Definition: G4ErrorMatrix.cc:684
G4ErrorMatrix::operator/=
G4ErrorMatrix & operator/=(G4double t)
Definition: G4ErrorMatrix.cc:339
G4ErrorMatrix::T
G4ErrorMatrix T() const
Definition: G4ErrorMatrix.cc:399
G4ErrorMatrix::~G4ErrorMatrix
virtual ~G4ErrorMatrix()
Definition: G4ErrorMatrix.cc:124
G4ErrorMatrix::determinant
G4double determinant() const
Definition: G4ErrorMatrix.cc:821
G4ErrorMatrix::invertHaywood5
virtual void invertHaywood5(G4int &ierr)
Definition: G4ErrorMatrix.cc:1056
G4ErrorMatrix::operator=
G4ErrorMatrix & operator=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:351
G4ErrorMatrix::num_col
virtual G4int num_col() const
G4ErrorMatrix::invertHaywood6
virtual void invertHaywood6(G4int &ierr)
Definition: G4ErrorMatrix.cc:1277
G4ErrorMatrix::operator*=
G4ErrorMatrix & operator*=(G4double t)
Definition: G4ErrorMatrix.cc:345
G4ErrorMatrix::operator-=
G4ErrorMatrix & operator-=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:332
G4ErrorMatrix::num_row
virtual G4int num_row() const
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:855
G4ErrorMatrix::operator+=
G4ErrorMatrix & operator+=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:325
G4ErrorMatrix::trace
G4double trace() const
Definition: G4ErrorMatrix.cc:845
G4ErrorMatrix::sub
G4ErrorMatrix sub(G4int min_row, G4int max_row, G4int min_col, G4int max_col) const
Definition: G4ErrorMatrix.cc:170
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:66
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77