ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ErrorMatrix.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ErrorMatrix.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // ------------------------------------------------------------
28 // GEANT 4 class implementation file
29 // ------------------------------------------------------------
30 
31 #include "globals.hh"
32 
33 #include <cmath>
34 #include <iostream>
35 
36 #include "G4ErrorMatrix.hh"
37 #include "G4ErrorSymMatrix.hh"
38 
39 // Simple operation for all elements
40 
41 #define SIMPLE_UOP(OPER) \
42  G4ErrorMatrixIter a=m.begin(); \
43  G4ErrorMatrixIter e=m.end(); \
44  for(;a!=e; a++) (*a) OPER t;
45 
46 #define SIMPLE_BOP(OPER) \
47  G4ErrorMatrixIter a=m.begin(); \
48  G4ErrorMatrixConstIter b=mat2.m.begin(); \
49  G4ErrorMatrixIter e=m.end(); \
50  for(;a!=e; a++, b++) (*a) OPER (*b);
51 
52 #define SIMPLE_TOP(OPER) \
53  G4ErrorMatrixConstIter a=mat1.m.begin(); \
54  G4ErrorMatrixConstIter b=mat2.m.begin(); \
55  G4ErrorMatrixIter t=mret.m.begin(); \
56  G4ErrorMatrixConstIter e=mat1.m.end(); \
57  for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
58 
59 // Static functions.
60 
61 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
62  if (r1!=r2 || c1!=c2) { \
63  G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
64  }
65 
66 #define CHK_DIM_1(c1,r2,fun) \
67  if (c1!=r2) { \
68  G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
69  }
70 
71 // Constructors. (Default constructors are inlined and in .icc file)
72 
73 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q)
74  : m(p*q), nrow(p), ncol(q)
75 {
76  size = nrow * ncol;
77 }
78 
79 G4ErrorMatrix::G4ErrorMatrix(G4int p,G4int q,G4int init)
80  : m(p*q), nrow(p), ncol(q)
81 {
82  size = nrow * ncol;
83 
84  if (size > 0)
85  {
86  switch(init)
87  {
88  case 0:
89  break;
90 
91  case 1:
92  {
93  if ( ncol == nrow )
94  {
95  G4ErrorMatrixIter a = m.begin();
96  G4ErrorMatrixIter b = m.end();
97  for( ; a<b; a+=(ncol+1)) *a = 1.0;
98  } else {
99  error("Invalid dimension in G4ErrorMatrix(G4int,G4int,1).");
100  }
101  break;
102  }
103  default:
104  error("G4ErrorMatrix: initialization must be either 0 or 1.");
105  }
106  }
107 }
108 
109 //
110 // Destructor
111 //
112 G4ErrorMatrix::~G4ErrorMatrix()
113 {
114 }
115 
116 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorMatrix &mat1)
117  : m(mat1.size), nrow(mat1.nrow), ncol(mat1.ncol), size(mat1.size)
118 {
119  m = mat1.m;
120 }
121 
122 G4ErrorMatrix::G4ErrorMatrix(const G4ErrorSymMatrix &mat1)
123  : m(mat1.nrow*mat1.nrow), nrow(mat1.nrow), ncol(mat1.nrow)
124 {
125  size = nrow * ncol;
126 
127  G4int n = ncol;
128  G4ErrorMatrixConstIter sjk = mat1.m.begin();
129  G4ErrorMatrixIter m1j = m.begin();
130  G4ErrorMatrixIter mj = m.begin();
131  // j >= k
132  for(G4int j=1;j<=nrow;j++)
133  {
134  G4ErrorMatrixIter mjk = mj;
135  G4ErrorMatrixIter mkj = m1j;
136  for(G4int k=1;k<=j;k++)
137  {
138  *(mjk++) = *sjk;
139  if(j!=k) *mkj = *sjk;
140  sjk++;
141  mkj += n;
142  }
143  mj += n;
144  m1j++;
145  }
146 }
147 
148 //
149 //
150 // Sub matrix
151 //
152 //
153 
154 G4ErrorMatrix G4ErrorMatrix::sub(G4int min_row, G4int max_row,
155  G4int min_col,G4int max_col) const
156 {
157  G4ErrorMatrix mret(max_row-min_row+1,max_col-min_col+1);
158  if(max_row > num_row() || max_col >num_col())
159  { error("G4ErrorMatrix::sub: Index out of range"); }
160  G4ErrorMatrixIter a = mret.m.begin();
161  G4int nc = num_col();
162  G4ErrorMatrixConstIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
163 
164  for(G4int irow=1; irow<=mret.num_row(); irow++)
165  {
166  G4ErrorMatrixConstIter brc = b1;
167  for(G4int icol=1; icol<=mret.num_col(); icol++)
168  {
169  *(a++) = *(brc++);
170  }
171  b1 += nc;
172  }
173  return mret;
174 }
175 
176 void G4ErrorMatrix::sub(G4int row,G4int col,const G4ErrorMatrix &mat1)
177 {
178  if(row <1 || row+mat1.num_row()-1 > num_row() ||
179  col <1 || col+mat1.num_col()-1 > num_col() )
180  { error("G4ErrorMatrix::sub: Index out of range"); }
181  G4ErrorMatrixConstIter a = mat1.m.begin();
182  G4int nc = num_col();
183  G4ErrorMatrixIter b1 = m.begin() + (row - 1) * nc + col - 1;
184 
185  for(G4int irow=1; irow<=mat1.num_row(); irow++)
186  {
187  G4ErrorMatrixIter brc = b1;
188  for(G4int icol=1; icol<=mat1.num_col(); icol++)
189  {
190  *(brc++) = *(a++);
191  }
192  b1 += nc;
193  }
194 }
195 
196 //
197 // Direct sum of two matricies
198 //
199 
200 G4ErrorMatrix dsum(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
201 {
202  G4ErrorMatrix mret(mat1.num_row() + mat2.num_row(),
203  mat1.num_col() + mat2.num_col(), 0);
204  mret.sub(1,1,mat1);
205  mret.sub(mat1.num_row()+1,mat1.num_col()+1,mat2);
206  return mret;
207 }
208 
209 /* -----------------------------------------------------------------------
210  This section contains support routines for matrix.h. This section contains
211  The two argument functions +,-. They call the copy constructor and +=,-=.
212  ----------------------------------------------------------------------- */
213 
214 G4ErrorMatrix G4ErrorMatrix::operator- () const
215 {
216  G4ErrorMatrix mat2(nrow, ncol);
217  G4ErrorMatrixConstIter a=m.begin();
218  G4ErrorMatrixIter b=mat2.m.begin();
219  G4ErrorMatrixConstIter e=m.end();
220  for(;a<e; a++, b++) (*b) = -(*a);
221  return mat2;
222 }
223 
224 G4ErrorMatrix operator+(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
225 {
226  G4ErrorMatrix mret(mat1.nrow, mat1.ncol);
227  CHK_DIM_2(mat1.num_row(),mat2.num_row(), mat1.num_col(),mat2.num_col(),+);
228  SIMPLE_TOP(+)
229  return mret;
230 }
231 
232 //
233 // operator -
234 //
235 
236 G4ErrorMatrix operator-(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
237 {
238  G4ErrorMatrix mret(mat1.num_row(), mat1.num_col());
239  CHK_DIM_2(mat1.num_row(),mat2.num_row(),
240  mat1.num_col(),mat2.num_col(),-);
241  SIMPLE_TOP(-)
242  return mret;
243 }
244 
245 /* -----------------------------------------------------------------------
246  This section contains support routines for matrix.h. This file contains
247  The two argument functions *,/. They call copy constructor and then /=,*=.
248  ----------------------------------------------------------------------- */
249 
250 G4ErrorMatrix operator/(const G4ErrorMatrix &mat1,G4double t)
251 {
252  G4ErrorMatrix mret(mat1);
253  mret /= t;
254  return mret;
255 }
256 
257 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,G4double t)
258 {
259  G4ErrorMatrix mret(mat1);
260  mret *= t;
261  return mret;
262 }
263 
264 G4ErrorMatrix operator*(G4double t,const G4ErrorMatrix &mat1)
265 {
266  G4ErrorMatrix mret(mat1);
267  mret *= t;
268  return mret;
269 }
270 
271 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1,const G4ErrorMatrix &mat2)
272 {
273  // initialize matrix to 0.0
274  G4ErrorMatrix mret(mat1.nrow,mat2.ncol,0);
275  CHK_DIM_1(mat1.ncol,mat2.nrow,*);
276 
277  G4int m1cols = mat1.ncol;
278  G4int m2cols = mat2.ncol;
279 
280  for (G4int i=0; i<mat1.nrow; i++)
281  {
282  for (G4int j=0; j<m1cols; j++)
283  {
284  G4double temp = mat1.m[i*m1cols+j];
285  G4ErrorMatrixIter pt = mret.m.begin() + i*m2cols;
286 
287  // Loop over k (the column index in matrix mat2)
288  G4ErrorMatrixConstIter pb = mat2.m.begin() + m2cols*j;
289  const G4ErrorMatrixConstIter pblast = pb + m2cols;
290  while (pb < pblast) // Loop checking, 06.08.2015, G.Cosmo
291  {
292  (*pt) += temp * (*pb);
293  pb++;
294  pt++;
295  }
296  }
297  }
298  return mret;
299 }
300 
301 /* -----------------------------------------------------------------------
302  This section contains the assignment and inplace operators =,+=,-=,*=,/=.
303  ----------------------------------------------------------------------- */
304 
305 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorMatrix &mat2)
306 {
307  CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
308  SIMPLE_BOP(+=)
309  return (*this);
310 }
311 
312 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorMatrix &mat2)
313 {
314  CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
315  SIMPLE_BOP(-=)
316  return (*this);
317 }
318 
319 G4ErrorMatrix & G4ErrorMatrix::operator/=(G4double t)
320 {
321  SIMPLE_UOP(/=)
322  return (*this);
323 }
324 
325 G4ErrorMatrix & G4ErrorMatrix::operator*=(G4double t)
326 {
327  SIMPLE_UOP(*=)
328  return (*this);
329 }
330 
331 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorMatrix &mat1)
332 {
333  if (&mat1 == this) { return *this; }
334 
335  if(mat1.nrow*mat1.ncol != size) //??fixme?? mat1.size != size
336  {
337  size = mat1.nrow * mat1.ncol;
338  m.resize(size); //??fixme?? if (size < mat1.size) m.resize(mat1.size);
339  }
340  nrow = mat1.nrow;
341  ncol = mat1.ncol;
342  m = mat1.m;
343  return (*this);
344 }
345 
346 
347 // Print the G4ErrorMatrix.
348 
349 std::ostream& operator<<(std::ostream &os, const G4ErrorMatrix &q)
350 {
351  os << "\n";
352 
353  // Fixed format needs 3 extra characters for field,
354  // while scientific needs 7
355 
356  G4int width;
357  if(os.flags() & std::ios::fixed)
358  { width = os.precision()+3; }
359  else
360  { width = os.precision()+7; }
361  for(G4int irow = 1; irow<= q.num_row(); irow++)
362  {
363  for(G4int icol = 1; icol <= q.num_col(); icol++)
364  {
365  os.width(width);
366  os << q(irow,icol) << " ";
367  }
368  os << G4endl;
369  }
370  return os;
371 }
372 
373 G4ErrorMatrix G4ErrorMatrix::T() const
374 {
375  G4ErrorMatrix mret(ncol,nrow);
376  G4ErrorMatrixConstIter pl = m.end();
377  G4ErrorMatrixConstIter pme = m.begin();
378  G4ErrorMatrixIter pt = mret.m.begin();
379  G4ErrorMatrixIter ptl = mret.m.end();
380  for (; pme < pl; pme++, pt+=nrow)
381  {
382  if (pt >= ptl) { pt -= (size-1); }
383  (*pt) = (*pme);
384  }
385  return mret;
386 }
387 
388 G4ErrorMatrix
389 G4ErrorMatrix::apply(G4double (*f)(G4double, G4int, G4int)) const
390 {
391  G4ErrorMatrix mret(num_row(),num_col());
392  G4ErrorMatrixConstIter a = m.begin();
393  G4ErrorMatrixIter b = mret.m.begin();
394  for(G4int ir=1;ir<=num_row();ir++)
395  {
396  for(G4int ic=1;ic<=num_col();ic++)
397  {
398  *(b++) = (*f)(*(a++), ir, ic);
399  }
400  }
401  return mret;
402 }
403 
404 G4int G4ErrorMatrix::dfinv_matrix(G4int *ir) {
405  if (num_col()!=num_row())
406  { error("dfinv_matrix: G4ErrorMatrix is not NxN"); }
407  G4int n = num_col();
408  if (n==1) { return 0; }
409 
410  G4double s31, s32;
411  G4double s33, s34;
412 
413  G4ErrorMatrixIter m11 = m.begin();
414  G4ErrorMatrixIter m12 = m11 + 1;
415  G4ErrorMatrixIter m21 = m11 + n;
416  G4ErrorMatrixIter m22 = m12 + n;
417  *m21 = -(*m22) * (*m11) * (*m21);
418  *m12 = -(*m12);
419  if (n>2)
420  {
421  G4ErrorMatrixIter mi = m.begin() + 2 * n;
422  G4ErrorMatrixIter mii= m.begin() + 2 * n + 2;
423  G4ErrorMatrixIter mimim = m.begin() + n + 1;
424  for (G4int i=3;i<=n;i++)
425  {
426  G4int im2 = i - 2;
427  G4ErrorMatrixIter mj = m.begin();
428  G4ErrorMatrixIter mji = mj + i - 1;
429  G4ErrorMatrixIter mij = mi;
430  for (G4int j=1;j<=im2;j++)
431  {
432  s31 = 0.0;
433  s32 = *mji;
434  G4ErrorMatrixIter mkj = mj + j - 1;
435  G4ErrorMatrixIter mik = mi + j - 1;
436  G4ErrorMatrixIter mjkp = mj + j;
437  G4ErrorMatrixIter mkpi = mj + n + i - 1;
438  for (G4int k=j;k<=im2;k++)
439  {
440  s31 += (*mkj) * (*(mik++));
441  s32 += (*(mjkp++)) * (*mkpi);
442  mkj += n;
443  mkpi += n;
444  }
445  *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
446  *mji = -s32;
447  mj += n;
448  mji += n;
449  mij++;
450  }
451  *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
452  *(mimim+1) = -(*(mimim+1));
453  mi += n;
454  mimim += (n+1);
455  mii += (n+1);
456  }
457  }
458  G4ErrorMatrixIter mi = m.begin();
459  G4ErrorMatrixIter mii = m.begin();
460  for (G4int i=1;i<n;i++)
461  {
462  G4int ni = n - i;
463  G4ErrorMatrixIter mij = mi;
464  G4int j;
465  for (j=1; j<=i;j++)
466  {
467  s33 = *mij;
468  G4ErrorMatrixIter mikj = mi + n + j - 1;
469  G4ErrorMatrixIter miik = mii + 1;
470  G4ErrorMatrixIter min_end = mi + n;
471  for (;miik<min_end;)
472  {
473  s33 += (*mikj) * (*(miik++));
474  mikj += n;
475  }
476  *(mij++) = s33;
477  }
478  for (j=1;j<=ni;j++)
479  {
480  s34 = 0.0;
481  G4ErrorMatrixIter miik = mii + j;
482  G4ErrorMatrixIter mikij = mii + j * n + j;
483  for (G4int k=j;k<=ni;k++)
484  {
485  s34 += *mikij * (*(miik++));
486  mikij += n;
487  }
488  *(mii+j) = s34;
489  }
490  mi += n;
491  mii += (n+1);
492  }
493  G4int nxch = ir[n];
494  if (nxch==0) return 0;
495  for (G4int mq=1;mq<=nxch;mq++)
496  {
497  G4int k = nxch - mq + 1;
498  G4int ij = ir[k];
499  G4int i = ij >> 12;
500  G4int j = ij%4096;
501  G4ErrorMatrixIter mki = m.begin() + i - 1;
502  G4ErrorMatrixIter mkj = m.begin() + j - 1;
503  for (k=1; k<=n;k++)
504  {
505  // 2/24/05 David Sachs fix of improper swap bug that was present
506  // for many years:
507  G4double ti = *mki; // 2/24/05
508  *mki = *mkj;
509  *mkj = ti; // 2/24/05
510  mki += n;
511  mkj += n;
512  }
513  }
514  return 0;
515 }
516 
517 G4int G4ErrorMatrix::dfact_matrix(G4double &det, G4int *ir)
518 {
519  if (ncol!=nrow)
520  error("dfact_matrix: G4ErrorMatrix is not NxN");
521 
522  G4int ifail, jfail;
523  G4int n = ncol;
524 
525  G4double tf;
526  G4double g1 = 1.0e-19, g2 = 1.0e19;
527 
528  G4double p, q, t;
529  G4double s11, s12;
530 
531  G4double epsilon = 8*DBL_EPSILON;
532  // could be set to zero (like it was before)
533  // but then the algorithm often doesn't detect
534  // that a matrix is singular
535 
536  G4int normal = 0, imposs = -1;
537  G4int jrange = 0, jover = 1, junder = -1;
538  ifail = normal;
539  jfail = jrange;
540  G4int nxch = 0;
541  det = 1.0;
542  G4ErrorMatrixIter mj = m.begin();
543  G4ErrorMatrixIter mjj = mj;
544  for (G4int j=1;j<=n;j++)
545  {
546  G4int k = j;
547  p = (std::fabs(*mjj));
548  if (j!=n) {
549  G4ErrorMatrixIter mij = mj + n + j - 1;
550  for (G4int i=j+1;i<=n;i++)
551  {
552  q = (std::fabs(*(mij)));
553  if (q > p)
554  {
555  k = i;
556  p = q;
557  }
558  mij += n;
559  }
560  if (k==j)
561  {
562  if (p <= epsilon)
563  {
564  det = 0;
565  ifail = imposs;
566  jfail = jrange;
567  return ifail;
568  }
569  det = -det; // in this case the sign of the determinant
570  // must not change. So I change it twice.
571  }
572  G4ErrorMatrixIter mjl = mj;
573  G4ErrorMatrixIter mkl = m.begin() + (k-1)*n;
574  for (G4int l=1;l<=n;l++)
575  {
576  tf = *mjl;
577  *(mjl++) = *mkl;
578  *(mkl++) = tf;
579  }
580  nxch = nxch + 1; // this makes the determinant change its sign
581  ir[nxch] = (((j)<<12)+(k));
582  }
583  else
584  {
585  if (p <= epsilon)
586  {
587  det = 0.0;
588  ifail = imposs;
589  jfail = jrange;
590  return ifail;
591  }
592  }
593  det *= *mjj;
594  *mjj = 1.0 / *mjj;
595  t = (std::fabs(det));
596  if (t < g1)
597  {
598  det = 0.0;
599  if (jfail == jrange) jfail = junder;
600  }
601  else if (t > g2)
602  {
603  det = 1.0;
604  if (jfail==jrange) jfail = jover;
605  }
606  if (j!=n)
607  {
608  G4ErrorMatrixIter mk = mj + n;
609  G4ErrorMatrixIter mkjp = mk + j;
610  G4ErrorMatrixIter mjk = mj + j;
611  for (k=j+1;k<=n;k++)
612  {
613  s11 = - (*mjk);
614  s12 = - (*mkjp);
615  if (j!=1)
616  {
617  G4ErrorMatrixIter mik = m.begin() + k - 1;
618  G4ErrorMatrixIter mijp = m.begin() + j;
619  G4ErrorMatrixIter mki = mk;
620  G4ErrorMatrixIter mji = mj;
621  for (G4int i=1;i<j;i++)
622  {
623  s11 += (*mik) * (*(mji++));
624  s12 += (*mijp) * (*(mki++));
625  mik += n;
626  mijp += n;
627  }
628  }
629  *(mjk++) = -s11 * (*mjj);
630  *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
631  mk += n;
632  mkjp += n;
633  }
634  }
635  mj += n;
636  mjj += (n+1);
637  }
638  if (nxch%2==1) det = -det;
639  if (jfail !=jrange) det = 0.0;
640  ir[n] = nxch;
641  return 0;
642 }
643 
644 void G4ErrorMatrix::invert(G4int &ierr)
645 {
646  if(ncol != nrow)
647  { error("G4ErrorMatrix::invert: G4ErrorMatrix is not NxN"); }
648 
649  static G4ThreadLocal G4int max_array = 20;
650  static G4ThreadLocal G4int *ir = 0 ; if (!ir) ir= new G4int [max_array+1];
651 
652  if (ncol > max_array)
653  {
654  delete [] ir;
655  max_array = nrow;
656  ir = new G4int [max_array+1];
657  }
658  G4double t1, t2, t3;
659  G4double det, temp, ss;
660  G4int ifail;
661  switch(nrow)
662  {
663  case 3:
664  G4double c11,c12,c13,c21,c22,c23,c31,c32,c33;
665  ifail = 0;
666  c11 = (*(m.begin()+4)) * (*(m.begin()+8))
667  - (*(m.begin()+5)) * (*(m.begin()+7));
668  c12 = (*(m.begin()+5)) * (*(m.begin()+6))
669  - (*(m.begin()+3)) * (*(m.begin()+8));
670  c13 = (*(m.begin()+3)) * (*(m.begin()+7))
671  - (*(m.begin()+4)) * (*(m.begin()+6));
672  c21 = (*(m.begin()+7)) * (*(m.begin()+2))
673  - (*(m.begin()+8)) * (*(m.begin()+1));
674  c22 = (*(m.begin()+8)) * (*m.begin())
675  - (*(m.begin()+6)) * (*(m.begin()+2));
676  c23 = (*(m.begin()+6)) * (*(m.begin()+1))
677  - (*(m.begin()+7)) * (*m.begin());
678  c31 = (*(m.begin()+1)) * (*(m.begin()+5))
679  - (*(m.begin()+2)) * (*(m.begin()+4));
680  c32 = (*(m.begin()+2)) * (*(m.begin()+3))
681  - (*m.begin()) * (*(m.begin()+5));
682  c33 = (*m.begin()) * (*(m.begin()+4))
683  - (*(m.begin()+1)) * (*(m.begin()+3));
684  t1 = std::fabs(*m.begin());
685  t2 = std::fabs(*(m.begin()+3));
686  t3 = std::fabs(*(m.begin()+6));
687  if (t1 >= t2)
688  {
689  if (t3 >= t1)
690  {
691  temp = *(m.begin()+6);
692  det = c23*c12-c22*c13;
693  }
694  else
695  {
696  temp = *(m.begin());
697  det = c22*c33-c23*c32;
698  }
699  }
700  else if (t3 >= t2)
701  {
702  temp = *(m.begin()+6);
703  det = c23*c12-c22*c13;
704  }
705  else
706  {
707  temp = *(m.begin()+3);
708  det = c13*c32-c12*c33;
709  }
710  if (det==0)
711  {
712  ierr = 1;
713  return;
714  }
715  {
716  ss = temp/det;
717  G4ErrorMatrixIter mq = m.begin();
718  *(mq++) = ss*c11;
719  *(mq++) = ss*c21;
720  *(mq++) = ss*c31;
721  *(mq++) = ss*c12;
722  *(mq++) = ss*c22;
723  *(mq++) = ss*c32;
724  *(mq++) = ss*c13;
725  *(mq++) = ss*c23;
726  *(mq) = ss*c33;
727  }
728  break;
729  case 2:
730  ifail = 0;
731  det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
732  if (det==0)
733  {
734  ierr = 1;
735  return;
736  }
737  ss = 1.0/det;
738  temp = ss*(*(m.begin()+3));
739  *(m.begin()+1) *= -ss;
740  *(m.begin()+2) *= -ss;
741  *(m.begin()+3) = ss*(*m.begin());
742  *(m.begin()) = temp;
743  break;
744  case 1:
745  ifail = 0;
746  if ((*(m.begin()))==0)
747  {
748  ierr = 1;
749  return;
750  }
751  *(m.begin()) = 1.0/(*(m.begin()));
752  break;
753  case 4:
754  invertHaywood4(ierr);
755  return;
756  case 5:
757  invertHaywood5(ierr);
758  return;
759  case 6:
760  invertHaywood6(ierr);
761  return;
762  default:
763  ifail = dfact_matrix(det, ir);
764  if(ifail) {
765  ierr = 1;
766  return;
767  }
768  dfinv_matrix(ir);
769  break;
770  }
771  ierr = 0;
772  return;
773 }
774 
775 G4double G4ErrorMatrix::determinant() const
776 {
777  static G4ThreadLocal G4int max_array = 20;
778  static G4ThreadLocal G4int *ir = 0 ; if (!ir) ir= new G4int [max_array+1];
779  if(ncol != nrow)
780  { error("G4ErrorMatrix::determinant: G4ErrorMatrix is not NxN"); }
781  if (ncol > max_array)
782  {
783  delete [] ir;
784  max_array = nrow;
785  ir = new G4int [max_array+1];
786  }
787  G4double det;
788  G4ErrorMatrix mt(*this);
789  G4int i = mt.dfact_matrix(det, ir);
790  if(i==0) return det;
791  return 0;
792 }
793 
794 G4double G4ErrorMatrix::trace() const
795 {
796  G4double t = 0.0;
797  for (G4ErrorMatrixConstIter d = m.begin(); d < m.end(); d += (ncol+1) )
798  {
799  t += *d;
800  }
801  return t;
802 }
803 
804 void G4ErrorMatrix::error(const char *msg)
805 {
806  std::ostringstream message;
807  message << msg;
808  G4Exception("G4ErrorMatrix::error()", "GEANT4e-Error",
809  FatalException, message, "Exiting to System.");
810 }
811 
812 
813 // Aij are indices for a 6x6 matrix.
814 // Mij are indices for a 5x5 matrix.
815 // Fij are indices for a 4x4 matrix.
816 
817 #define A00 0
818 #define A01 1
819 #define A02 2
820 #define A03 3
821 #define A04 4
822 #define A05 5
823 
824 #define A10 6
825 #define A11 7
826 #define A12 8
827 #define A13 9
828 #define A14 10
829 #define A15 11
830 
831 #define A20 12
832 #define A21 13
833 #define A22 14
834 #define A23 15
835 #define A24 16
836 #define A25 17
837 
838 #define A30 18
839 #define A31 19
840 #define A32 20
841 #define A33 21
842 #define A34 22
843 #define A35 23
844 
845 #define A40 24
846 #define A41 25
847 #define A42 26
848 #define A43 27
849 #define A44 28
850 #define A45 29
851 
852 #define A50 30
853 #define A51 31
854 #define A52 32
855 #define A53 33
856 #define A54 34
857 #define A55 35
858 
859 #define M00 0
860 #define M01 1
861 #define M02 2
862 #define M03 3
863 #define M04 4
864 
865 #define M10 5
866 #define M11 6
867 #define M12 7
868 #define M13 8
869 #define M14 9
870 
871 #define M20 10
872 #define M21 11
873 #define M22 12
874 #define M23 13
875 #define M24 14
876 
877 #define M30 15
878 #define M31 16
879 #define M32 17
880 #define M33 18
881 #define M34 19
882 
883 #define M40 20
884 #define M41 21
885 #define M42 22
886 #define M43 23
887 #define M44 24
888 
889 #define F00 0
890 #define F01 1
891 #define F02 2
892 #define F03 3
893 
894 #define F10 4
895 #define F11 5
896 #define F12 6
897 #define F13 7
898 
899 #define F20 8
900 #define F21 9
901 #define F22 10
902 #define F23 11
903 
904 #define F30 12
905 #define F31 13
906 #define F32 14
907 #define F33 15
908 
909 
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 }