Geant4 10.7.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//
27// ------------------------------------------------------------
28// GEANT 4 class implementation file
29// ------------------------------------------------------------
30
31#include "globals.hh"
32
33#include <cmath>
34#include <iostream>
35
36#include "G4ErrorMatrix.hh"
37#include "G4ErrorSymMatrix.hh"
38
39// Simple operation for all elements
40
41#define SIMPLE_UOP(OPER) \
42 G4ErrorMatrixIter a=m.begin(); \
43 G4ErrorMatrixIter e=m.end(); \
44 for(;a!=e; a++) (*a) OPER t;
45
46#define SIMPLE_BOP(OPER) \
47 G4ErrorMatrixIter a=m.begin(); \
48 G4ErrorMatrixConstIter b=mat2.m.begin(); \
49 G4ErrorMatrixIter e=m.end(); \
50 for(;a!=e; a++, b++) (*a) OPER (*b);
51
52#define SIMPLE_TOP(OPER) \
53 G4ErrorMatrixConstIter a=mat1.m.begin(); \
54 G4ErrorMatrixConstIter b=mat2.m.begin(); \
55 G4ErrorMatrixIter t=mret.m.begin(); \
56 G4ErrorMatrixConstIter e=mat1.m.end(); \
57 for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
58
59// Static functions.
60
61#define CHK_DIM_2(r1,r2,c1,c2,fun) \
62 if (r1!=r2 || c1!=c2) { \
63 G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
64 }
65
66#define CHK_DIM_1(c1,r2,fun) \
67 if (c1!=r2) { \
68 G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
69 }
70
71// Constructors. (Default constructors are inlined and in .icc file)
72
73G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q)
74 : m(p*q), nrow(p), ncol(q)
75{
76 size = nrow * ncol;
77}
78
79G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q,G4int init)
80 : m(p*q), nrow(p), ncol(q)
81{
82 size = nrow * ncol;
83
84 if (size > 0)
85 {
86 switch(init)
87 {
88 case 0:
89 break;
90
91 case 1:
92 {
93 if ( ncol == nrow )
94 {
95 G4ErrorMatrixIter a = m.begin();
96 G4ErrorMatrixIter b = m.end();
97 for( ; a<b; a+=(ncol+1)) *a = 1.0;
98 } else {
99 error("Invalid dimension in G4ErrorMatrix(G4int,G4int,1).");
100 }
101 break;
102 }
103 default:
104 error("G4ErrorMatrix: initialization must be either 0 or 1.");
105 }
106 }
107}
108
109//
110// Destructor
111//
112G4ErrorMatrix::~G4ErrorMatrix()
113{
114}
115
116G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatrix &mat1)
117 : m(mat1.size), nrow(mat1.nrow), ncol(mat1.ncol), size(mat1.size)
118{
119 m = mat1.m;
120}
121
122G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymMatrix &mat1)
123 : m(mat1.nrow*mat1.nrow), nrow(mat1.nrow), ncol(mat1.nrow)
124{
125 size = nrow * ncol;
126
127 G4int n = ncol;
128 G4ErrorMatrixConstIter sjk = mat1.m.begin();
129 G4ErrorMatrixIter m1j = m.begin();
130 G4ErrorMatrixIter mj = m.begin();
131 // j >= k
132 for(G4int j=1;j<=nrow;j++)
133 {
134 G4ErrorMatrixIter mjk = mj;
135 G4ErrorMatrixIter mkj = m1j;
136 for(G4int k=1;k<=j;k++)
137 {
138 *(mjk++) = *sjk;
139 if(j!=k) *mkj = *sjk;
140 sjk++;
141 mkj += n;
142 }
143 mj += n;
144 m1j++;
145 }
146}
147
148//
149//
150// Sub matrix
151//
152//
153
154G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row, G4int max_row,
155 G4int min_col,G4int max_col) const
156{
157 G4ErrorMatrix mret(max_row-min_row+1,max_col-min_col+1);
158 if(max_row > num_row() || max_col >num_col())
159 { error("G4ErrorMatrix::sub: Index out of range"); }
160 G4ErrorMatrixIter a = mret.m.begin();
161 G4int nc = num_col();
162 G4ErrorMatrixConstIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
163
164 for(G4int irow=1; irow<=mret.num_row(); irow++)
165 {
166 G4ErrorMatrixConstIter brc = b1;
167 for(G4int icol=1; icol<=mret.num_col(); icol++)
168 {
169 *(a++) = *(brc++);
170 }
171 b1 += nc;
172 }
173 return mret;
174}
175
176void G4ErrorMatrix::sub(G4int row,G4int col,const G4ErrorMatrix &mat1)
177{
178 if(row <1 || row+mat1.num_row()-1 > num_row() ||
179 col <1 || col+mat1.num_col()-1 > num_col() )
180 { error("G4ErrorMatrix::sub: Index out of range"); }
181 G4ErrorMatrixConstIter a = mat1.m.begin();
182 G4int nc = num_col();
183 G4ErrorMatrixIter b1 = m.begin() + (row - 1) * nc + col - 1;
184
185 for(G4int irow=1; irow<=mat1.num_row(); irow++)
186 {
187 G4ErrorMatrixIter brc = b1;
188 for(G4int icol=1; icol<=mat1.num_col(); icol++)
189 {
190 *(brc++) = *(a++);
191 }
192 b1 += nc;
193 }
194}
195
196//
197// Direct sum of two matricies
198//
199
200G4ErrorMatrix dsum(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
201{
202 G4ErrorMatrix mret(mat1.num_row() + mat2.num_row(),
203 mat1.num_col() + mat2.num_col(), 0);
204 mret.sub(1,1,mat1);
205 mret.sub(mat1.num_row()+1,mat1.num_col()+1,mat2);
206 return mret;
207}
208
209/* -----------------------------------------------------------------------
210 This section contains support routines for matrix.h. This section contains
211 The two argument functions +,-. They call the copy constructor and +=,-=.
212 ----------------------------------------------------------------------- */
213
214G4ErrorMatrix G4ErrorMatrix::operator- () const
215{
216 G4ErrorMatrix mat2(nrow, ncol);
217 G4ErrorMatrixConstIter a=m.begin();
218 G4ErrorMatrixIter b=mat2.m.begin();
219 G4ErrorMatrixConstIter e=m.end();
220 for(;a<e; a++, b++) (*b) = -(*a);
221 return mat2;
222}
223
224G4ErrorMatrix operator+(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
225{
226 G4ErrorMatrix mret(mat1.nrow, mat1.ncol);
227 CHK_DIM_2(mat1.num_row(),mat2.num_row(), mat1.num_col(),mat2.num_col(),+);
228 SIMPLE_TOP(+)
229 return mret;
230}
231
232//
233// operator -
234//
235
236G4ErrorMatrix operator-(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
237{
238 G4ErrorMatrix mret(mat1.num_row(), mat1.num_col());
239 CHK_DIM_2(mat1.num_row(),mat2.num_row(),
240 mat1.num_col(),mat2.num_col(),-);
241 SIMPLE_TOP(-)
242 return mret;
243}
244
245/* -----------------------------------------------------------------------
246 This section contains support routines for matrix.h. This file contains
247 The two argument functions *,/. They call copy constructor and then /=,*=.
248 ----------------------------------------------------------------------- */
249
250G4ErrorMatrix operator/(const G4ErrorMatrix &mat1,G4double t)
251{
252 G4ErrorMatrix mret(mat1);
253 mret /= t;
254 return mret;
255}
256
257G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,G4double t)
258{
259 G4ErrorMatrix mret(mat1);
260 mret *= t;
261 return mret;
262}
263
264G4ErrorMatrix operator*(G4double t,const G4ErrorMatrix &mat1)
265{
266 G4ErrorMatrix mret(mat1);
267 mret *= t;
268 return mret;
269}
270
271G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
272{
273 // initialize matrix to 0.0
274 G4ErrorMatrix mret(mat1.nrow,mat2.ncol,0);
275 CHK_DIM_1(mat1.ncol,mat2.nrow,*);
276
277 G4int m1cols = mat1.ncol;
278 G4int m2cols = mat2.ncol;
279
280 for (G4int i=0; i<mat1.nrow; i++)
281 {
282 for (G4int j=0; j<m1cols; j++)
283 {
284 G4double temp = mat1.m[i*m1cols+j];
285 G4ErrorMatrixIter pt = mret.m.begin() + i*m2cols;
286
287 // Loop over k (the column index in matrix mat2)
288 G4ErrorMatrixConstIter pb = mat2.m.begin() + m2cols*j;
289 const G4ErrorMatrixConstIter pblast = pb + m2cols;
290 while (pb < pblast) // Loop checking, 06.08.2015, G.Cosmo
291 {
292 (*pt) += temp * (*pb);
293 pb++;
294 pt++;
295 }
296 }
297 }
298 return mret;
299}
300
301/* -----------------------------------------------------------------------
302 This section contains the assignment and inplace operators =,+=,-=,*=,/=.
303 ----------------------------------------------------------------------- */
304
305G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorMatrix &mat2)
306{
307 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
308 SIMPLE_BOP(+=)
309 return (*this);
310}
311
312G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorMatrix &mat2)
313{
314 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
315 SIMPLE_BOP(-=)
316 return (*this);
317}
318
319G4ErrorMatrix & G4ErrorMatrix::operator/=(G4double t)
320{
321 SIMPLE_UOP(/=)
322 return (*this);
323}
324
325G4ErrorMatrix & G4ErrorMatrix::operator*=(G4double t)
326{
327 SIMPLE_UOP(*=)
328 return (*this);
329}
330
331G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorMatrix &mat1)
332{
333 if (&mat1 == this) { return *this; }
334
335 if(mat1.nrow*mat1.ncol != size) //??fixme?? mat1.size != size
336 {
337 size = mat1.nrow * mat1.ncol;
338 m.resize(size); //??fixme?? if (size < mat1.size) m.resize(mat1.size);
339 }
340 nrow = mat1.nrow;
341 ncol = mat1.ncol;
342 m = mat1.m;
343 return (*this);
344}
345
346
347// Print the G4ErrorMatrix.
348
349std::ostream& operator<<(std::ostream &os, const G4ErrorMatrix &q)
350{
351 os << "\n";
352
353 // Fixed format needs 3 extra characters for field,
354 // while scientific needs 7
355
356 G4int width;
357 if(os.flags() & std::ios::fixed)
358 { width = os.precision()+3; }
359 else
360 { width = os.precision()+7; }
361 for(G4int irow = 1; irow<= q.num_row(); irow++)
362 {
363 for(G4int icol = 1; icol <= q.num_col(); icol++)
364 {
365 os.width(width);
366 os << q(irow,icol) << " ";
367 }
368 os << G4endl;
369 }
370 return os;
371}
372
373G4ErrorMatrix G4ErrorMatrix::T() const
374{
375 G4ErrorMatrix mret(ncol,nrow);
376 G4ErrorMatrixConstIter pl = m.end();
377 G4ErrorMatrixConstIter pme = m.begin();
378 G4ErrorMatrixIter pt = mret.m.begin();
379 G4ErrorMatrixIter ptl = mret.m.end();
380 for (; pme < pl; pme++, pt+=nrow)
381 {
382 if (pt >= ptl) { pt -= (size-1); }
383 (*pt) = (*pme);
384 }
385 return mret;
386}
387
388G4ErrorMatrix
389G4ErrorMatrix::apply(G4double (*f)(G4double, G4int, G4int)) const
390{
391 G4ErrorMatrix mret(num_row(),num_col());
392 G4ErrorMatrixConstIter a = m.begin();
393 G4ErrorMatrixIter b = mret.m.begin();
394 for(G4int ir=1;ir<=num_row();ir++)
395 {
396 for(G4int ic=1;ic<=num_col();ic++)
397 {
398 *(b++) = (*f)(*(a++), ir, ic);
399 }
400 }
401 return mret;
402}
403
404G4int G4ErrorMatrix::dfinv_matrix(G4int *ir) {
405 if (num_col()!=num_row())
406 { error("dfinv_matrix: G4ErrorMatrix is not NxN"); }
407 G4int n = num_col();
408 if (n==1) { return 0; }
409
410 G4double s31, s32;
411 G4double s33, s34;
412
413 G4ErrorMatrixIter m11 = m.begin();
414 G4ErrorMatrixIter m12 = m11 + 1;
415 G4ErrorMatrixIter m21 = m11 + n;
416 G4ErrorMatrixIter m22 = m12 + n;
417 *m21 = -(*m22) * (*m11) * (*m21);
418 *m12 = -(*m12);
419 if (n>2)
420 {
421 G4ErrorMatrixIter mi = m.begin() + 2 * n;
422 G4ErrorMatrixIter mii= m.begin() + 2 * n + 2;
423 G4ErrorMatrixIter mimim = m.begin() + n + 1;
424 for (G4int i=3;i<=n;i++)
425 {
426 G4int im2 = i - 2;
427 G4ErrorMatrixIter mj = m.begin();
428 G4ErrorMatrixIter mji = mj + i - 1;
429 G4ErrorMatrixIter mij = mi;
430 for (G4int j=1;j<=im2;j++)
431 {
432 s31 = 0.0;
433 s32 = *mji;
434 G4ErrorMatrixIter mkj = mj + j - 1;
435 G4ErrorMatrixIter mik = mi + j - 1;
436 G4ErrorMatrixIter mjkp = mj + j;
437 G4ErrorMatrixIter mkpi = mj + n + i - 1;
438 for (G4int k=j;k<=im2;k++)
439 {
440 s31 += (*mkj) * (*(mik++));
441 s32 += (*(mjkp++)) * (*mkpi);
442 mkj += n;
443 mkpi += n;
444 }
445 *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
446 *mji = -s32;
447 mj += n;
448 mji += n;
449 mij++;
450 }
451 *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
452 *(mimim+1) = -(*(mimim+1));
453 mi += n;
454 mimim += (n+1);
455 mii += (n+1);
456 }
457 }
458 G4ErrorMatrixIter mi = m.begin();
459 G4ErrorMatrixIter mii = m.begin();
460 for (G4int i=1;i<n;i++)
461 {
462 G4int ni = n - i;
463 G4ErrorMatrixIter mij = mi;
464 G4int j;
465 for (j=1; j<=i;j++)
466 {
467 s33 = *mij;
468 G4ErrorMatrixIter mikj = mi + n + j - 1;
469 G4ErrorMatrixIter miik = mii + 1;
470 G4ErrorMatrixIter min_end = mi + n;
471 for (;miik<min_end;)
472 {
473 s33 += (*mikj) * (*(miik++));
474 mikj += n;
475 }
476 *(mij++) = s33;
477 }
478 for (j=1;j<=ni;j++)
479 {
480 s34 = 0.0;
481 G4ErrorMatrixIter miik = mii + j;
482 G4ErrorMatrixIter mikij = mii + j * n + j;
483 for (G4int k=j;k<=ni;k++)
484 {
485 s34 += *mikij * (*(miik++));
486 mikij += n;
487 }
488 *(mii+j) = s34;
489 }
490 mi += n;
491 mii += (n+1);
492 }
493 G4int nxch = ir[n];
494 if (nxch==0) return 0;
495 for (G4int mq=1;mq<=nxch;mq++)
496 {
497 G4int k = nxch - mq + 1;
498 G4int ij = ir[k];
499 G4int i = ij >> 12;
500 G4int j = ij%4096;
501 G4ErrorMatrixIter mki = m.begin() + i - 1;
502 G4ErrorMatrixIter mkj = m.begin() + j - 1;
503 for (k=1; k<=n;k++)
504 {
505 // 2/24/05 David Sachs fix of improper swap bug that was present
506 // for many years:
507 G4double ti = *mki; // 2/24/05
508 *mki = *mkj;
509 *mkj = ti; // 2/24/05
510 mki += n;
511 mkj += n;
512 }
513 }
514 return 0;
515}
516
517G4int G4ErrorMatrix::dfact_matrix(G4double &det, G4int *ir)
518{
519 if (ncol!=nrow)
520 error("dfact_matrix: G4ErrorMatrix is not NxN");
521
522 G4int ifail, jfail;
523 G4int n = ncol;
524
525 G4double tf;
526 G4double g1 = 1.0e-19, g2 = 1.0e19;
527
528 G4double p, q, t;
529 G4double s11, s12;
530
531 G4double epsilon = 8*DBL_EPSILON;
532 // could be set to zero (like it was before)
533 // but then the algorithm often doesn't detect
534 // that a matrix is singular
535
536 G4int normal = 0, imposs = -1;
537 G4int jrange = 0, jover = 1, junder = -1;
538 ifail = normal;
539 jfail = jrange;
540 G4int nxch = 0;
541 det = 1.0;
542 G4ErrorMatrixIter mj = m.begin();
543 G4ErrorMatrixIter mjj = mj;
544 for (G4int j=1;j<=n;j++)
545 {
546 G4int k = j;
547 p = (std::fabs(*mjj));
548 if (j!=n) {
549 G4ErrorMatrixIter mij = mj + n + j - 1;
550 for (G4int i=j+1;i<=n;i++)
551 {
552 q = (std::fabs(*(mij)));
553 if (q > p)
554 {
555 k = i;
556 p = q;
557 }
558 mij += n;
559 }
560 if (k==j)
561 {
562 if (p <= epsilon)
563 {
564 det = 0;
565 ifail = imposs;
566 jfail = jrange;
567 return ifail;
568 }
569 det = -det; // in this case the sign of the determinant
570 // must not change. So I change it twice.
571 }
572 G4ErrorMatrixIter mjl = mj;
573 G4ErrorMatrixIter mkl = m.begin() + (k-1)*n;
574 for (G4int l=1;l<=n;l++)
575 {
576 tf = *mjl;
577 *(mjl++) = *mkl;
578 *(mkl++) = tf;
579 }
580 nxch = nxch + 1; // this makes the determinant change its sign
581 ir[nxch] = (((j)<<12)+(k));
582 }
583 else
584 {
585 if (p <= epsilon)
586 {
587 det = 0.0;
588 ifail = imposs;
589 jfail = jrange;
590 return ifail;
591 }
592 }
593 det *= *mjj;
594 *mjj = 1.0 / *mjj;
595 t = (std::fabs(det));
596 if (t < g1)
597 {
598 det = 0.0;
599 if (jfail == jrange) jfail = junder;
600 }
601 else if (t > g2)
602 {
603 det = 1.0;
604 if (jfail==jrange) jfail = jover;
605 }
606 if (j!=n)
607 {
608 G4ErrorMatrixIter mk = mj + n;
609 G4ErrorMatrixIter mkjp = mk + j;
610 G4ErrorMatrixIter mjk = mj + j;
611 for (k=j+1;k<=n;k++)
612 {
613 s11 = - (*mjk);
614 s12 = - (*mkjp);
615 if (j!=1)
616 {
617 G4ErrorMatrixIter mik = m.begin() + k - 1;
618 G4ErrorMatrixIter mijp = m.begin() + j;
619 G4ErrorMatrixIter mki = mk;
620 G4ErrorMatrixIter mji = mj;
621 for (G4int i=1;i<j;i++)
622 {
623 s11 += (*mik) * (*(mji++));
624 s12 += (*mijp) * (*(mki++));
625 mik += n;
626 mijp += n;
627 }
628 }
629 *(mjk++) = -s11 * (*mjj);
630 *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
631 mk += n;
632 mkjp += n;
633 }
634 }
635 mj += n;
636 mjj += (n+1);
637 }
638 if (nxch%2==1) det = -det;
639 if (jfail !=jrange) det = 0.0;
640 ir[n] = nxch;
641 return 0;
642}
643
644void G4ErrorMatrix::invert(G4int &ierr)
645{
646 if(ncol != nrow)
647 { error("G4ErrorMatrix::invert: G4ErrorMatrix is not NxN"); }
648
649 static G4ThreadLocal G4int max_array = 20;
650 static G4ThreadLocal G4int *ir = 0 ; if (!ir) ir= new G4int [max_array+1];
651
652 if (ncol > max_array)
653 {
654 delete [] ir;
655 max_array = nrow;
656 ir = new G4int [max_array+1];
657 }
658 G4double t1, t2, t3;
659 G4double det, temp, ss;
660 G4int ifail;
661 switch(nrow)
662 {
663 case 3:
664 G4double c11,c12,c13,c21,c22,c23,c31,c32,c33;
665 ifail = 0;
666 c11 = (*(m.begin()+4)) * (*(m.begin()+8))
667 - (*(m.begin()+5)) * (*(m.begin()+7));
668 c12 = (*(m.begin()+5)) * (*(m.begin()+6))
669 - (*(m.begin()+3)) * (*(m.begin()+8));
670 c13 = (*(m.begin()+3)) * (*(m.begin()+7))
671 - (*(m.begin()+4)) * (*(m.begin()+6));
672 c21 = (*(m.begin()+7)) * (*(m.begin()+2))
673 - (*(m.begin()+8)) * (*(m.begin()+1));
674 c22 = (*(m.begin()+8)) * (*m.begin())
675 - (*(m.begin()+6)) * (*(m.begin()+2));
676 c23 = (*(m.begin()+6)) * (*(m.begin()+1))
677 - (*(m.begin()+7)) * (*m.begin());
678 c31 = (*(m.begin()+1)) * (*(m.begin()+5))
679 - (*(m.begin()+2)) * (*(m.begin()+4));
680 c32 = (*(m.begin()+2)) * (*(m.begin()+3))
681 - (*m.begin()) * (*(m.begin()+5));
682 c33 = (*m.begin()) * (*(m.begin()+4))
683 - (*(m.begin()+1)) * (*(m.begin()+3));
684 t1 = std::fabs(*m.begin());
685 t2 = std::fabs(*(m.begin()+3));
686 t3 = std::fabs(*(m.begin()+6));
687 if (t1 >= t2)
688 {
689 if (t3 >= t1)
690 {
691 temp = *(m.begin()+6);
692 det = c23*c12-c22*c13;
693 }
694 else
695 {
696 temp = *(m.begin());
697 det = c22*c33-c23*c32;
698 }
699 }
700 else if (t3 >= t2)
701 {
702 temp = *(m.begin()+6);
703 det = c23*c12-c22*c13;
704 }
705 else
706 {
707 temp = *(m.begin()+3);
708 det = c13*c32-c12*c33;
709 }
710 if (det==0)
711 {
712 ierr = 1;
713 return;
714 }
715 {
716 ss = temp/det;
717 G4ErrorMatrixIter mq = m.begin();
718 *(mq++) = ss*c11;
719 *(mq++) = ss*c21;
720 *(mq++) = ss*c31;
721 *(mq++) = ss*c12;
722 *(mq++) = ss*c22;
723 *(mq++) = ss*c32;
724 *(mq++) = ss*c13;
725 *(mq++) = ss*c23;
726 *(mq) = ss*c33;
727 }
728 break;
729 case 2:
730 ifail = 0;
731 det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
732 if (det==0)
733 {
734 ierr = 1;
735 return;
736 }
737 ss = 1.0/det;
738 temp = ss*(*(m.begin()+3));
739 *(m.begin()+1) *= -ss;
740 *(m.begin()+2) *= -ss;
741 *(m.begin()+3) = ss*(*m.begin());
742 *(m.begin()) = temp;
743 break;
744 case 1:
745 ifail = 0;
746 if ((*(m.begin()))==0)
747 {
748 ierr = 1;
749 return;
750 }
751 *(m.begin()) = 1.0/(*(m.begin()));
752 break;
753 case 4:
754 invertHaywood4(ierr);
755 return;
756 case 5:
757 invertHaywood5(ierr);
758 return;
759 case 6:
760 invertHaywood6(ierr);
761 return;
762 default:
763 ifail = dfact_matrix(det, ir);
764 if(ifail) {
765 ierr = 1;
766 return;
767 }
768 dfinv_matrix(ir);
769 break;
770 }
771 ierr = 0;
772 return;
773}
774
775G4double G4ErrorMatrix::determinant() const
776{
777 static G4ThreadLocal G4int max_array = 20;
778 static G4ThreadLocal G4int *ir = 0 ; if (!ir) ir= new G4int [max_array+1];
779 if(ncol != nrow)
780 { error("G4ErrorMatrix::determinant: G4ErrorMatrix is not NxN"); }
781 if (ncol > max_array)
782 {
783 delete [] ir;
784 max_array = nrow;
785 ir = new G4int [max_array+1];
786 }
787 G4double det;
788 G4ErrorMatrix mt(*this);
789 G4int i = mt.dfact_matrix(det, ir);
790 if(i==0) return det;
791 return 0;
792}
793
794G4double G4ErrorMatrix::trace() const
795{
796 G4double t = 0.0;
797 for (G4ErrorMatrixConstIter d = m.begin(); d < m.end(); d += (ncol+1) )
798 {
799 t += *d;
800 }
801 return t;
802}
803
804void G4ErrorMatrix::error(const char *msg)
805{
806 std::ostringstream message;
807 message << msg;
808 G4Exception("G4ErrorMatrix::error()", "GEANT4e-Error",
809 FatalException, message, "Exiting to System.");
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}
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:96
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:61
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:200
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:257
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:349
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:46
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:236
A30
#define A30
Definition: G4ErrorMatrix.cc:838
M23
#define M23
Definition: G4ErrorMatrix.cc:874
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorMatrix.cc:41
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:250
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:224
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:52
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:66
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
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4ErrorMatrix::apply
G4ErrorMatrix apply(G4double(*f)(G4double, G4int, G4int)) const
Definition: G4ErrorMatrix.cc:389
G4ErrorMatrix::invertHaywood4
virtual void invertHaywood4(G4int &ierr)
Definition: G4ErrorMatrix.cc:910
G4ErrorMatrix::operator-
G4ErrorMatrix operator-() const
Definition: G4ErrorMatrix.cc:214
G4ErrorMatrix::invert
virtual void invert(G4int &ierr)
Definition: G4ErrorMatrix.cc:644
G4ErrorMatrix::operator/=
G4ErrorMatrix & operator/=(G4double t)
Definition: G4ErrorMatrix.cc:319
G4ErrorMatrix::T
G4ErrorMatrix T() const
Definition: G4ErrorMatrix.cc:373
G4ErrorMatrix::~G4ErrorMatrix
virtual ~G4ErrorMatrix()
Definition: G4ErrorMatrix.cc:112
G4ErrorMatrix::determinant
G4double determinant() const
Definition: G4ErrorMatrix.cc:775
G4ErrorMatrix::invertHaywood5
virtual void invertHaywood5(G4int &ierr)
Definition: G4ErrorMatrix.cc:1009
G4ErrorMatrix::operator=
G4ErrorMatrix & operator=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:331
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:325
G4ErrorMatrix::operator-=
G4ErrorMatrix & operator-=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:312
G4ErrorMatrix::num_row
virtual G4int num_row() const
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:804
G4ErrorMatrix::operator+=
G4ErrorMatrix & operator+=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:305
G4ErrorMatrix::trace
G4double trace() const
Definition: G4ErrorMatrix.cc:794
G4ErrorMatrix::sub
G4ErrorMatrix sub(G4int min_row, G4int max_row, G4int min_col, G4int max_col) const
Definition: G4ErrorMatrix.cc:154
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:66
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77