Geant4-11
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 G4int 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.begin();
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
A02
#define A02
Definition: G4ErrorMatrix.cc:869
F30
#define F30
Definition: G4ErrorMatrix.cc:954
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
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
M22
#define M22
Definition: G4ErrorMatrix.cc:923
A42
#define A42
Definition: G4ErrorMatrix.cc:897
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:41
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:35
m
static constexpr double m
Definition: G4SIunits.hh:109
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::m
std::vector< G4double > m
Definition: G4ErrorMatrix.hh:222
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::dfinv_matrix
G4int dfinv_matrix(G4int *ir)
Definition: G4ErrorMatrix.cc:432
G4ErrorMatrix::~G4ErrorMatrix
virtual ~G4ErrorMatrix()
Definition: G4ErrorMatrix.cc:124
G4ErrorMatrix::dfact_matrix
G4int dfact_matrix(G4double &det, G4int *ir)
Definition: G4ErrorMatrix.cc:552
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
G4ErrorSymMatrix::m
std::vector< G4double > m
Definition: G4ErrorSymMatrix.hh:219
HepGeom::BasicVector3D::operator*
BasicVector3D< float > operator*(const BasicVector3D< float > &v, double a)
Definition: BasicVector3D.h:408
HepGeom::BasicVector3D::operator/
BasicVector3D< float > operator/(const BasicVector3D< float > &v, double a)
Definition: BasicVector3D.h:435
HepGeom::BasicVector3D::operator-
BasicVector3D< float > operator-(const BasicVector3D< float > &v)
Definition: BasicVector3D.h:390
HepGeom::BasicVector3D::operator<<
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
Definition: BasicVector3D.cc:106
HepGeom::BasicVector3D::operator+
BasicVector3D< float > operator+(const BasicVector3D< float > &v)
Definition: BasicVector3D.h:374
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:79
anonymous_namespace{G4QuasiElRatios.cc}::mi
const G4double mi
Definition: G4QuasiElRatios.cc:65
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:66
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77
Generated by doxygen 1.9.3