Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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