source: trunk/source/error_propagation/src/G4ErrorMatrix.cc@ 971

Last change on this file since 971 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 67.1 KB
Line 
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: G4ErrorMatrix.cc,v 1.3 2007/06/21 15:04:06 gunter Exp $
27// GEANT4 tag $Name: HEAD $
28//
29// ------------------------------------------------------------
30// GEANT 4 class implementation file
31// ------------------------------------------------------------
32
33#include "globals.hh"
34
35#include <cmath>
36#include <iostream>
37
38#include "G4ErrorMatrix.hh"
39#include "G4ErrorSymMatrix.hh"
40
41// Simple operation for all elements
42
43#define SIMPLE_UOP(OPER) \
44 G4ErrorMatrixIter a=m.begin(); \
45 G4ErrorMatrixIter e=m.end(); \
46 for(;a!=e; a++) (*a) OPER t;
47
48#define SIMPLE_BOP(OPER) \
49 G4ErrorMatrixIter a=m.begin(); \
50 G4ErrorMatrixConstIter b=m2.m.begin(); \
51 G4ErrorMatrixIter e=m.end(); \
52 for(;a!=e; a++, b++) (*a) OPER (*b);
53
54#define SIMPLE_TOP(OPER) \
55 G4ErrorMatrixConstIter a=m1.m.begin(); \
56 G4ErrorMatrixConstIter b=m2.m.begin(); \
57 G4ErrorMatrixIter t=mret.m.begin(); \
58 G4ErrorMatrixConstIter e=m1.m.end(); \
59 for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
60
61// Static functions.
62
63#define CHK_DIM_2(r1,r2,c1,c2,fun) \
64 if (r1!=r2 || c1!=c2) { \
65 G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
66 }
67
68#define CHK_DIM_1(c1,r2,fun) \
69 if (c1!=r2) { \
70 G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
71 }
72
73// Constructors. (Default constructors are inlined and in .icc file)
74
75G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q)
76 : m(p*q), nrow(p), ncol(q)
77{
78 size = nrow * ncol;
79}
80
81G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q,G4int init)
82 : m(p*q), nrow(p), ncol(q)
83{
84 size = nrow * ncol;
85
86 if (size > 0)
87 {
88 switch(init)
89 {
90 case 0:
91 break;
92
93 case 1:
94 {
95 if ( ncol == nrow )
96 {
97 G4ErrorMatrixIter a = m.begin();
98 G4ErrorMatrixIter b = m.end();
99 for( ; a<b; a+=(ncol+1)) *a = 1.0;
100 } else {
101 error("Invalid dimension in G4ErrorMatrix(G4int,G4int,1).");
102 }
103 break;
104 }
105 default:
106 error("G4ErrorMatrix: initialization must be either 0 or 1.");
107 }
108 }
109}
110
111//
112// Destructor
113//
114G4ErrorMatrix::~G4ErrorMatrix()
115{
116}
117
118G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatrix &m1)
119 : m(m1.size), nrow(m1.nrow), ncol(m1.ncol), size(m1.size)
120{
121 m = m1.m;
122}
123
124G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymMatrix &m1)
125 : m(m1.nrow*m1.nrow), nrow(m1.nrow), ncol(m1.nrow)
126{
127 size = nrow * ncol;
128
129 G4int n = ncol;
130 G4ErrorMatrixConstIter sjk = m1.m.begin();
131 G4ErrorMatrixIter m1j = m.begin();
132 G4ErrorMatrixIter mj = m.begin();
133 // j >= k
134 for(G4int j=1;j<=nrow;j++)
135 {
136 G4ErrorMatrixIter mjk = mj;
137 G4ErrorMatrixIter mkj = m1j;
138 for(G4int k=1;k<=j;k++)
139 {
140 *(mjk++) = *sjk;
141 if(j!=k) *mkj = *sjk;
142 sjk++;
143 mkj += n;
144 }
145 mj += n;
146 m1j++;
147 }
148}
149
150//
151//
152// Sub matrix
153//
154//
155
156G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row, G4int max_row,
157 G4int min_col,G4int max_col) const
158{
159 G4ErrorMatrix mret(max_row-min_row+1,max_col-min_col+1);
160 if(max_row > num_row() || max_col >num_col())
161 { error("G4ErrorMatrix::sub: Index out of range"); }
162 G4ErrorMatrixIter a = mret.m.begin();
163 G4int nc = num_col();
164 G4ErrorMatrixConstIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
165
166 for(G4int irow=1; irow<=mret.num_row(); irow++)
167 {
168 G4ErrorMatrixConstIter brc = b1;
169 for(G4int icol=1; icol<=mret.num_col(); icol++)
170 {
171 *(a++) = *(brc++);
172 }
173 b1 += nc;
174 }
175 return mret;
176}
177
178void G4ErrorMatrix::sub(G4int row,G4int col,const G4ErrorMatrix &m1)
179{
180 if(row <1 || row+m1.num_row()-1 > num_row() ||
181 col <1 || col+m1.num_col()-1 > num_col() )
182 { error("G4ErrorMatrix::sub: Index out of range"); }
183 G4ErrorMatrixConstIter a = m1.m.begin();
184 G4int nc = num_col();
185 G4ErrorMatrixIter b1 = m.begin() + (row - 1) * nc + col - 1;
186
187 for(G4int irow=1; irow<=m1.num_row(); irow++)
188 {
189 G4ErrorMatrixIter brc = b1;
190 for(G4int icol=1; icol<=m1.num_col(); icol++)
191 {
192 *(brc++) = *(a++);
193 }
194 b1 += nc;
195 }
196}
197
198//
199// Direct sum of two matricies
200//
201
202G4ErrorMatrix dsum(const G4ErrorMatrix &m1, const G4ErrorMatrix &m2)
203{
204 G4ErrorMatrix mret(m1.num_row() + m2.num_row(), m1.num_col() + m2.num_col(),
205 0);
206 mret.sub(1,1,m1);
207 mret.sub(m1.num_row()+1,m1.num_col()+1,m2);
208 return mret;
209}
210
211/* -----------------------------------------------------------------------
212 This section contains support routines for matrix.h. This section contains
213 The two argument functions +,-. They call the copy constructor and +=,-=.
214 ----------------------------------------------------------------------- */
215
216G4ErrorMatrix G4ErrorMatrix::operator- () const
217{
218 G4ErrorMatrix m2(nrow, ncol);
219 G4ErrorMatrixConstIter a=m.begin();
220 G4ErrorMatrixIter b=m2.m.begin();
221 G4ErrorMatrixConstIter e=m.end();
222 for(;a<e; a++, b++) (*b) = -(*a);
223 return m2;
224}
225
226G4ErrorMatrix operator+(const G4ErrorMatrix &m1,const G4ErrorMatrix &m2)
227{
228 G4ErrorMatrix mret(m1.nrow, m1.ncol);
229 CHK_DIM_2(m1.num_row(),m2.num_row(), m1.num_col(),m2.num_col(),+);
230 SIMPLE_TOP(+)
231 return mret;
232}
233
234//
235// operator -
236//
237
238G4ErrorMatrix operator-(const G4ErrorMatrix &m1,const G4ErrorMatrix &m2)
239{
240 G4ErrorMatrix mret(m1.num_row(), m1.num_col());
241 CHK_DIM_2(m1.num_row(),m2.num_row(),
242 m1.num_col(),m2.num_col(),-);
243 SIMPLE_TOP(-)
244 return mret;
245}
246
247/* -----------------------------------------------------------------------
248 This section contains support routines for matrix.h. This file contains
249 The two argument functions *,/. They call copy constructor and then /=,*=.
250 ----------------------------------------------------------------------- */
251
252G4ErrorMatrix operator/(const G4ErrorMatrix &m1,G4double t)
253{
254 G4ErrorMatrix mret(m1);
255 mret /= t;
256 return mret;
257}
258
259G4ErrorMatrix operator*(const G4ErrorMatrix &m1,G4double t)
260{
261 G4ErrorMatrix mret(m1);
262 mret *= t;
263 return mret;
264}
265
266G4ErrorMatrix operator*(G4double t,const G4ErrorMatrix &m1)
267{
268 G4ErrorMatrix mret(m1);
269 mret *= t;
270 return mret;
271}
272
273G4ErrorMatrix operator*(const G4ErrorMatrix &m1,const G4ErrorMatrix &m2)
274{
275 // initialize matrix to 0.0
276 G4ErrorMatrix mret(m1.nrow,m2.ncol,0);
277 CHK_DIM_1(m1.ncol,m2.nrow,*);
278
279 G4int m1cols = m1.ncol;
280 G4int m2cols = m2.ncol;
281
282 for (G4int i=0; i<m1.nrow; i++)
283 {
284 for (G4int j=0; j<m1cols; j++)
285 {
286 G4double temp = m1.m[i*m1cols+j];
287 G4ErrorMatrixIter pt = mret.m.begin() + i*m2cols;
288
289 // Loop over k (the column index in matrix m2)
290 G4ErrorMatrixConstIter pb = m2.m.begin() + m2cols*j;
291 const G4ErrorMatrixConstIter pblast = pb + m2cols;
292 while (pb < pblast)
293 {
294 (*pt) += temp * (*pb);
295 pb++;
296 pt++;
297 }
298 }
299 }
300 return mret;
301}
302
303/* -----------------------------------------------------------------------
304 This section contains the assignment and inplace operators =,+=,-=,*=,/=.
305 ----------------------------------------------------------------------- */
306
307G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorMatrix &m2)
308{
309 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
310 SIMPLE_BOP(+=)
311 return (*this);
312}
313
314G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorMatrix &m2)
315{
316 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
317 SIMPLE_BOP(-=)
318 return (*this);
319}
320
321G4ErrorMatrix & G4ErrorMatrix::operator/=(G4double t)
322{
323 SIMPLE_UOP(/=)
324 return (*this);
325}
326
327G4ErrorMatrix & G4ErrorMatrix::operator*=(G4double t)
328{
329 SIMPLE_UOP(*=)
330 return (*this);
331}
332
333G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorMatrix &m1)
334{
335 if(m1.nrow*m1.ncol != size) //??fixme?? m1.size != size
336 {
337 size = m1.nrow * m1.ncol;
338 m.resize(size); //??fixme?? if (size < m1.size) m.resize(m1.size);
339 }
340 nrow = m1.nrow;
341 ncol = m1.ncol;
342 m = m1.m;
343 return (*this);
344}
345
346
347// Print the G4ErrorMatrix.
348
349std::ostream& operator<<(std::ostream &s, const G4ErrorMatrix &q)
350{
351 s << "\n";
352
353 // Fixed format needs 3 extra characters for field,
354 // while scientific needs 7
355
356 G4int width;
357 if(s.flags() & std::ios::fixed)
358 { width = s.precision()+3; }
359 else
360 { width = s.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 s.width(width);
366 s << q(irow,icol) << " ";
367 }
368 s << G4endl;
369 }
370 return s;
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 mm=1;mm<=nxch;mm++)
496 {
497 G4int k = nxch - mm + 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 G4int max_array = 20;
650 static G4int *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, s;
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 G4double s = temp/det;
717 G4ErrorMatrixIter mm = m.begin();
718 *(mm++) = s*c11;
719 *(mm++) = s*c21;
720 *(mm++) = s*c31;
721 *(mm++) = s*c12;
722 *(mm++) = s*c22;
723 *(mm++) = s*c32;
724 *(mm++) = s*c13;
725 *(mm++) = s*c23;
726 *(mm) = s*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 s = 1.0/det;
738 temp = s*(*(m.begin()+3));
739 *(m.begin()+1) *= -s;
740 *(m.begin()+2) *= -s;
741 *(m.begin()+3) = s*(*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 G4int max_array = 20;
778 static G4int *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 *s)
805{
806 G4cerr << s << G4endl;
807 G4cerr << "---Exiting to System." << G4endl;
808 abort();
809}
810
811
812// Aij are indices for a 6x6 matrix.
813// Mij are indices for a 5x5 matrix.
814// Fij are indices for a 4x4 matrix.
815
816#define A00 0
817#define A01 1
818#define A02 2
819#define A03 3
820#define A04 4
821#define A05 5
822
823#define A10 6
824#define A11 7
825#define A12 8
826#define A13 9
827#define A14 10
828#define A15 11
829
830#define A20 12
831#define A21 13
832#define A22 14
833#define A23 15
834#define A24 16
835#define A25 17
836
837#define A30 18
838#define A31 19
839#define A32 20
840#define A33 21
841#define A34 22
842#define A35 23
843
844#define A40 24
845#define A41 25
846#define A42 26
847#define A43 27
848#define A44 28
849#define A45 29
850
851#define A50 30
852#define A51 31
853#define A52 32
854#define A53 33
855#define A54 34
856#define A55 35
857
858#define M00 0
859#define M01 1
860#define M02 2
861#define M03 3
862#define M04 4
863
864#define M10 5
865#define M11 6
866#define M12 7
867#define M13 8
868#define M14 9
869
870#define M20 10
871#define M21 11
872#define M22 12
873#define M23 13
874#define M24 14
875
876#define M30 15
877#define M31 16
878#define M32 17
879#define M33 18
880#define M34 19
881
882#define M40 20
883#define M41 21
884#define M42 22
885#define M43 23
886#define M44 24
887
888#define F00 0
889#define F01 1
890#define F02 2
891#define F03 3
892
893#define F10 4
894#define F11 5
895#define F12 6
896#define F13 7
897
898#define F20 8
899#define F21 9
900#define F22 10
901#define F23 11
902
903#define F30 12
904#define F31 13
905#define F32 14
906#define F33 15
907
908
909void G4ErrorMatrix::invertHaywood4 (G4int & ifail)
910{
911 ifail = 0;
912
913 // Find all NECESSARY 2x2 dets: (18 of them)
914
915 G4double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
916 G4double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
917 G4double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20];
918 G4double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21];
919 G4double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22];
920 G4double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
921 G4double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
922 G4double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
923 G4double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
924 G4double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
925 G4double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
926 G4double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32];
927 G4double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
928 G4double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
929 G4double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
930 G4double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
931 G4double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
932 G4double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
933
934 // Find all NECESSARY 3x3 dets: (16 of them)
935
936 G4double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02
937 + m[F02]*Det2_12_01;
938 G4double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03
939 + m[F03]*Det2_12_01;
940 G4double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03
941 + m[F03]*Det2_12_02;
942 G4double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13
943 + m[F03]*Det2_12_12;
944 G4double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02
945 + m[F02]*Det2_13_01;
946 G4double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
947 + m[F03]*Det2_13_01;
948 G4double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
949 + m[F03]*Det2_13_02;
950 G4double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
951 + m[F03]*Det2_13_12;
952 G4double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02
953 + m[F02]*Det2_23_01;
954 G4double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
955 + m[F03]*Det2_23_01;
956 G4double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
957 + m[F03]*Det2_23_02;
958 G4double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
959 + m[F03]*Det2_23_12;
960 G4double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02
961 + m[F12]*Det2_23_01;
962 G4double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03
963 + m[F13]*Det2_23_01;
964 G4double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03
965 + m[F13]*Det2_23_02;
966 G4double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13
967 + m[F13]*Det2_23_12;
968
969 // Find the 4x4 det:
970
971 G4double det = m[F00]*Det3_123_123
972 - m[F01]*Det3_123_023
973 + m[F02]*Det3_123_013
974 - m[F03]*Det3_123_012;
975
976 if ( det == 0 )
977 {
978 ifail = 1;
979 return;
980 }
981
982 G4double oneOverDet = 1.0/det;
983 G4double mn1OverDet = - oneOverDet;
984
985 m[F00] = Det3_123_123 * oneOverDet;
986 m[F01] = Det3_023_123 * mn1OverDet;
987 m[F02] = Det3_013_123 * oneOverDet;
988 m[F03] = Det3_012_123 * mn1OverDet;
989
990 m[F10] = Det3_123_023 * mn1OverDet;
991 m[F11] = Det3_023_023 * oneOverDet;
992 m[F12] = Det3_013_023 * mn1OverDet;
993 m[F13] = Det3_012_023 * oneOverDet;
994
995 m[F20] = Det3_123_013 * oneOverDet;
996 m[F21] = Det3_023_013 * mn1OverDet;
997 m[F22] = Det3_013_013 * oneOverDet;
998 m[F23] = Det3_012_013 * mn1OverDet;
999
1000 m[F30] = Det3_123_012 * mn1OverDet;
1001 m[F31] = Det3_023_012 * oneOverDet;
1002 m[F32] = Det3_013_012 * mn1OverDet;
1003 m[F33] = Det3_012_012 * oneOverDet;
1004
1005 return;
1006}
1007
1008void G4ErrorMatrix::invertHaywood5 (G4int & ifail)
1009{
1010 ifail = 0;
1011
1012 // Find all NECESSARY 2x2 dets: (30 of them)
1013
1014 G4double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
1015 G4double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
1016 G4double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
1017 G4double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
1018 G4double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31];
1019 G4double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
1020 G4double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31];
1021 G4double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
1022 G4double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32];
1023 G4double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33];
1024 G4double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
1025 G4double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
1026 G4double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
1027 G4double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
1028 G4double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
1029 G4double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
1030 G4double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
1031 G4double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
1032 G4double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
1033 G4double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43];
1034 G4double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
1035 G4double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
1036 G4double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
1037 G4double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
1038 G4double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
1039 G4double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
1040 G4double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
1041 G4double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
1042 G4double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
1043 G4double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
1044
1045 // Find all NECESSARY 3x3 dets: (40 of them)
1046
1047 G4double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02
1048 + m[M12]*Det2_23_01;
1049 G4double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03
1050 + m[M13]*Det2_23_01;
1051 G4double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04
1052 + m[M14]*Det2_23_01;
1053 G4double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03
1054 + m[M13]*Det2_23_02;
1055 G4double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04
1056 + m[M14]*Det2_23_02;
1057 G4double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04
1058 + m[M14]*Det2_23_03;
1059 G4double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13
1060 + m[M13]*Det2_23_12;
1061 G4double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14
1062 + m[M14]*Det2_23_12;
1063 G4double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14
1064 + m[M14]*Det2_23_13;
1065 G4double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24
1066 + m[M14]*Det2_23_23;
1067 G4double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02
1068 + m[M12]*Det2_24_01;
1069 G4double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03
1070 + m[M13]*Det2_24_01;
1071 G4double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04
1072 + m[M14]*Det2_24_01;
1073 G4double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03
1074 + m[M13]*Det2_24_02;
1075 G4double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04
1076 + m[M14]*Det2_24_02;
1077 G4double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04
1078 + m[M14]*Det2_24_03;
1079 G4double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13
1080 + m[M13]*Det2_24_12;
1081 G4double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14
1082 + m[M14]*Det2_24_12;
1083 G4double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14
1084 + m[M14]*Det2_24_13;
1085 G4double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24
1086 + m[M14]*Det2_24_23;
1087 G4double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02
1088 + m[M12]*Det2_34_01;
1089 G4double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03
1090 + m[M13]*Det2_34_01;
1091 G4double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04
1092 + m[M14]*Det2_34_01;
1093 G4double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03
1094 + m[M13]*Det2_34_02;
1095 G4double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04
1096 + m[M14]*Det2_34_02;
1097 G4double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04
1098 + m[M14]*Det2_34_03;
1099 G4double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13
1100 + m[M13]*Det2_34_12;
1101 G4double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14
1102 + m[M14]*Det2_34_12;
1103 G4double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14
1104 + m[M14]*Det2_34_13;
1105 G4double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24
1106 + m[M14]*Det2_34_23;
1107 G4double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02
1108 + m[M22]*Det2_34_01;
1109 G4double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03
1110 + m[M23]*Det2_34_01;
1111 G4double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04
1112 + m[M24]*Det2_34_01;
1113 G4double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03
1114 + m[M23]*Det2_34_02;
1115 G4double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04
1116 + m[M24]*Det2_34_02;
1117 G4double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04
1118 + m[M24]*Det2_34_03;
1119 G4double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13
1120 + m[M23]*Det2_34_12;
1121 G4double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14
1122 + m[M24]*Det2_34_12;
1123 G4double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14
1124 + m[M24]*Det2_34_13;
1125 G4double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24
1126 + m[M24]*Det2_34_23;
1127
1128 // Find all NECESSARY 4x4 dets: (25 of them)
1129
1130 G4double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023
1131 + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
1132 G4double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024
1133 + m[M02]*Det3_123_014 - m[M04]*Det3_123_012;
1134 G4double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034
1135 + m[M03]*Det3_123_014 - m[M04]*Det3_123_013;
1136 G4double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034
1137 + m[M03]*Det3_123_024 - m[M04]*Det3_123_023;
1138 G4double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134
1139 + m[M03]*Det3_123_124 - m[M04]*Det3_123_123;
1140 G4double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023
1141 + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
1142 G4double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024
1143 + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
1144 G4double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034
1145 + m[M03]*Det3_124_014 - m[M04]*Det3_124_013;
1146 G4double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034
1147 + m[M03]*Det3_124_024 - m[M04]*Det3_124_023;
1148 G4double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134
1149 + m[M03]*Det3_124_124 - m[M04]*Det3_124_123;
1150 G4double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023
1151 + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
1152 G4double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024
1153 + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
1154 G4double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034
1155 + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
1156 G4double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034
1157 + m[M03]*Det3_134_024 - m[M04]*Det3_134_023;
1158 G4double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134
1159 + m[M03]*Det3_134_124 - m[M04]*Det3_134_123;
1160 G4double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023
1161 + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
1162 G4double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024
1163 + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
1164 G4double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034
1165 + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
1166 G4double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034
1167 + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
1168 G4double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134
1169 + m[M03]*Det3_234_124 - m[M04]*Det3_234_123;
1170 G4double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023
1171 + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
1172 G4double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024
1173 + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
1174 G4double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034
1175 + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
1176 G4double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034
1177 + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
1178 G4double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134
1179 + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
1180
1181 // Find the 5x5 det:
1182
1183 G4double det = m[M00]*Det4_1234_1234
1184 - m[M01]*Det4_1234_0234
1185 + m[M02]*Det4_1234_0134
1186 - m[M03]*Det4_1234_0124
1187 + m[M04]*Det4_1234_0123;
1188
1189 if ( det == 0 )
1190 {
1191 ifail = 1;
1192 return;
1193 }
1194
1195 G4double oneOverDet = 1.0/det;
1196 G4double mn1OverDet = - oneOverDet;
1197
1198 m[M00] = Det4_1234_1234 * oneOverDet;
1199 m[M01] = Det4_0234_1234 * mn1OverDet;
1200 m[M02] = Det4_0134_1234 * oneOverDet;
1201 m[M03] = Det4_0124_1234 * mn1OverDet;
1202 m[M04] = Det4_0123_1234 * oneOverDet;
1203
1204 m[M10] = Det4_1234_0234 * mn1OverDet;
1205 m[M11] = Det4_0234_0234 * oneOverDet;
1206 m[M12] = Det4_0134_0234 * mn1OverDet;
1207 m[M13] = Det4_0124_0234 * oneOverDet;
1208 m[M14] = Det4_0123_0234 * mn1OverDet;
1209
1210 m[M20] = Det4_1234_0134 * oneOverDet;
1211 m[M21] = Det4_0234_0134 * mn1OverDet;
1212 m[M22] = Det4_0134_0134 * oneOverDet;
1213 m[M23] = Det4_0124_0134 * mn1OverDet;
1214 m[M24] = Det4_0123_0134 * oneOverDet;
1215
1216 m[M30] = Det4_1234_0124 * mn1OverDet;
1217 m[M31] = Det4_0234_0124 * oneOverDet;
1218 m[M32] = Det4_0134_0124 * mn1OverDet;
1219 m[M33] = Det4_0124_0124 * oneOverDet;
1220 m[M34] = Det4_0123_0124 * mn1OverDet;
1221
1222 m[M40] = Det4_1234_0123 * oneOverDet;
1223 m[M41] = Det4_0234_0123 * mn1OverDet;
1224 m[M42] = Det4_0134_0123 * oneOverDet;
1225 m[M43] = Det4_0124_0123 * mn1OverDet;
1226 m[M44] = Det4_0123_0123 * oneOverDet;
1227
1228 return;
1229}
1230
1231void G4ErrorMatrix::invertHaywood6 (G4int & ifail)
1232{
1233 ifail = 0;
1234
1235 // Find all NECESSARY 2x2 dets: (45 of them)
1236
1237 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1238 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1239 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1240 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1241 G4double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40];
1242 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1243 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1244 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1245 G4double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41];
1246 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1247 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1248 G4double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42];
1249 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1250 G4double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43];
1251 G4double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44];
1252 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1253 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1254 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1255 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1256 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1257 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1258 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1259 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1260 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1261 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1262 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1263 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1264 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1265 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1266 G4double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54];
1267 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1268 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1269 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1270 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1271 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1272 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1273 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1274 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1275 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1276 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1277 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1278 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1279 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1280 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1281 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1282
1283 // Find all NECESSARY 3x3 dets: (80 of them)
1284
1285 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1286 + m[A22]*Det2_34_01;
1287 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1288 + m[A23]*Det2_34_01;
1289 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1290 + m[A24]*Det2_34_01;
1291 G4double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
1292 + m[A25]*Det2_34_01;
1293 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1294 + m[A23]*Det2_34_02;
1295 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1296 + m[A24]*Det2_34_02;
1297 G4double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05
1298 + m[A25]*Det2_34_02;
1299 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1300 + m[A24]*Det2_34_03;
1301 G4double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05
1302 + m[A25]*Det2_34_03;
1303 G4double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05
1304 + m[A25]*Det2_34_04;
1305 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1306 + m[A23]*Det2_34_12;
1307 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1308 + m[A24]*Det2_34_12;
1309 G4double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15
1310 + m[A25]*Det2_34_12;
1311 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1312 + m[A24]*Det2_34_13;
1313 G4double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15
1314 + m[A25]*Det2_34_13;
1315 G4double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15
1316 + m[A25]*Det2_34_14;
1317 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1318 + m[A24]*Det2_34_23;
1319 G4double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
1320 + m[A25]*Det2_34_23;
1321 G4double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
1322 + m[A25]*Det2_34_24;
1323 G4double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
1324 + m[A25]*Det2_34_34;
1325 G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
1326 + m[A22]*Det2_35_01;
1327 G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
1328 + m[A23]*Det2_35_01;
1329 G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
1330 + m[A24]*Det2_35_01;
1331 G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
1332 + m[A25]*Det2_35_01;
1333 G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
1334 + m[A23]*Det2_35_02;
1335 G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
1336 + m[A24]*Det2_35_02;
1337 G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
1338 + m[A25]*Det2_35_02;
1339 G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
1340 + m[A24]*Det2_35_03;
1341 G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
1342 + m[A25]*Det2_35_03;
1343 G4double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05
1344 + m[A25]*Det2_35_04;
1345 G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
1346 + m[A23]*Det2_35_12;
1347 G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
1348 + m[A24]*Det2_35_12;
1349 G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
1350 + m[A25]*Det2_35_12;
1351 G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
1352 + m[A24]*Det2_35_13;
1353 G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
1354 + m[A25]*Det2_35_13;
1355 G4double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15
1356 + m[A25]*Det2_35_14;
1357 G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
1358 + m[A24]*Det2_35_23;
1359 G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
1360 + m[A25]*Det2_35_23;
1361 G4double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25
1362 + m[A25]*Det2_35_24;
1363 G4double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35
1364 + m[A25]*Det2_35_34;
1365 G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
1366 + m[A22]*Det2_45_01;
1367 G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
1368 + m[A23]*Det2_45_01;
1369 G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
1370 + m[A24]*Det2_45_01;
1371 G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
1372 + m[A25]*Det2_45_01;
1373 G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
1374 + m[A23]*Det2_45_02;
1375 G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
1376 + m[A24]*Det2_45_02;
1377 G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
1378 + m[A25]*Det2_45_02;
1379 G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
1380 + m[A24]*Det2_45_03;
1381 G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
1382 + m[A25]*Det2_45_03;
1383 G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
1384 + m[A25]*Det2_45_04;
1385 G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
1386 + m[A23]*Det2_45_12;
1387 G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
1388 + m[A24]*Det2_45_12;
1389 G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
1390 + m[A25]*Det2_45_12;
1391 G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
1392 + m[A24]*Det2_45_13;
1393 G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
1394 + m[A25]*Det2_45_13;
1395 G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
1396 + m[A25]*Det2_45_14;
1397 G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
1398 + m[A24]*Det2_45_23;
1399 G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
1400 + m[A25]*Det2_45_23;
1401 G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
1402 + m[A25]*Det2_45_24;
1403 G4double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35
1404 + m[A25]*Det2_45_34;
1405 G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
1406 + m[A32]*Det2_45_01;
1407 G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
1408 + m[A33]*Det2_45_01;
1409 G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
1410 + m[A34]*Det2_45_01;
1411 G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
1412 + m[A35]*Det2_45_01;
1413 G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
1414 + m[A33]*Det2_45_02;
1415 G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
1416 + m[A34]*Det2_45_02;
1417 G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
1418 + m[A35]*Det2_45_02;
1419 G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
1420 + m[A34]*Det2_45_03;
1421 G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
1422 + m[A35]*Det2_45_03;
1423 G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
1424 + m[A35]*Det2_45_04;
1425 G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
1426 + m[A33]*Det2_45_12;
1427 G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
1428 + m[A34]*Det2_45_12;
1429 G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
1430 + m[A35]*Det2_45_12;
1431 G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
1432 + m[A34]*Det2_45_13;
1433 G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
1434 + m[A35]*Det2_45_13;
1435 G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
1436 + m[A35]*Det2_45_14;
1437 G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
1438 + m[A34]*Det2_45_23;
1439 G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
1440 + m[A35]*Det2_45_23;
1441 G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
1442 + m[A35]*Det2_45_24;
1443 G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
1444 + m[A35]*Det2_45_34;
1445
1446 // Find all NECESSARY 4x4 dets: (75 of them)
1447
1448 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1449 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1450 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1451 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1452 G4double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025
1453 + m[A12]*Det3_234_015 - m[A15]*Det3_234_012;
1454 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1455 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1456 G4double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
1457 + m[A13]*Det3_234_015 - m[A15]*Det3_234_013;
1458 G4double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
1459 + m[A14]*Det3_234_015 - m[A15]*Det3_234_014;
1460 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1461 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1462 G4double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035
1463 + m[A13]*Det3_234_025 - m[A15]*Det3_234_023;
1464 G4double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045
1465 + m[A14]*Det3_234_025 - m[A15]*Det3_234_024;
1466 G4double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045
1467 + m[A14]*Det3_234_035 - m[A15]*Det3_234_034;
1468 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1469 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1470 G4double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135
1471 + m[A13]*Det3_234_125 - m[A15]*Det3_234_123;
1472 G4double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145
1473 + m[A14]*Det3_234_125 - m[A15]*Det3_234_124;
1474 G4double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145
1475 + m[A14]*Det3_234_135 - m[A15]*Det3_234_134;
1476 G4double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245
1477 + m[A14]*Det3_234_235 - m[A15]*Det3_234_234;
1478 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
1479 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1480 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
1481 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1482 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
1483 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1484 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
1485 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1486 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
1487 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1488 G4double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045
1489 + m[A14]*Det3_235_015 - m[A15]*Det3_235_014;
1490 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
1491 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1492 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
1493 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1494 G4double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045
1495 + m[A14]*Det3_235_025 - m[A15]*Det3_235_024;
1496 G4double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045
1497 + m[A14]*Det3_235_035 - m[A15]*Det3_235_034;
1498 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
1499 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1500 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
1501 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1502 G4double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145
1503 + m[A14]*Det3_235_125 - m[A15]*Det3_235_124;
1504 G4double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145
1505 + m[A14]*Det3_235_135 - m[A15]*Det3_235_134;
1506 G4double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245
1507 + m[A14]*Det3_235_235 - m[A15]*Det3_235_234;
1508 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
1509 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1510 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
1511 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1512 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
1513 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1514 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
1515 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1516 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
1517 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1518 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
1519 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1520 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
1521 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1522 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
1523 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1524 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
1525 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1526 G4double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045
1527 + m[A14]*Det3_245_035 - m[A15]*Det3_245_034;
1528 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
1529 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1530 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
1531 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1532 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
1533 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1534 G4double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145
1535 + m[A14]*Det3_245_135 - m[A15]*Det3_245_134;
1536 G4double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245
1537 + m[A14]*Det3_245_235 - m[A15]*Det3_245_234;
1538 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
1539 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1540 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
1541 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1542 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
1543 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1544 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
1545 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1546 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
1547 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1548 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
1549 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1550 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
1551 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1552 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
1553 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1554 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
1555 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1556 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
1557 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1558 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
1559 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1560 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
1561 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1562 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
1563 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1564 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
1565 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1566 G4double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245
1567 + m[A14]*Det3_345_235 - m[A15]*Det3_345_234;
1568 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
1569 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1570 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
1571 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1572 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
1573 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1574 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
1575 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1576 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
1577 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1578 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
1579 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1580 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
1581 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1582 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
1583 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1584 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
1585 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1586 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
1587 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1588 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
1589 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1590 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
1591 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1592 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
1593 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1594 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
1595 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1596 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
1597 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1598
1599 // Find all NECESSARY 5x5 dets: (36 of them)
1600
1601 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
1602 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1603 G4double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
1604 + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123;
1605 G4double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
1606 + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124;
1607 G4double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
1608 + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134;
1609 G4double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
1610 + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234;
1611 G4double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
1612 + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234;
1613 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
1614 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1615 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
1616 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1617 G4double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245
1618 + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124;
1619 G4double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345
1620 + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134;
1621 G4double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345
1622 + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234;
1623 G4double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345
1624 + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234;
1625 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
1626 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1627 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
1628 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1629 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
1630 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1631 G4double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345
1632 + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134;
1633 G4double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345
1634 + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234;
1635 G4double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345
1636 + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234;
1637 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
1638 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1639 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
1640 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1641 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
1642 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1643 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
1644 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1645 G4double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345
1646 + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234;
1647 G4double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345
1648 + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234; //
1649 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
1650 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1651 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
1652 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1653 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
1654 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1655 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
1656 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1657 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
1658 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1659 G4double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345
1660 + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234;
1661 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
1662 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1663 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
1664 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1665 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
1666 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1667 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
1668 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1669 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
1670 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1671 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
1672 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1673
1674 // Find the determinant
1675
1676 G4double det = m[A00]*Det5_12345_12345
1677 - m[A01]*Det5_12345_02345
1678 + m[A02]*Det5_12345_01345
1679 - m[A03]*Det5_12345_01245
1680 + m[A04]*Det5_12345_01235
1681 - m[A05]*Det5_12345_01234;
1682
1683 if ( det == 0 )
1684 {
1685 ifail = 1;
1686 return;
1687 }
1688
1689 G4double oneOverDet = 1.0/det;
1690 G4double mn1OverDet = - oneOverDet;
1691
1692 m[A00] = Det5_12345_12345*oneOverDet;
1693 m[A01] = Det5_02345_12345*mn1OverDet;
1694 m[A02] = Det5_01345_12345*oneOverDet;
1695 m[A03] = Det5_01245_12345*mn1OverDet;
1696 m[A04] = Det5_01235_12345*oneOverDet;
1697 m[A05] = Det5_01234_12345*mn1OverDet;
1698
1699 m[A10] = Det5_12345_02345*mn1OverDet;
1700 m[A11] = Det5_02345_02345*oneOverDet;
1701 m[A12] = Det5_01345_02345*mn1OverDet;
1702 m[A13] = Det5_01245_02345*oneOverDet;
1703 m[A14] = Det5_01235_02345*mn1OverDet;
1704 m[A15] = Det5_01234_02345*oneOverDet;
1705
1706 m[A20] = Det5_12345_01345*oneOverDet;
1707 m[A21] = Det5_02345_01345*mn1OverDet;
1708 m[A22] = Det5_01345_01345*oneOverDet;
1709 m[A23] = Det5_01245_01345*mn1OverDet;
1710 m[A24] = Det5_01235_01345*oneOverDet;
1711 m[A25] = Det5_01234_01345*mn1OverDet;
1712
1713 m[A30] = Det5_12345_01245*mn1OverDet;
1714 m[A31] = Det5_02345_01245*oneOverDet;
1715 m[A32] = Det5_01345_01245*mn1OverDet;
1716 m[A33] = Det5_01245_01245*oneOverDet;
1717 m[A34] = Det5_01235_01245*mn1OverDet;
1718 m[A35] = Det5_01234_01245*oneOverDet;
1719
1720 m[A40] = Det5_12345_01235*oneOverDet;
1721 m[A41] = Det5_02345_01235*mn1OverDet;
1722 m[A42] = Det5_01345_01235*oneOverDet;
1723 m[A43] = Det5_01245_01235*mn1OverDet;
1724 m[A44] = Det5_01235_01235*oneOverDet;
1725 m[A45] = Det5_01234_01235*mn1OverDet;
1726
1727 m[A50] = Det5_12345_01234*mn1OverDet;
1728 m[A51] = Det5_02345_01234*oneOverDet;
1729 m[A52] = Det5_01345_01234*mn1OverDet;
1730 m[A53] = Det5_01245_01234*oneOverDet;
1731 m[A54] = Det5_01235_01234*mn1OverDet;
1732 m[A55] = Det5_01234_01234*oneOverDet;
1733
1734 return;
1735}
Note: See TracBrowser for help on using the repository browser.