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