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