00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "globals.hh"
00033
00034 #include <cmath>
00035 #include <iostream>
00036
00037 #include "G4ErrorMatrix.hh"
00038 #include "G4ErrorSymMatrix.hh"
00039
00040
00041
00042 #define SIMPLE_UOP(OPER) \
00043 G4ErrorMatrixIter a=m.begin(); \
00044 G4ErrorMatrixIter e=m.end(); \
00045 for(;a!=e; a++) (*a) OPER t;
00046
00047 #define SIMPLE_BOP(OPER) \
00048 G4ErrorMatrixIter a=m.begin(); \
00049 G4ErrorMatrixConstIter b=mat2.m.begin(); \
00050 G4ErrorMatrixIter e=m.end(); \
00051 for(;a!=e; a++, b++) (*a) OPER (*b);
00052
00053 #define SIMPLE_TOP(OPER) \
00054 G4ErrorMatrixConstIter a=mat1.m.begin(); \
00055 G4ErrorMatrixConstIter b=mat2.m.begin(); \
00056 G4ErrorMatrixIter t=mret.m.begin(); \
00057 G4ErrorMatrixConstIter e=mat1.m.end(); \
00058 for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
00059
00060
00061
00062 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
00063 if (r1!=r2 || c1!=c2) { \
00064 G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
00065 }
00066
00067 #define CHK_DIM_1(c1,r2,fun) \
00068 if (c1!=r2) { \
00069 G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
00070 }
00071
00072
00073
00074 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q)
00075 : m(p*q), nrow(p), ncol(q)
00076 {
00077 size = nrow * ncol;
00078 }
00079
00080 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q,G4int init)
00081 : m(p*q), nrow(p), ncol(q)
00082 {
00083 size = nrow * ncol;
00084
00085 if (size > 0)
00086 {
00087 switch(init)
00088 {
00089 case 0:
00090 break;
00091
00092 case 1:
00093 {
00094 if ( ncol == nrow )
00095 {
00096 G4ErrorMatrixIter a = m.begin();
00097 G4ErrorMatrixIter b = m.end();
00098 for( ; a<b; a+=(ncol+1)) *a = 1.0;
00099 } else {
00100 error("Invalid dimension in G4ErrorMatrix(G4int,G4int,1).");
00101 }
00102 break;
00103 }
00104 default:
00105 error("G4ErrorMatrix: initialization must be either 0 or 1.");
00106 }
00107 }
00108 }
00109
00110
00111
00112
00113 G4ErrorMatrix::~G4ErrorMatrix()
00114 {
00115 }
00116
00117 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatrix &mat1)
00118 : m(mat1.size), nrow(mat1.nrow), ncol(mat1.ncol), size(mat1.size)
00119 {
00120 m = mat1.m;
00121 }
00122
00123 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymMatrix &mat1)
00124 : m(mat1.nrow*mat1.nrow), nrow(mat1.nrow), ncol(mat1.nrow)
00125 {
00126 size = nrow * ncol;
00127
00128 G4int n = ncol;
00129 G4ErrorMatrixConstIter sjk = mat1.m.begin();
00130 G4ErrorMatrixIter m1j = m.begin();
00131 G4ErrorMatrixIter mj = m.begin();
00132
00133 for(G4int j=1;j<=nrow;j++)
00134 {
00135 G4ErrorMatrixIter mjk = mj;
00136 G4ErrorMatrixIter mkj = m1j;
00137 for(G4int k=1;k<=j;k++)
00138 {
00139 *(mjk++) = *sjk;
00140 if(j!=k) *mkj = *sjk;
00141 sjk++;
00142 mkj += n;
00143 }
00144 mj += n;
00145 m1j++;
00146 }
00147 }
00148
00149
00150
00151
00152
00153
00154
00155 G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row, G4int max_row,
00156 G4int min_col,G4int max_col) const
00157 {
00158 G4ErrorMatrix mret(max_row-min_row+1,max_col-min_col+1);
00159 if(max_row > num_row() || max_col >num_col())
00160 { error("G4ErrorMatrix::sub: Index out of range"); }
00161 G4ErrorMatrixIter a = mret.m.begin();
00162 G4int nc = num_col();
00163 G4ErrorMatrixConstIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
00164
00165 for(G4int irow=1; irow<=mret.num_row(); irow++)
00166 {
00167 G4ErrorMatrixConstIter brc = b1;
00168 for(G4int icol=1; icol<=mret.num_col(); icol++)
00169 {
00170 *(a++) = *(brc++);
00171 }
00172 b1 += nc;
00173 }
00174 return mret;
00175 }
00176
00177 void G4ErrorMatrix::sub(G4int row,G4int col,const G4ErrorMatrix &mat1)
00178 {
00179 if(row <1 || row+mat1.num_row()-1 > num_row() ||
00180 col <1 || col+mat1.num_col()-1 > num_col() )
00181 { error("G4ErrorMatrix::sub: Index out of range"); }
00182 G4ErrorMatrixConstIter a = mat1.m.begin();
00183 G4int nc = num_col();
00184 G4ErrorMatrixIter b1 = m.begin() + (row - 1) * nc + col - 1;
00185
00186 for(G4int irow=1; irow<=mat1.num_row(); irow++)
00187 {
00188 G4ErrorMatrixIter brc = b1;
00189 for(G4int icol=1; icol<=mat1.num_col(); icol++)
00190 {
00191 *(brc++) = *(a++);
00192 }
00193 b1 += nc;
00194 }
00195 }
00196
00197
00198
00199
00200
00201 G4ErrorMatrix dsum(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
00202 {
00203 G4ErrorMatrix mret(mat1.num_row() + mat2.num_row(),
00204 mat1.num_col() + mat2.num_col(), 0);
00205 mret.sub(1,1,mat1);
00206 mret.sub(mat1.num_row()+1,mat1.num_col()+1,mat2);
00207 return mret;
00208 }
00209
00210
00211
00212
00213
00214
00215 G4ErrorMatrix G4ErrorMatrix::operator- () const
00216 {
00217 G4ErrorMatrix mat2(nrow, ncol);
00218 G4ErrorMatrixConstIter a=m.begin();
00219 G4ErrorMatrixIter b=mat2.m.begin();
00220 G4ErrorMatrixConstIter e=m.end();
00221 for(;a<e; a++, b++) (*b) = -(*a);
00222 return mat2;
00223 }
00224
00225 G4ErrorMatrix operator+(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
00226 {
00227 G4ErrorMatrix mret(mat1.nrow, mat1.ncol);
00228 CHK_DIM_2(mat1.num_row(),mat2.num_row(), mat1.num_col(),mat2.num_col(),+);
00229 SIMPLE_TOP(+)
00230 return mret;
00231 }
00232
00233
00234
00235
00236
00237 G4ErrorMatrix operator-(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
00238 {
00239 G4ErrorMatrix mret(mat1.num_row(), mat1.num_col());
00240 CHK_DIM_2(mat1.num_row(),mat2.num_row(),
00241 mat1.num_col(),mat2.num_col(),-);
00242 SIMPLE_TOP(-)
00243 return mret;
00244 }
00245
00246
00247
00248
00249
00250
00251 G4ErrorMatrix operator/(const G4ErrorMatrix &mat1,G4double t)
00252 {
00253 G4ErrorMatrix mret(mat1);
00254 mret /= t;
00255 return mret;
00256 }
00257
00258 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,G4double t)
00259 {
00260 G4ErrorMatrix mret(mat1);
00261 mret *= t;
00262 return mret;
00263 }
00264
00265 G4ErrorMatrix operator*(G4double t,const G4ErrorMatrix &mat1)
00266 {
00267 G4ErrorMatrix mret(mat1);
00268 mret *= t;
00269 return mret;
00270 }
00271
00272 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
00273 {
00274
00275 G4ErrorMatrix mret(mat1.nrow,mat2.ncol,0);
00276 CHK_DIM_1(mat1.ncol,mat2.nrow,*);
00277
00278 G4int m1cols = mat1.ncol;
00279 G4int m2cols = mat2.ncol;
00280
00281 for (G4int i=0; i<mat1.nrow; i++)
00282 {
00283 for (G4int j=0; j<m1cols; j++)
00284 {
00285 G4double temp = mat1.m[i*m1cols+j];
00286 G4ErrorMatrixIter pt = mret.m.begin() + i*m2cols;
00287
00288
00289 G4ErrorMatrixConstIter pb = mat2.m.begin() + m2cols*j;
00290 const G4ErrorMatrixConstIter pblast = pb + m2cols;
00291 while (pb < pblast)
00292 {
00293 (*pt) += temp * (*pb);
00294 pb++;
00295 pt++;
00296 }
00297 }
00298 }
00299 return mret;
00300 }
00301
00302
00303
00304
00305
00306 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorMatrix &mat2)
00307 {
00308 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
00309 SIMPLE_BOP(+=)
00310 return (*this);
00311 }
00312
00313 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorMatrix &mat2)
00314 {
00315 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
00316 SIMPLE_BOP(-=)
00317 return (*this);
00318 }
00319
00320 G4ErrorMatrix & G4ErrorMatrix::operator/=(G4double t)
00321 {
00322 SIMPLE_UOP(/=)
00323 return (*this);
00324 }
00325
00326 G4ErrorMatrix & G4ErrorMatrix::operator*=(G4double t)
00327 {
00328 SIMPLE_UOP(*=)
00329 return (*this);
00330 }
00331
00332 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorMatrix &mat1)
00333 {
00334 if (&mat1 == this) { return *this; }
00335
00336 if(mat1.nrow*mat1.ncol != size)
00337 {
00338 size = mat1.nrow * mat1.ncol;
00339 m.resize(size);
00340 }
00341 nrow = mat1.nrow;
00342 ncol = mat1.ncol;
00343 m = mat1.m;
00344 return (*this);
00345 }
00346
00347
00348
00349
00350 std::ostream& operator<<(std::ostream &os, const G4ErrorMatrix &q)
00351 {
00352 os << "\n";
00353
00354
00355
00356
00357 G4int width;
00358 if(os.flags() & std::ios::fixed)
00359 { width = os.precision()+3; }
00360 else
00361 { width = os.precision()+7; }
00362 for(G4int irow = 1; irow<= q.num_row(); irow++)
00363 {
00364 for(G4int icol = 1; icol <= q.num_col(); icol++)
00365 {
00366 os.width(width);
00367 os << q(irow,icol) << " ";
00368 }
00369 os << G4endl;
00370 }
00371 return os;
00372 }
00373
00374 G4ErrorMatrix G4ErrorMatrix::T() const
00375 {
00376 G4ErrorMatrix mret(ncol,nrow);
00377 G4ErrorMatrixConstIter pl = m.end();
00378 G4ErrorMatrixConstIter pme = m.begin();
00379 G4ErrorMatrixIter pt = mret.m.begin();
00380 G4ErrorMatrixIter ptl = mret.m.end();
00381 for (; pme < pl; pme++, pt+=nrow)
00382 {
00383 if (pt >= ptl) { pt -= (size-1); }
00384 (*pt) = (*pme);
00385 }
00386 return mret;
00387 }
00388
00389 G4ErrorMatrix
00390 G4ErrorMatrix::apply(G4double (*f)(G4double, G4int, G4int)) const
00391 {
00392 G4ErrorMatrix mret(num_row(),num_col());
00393 G4ErrorMatrixConstIter a = m.begin();
00394 G4ErrorMatrixIter b = mret.m.begin();
00395 for(G4int ir=1;ir<=num_row();ir++)
00396 {
00397 for(G4int ic=1;ic<=num_col();ic++)
00398 {
00399 *(b++) = (*f)(*(a++), ir, ic);
00400 }
00401 }
00402 return mret;
00403 }
00404
00405 G4int G4ErrorMatrix::dfinv_matrix(G4int *ir) {
00406 if (num_col()!=num_row())
00407 { error("dfinv_matrix: G4ErrorMatrix is not NxN"); }
00408 G4int n = num_col();
00409 if (n==1) { return 0; }
00410
00411 G4double s31, s32;
00412 G4double s33, s34;
00413
00414 G4ErrorMatrixIter m11 = m.begin();
00415 G4ErrorMatrixIter m12 = m11 + 1;
00416 G4ErrorMatrixIter m21 = m11 + n;
00417 G4ErrorMatrixIter m22 = m12 + n;
00418 *m21 = -(*m22) * (*m11) * (*m21);
00419 *m12 = -(*m12);
00420 if (n>2)
00421 {
00422 G4ErrorMatrixIter mi = m.begin() + 2 * n;
00423 G4ErrorMatrixIter mii= m.begin() + 2 * n + 2;
00424 G4ErrorMatrixIter mimim = m.begin() + n + 1;
00425 for (G4int i=3;i<=n;i++)
00426 {
00427 G4int im2 = i - 2;
00428 G4ErrorMatrixIter mj = m.begin();
00429 G4ErrorMatrixIter mji = mj + i - 1;
00430 G4ErrorMatrixIter mij = mi;
00431 for (G4int j=1;j<=im2;j++)
00432 {
00433 s31 = 0.0;
00434 s32 = *mji;
00435 G4ErrorMatrixIter mkj = mj + j - 1;
00436 G4ErrorMatrixIter mik = mi + j - 1;
00437 G4ErrorMatrixIter mjkp = mj + j;
00438 G4ErrorMatrixIter mkpi = mj + n + i - 1;
00439 for (G4int k=j;k<=im2;k++)
00440 {
00441 s31 += (*mkj) * (*(mik++));
00442 s32 += (*(mjkp++)) * (*mkpi);
00443 mkj += n;
00444 mkpi += n;
00445 }
00446 *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
00447 *mji = -s32;
00448 mj += n;
00449 mji += n;
00450 mij++;
00451 }
00452 *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
00453 *(mimim+1) = -(*(mimim+1));
00454 mi += n;
00455 mimim += (n+1);
00456 mii += (n+1);
00457 }
00458 }
00459 G4ErrorMatrixIter mi = m.begin();
00460 G4ErrorMatrixIter mii = m.begin();
00461 for (G4int i=1;i<n;i++)
00462 {
00463 G4int ni = n - i;
00464 G4ErrorMatrixIter mij = mi;
00465 G4int j;
00466 for (j=1; j<=i;j++)
00467 {
00468 s33 = *mij;
00469 G4ErrorMatrixIter mikj = mi + n + j - 1;
00470 G4ErrorMatrixIter miik = mii + 1;
00471 G4ErrorMatrixIter min_end = mi + n;
00472 for (;miik<min_end;)
00473 {
00474 s33 += (*mikj) * (*(miik++));
00475 mikj += n;
00476 }
00477 *(mij++) = s33;
00478 }
00479 for (j=1;j<=ni;j++)
00480 {
00481 s34 = 0.0;
00482 G4ErrorMatrixIter miik = mii + j;
00483 G4ErrorMatrixIter mikij = mii + j * n + j;
00484 for (G4int k=j;k<=ni;k++)
00485 {
00486 s34 += *mikij * (*(miik++));
00487 mikij += n;
00488 }
00489 *(mii+j) = s34;
00490 }
00491 mi += n;
00492 mii += (n+1);
00493 }
00494 G4int nxch = ir[n];
00495 if (nxch==0) return 0;
00496 for (G4int mq=1;mq<=nxch;mq++)
00497 {
00498 G4int k = nxch - mq + 1;
00499 G4int ij = ir[k];
00500 G4int i = ij >> 12;
00501 G4int j = ij%4096;
00502 G4ErrorMatrixIter mki = m.begin() + i - 1;
00503 G4ErrorMatrixIter mkj = m.begin() + j - 1;
00504 for (k=1; k<=n;k++)
00505 {
00506
00507
00508 G4double ti = *mki;
00509 *mki = *mkj;
00510 *mkj = ti;
00511 mki += n;
00512 mkj += n;
00513 }
00514 }
00515 return 0;
00516 }
00517
00518 G4int G4ErrorMatrix::dfact_matrix(G4double &det, G4int *ir)
00519 {
00520 if (ncol!=nrow)
00521 error("dfact_matrix: G4ErrorMatrix is not NxN");
00522
00523 G4int ifail, jfail;
00524 G4int n = ncol;
00525
00526 G4double tf;
00527 G4double g1 = 1.0e-19, g2 = 1.0e19;
00528
00529 G4double p, q, t;
00530 G4double s11, s12;
00531
00532 G4double epsilon = 8*DBL_EPSILON;
00533
00534
00535
00536
00537 G4int normal = 0, imposs = -1;
00538 G4int jrange = 0, jover = 1, junder = -1;
00539 ifail = normal;
00540 jfail = jrange;
00541 G4int nxch = 0;
00542 det = 1.0;
00543 G4ErrorMatrixIter mj = m.begin();
00544 G4ErrorMatrixIter mjj = mj;
00545 for (G4int j=1;j<=n;j++)
00546 {
00547 G4int k = j;
00548 p = (std::fabs(*mjj));
00549 if (j!=n) {
00550 G4ErrorMatrixIter mij = mj + n + j - 1;
00551 for (G4int i=j+1;i<=n;i++)
00552 {
00553 q = (std::fabs(*(mij)));
00554 if (q > p)
00555 {
00556 k = i;
00557 p = q;
00558 }
00559 mij += n;
00560 }
00561 if (k==j)
00562 {
00563 if (p <= epsilon)
00564 {
00565 det = 0;
00566 ifail = imposs;
00567 jfail = jrange;
00568 return ifail;
00569 }
00570 det = -det;
00571
00572 }
00573 G4ErrorMatrixIter mjl = mj;
00574 G4ErrorMatrixIter mkl = m.begin() + (k-1)*n;
00575 for (G4int l=1;l<=n;l++)
00576 {
00577 tf = *mjl;
00578 *(mjl++) = *mkl;
00579 *(mkl++) = tf;
00580 }
00581 nxch = nxch + 1;
00582 ir[nxch] = (((j)<<12)+(k));
00583 }
00584 else
00585 {
00586 if (p <= epsilon)
00587 {
00588 det = 0.0;
00589 ifail = imposs;
00590 jfail = jrange;
00591 return ifail;
00592 }
00593 }
00594 det *= *mjj;
00595 *mjj = 1.0 / *mjj;
00596 t = (std::fabs(det));
00597 if (t < g1)
00598 {
00599 det = 0.0;
00600 if (jfail == jrange) jfail = junder;
00601 }
00602 else if (t > g2)
00603 {
00604 det = 1.0;
00605 if (jfail==jrange) jfail = jover;
00606 }
00607 if (j!=n)
00608 {
00609 G4ErrorMatrixIter mk = mj + n;
00610 G4ErrorMatrixIter mkjp = mk + j;
00611 G4ErrorMatrixIter mjk = mj + j;
00612 for (k=j+1;k<=n;k++)
00613 {
00614 s11 = - (*mjk);
00615 s12 = - (*mkjp);
00616 if (j!=1)
00617 {
00618 G4ErrorMatrixIter mik = m.begin() + k - 1;
00619 G4ErrorMatrixIter mijp = m.begin() + j;
00620 G4ErrorMatrixIter mki = mk;
00621 G4ErrorMatrixIter mji = mj;
00622 for (G4int i=1;i<j;i++)
00623 {
00624 s11 += (*mik) * (*(mji++));
00625 s12 += (*mijp) * (*(mki++));
00626 mik += n;
00627 mijp += n;
00628 }
00629 }
00630 *(mjk++) = -s11 * (*mjj);
00631 *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
00632 mk += n;
00633 mkjp += n;
00634 }
00635 }
00636 mj += n;
00637 mjj += (n+1);
00638 }
00639 if (nxch%2==1) det = -det;
00640 if (jfail !=jrange) det = 0.0;
00641 ir[n] = nxch;
00642 return 0;
00643 }
00644
00645 void G4ErrorMatrix::invert(G4int &ierr)
00646 {
00647 if(ncol != nrow)
00648 { error("G4ErrorMatrix::invert: G4ErrorMatrix is not NxN"); }
00649
00650 static G4int max_array = 20;
00651 static G4int *ir = new G4int [max_array+1];
00652
00653 if (ncol > max_array)
00654 {
00655 delete [] ir;
00656 max_array = nrow;
00657 ir = new G4int [max_array+1];
00658 }
00659 G4double t1, t2, t3;
00660 G4double det, temp, ss;
00661 G4int ifail;
00662 switch(nrow)
00663 {
00664 case 3:
00665 G4double c11,c12,c13,c21,c22,c23,c31,c32,c33;
00666 ifail = 0;
00667 c11 = (*(m.begin()+4)) * (*(m.begin()+8))
00668 - (*(m.begin()+5)) * (*(m.begin()+7));
00669 c12 = (*(m.begin()+5)) * (*(m.begin()+6))
00670 - (*(m.begin()+3)) * (*(m.begin()+8));
00671 c13 = (*(m.begin()+3)) * (*(m.begin()+7))
00672 - (*(m.begin()+4)) * (*(m.begin()+6));
00673 c21 = (*(m.begin()+7)) * (*(m.begin()+2))
00674 - (*(m.begin()+8)) * (*(m.begin()+1));
00675 c22 = (*(m.begin()+8)) * (*m.begin())
00676 - (*(m.begin()+6)) * (*(m.begin()+2));
00677 c23 = (*(m.begin()+6)) * (*(m.begin()+1))
00678 - (*(m.begin()+7)) * (*m.begin());
00679 c31 = (*(m.begin()+1)) * (*(m.begin()+5))
00680 - (*(m.begin()+2)) * (*(m.begin()+4));
00681 c32 = (*(m.begin()+2)) * (*(m.begin()+3))
00682 - (*m.begin()) * (*(m.begin()+5));
00683 c33 = (*m.begin()) * (*(m.begin()+4))
00684 - (*(m.begin()+1)) * (*(m.begin()+3));
00685 t1 = std::fabs(*m.begin());
00686 t2 = std::fabs(*(m.begin()+3));
00687 t3 = std::fabs(*(m.begin()+6));
00688 if (t1 >= t2)
00689 {
00690 if (t3 >= t1)
00691 {
00692 temp = *(m.begin()+6);
00693 det = c23*c12-c22*c13;
00694 }
00695 else
00696 {
00697 temp = *(m.begin());
00698 det = c22*c33-c23*c32;
00699 }
00700 }
00701 else if (t3 >= t2)
00702 {
00703 temp = *(m.begin()+6);
00704 det = c23*c12-c22*c13;
00705 }
00706 else
00707 {
00708 temp = *(m.begin()+3);
00709 det = c13*c32-c12*c33;
00710 }
00711 if (det==0)
00712 {
00713 ierr = 1;
00714 return;
00715 }
00716 {
00717 ss = temp/det;
00718 G4ErrorMatrixIter mq = m.begin();
00719 *(mq++) = ss*c11;
00720 *(mq++) = ss*c21;
00721 *(mq++) = ss*c31;
00722 *(mq++) = ss*c12;
00723 *(mq++) = ss*c22;
00724 *(mq++) = ss*c32;
00725 *(mq++) = ss*c13;
00726 *(mq++) = ss*c23;
00727 *(mq) = ss*c33;
00728 }
00729 break;
00730 case 2:
00731 ifail = 0;
00732 det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
00733 if (det==0)
00734 {
00735 ierr = 1;
00736 return;
00737 }
00738 ss = 1.0/det;
00739 temp = ss*(*(m.begin()+3));
00740 *(m.begin()+1) *= -ss;
00741 *(m.begin()+2) *= -ss;
00742 *(m.begin()+3) = ss*(*m.begin());
00743 *(m.begin()) = temp;
00744 break;
00745 case 1:
00746 ifail = 0;
00747 if ((*(m.begin()))==0)
00748 {
00749 ierr = 1;
00750 return;
00751 }
00752 *(m.begin()) = 1.0/(*(m.begin()));
00753 break;
00754 case 4:
00755 invertHaywood4(ierr);
00756 return;
00757 case 5:
00758 invertHaywood5(ierr);
00759 return;
00760 case 6:
00761 invertHaywood6(ierr);
00762 return;
00763 default:
00764 ifail = dfact_matrix(det, ir);
00765 if(ifail) {
00766 ierr = 1;
00767 return;
00768 }
00769 dfinv_matrix(ir);
00770 break;
00771 }
00772 ierr = 0;
00773 return;
00774 }
00775
00776 G4double G4ErrorMatrix::determinant() const
00777 {
00778 static G4int max_array = 20;
00779 static G4int *ir = new G4int [max_array+1];
00780 if(ncol != nrow)
00781 { error("G4ErrorMatrix::determinant: G4ErrorMatrix is not NxN"); }
00782 if (ncol > max_array)
00783 {
00784 delete [] ir;
00785 max_array = nrow;
00786 ir = new G4int [max_array+1];
00787 }
00788 G4double det;
00789 G4ErrorMatrix mt(*this);
00790 G4int i = mt.dfact_matrix(det, ir);
00791 if(i==0) return det;
00792 return 0;
00793 }
00794
00795 G4double G4ErrorMatrix::trace() const
00796 {
00797 G4double t = 0.0;
00798 for (G4ErrorMatrixConstIter d = m.begin(); d < m.end(); d += (ncol+1) )
00799 {
00800 t += *d;
00801 }
00802 return t;
00803 }
00804
00805 void G4ErrorMatrix::error(const char *msg)
00806 {
00807 G4cerr << msg << G4endl;
00808 G4cerr << "---Exiting to System." << G4endl;
00809 abort();
00810 }
00811
00812
00813
00814
00815
00816
00817 #define A00 0
00818 #define A01 1
00819 #define A02 2
00820 #define A03 3
00821 #define A04 4
00822 #define A05 5
00823
00824 #define A10 6
00825 #define A11 7
00826 #define A12 8
00827 #define A13 9
00828 #define A14 10
00829 #define A15 11
00830
00831 #define A20 12
00832 #define A21 13
00833 #define A22 14
00834 #define A23 15
00835 #define A24 16
00836 #define A25 17
00837
00838 #define A30 18
00839 #define A31 19
00840 #define A32 20
00841 #define A33 21
00842 #define A34 22
00843 #define A35 23
00844
00845 #define A40 24
00846 #define A41 25
00847 #define A42 26
00848 #define A43 27
00849 #define A44 28
00850 #define A45 29
00851
00852 #define A50 30
00853 #define A51 31
00854 #define A52 32
00855 #define A53 33
00856 #define A54 34
00857 #define A55 35
00858
00859 #define M00 0
00860 #define M01 1
00861 #define M02 2
00862 #define M03 3
00863 #define M04 4
00864
00865 #define M10 5
00866 #define M11 6
00867 #define M12 7
00868 #define M13 8
00869 #define M14 9
00870
00871 #define M20 10
00872 #define M21 11
00873 #define M22 12
00874 #define M23 13
00875 #define M24 14
00876
00877 #define M30 15
00878 #define M31 16
00879 #define M32 17
00880 #define M33 18
00881 #define M34 19
00882
00883 #define M40 20
00884 #define M41 21
00885 #define M42 22
00886 #define M43 23
00887 #define M44 24
00888
00889 #define F00 0
00890 #define F01 1
00891 #define F02 2
00892 #define F03 3
00893
00894 #define F10 4
00895 #define F11 5
00896 #define F12 6
00897 #define F13 7
00898
00899 #define F20 8
00900 #define F21 9
00901 #define F22 10
00902 #define F23 11
00903
00904 #define F30 12
00905 #define F31 13
00906 #define F32 14
00907 #define F33 15
00908
00909
00910 void G4ErrorMatrix::invertHaywood4 (G4int & ifail)
00911 {
00912 ifail = 0;
00913
00914
00915
00916 G4double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
00917 G4double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
00918 G4double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20];
00919 G4double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21];
00920 G4double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22];
00921 G4double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
00922 G4double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
00923 G4double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
00924 G4double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
00925 G4double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
00926 G4double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
00927 G4double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32];
00928 G4double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
00929 G4double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
00930 G4double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
00931 G4double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
00932 G4double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
00933 G4double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
00934
00935
00936
00937 G4double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02
00938 + m[F02]*Det2_12_01;
00939 G4double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03
00940 + m[F03]*Det2_12_01;
00941 G4double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03
00942 + m[F03]*Det2_12_02;
00943 G4double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13
00944 + m[F03]*Det2_12_12;
00945 G4double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02
00946 + m[F02]*Det2_13_01;
00947 G4double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
00948 + m[F03]*Det2_13_01;
00949 G4double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
00950 + m[F03]*Det2_13_02;
00951 G4double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
00952 + m[F03]*Det2_13_12;
00953 G4double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02
00954 + m[F02]*Det2_23_01;
00955 G4double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
00956 + m[F03]*Det2_23_01;
00957 G4double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
00958 + m[F03]*Det2_23_02;
00959 G4double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
00960 + m[F03]*Det2_23_12;
00961 G4double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02
00962 + m[F12]*Det2_23_01;
00963 G4double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03
00964 + m[F13]*Det2_23_01;
00965 G4double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03
00966 + m[F13]*Det2_23_02;
00967 G4double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13
00968 + m[F13]*Det2_23_12;
00969
00970
00971
00972 G4double det = m[F00]*Det3_123_123
00973 - m[F01]*Det3_123_023
00974 + m[F02]*Det3_123_013
00975 - m[F03]*Det3_123_012;
00976
00977 if ( det == 0 )
00978 {
00979 ifail = 1;
00980 return;
00981 }
00982
00983 G4double oneOverDet = 1.0/det;
00984 G4double mn1OverDet = - oneOverDet;
00985
00986 m[F00] = Det3_123_123 * oneOverDet;
00987 m[F01] = Det3_023_123 * mn1OverDet;
00988 m[F02] = Det3_013_123 * oneOverDet;
00989 m[F03] = Det3_012_123 * mn1OverDet;
00990
00991 m[F10] = Det3_123_023 * mn1OverDet;
00992 m[F11] = Det3_023_023 * oneOverDet;
00993 m[F12] = Det3_013_023 * mn1OverDet;
00994 m[F13] = Det3_012_023 * oneOverDet;
00995
00996 m[F20] = Det3_123_013 * oneOverDet;
00997 m[F21] = Det3_023_013 * mn1OverDet;
00998 m[F22] = Det3_013_013 * oneOverDet;
00999 m[F23] = Det3_012_013 * mn1OverDet;
01000
01001 m[F30] = Det3_123_012 * mn1OverDet;
01002 m[F31] = Det3_023_012 * oneOverDet;
01003 m[F32] = Det3_013_012 * mn1OverDet;
01004 m[F33] = Det3_012_012 * oneOverDet;
01005
01006 return;
01007 }
01008
01009 void G4ErrorMatrix::invertHaywood5 (G4int & ifail)
01010 {
01011 ifail = 0;
01012
01013
01014
01015 G4double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
01016 G4double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
01017 G4double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
01018 G4double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
01019 G4double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31];
01020 G4double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
01021 G4double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31];
01022 G4double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
01023 G4double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32];
01024 G4double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33];
01025 G4double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
01026 G4double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
01027 G4double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
01028 G4double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
01029 G4double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
01030 G4double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
01031 G4double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
01032 G4double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
01033 G4double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
01034 G4double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43];
01035 G4double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
01036 G4double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
01037 G4double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
01038 G4double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
01039 G4double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
01040 G4double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
01041 G4double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
01042 G4double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
01043 G4double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
01044 G4double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
01045
01046
01047
01048 G4double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02
01049 + m[M12]*Det2_23_01;
01050 G4double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03
01051 + m[M13]*Det2_23_01;
01052 G4double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04
01053 + m[M14]*Det2_23_01;
01054 G4double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03
01055 + m[M13]*Det2_23_02;
01056 G4double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04
01057 + m[M14]*Det2_23_02;
01058 G4double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04
01059 + m[M14]*Det2_23_03;
01060 G4double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13
01061 + m[M13]*Det2_23_12;
01062 G4double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14
01063 + m[M14]*Det2_23_12;
01064 G4double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14
01065 + m[M14]*Det2_23_13;
01066 G4double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24
01067 + m[M14]*Det2_23_23;
01068 G4double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02
01069 + m[M12]*Det2_24_01;
01070 G4double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03
01071 + m[M13]*Det2_24_01;
01072 G4double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04
01073 + m[M14]*Det2_24_01;
01074 G4double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03
01075 + m[M13]*Det2_24_02;
01076 G4double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04
01077 + m[M14]*Det2_24_02;
01078 G4double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04
01079 + m[M14]*Det2_24_03;
01080 G4double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13
01081 + m[M13]*Det2_24_12;
01082 G4double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14
01083 + m[M14]*Det2_24_12;
01084 G4double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14
01085 + m[M14]*Det2_24_13;
01086 G4double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24
01087 + m[M14]*Det2_24_23;
01088 G4double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02
01089 + m[M12]*Det2_34_01;
01090 G4double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03
01091 + m[M13]*Det2_34_01;
01092 G4double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04
01093 + m[M14]*Det2_34_01;
01094 G4double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03
01095 + m[M13]*Det2_34_02;
01096 G4double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04
01097 + m[M14]*Det2_34_02;
01098 G4double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04
01099 + m[M14]*Det2_34_03;
01100 G4double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13
01101 + m[M13]*Det2_34_12;
01102 G4double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14
01103 + m[M14]*Det2_34_12;
01104 G4double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14
01105 + m[M14]*Det2_34_13;
01106 G4double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24
01107 + m[M14]*Det2_34_23;
01108 G4double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02
01109 + m[M22]*Det2_34_01;
01110 G4double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03
01111 + m[M23]*Det2_34_01;
01112 G4double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04
01113 + m[M24]*Det2_34_01;
01114 G4double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03
01115 + m[M23]*Det2_34_02;
01116 G4double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04
01117 + m[M24]*Det2_34_02;
01118 G4double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04
01119 + m[M24]*Det2_34_03;
01120 G4double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13
01121 + m[M23]*Det2_34_12;
01122 G4double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14
01123 + m[M24]*Det2_34_12;
01124 G4double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14
01125 + m[M24]*Det2_34_13;
01126 G4double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24
01127 + m[M24]*Det2_34_23;
01128
01129
01130
01131 G4double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023
01132 + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
01133 G4double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024
01134 + m[M02]*Det3_123_014 - m[M04]*Det3_123_012;
01135 G4double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034
01136 + m[M03]*Det3_123_014 - m[M04]*Det3_123_013;
01137 G4double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034
01138 + m[M03]*Det3_123_024 - m[M04]*Det3_123_023;
01139 G4double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134
01140 + m[M03]*Det3_123_124 - m[M04]*Det3_123_123;
01141 G4double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023
01142 + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
01143 G4double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024
01144 + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
01145 G4double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034
01146 + m[M03]*Det3_124_014 - m[M04]*Det3_124_013;
01147 G4double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034
01148 + m[M03]*Det3_124_024 - m[M04]*Det3_124_023;
01149 G4double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134
01150 + m[M03]*Det3_124_124 - m[M04]*Det3_124_123;
01151 G4double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023
01152 + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
01153 G4double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024
01154 + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
01155 G4double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034
01156 + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
01157 G4double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034
01158 + m[M03]*Det3_134_024 - m[M04]*Det3_134_023;
01159 G4double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134
01160 + m[M03]*Det3_134_124 - m[M04]*Det3_134_123;
01161 G4double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023
01162 + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
01163 G4double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024
01164 + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
01165 G4double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034
01166 + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
01167 G4double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034
01168 + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
01169 G4double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134
01170 + m[M03]*Det3_234_124 - m[M04]*Det3_234_123;
01171 G4double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023
01172 + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
01173 G4double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024
01174 + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
01175 G4double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034
01176 + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
01177 G4double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034
01178 + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
01179 G4double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134
01180 + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
01181
01182
01183
01184 G4double det = m[M00]*Det4_1234_1234
01185 - m[M01]*Det4_1234_0234
01186 + m[M02]*Det4_1234_0134
01187 - m[M03]*Det4_1234_0124
01188 + m[M04]*Det4_1234_0123;
01189
01190 if ( det == 0 )
01191 {
01192 ifail = 1;
01193 return;
01194 }
01195
01196 G4double oneOverDet = 1.0/det;
01197 G4double mn1OverDet = - oneOverDet;
01198
01199 m[M00] = Det4_1234_1234 * oneOverDet;
01200 m[M01] = Det4_0234_1234 * mn1OverDet;
01201 m[M02] = Det4_0134_1234 * oneOverDet;
01202 m[M03] = Det4_0124_1234 * mn1OverDet;
01203 m[M04] = Det4_0123_1234 * oneOverDet;
01204
01205 m[M10] = Det4_1234_0234 * mn1OverDet;
01206 m[M11] = Det4_0234_0234 * oneOverDet;
01207 m[M12] = Det4_0134_0234 * mn1OverDet;
01208 m[M13] = Det4_0124_0234 * oneOverDet;
01209 m[M14] = Det4_0123_0234 * mn1OverDet;
01210
01211 m[M20] = Det4_1234_0134 * oneOverDet;
01212 m[M21] = Det4_0234_0134 * mn1OverDet;
01213 m[M22] = Det4_0134_0134 * oneOverDet;
01214 m[M23] = Det4_0124_0134 * mn1OverDet;
01215 m[M24] = Det4_0123_0134 * oneOverDet;
01216
01217 m[M30] = Det4_1234_0124 * mn1OverDet;
01218 m[M31] = Det4_0234_0124 * oneOverDet;
01219 m[M32] = Det4_0134_0124 * mn1OverDet;
01220 m[M33] = Det4_0124_0124 * oneOverDet;
01221 m[M34] = Det4_0123_0124 * mn1OverDet;
01222
01223 m[M40] = Det4_1234_0123 * oneOverDet;
01224 m[M41] = Det4_0234_0123 * mn1OverDet;
01225 m[M42] = Det4_0134_0123 * oneOverDet;
01226 m[M43] = Det4_0124_0123 * mn1OverDet;
01227 m[M44] = Det4_0123_0123 * oneOverDet;
01228
01229 return;
01230 }
01231
01232 void G4ErrorMatrix::invertHaywood6 (G4int & ifail)
01233 {
01234 ifail = 0;
01235
01236
01237
01238 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
01239 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
01240 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
01241 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
01242 G4double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40];
01243 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
01244 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
01245 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
01246 G4double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41];
01247 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
01248 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
01249 G4double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42];
01250 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
01251 G4double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43];
01252 G4double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44];
01253 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
01254 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
01255 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
01256 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
01257 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
01258 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
01259 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
01260 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
01261 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
01262 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
01263 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
01264 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
01265 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
01266 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
01267 G4double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54];
01268 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
01269 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
01270 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
01271 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
01272 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
01273 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
01274 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
01275 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
01276 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
01277 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
01278 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
01279 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
01280 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
01281 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
01282 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
01283
01284
01285
01286 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
01287 + m[A22]*Det2_34_01;
01288 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
01289 + m[A23]*Det2_34_01;
01290 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
01291 + m[A24]*Det2_34_01;
01292 G4double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
01293 + m[A25]*Det2_34_01;
01294 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
01295 + m[A23]*Det2_34_02;
01296 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
01297 + m[A24]*Det2_34_02;
01298 G4double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05
01299 + m[A25]*Det2_34_02;
01300 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
01301 + m[A24]*Det2_34_03;
01302 G4double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05
01303 + m[A25]*Det2_34_03;
01304 G4double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05
01305 + m[A25]*Det2_34_04;
01306 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
01307 + m[A23]*Det2_34_12;
01308 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
01309 + m[A24]*Det2_34_12;
01310 G4double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15
01311 + m[A25]*Det2_34_12;
01312 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
01313 + m[A24]*Det2_34_13;
01314 G4double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15
01315 + m[A25]*Det2_34_13;
01316 G4double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15
01317 + m[A25]*Det2_34_14;
01318 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
01319 + m[A24]*Det2_34_23;
01320 G4double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
01321 + m[A25]*Det2_34_23;
01322 G4double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
01323 + m[A25]*Det2_34_24;
01324 G4double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
01325 + m[A25]*Det2_34_34;
01326 G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
01327 + m[A22]*Det2_35_01;
01328 G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
01329 + m[A23]*Det2_35_01;
01330 G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
01331 + m[A24]*Det2_35_01;
01332 G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
01333 + m[A25]*Det2_35_01;
01334 G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
01335 + m[A23]*Det2_35_02;
01336 G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
01337 + m[A24]*Det2_35_02;
01338 G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
01339 + m[A25]*Det2_35_02;
01340 G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
01341 + m[A24]*Det2_35_03;
01342 G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
01343 + m[A25]*Det2_35_03;
01344 G4double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05
01345 + m[A25]*Det2_35_04;
01346 G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
01347 + m[A23]*Det2_35_12;
01348 G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
01349 + m[A24]*Det2_35_12;
01350 G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
01351 + m[A25]*Det2_35_12;
01352 G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
01353 + m[A24]*Det2_35_13;
01354 G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
01355 + m[A25]*Det2_35_13;
01356 G4double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15
01357 + m[A25]*Det2_35_14;
01358 G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
01359 + m[A24]*Det2_35_23;
01360 G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
01361 + m[A25]*Det2_35_23;
01362 G4double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25
01363 + m[A25]*Det2_35_24;
01364 G4double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35
01365 + m[A25]*Det2_35_34;
01366 G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
01367 + m[A22]*Det2_45_01;
01368 G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
01369 + m[A23]*Det2_45_01;
01370 G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
01371 + m[A24]*Det2_45_01;
01372 G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
01373 + m[A25]*Det2_45_01;
01374 G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
01375 + m[A23]*Det2_45_02;
01376 G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
01377 + m[A24]*Det2_45_02;
01378 G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
01379 + m[A25]*Det2_45_02;
01380 G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
01381 + m[A24]*Det2_45_03;
01382 G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
01383 + m[A25]*Det2_45_03;
01384 G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
01385 + m[A25]*Det2_45_04;
01386 G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
01387 + m[A23]*Det2_45_12;
01388 G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
01389 + m[A24]*Det2_45_12;
01390 G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
01391 + m[A25]*Det2_45_12;
01392 G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
01393 + m[A24]*Det2_45_13;
01394 G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
01395 + m[A25]*Det2_45_13;
01396 G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
01397 + m[A25]*Det2_45_14;
01398 G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
01399 + m[A24]*Det2_45_23;
01400 G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
01401 + m[A25]*Det2_45_23;
01402 G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
01403 + m[A25]*Det2_45_24;
01404 G4double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35
01405 + m[A25]*Det2_45_34;
01406 G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
01407 + m[A32]*Det2_45_01;
01408 G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
01409 + m[A33]*Det2_45_01;
01410 G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
01411 + m[A34]*Det2_45_01;
01412 G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
01413 + m[A35]*Det2_45_01;
01414 G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
01415 + m[A33]*Det2_45_02;
01416 G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
01417 + m[A34]*Det2_45_02;
01418 G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
01419 + m[A35]*Det2_45_02;
01420 G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
01421 + m[A34]*Det2_45_03;
01422 G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
01423 + m[A35]*Det2_45_03;
01424 G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
01425 + m[A35]*Det2_45_04;
01426 G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
01427 + m[A33]*Det2_45_12;
01428 G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
01429 + m[A34]*Det2_45_12;
01430 G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
01431 + m[A35]*Det2_45_12;
01432 G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
01433 + m[A34]*Det2_45_13;
01434 G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
01435 + m[A35]*Det2_45_13;
01436 G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
01437 + m[A35]*Det2_45_14;
01438 G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
01439 + m[A34]*Det2_45_23;
01440 G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
01441 + m[A35]*Det2_45_23;
01442 G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
01443 + m[A35]*Det2_45_24;
01444 G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
01445 + m[A35]*Det2_45_34;
01446
01447
01448
01449 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
01450 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
01451 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
01452 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
01453 G4double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025
01454 + m[A12]*Det3_234_015 - m[A15]*Det3_234_012;
01455 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
01456 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
01457 G4double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
01458 + m[A13]*Det3_234_015 - m[A15]*Det3_234_013;
01459 G4double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
01460 + m[A14]*Det3_234_015 - m[A15]*Det3_234_014;
01461 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
01462 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
01463 G4double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035
01464 + m[A13]*Det3_234_025 - m[A15]*Det3_234_023;
01465 G4double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045
01466 + m[A14]*Det3_234_025 - m[A15]*Det3_234_024;
01467 G4double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045
01468 + m[A14]*Det3_234_035 - m[A15]*Det3_234_034;
01469 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
01470 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
01471 G4double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135
01472 + m[A13]*Det3_234_125 - m[A15]*Det3_234_123;
01473 G4double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145
01474 + m[A14]*Det3_234_125 - m[A15]*Det3_234_124;
01475 G4double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145
01476 + m[A14]*Det3_234_135 - m[A15]*Det3_234_134;
01477 G4double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245
01478 + m[A14]*Det3_234_235 - m[A15]*Det3_234_234;
01479 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
01480 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
01481 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
01482 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
01483 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
01484 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
01485 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
01486 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
01487 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
01488 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
01489 G4double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045
01490 + m[A14]*Det3_235_015 - m[A15]*Det3_235_014;
01491 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
01492 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
01493 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
01494 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
01495 G4double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045
01496 + m[A14]*Det3_235_025 - m[A15]*Det3_235_024;
01497 G4double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045
01498 + m[A14]*Det3_235_035 - m[A15]*Det3_235_034;
01499 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
01500 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
01501 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
01502 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
01503 G4double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145
01504 + m[A14]*Det3_235_125 - m[A15]*Det3_235_124;
01505 G4double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145
01506 + m[A14]*Det3_235_135 - m[A15]*Det3_235_134;
01507 G4double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245
01508 + m[A14]*Det3_235_235 - m[A15]*Det3_235_234;
01509 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
01510 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
01511 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
01512 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
01513 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
01514 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
01515 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
01516 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
01517 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
01518 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
01519 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
01520 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
01521 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
01522 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
01523 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
01524 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
01525 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
01526 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
01527 G4double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045
01528 + m[A14]*Det3_245_035 - m[A15]*Det3_245_034;
01529 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
01530 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
01531 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
01532 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
01533 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
01534 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
01535 G4double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145
01536 + m[A14]*Det3_245_135 - m[A15]*Det3_245_134;
01537 G4double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245
01538 + m[A14]*Det3_245_235 - m[A15]*Det3_245_234;
01539 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
01540 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
01541 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
01542 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
01543 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
01544 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
01545 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
01546 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
01547 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
01548 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
01549 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
01550 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
01551 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
01552 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
01553 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
01554 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
01555 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
01556 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
01557 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
01558 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
01559 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
01560 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
01561 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
01562 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
01563 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
01564 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
01565 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
01566 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
01567 G4double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245
01568 + m[A14]*Det3_345_235 - m[A15]*Det3_345_234;
01569 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
01570 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
01571 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
01572 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
01573 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
01574 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
01575 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
01576 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
01577 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
01578 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
01579 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
01580 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
01581 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
01582 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
01583 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
01584 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
01585 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
01586 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
01587 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
01588 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
01589 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
01590 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
01591 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
01592 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
01593 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
01594 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
01595 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
01596 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
01597 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
01598 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
01599
01600
01601
01602 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
01603 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
01604 G4double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
01605 + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123;
01606 G4double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
01607 + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124;
01608 G4double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
01609 + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134;
01610 G4double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
01611 + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234;
01612 G4double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
01613 + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234;
01614 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
01615 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
01616 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
01617 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
01618 G4double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245
01619 + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124;
01620 G4double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345
01621 + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134;
01622 G4double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345
01623 + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234;
01624 G4double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345
01625 + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234;
01626 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
01627 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
01628 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
01629 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
01630 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
01631 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
01632 G4double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345
01633 + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134;
01634 G4double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345
01635 + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234;
01636 G4double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345
01637 + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234;
01638 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
01639 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
01640 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
01641 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
01642 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
01643 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
01644 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
01645 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
01646 G4double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345
01647 + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234;
01648 G4double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345
01649 + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234;
01650 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
01651 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
01652 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
01653 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
01654 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
01655 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
01656 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
01657 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
01658 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
01659 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
01660 G4double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345
01661 + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234;
01662 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
01663 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
01664 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
01665 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
01666 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
01667 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
01668 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
01669 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
01670 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
01671 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
01672 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
01673 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
01674
01675
01676
01677 G4double det = m[A00]*Det5_12345_12345
01678 - m[A01]*Det5_12345_02345
01679 + m[A02]*Det5_12345_01345
01680 - m[A03]*Det5_12345_01245
01681 + m[A04]*Det5_12345_01235
01682 - m[A05]*Det5_12345_01234;
01683
01684 if ( det == 0 )
01685 {
01686 ifail = 1;
01687 return;
01688 }
01689
01690 G4double oneOverDet = 1.0/det;
01691 G4double mn1OverDet = - oneOverDet;
01692
01693 m[A00] = Det5_12345_12345*oneOverDet;
01694 m[A01] = Det5_02345_12345*mn1OverDet;
01695 m[A02] = Det5_01345_12345*oneOverDet;
01696 m[A03] = Det5_01245_12345*mn1OverDet;
01697 m[A04] = Det5_01235_12345*oneOverDet;
01698 m[A05] = Det5_01234_12345*mn1OverDet;
01699
01700 m[A10] = Det5_12345_02345*mn1OverDet;
01701 m[A11] = Det5_02345_02345*oneOverDet;
01702 m[A12] = Det5_01345_02345*mn1OverDet;
01703 m[A13] = Det5_01245_02345*oneOverDet;
01704 m[A14] = Det5_01235_02345*mn1OverDet;
01705 m[A15] = Det5_01234_02345*oneOverDet;
01706
01707 m[A20] = Det5_12345_01345*oneOverDet;
01708 m[A21] = Det5_02345_01345*mn1OverDet;
01709 m[A22] = Det5_01345_01345*oneOverDet;
01710 m[A23] = Det5_01245_01345*mn1OverDet;
01711 m[A24] = Det5_01235_01345*oneOverDet;
01712 m[A25] = Det5_01234_01345*mn1OverDet;
01713
01714 m[A30] = Det5_12345_01245*mn1OverDet;
01715 m[A31] = Det5_02345_01245*oneOverDet;
01716 m[A32] = Det5_01345_01245*mn1OverDet;
01717 m[A33] = Det5_01245_01245*oneOverDet;
01718 m[A34] = Det5_01235_01245*mn1OverDet;
01719 m[A35] = Det5_01234_01245*oneOverDet;
01720
01721 m[A40] = Det5_12345_01235*oneOverDet;
01722 m[A41] = Det5_02345_01235*mn1OverDet;
01723 m[A42] = Det5_01345_01235*oneOverDet;
01724 m[A43] = Det5_01245_01235*mn1OverDet;
01725 m[A44] = Det5_01235_01235*oneOverDet;
01726 m[A45] = Det5_01234_01235*mn1OverDet;
01727
01728 m[A50] = Det5_12345_01234*mn1OverDet;
01729 m[A51] = Det5_02345_01234*oneOverDet;
01730 m[A52] = Det5_01345_01234*mn1OverDet;
01731 m[A53] = Det5_01245_01234*oneOverDet;
01732 m[A54] = Det5_01235_01234*mn1OverDet;
01733 m[A55] = Det5_01234_01234*oneOverDet;
01734
01735 return;
01736 }