Geant4  10.02.p02
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 91809 2015-08-06 13:00:09Z 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) // Loop checking, 06.08.2015, G.Cosmo
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 G4ThreadLocal G4int max_array = 20;
651  static G4ThreadLocal G4int *ir = 0 ; if (!ir) 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 G4ThreadLocal G4int max_array = 20;
779  static G4ThreadLocal G4int *ir = 0 ; if (!ir) 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  std::ostringstream message;
808  message << msg;
809  G4Exception("G4ErrorMatrix::error()", "GEANT4e-Error",
810  FatalException, message, "Exiting to System.");
811 }
812 
813 
814 // Aij are indices for a 6x6 matrix.
815 // Mij are indices for a 5x5 matrix.
816 // Fij are indices for a 4x4 matrix.
817 
818 #define A00 0
819 #define A01 1
820 #define A02 2
821 #define A03 3
822 #define A04 4
823 #define A05 5
824 
825 #define A10 6
826 #define A11 7
827 #define A12 8
828 #define A13 9
829 #define A14 10
830 #define A15 11
831 
832 #define A20 12
833 #define A21 13
834 #define A22 14
835 #define A23 15
836 #define A24 16
837 #define A25 17
838 
839 #define A30 18
840 #define A31 19
841 #define A32 20
842 #define A33 21
843 #define A34 22
844 #define A35 23
845 
846 #define A40 24
847 #define A41 25
848 #define A42 26
849 #define A43 27
850 #define A44 28
851 #define A45 29
852 
853 #define A50 30
854 #define A51 31
855 #define A52 32
856 #define A53 33
857 #define A54 34
858 #define A55 35
859 
860 #define M00 0
861 #define M01 1
862 #define M02 2
863 #define M03 3
864 #define M04 4
865 
866 #define M10 5
867 #define M11 6
868 #define M12 7
869 #define M13 8
870 #define M14 9
871 
872 #define M20 10
873 #define M21 11
874 #define M22 12
875 #define M23 13
876 #define M24 14
877 
878 #define M30 15
879 #define M31 16
880 #define M32 17
881 #define M33 18
882 #define M34 19
883 
884 #define M40 20
885 #define M41 21
886 #define M42 22
887 #define M43 23
888 #define M44 24
889 
890 #define F00 0
891 #define F01 1
892 #define F02 2
893 #define F03 3
894 
895 #define F10 4
896 #define F11 5
897 #define F12 6
898 #define F13 7
899 
900 #define F20 8
901 #define F21 9
902 #define F22 10
903 #define F23 11
904 
905 #define F30 12
906 #define F31 13
907 #define F32 14
908 #define F33 15
909 
910 
911 void G4ErrorMatrix::invertHaywood4 (G4int & ifail)
912 {
913  ifail = 0;
914 
915  // Find all NECESSARY 2x2 dets: (18 of them)
916 
917  G4double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
918  G4double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
919  G4double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20];
920  G4double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21];
921  G4double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22];
922  G4double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
923  G4double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
924  G4double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
925  G4double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
926  G4double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
927  G4double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
928  G4double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32];
929  G4double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
930  G4double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
931  G4double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
932  G4double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
933  G4double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
934  G4double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
935 
936  // Find all NECESSARY 3x3 dets: (16 of them)
937 
938  G4double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02
939  + m[F02]*Det2_12_01;
940  G4double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03
941  + m[F03]*Det2_12_01;
942  G4double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03
943  + m[F03]*Det2_12_02;
944  G4double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13
945  + m[F03]*Det2_12_12;
946  G4double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02
947  + m[F02]*Det2_13_01;
948  G4double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
949  + m[F03]*Det2_13_01;
950  G4double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
951  + m[F03]*Det2_13_02;
952  G4double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
953  + m[F03]*Det2_13_12;
954  G4double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02
955  + m[F02]*Det2_23_01;
956  G4double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
957  + m[F03]*Det2_23_01;
958  G4double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
959  + m[F03]*Det2_23_02;
960  G4double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
961  + m[F03]*Det2_23_12;
962  G4double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02
963  + m[F12]*Det2_23_01;
964  G4double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03
965  + m[F13]*Det2_23_01;
966  G4double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03
967  + m[F13]*Det2_23_02;
968  G4double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13
969  + m[F13]*Det2_23_12;
970 
971  // Find the 4x4 det:
972 
973  G4double det = m[F00]*Det3_123_123
974  - m[F01]*Det3_123_023
975  + m[F02]*Det3_123_013
976  - m[F03]*Det3_123_012;
977 
978  if ( det == 0 )
979  {
980  ifail = 1;
981  return;
982  }
983 
984  G4double oneOverDet = 1.0/det;
985  G4double mn1OverDet = - oneOverDet;
986 
987  m[F00] = Det3_123_123 * oneOverDet;
988  m[F01] = Det3_023_123 * mn1OverDet;
989  m[F02] = Det3_013_123 * oneOverDet;
990  m[F03] = Det3_012_123 * mn1OverDet;
991 
992  m[F10] = Det3_123_023 * mn1OverDet;
993  m[F11] = Det3_023_023 * oneOverDet;
994  m[F12] = Det3_013_023 * mn1OverDet;
995  m[F13] = Det3_012_023 * oneOverDet;
996 
997  m[F20] = Det3_123_013 * oneOverDet;
998  m[F21] = Det3_023_013 * mn1OverDet;
999  m[F22] = Det3_013_013 * oneOverDet;
1000  m[F23] = Det3_012_013 * mn1OverDet;
1001 
1002  m[F30] = Det3_123_012 * mn1OverDet;
1003  m[F31] = Det3_023_012 * oneOverDet;
1004  m[F32] = Det3_013_012 * mn1OverDet;
1005  m[F33] = Det3_012_012 * oneOverDet;
1006 
1007  return;
1008 }
1009 
1010 void G4ErrorMatrix::invertHaywood5 (G4int & ifail)
1011 {
1012  ifail = 0;
1013 
1014  // Find all NECESSARY 2x2 dets: (30 of them)
1015 
1016  G4double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
1017  G4double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
1018  G4double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
1019  G4double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
1020  G4double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31];
1021  G4double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
1022  G4double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31];
1023  G4double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
1024  G4double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32];
1025  G4double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33];
1026  G4double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
1027  G4double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
1028  G4double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
1029  G4double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
1030  G4double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
1031  G4double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
1032  G4double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
1033  G4double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
1034  G4double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
1035  G4double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43];
1036  G4double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
1037  G4double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
1038  G4double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
1039  G4double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
1040  G4double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
1041  G4double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
1042  G4double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
1043  G4double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
1044  G4double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
1045  G4double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
1046 
1047  // Find all NECESSARY 3x3 dets: (40 of them)
1048 
1049  G4double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02
1050  + m[M12]*Det2_23_01;
1051  G4double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03
1052  + m[M13]*Det2_23_01;
1053  G4double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04
1054  + m[M14]*Det2_23_01;
1055  G4double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03
1056  + m[M13]*Det2_23_02;
1057  G4double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04
1058  + m[M14]*Det2_23_02;
1059  G4double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04
1060  + m[M14]*Det2_23_03;
1061  G4double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13
1062  + m[M13]*Det2_23_12;
1063  G4double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14
1064  + m[M14]*Det2_23_12;
1065  G4double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14
1066  + m[M14]*Det2_23_13;
1067  G4double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24
1068  + m[M14]*Det2_23_23;
1069  G4double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02
1070  + m[M12]*Det2_24_01;
1071  G4double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03
1072  + m[M13]*Det2_24_01;
1073  G4double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04
1074  + m[M14]*Det2_24_01;
1075  G4double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03
1076  + m[M13]*Det2_24_02;
1077  G4double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04
1078  + m[M14]*Det2_24_02;
1079  G4double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04
1080  + m[M14]*Det2_24_03;
1081  G4double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13
1082  + m[M13]*Det2_24_12;
1083  G4double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14
1084  + m[M14]*Det2_24_12;
1085  G4double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14
1086  + m[M14]*Det2_24_13;
1087  G4double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24
1088  + m[M14]*Det2_24_23;
1089  G4double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02
1090  + m[M12]*Det2_34_01;
1091  G4double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03
1092  + m[M13]*Det2_34_01;
1093  G4double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04
1094  + m[M14]*Det2_34_01;
1095  G4double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03
1096  + m[M13]*Det2_34_02;
1097  G4double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04
1098  + m[M14]*Det2_34_02;
1099  G4double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04
1100  + m[M14]*Det2_34_03;
1101  G4double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13
1102  + m[M13]*Det2_34_12;
1103  G4double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14
1104  + m[M14]*Det2_34_12;
1105  G4double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14
1106  + m[M14]*Det2_34_13;
1107  G4double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24
1108  + m[M14]*Det2_34_23;
1109  G4double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02
1110  + m[M22]*Det2_34_01;
1111  G4double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03
1112  + m[M23]*Det2_34_01;
1113  G4double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04
1114  + m[M24]*Det2_34_01;
1115  G4double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03
1116  + m[M23]*Det2_34_02;
1117  G4double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04
1118  + m[M24]*Det2_34_02;
1119  G4double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04
1120  + m[M24]*Det2_34_03;
1121  G4double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13
1122  + m[M23]*Det2_34_12;
1123  G4double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14
1124  + m[M24]*Det2_34_12;
1125  G4double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14
1126  + m[M24]*Det2_34_13;
1127  G4double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24
1128  + m[M24]*Det2_34_23;
1129 
1130  // Find all NECESSARY 4x4 dets: (25 of them)
1131 
1132  G4double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023
1133  + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
1134  G4double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024
1135  + m[M02]*Det3_123_014 - m[M04]*Det3_123_012;
1136  G4double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034
1137  + m[M03]*Det3_123_014 - m[M04]*Det3_123_013;
1138  G4double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034
1139  + m[M03]*Det3_123_024 - m[M04]*Det3_123_023;
1140  G4double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134
1141  + m[M03]*Det3_123_124 - m[M04]*Det3_123_123;
1142  G4double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023
1143  + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
1144  G4double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024
1145  + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
1146  G4double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034
1147  + m[M03]*Det3_124_014 - m[M04]*Det3_124_013;
1148  G4double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034
1149  + m[M03]*Det3_124_024 - m[M04]*Det3_124_023;
1150  G4double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134
1151  + m[M03]*Det3_124_124 - m[M04]*Det3_124_123;
1152  G4double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023
1153  + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
1154  G4double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024
1155  + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
1156  G4double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034
1157  + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
1158  G4double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034
1159  + m[M03]*Det3_134_024 - m[M04]*Det3_134_023;
1160  G4double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134
1161  + m[M03]*Det3_134_124 - m[M04]*Det3_134_123;
1162  G4double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023
1163  + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
1164  G4double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024
1165  + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
1166  G4double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034
1167  + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
1168  G4double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034
1169  + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
1170  G4double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134
1171  + m[M03]*Det3_234_124 - m[M04]*Det3_234_123;
1172  G4double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023
1173  + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
1174  G4double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024
1175  + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
1176  G4double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034
1177  + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
1178  G4double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034
1179  + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
1180  G4double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134
1181  + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
1182 
1183  // Find the 5x5 det:
1184 
1185  G4double det = m[M00]*Det4_1234_1234
1186  - m[M01]*Det4_1234_0234
1187  + m[M02]*Det4_1234_0134
1188  - m[M03]*Det4_1234_0124
1189  + m[M04]*Det4_1234_0123;
1190 
1191  if ( det == 0 )
1192  {
1193  ifail = 1;
1194  return;
1195  }
1196 
1197  G4double oneOverDet = 1.0/det;
1198  G4double mn1OverDet = - oneOverDet;
1199 
1200  m[M00] = Det4_1234_1234 * oneOverDet;
1201  m[M01] = Det4_0234_1234 * mn1OverDet;
1202  m[M02] = Det4_0134_1234 * oneOverDet;
1203  m[M03] = Det4_0124_1234 * mn1OverDet;
1204  m[M04] = Det4_0123_1234 * oneOverDet;
1205 
1206  m[M10] = Det4_1234_0234 * mn1OverDet;
1207  m[M11] = Det4_0234_0234 * oneOverDet;
1208  m[M12] = Det4_0134_0234 * mn1OverDet;
1209  m[M13] = Det4_0124_0234 * oneOverDet;
1210  m[M14] = Det4_0123_0234 * mn1OverDet;
1211 
1212  m[M20] = Det4_1234_0134 * oneOverDet;
1213  m[M21] = Det4_0234_0134 * mn1OverDet;
1214  m[M22] = Det4_0134_0134 * oneOverDet;
1215  m[M23] = Det4_0124_0134 * mn1OverDet;
1216  m[M24] = Det4_0123_0134 * oneOverDet;
1217 
1218  m[M30] = Det4_1234_0124 * mn1OverDet;
1219  m[M31] = Det4_0234_0124 * oneOverDet;
1220  m[M32] = Det4_0134_0124 * mn1OverDet;
1221  m[M33] = Det4_0124_0124 * oneOverDet;
1222  m[M34] = Det4_0123_0124 * mn1OverDet;
1223 
1224  m[M40] = Det4_1234_0123 * oneOverDet;
1225  m[M41] = Det4_0234_0123 * mn1OverDet;
1226  m[M42] = Det4_0134_0123 * oneOverDet;
1227  m[M43] = Det4_0124_0123 * mn1OverDet;
1228  m[M44] = Det4_0123_0123 * oneOverDet;
1229 
1230  return;
1231 }
1232 
1233 void G4ErrorMatrix::invertHaywood6 (G4int & ifail)
1234 {
1235  ifail = 0;
1236 
1237  // Find all NECESSARY 2x2 dets: (45 of them)
1238 
1239  G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1240  G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1241  G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1242  G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1243  G4double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40];
1244  G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1245  G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1246  G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1247  G4double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41];
1248  G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1249  G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1250  G4double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42];
1251  G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1252  G4double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43];
1253  G4double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44];
1254  G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1255  G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1256  G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1257  G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1258  G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1259  G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1260  G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1261  G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1262  G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1263  G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1264  G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1265  G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1266  G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1267  G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1268  G4double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54];
1269  G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1270  G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1271  G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1272  G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1273  G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1274  G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1275  G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1276  G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1277  G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1278  G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1279  G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1280  G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1281  G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1282  G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1283  G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1284 
1285  // Find all NECESSARY 3x3 dets: (80 of them)
1286 
1287  G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1288  + m[A22]*Det2_34_01;
1289  G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1290  + m[A23]*Det2_34_01;
1291  G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1292  + m[A24]*Det2_34_01;
1293  G4double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
1294  + m[A25]*Det2_34_01;
1295  G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1296  + m[A23]*Det2_34_02;
1297  G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1298  + m[A24]*Det2_34_02;
1299  G4double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05
1300  + m[A25]*Det2_34_02;
1301  G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1302  + m[A24]*Det2_34_03;
1303  G4double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05
1304  + m[A25]*Det2_34_03;
1305  G4double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05
1306  + m[A25]*Det2_34_04;
1307  G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1308  + m[A23]*Det2_34_12;
1309  G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1310  + m[A24]*Det2_34_12;
1311  G4double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15
1312  + m[A25]*Det2_34_12;
1313  G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1314  + m[A24]*Det2_34_13;
1315  G4double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15
1316  + m[A25]*Det2_34_13;
1317  G4double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15
1318  + m[A25]*Det2_34_14;
1319  G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1320  + m[A24]*Det2_34_23;
1321  G4double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
1322  + m[A25]*Det2_34_23;
1323  G4double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
1324  + m[A25]*Det2_34_24;
1325  G4double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
1326  + m[A25]*Det2_34_34;
1327  G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
1328  + m[A22]*Det2_35_01;
1329  G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
1330  + m[A23]*Det2_35_01;
1331  G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
1332  + m[A24]*Det2_35_01;
1333  G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
1334  + m[A25]*Det2_35_01;
1335  G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
1336  + m[A23]*Det2_35_02;
1337  G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
1338  + m[A24]*Det2_35_02;
1339  G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
1340  + m[A25]*Det2_35_02;
1341  G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
1342  + m[A24]*Det2_35_03;
1343  G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
1344  + m[A25]*Det2_35_03;
1345  G4double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05
1346  + m[A25]*Det2_35_04;
1347  G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
1348  + m[A23]*Det2_35_12;
1349  G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
1350  + m[A24]*Det2_35_12;
1351  G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
1352  + m[A25]*Det2_35_12;
1353  G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
1354  + m[A24]*Det2_35_13;
1355  G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
1356  + m[A25]*Det2_35_13;
1357  G4double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15
1358  + m[A25]*Det2_35_14;
1359  G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
1360  + m[A24]*Det2_35_23;
1361  G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
1362  + m[A25]*Det2_35_23;
1363  G4double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25
1364  + m[A25]*Det2_35_24;
1365  G4double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35
1366  + m[A25]*Det2_35_34;
1367  G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
1368  + m[A22]*Det2_45_01;
1369  G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
1370  + m[A23]*Det2_45_01;
1371  G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
1372  + m[A24]*Det2_45_01;
1373  G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
1374  + m[A25]*Det2_45_01;
1375  G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
1376  + m[A23]*Det2_45_02;
1377  G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
1378  + m[A24]*Det2_45_02;
1379  G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
1380  + m[A25]*Det2_45_02;
1381  G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
1382  + m[A24]*Det2_45_03;
1383  G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
1384  + m[A25]*Det2_45_03;
1385  G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
1386  + m[A25]*Det2_45_04;
1387  G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
1388  + m[A23]*Det2_45_12;
1389  G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
1390  + m[A24]*Det2_45_12;
1391  G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
1392  + m[A25]*Det2_45_12;
1393  G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
1394  + m[A24]*Det2_45_13;
1395  G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
1396  + m[A25]*Det2_45_13;
1397  G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
1398  + m[A25]*Det2_45_14;
1399  G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
1400  + m[A24]*Det2_45_23;
1401  G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
1402  + m[A25]*Det2_45_23;
1403  G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
1404  + m[A25]*Det2_45_24;
1405  G4double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35
1406  + m[A25]*Det2_45_34;
1407  G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
1408  + m[A32]*Det2_45_01;
1409  G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
1410  + m[A33]*Det2_45_01;
1411  G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
1412  + m[A34]*Det2_45_01;
1413  G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
1414  + m[A35]*Det2_45_01;
1415  G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
1416  + m[A33]*Det2_45_02;
1417  G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
1418  + m[A34]*Det2_45_02;
1419  G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
1420  + m[A35]*Det2_45_02;
1421  G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
1422  + m[A34]*Det2_45_03;
1423  G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
1424  + m[A35]*Det2_45_03;
1425  G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
1426  + m[A35]*Det2_45_04;
1427  G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
1428  + m[A33]*Det2_45_12;
1429  G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
1430  + m[A34]*Det2_45_12;
1431  G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
1432  + m[A35]*Det2_45_12;
1433  G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
1434  + m[A34]*Det2_45_13;
1435  G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
1436  + m[A35]*Det2_45_13;
1437  G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
1438  + m[A35]*Det2_45_14;
1439  G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
1440  + m[A34]*Det2_45_23;
1441  G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
1442  + m[A35]*Det2_45_23;
1443  G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
1444  + m[A35]*Det2_45_24;
1445  G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
1446  + m[A35]*Det2_45_34;
1447 
1448  // Find all NECESSARY 4x4 dets: (75 of them)
1449 
1450  G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1451  + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1452  G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1453  + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1454  G4double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025
1455  + m[A12]*Det3_234_015 - m[A15]*Det3_234_012;
1456  G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1457  + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1458  G4double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
1459  + m[A13]*Det3_234_015 - m[A15]*Det3_234_013;
1460  G4double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
1461  + m[A14]*Det3_234_015 - m[A15]*Det3_234_014;
1462  G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1463  + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1464  G4double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035
1465  + m[A13]*Det3_234_025 - m[A15]*Det3_234_023;
1466  G4double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045
1467  + m[A14]*Det3_234_025 - m[A15]*Det3_234_024;
1468  G4double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045
1469  + m[A14]*Det3_234_035 - m[A15]*Det3_234_034;
1470  G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1471  + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1472  G4double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135
1473  + m[A13]*Det3_234_125 - m[A15]*Det3_234_123;
1474  G4double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145
1475  + m[A14]*Det3_234_125 - m[A15]*Det3_234_124;
1476  G4double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145
1477  + m[A14]*Det3_234_135 - m[A15]*Det3_234_134;
1478  G4double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245
1479  + m[A14]*Det3_234_235 - m[A15]*Det3_234_234;
1480  G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
1481  + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1482  G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
1483  + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1484  G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
1485  + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1486  G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
1487  + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1488  G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
1489  + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1490  G4double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045
1491  + m[A14]*Det3_235_015 - m[A15]*Det3_235_014;
1492  G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
1493  + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1494  G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
1495  + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1496  G4double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045
1497  + m[A14]*Det3_235_025 - m[A15]*Det3_235_024;
1498  G4double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045
1499  + m[A14]*Det3_235_035 - m[A15]*Det3_235_034;
1500  G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
1501  + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1502  G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
1503  + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1504  G4double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145
1505  + m[A14]*Det3_235_125 - m[A15]*Det3_235_124;
1506  G4double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145
1507  + m[A14]*Det3_235_135 - m[A15]*Det3_235_134;
1508  G4double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245
1509  + m[A14]*Det3_235_235 - m[A15]*Det3_235_234;
1510  G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
1511  + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1512  G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
1513  + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1514  G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
1515  + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1516  G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
1517  + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1518  G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
1519  + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1520  G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
1521  + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1522  G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
1523  + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1524  G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
1525  + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1526  G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
1527  + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1528  G4double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045
1529  + m[A14]*Det3_245_035 - m[A15]*Det3_245_034;
1530  G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
1531  + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1532  G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
1533  + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1534  G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
1535  + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1536  G4double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145
1537  + m[A14]*Det3_245_135 - m[A15]*Det3_245_134;
1538  G4double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245
1539  + m[A14]*Det3_245_235 - m[A15]*Det3_245_234;
1540  G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
1541  + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1542  G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
1543  + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1544  G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
1545  + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1546  G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
1547  + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1548  G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
1549  + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1550  G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
1551  + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1552  G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
1553  + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1554  G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
1555  + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1556  G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
1557  + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1558  G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
1559  + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1560  G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
1561  + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1562  G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
1563  + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1564  G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
1565  + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1566  G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
1567  + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1568  G4double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245
1569  + m[A14]*Det3_345_235 - m[A15]*Det3_345_234;
1570  G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
1571  + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1572  G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
1573  + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1574  G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
1575  + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1576  G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
1577  + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1578  G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
1579  + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1580  G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
1581  + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1582  G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
1583  + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1584  G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
1585  + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1586  G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
1587  + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1588  G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
1589  + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1590  G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
1591  + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1592  G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
1593  + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1594  G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
1595  + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1596  G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
1597  + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1598  G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
1599  + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1600 
1601  // Find all NECESSARY 5x5 dets: (36 of them)
1602 
1603  G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
1604  + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1605  G4double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
1606  + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123;
1607  G4double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
1608  + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124;
1609  G4double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
1610  + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134;
1611  G4double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
1612  + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234;
1613  G4double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
1614  + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234;
1615  G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
1616  + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1617  G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
1618  + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1619  G4double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245
1620  + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124;
1621  G4double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345
1622  + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134;
1623  G4double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345
1624  + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234;
1625  G4double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345
1626  + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234;
1627  G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
1628  + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1629  G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
1630  + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1631  G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
1632  + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1633  G4double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345
1634  + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134;
1635  G4double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345
1636  + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234;
1637  G4double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345
1638  + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234;
1639  G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
1640  + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1641  G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
1642  + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1643  G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
1644  + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1645  G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
1646  + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1647  G4double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345
1648  + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234;
1649  G4double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345
1650  + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234; //
1651  G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
1652  + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1653  G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
1654  + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1655  G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
1656  + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1657  G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
1658  + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1659  G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
1660  + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1661  G4double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345
1662  + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234;
1663  G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
1664  + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1665  G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
1666  + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1667  G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
1668  + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1669  G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
1670  + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1671  G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
1672  + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1673  G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
1674  + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1675 
1676  // Find the determinant
1677 
1678  G4double det = m[A00]*Det5_12345_12345
1679  - m[A01]*Det5_12345_02345
1680  + m[A02]*Det5_12345_01345
1681  - m[A03]*Det5_12345_01245
1682  + m[A04]*Det5_12345_01235
1683  - m[A05]*Det5_12345_01234;
1684 
1685  if ( det == 0 )
1686  {
1687  ifail = 1;
1688  return;
1689  }
1690 
1691  G4double oneOverDet = 1.0/det;
1692  G4double mn1OverDet = - oneOverDet;
1693 
1694  m[A00] = Det5_12345_12345*oneOverDet;
1695  m[A01] = Det5_02345_12345*mn1OverDet;
1696  m[A02] = Det5_01345_12345*oneOverDet;
1697  m[A03] = Det5_01245_12345*mn1OverDet;
1698  m[A04] = Det5_01235_12345*oneOverDet;
1699  m[A05] = Det5_01234_12345*mn1OverDet;
1700 
1701  m[A10] = Det5_12345_02345*mn1OverDet;
1702  m[A11] = Det5_02345_02345*oneOverDet;
1703  m[A12] = Det5_01345_02345*mn1OverDet;
1704  m[A13] = Det5_01245_02345*oneOverDet;
1705  m[A14] = Det5_01235_02345*mn1OverDet;
1706  m[A15] = Det5_01234_02345*oneOverDet;
1707 
1708  m[A20] = Det5_12345_01345*oneOverDet;
1709  m[A21] = Det5_02345_01345*mn1OverDet;
1710  m[A22] = Det5_01345_01345*oneOverDet;
1711  m[A23] = Det5_01245_01345*mn1OverDet;
1712  m[A24] = Det5_01235_01345*oneOverDet;
1713  m[A25] = Det5_01234_01345*mn1OverDet;
1714 
1715  m[A30] = Det5_12345_01245*mn1OverDet;
1716  m[A31] = Det5_02345_01245*oneOverDet;
1717  m[A32] = Det5_01345_01245*mn1OverDet;
1718  m[A33] = Det5_01245_01245*oneOverDet;
1719  m[A34] = Det5_01235_01245*mn1OverDet;
1720  m[A35] = Det5_01234_01245*oneOverDet;
1721 
1722  m[A40] = Det5_12345_01235*oneOverDet;
1723  m[A41] = Det5_02345_01235*mn1OverDet;
1724  m[A42] = Det5_01345_01235*oneOverDet;
1725  m[A43] = Det5_01245_01235*mn1OverDet;
1726  m[A44] = Det5_01235_01235*oneOverDet;
1727  m[A45] = Det5_01234_01235*mn1OverDet;
1728 
1729  m[A50] = Det5_12345_01234*mn1OverDet;
1730  m[A51] = Det5_02345_01234*oneOverDet;
1731  m[A52] = Det5_01345_01234*mn1OverDet;
1732  m[A53] = Det5_01245_01234*oneOverDet;
1733  m[A54] = Det5_01235_01234*mn1OverDet;
1734  m[A55] = Det5_01234_01234*oneOverDet;
1735 
1736  return;
1737 }
A21
#define A21
Definition: G4ErrorMatrix.cc:833
A14
#define A14
Definition: G4ErrorMatrix.cc:829
M30
#define M30
Definition: G4ErrorMatrix.cc:878
F33
#define F33
Definition: G4ErrorMatrix.cc:908
G4ErrorMatrix::num_row
virtual G4int num_row() const
operator<<
std::ostream & operator<<(std::ostream &os, const G4ErrorMatrix &q)
Definition: G4ErrorMatrix.cc:350
A42
#define A42
Definition: G4ErrorMatrix.cc:848
A13
#define A13
Definition: G4ErrorMatrix.cc:828
F13
#define F13
Definition: G4ErrorMatrix.cc:898
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: G4ErrorMatrix.cc:62
M04
#define M04
Definition: G4ErrorMatrix.cc:864
F00
#define F00
Definition: G4ErrorMatrix.cc:890
A01
#define A01
Definition: G4ErrorMatrix.cc:819
A00
#define A00
Definition: G4ErrorMatrix.cc:818
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorMatrix.cc:67
A31
#define A31
Definition: G4ErrorMatrix.cc:840
M43
#define M43
Definition: G4ErrorMatrix.cc:887
M44
#define M44
Definition: G4ErrorMatrix.cc:888
G4ErrorMatrix::num_col
virtual G4int num_col() const
A22
#define A22
Definition: G4ErrorMatrix.cc:834
A10
#define A10
Definition: G4ErrorMatrix.cc:825
M03
#define M03
Definition: G4ErrorMatrix.cc:863
G4ErrorMatrix.hh
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorMatrix.cc:42
F30
#define F30
Definition: G4ErrorMatrix.cc:905
A52
#define A52
Definition: G4ErrorMatrix.cc:855
A50
#define A50
Definition: G4ErrorMatrix.cc:853
F22
#define F22
Definition: G4ErrorMatrix.cc:902
width
#define width
F11
#define F11
Definition: G4ErrorMatrix.cc:896
M31
#define M31
Definition: G4ErrorMatrix.cc:879
G4ErrorMatrix::~G4ErrorMatrix
virtual ~G4ErrorMatrix()
Definition: G4ErrorMatrix.cc:113
a
G4double a
Definition: TRTMaterials.hh:39
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: G4ErrorMatrix.cc:53
dsum
G4ErrorMatrix dsum(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
Definition: G4ErrorMatrix.cc:201
F20
#define F20
Definition: G4ErrorMatrix.cc:900
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
F23
#define F23
Definition: G4ErrorMatrix.cc:903
G4int
int G4int
Definition: G4Types.hh:78
A12
#define A12
Definition: G4ErrorMatrix.cc:827
G4ErrorSymMatrix.hh
A54
#define A54
Definition: G4ErrorMatrix.cc:857
G4ErrorMatrix::operator*=
G4ErrorMatrix & operator*=(G4double t)
Definition: G4ErrorMatrix.cc:326
G4ErrorMatrix::operator=
G4ErrorMatrix & operator=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:332
A25
#define A25
Definition: G4ErrorMatrix.cc:837
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
CLHEP::normal
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4ErrorMatrix::sub
G4ErrorMatrix sub(G4int min_row, G4int max_row, G4int min_col, G4int max_col) const
Definition: G4ErrorMatrix.cc:155
operator/
G4ErrorMatrix operator/(const G4ErrorMatrix &mat1, G4double t)
Definition: G4ErrorMatrix.cc:251
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: G4ErrorMatrix.cc:47
M13
#define M13
Definition: G4ErrorMatrix.cc:869
G4ErrorMatrix::invert
virtual void invert(G4int &ierr)
Definition: G4ErrorMatrix.cc:645
M01
#define M01
Definition: G4ErrorMatrix.cc:861
A03
#define A03
Definition: G4ErrorMatrix.cc:821
M34
#define M34
Definition: G4ErrorMatrix.cc:882
M14
#define M14
Definition: G4ErrorMatrix.cc:870
M41
#define M41
Definition: G4ErrorMatrix.cc:885
A35
#define A35
Definition: G4ErrorMatrix.cc:844
A02
#define A02
Definition: G4ErrorMatrix.cc:820
M40
#define M40
Definition: G4ErrorMatrix.cc:884
A55
#define A55
Definition: G4ErrorMatrix.cc:858
G4ErrorSymMatrix::m
std::vector< G4double > m
Definition: G4ErrorSymMatrix.hh:217
A30
#define A30
Definition: G4ErrorMatrix.cc:839
M12
#define M12
Definition: G4ErrorMatrix.cc:868
operator+
G4ErrorMatrix operator+(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
Definition: G4ErrorMatrix.cc:225
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:87
G4ErrorMatrix::operator/=
G4ErrorMatrix & operator/=(G4double t)
Definition: G4ErrorMatrix.cc:320
G4ErrorMatrix::m
std::vector< G4double > m
Definition: G4ErrorMatrix.hh:224
M24
#define M24
Definition: G4ErrorMatrix.cc:876
M00
#define M00
Definition: G4ErrorMatrix.cc:860
A45
#define A45
Definition: G4ErrorMatrix.cc:851
F10
#define F10
Definition: G4ErrorMatrix.cc:895
A15
#define A15
Definition: G4ErrorMatrix.cc:830
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
M21
#define M21
Definition: G4ErrorMatrix.cc:873
M42
#define M42
Definition: G4ErrorMatrix.cc:886
G4ErrorMatrix::dfact_matrix
G4int dfact_matrix(G4double &det, G4int *ir)
Definition: G4ErrorMatrix.cc:518
A43
#define A43
Definition: G4ErrorMatrix.cc:849
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:805
A32
#define A32
Definition: G4ErrorMatrix.cc:841
A23
#define A23
Definition: G4ErrorMatrix.cc:835
G4ErrorMatrix::T
G4ErrorMatrix T() const
Definition: G4ErrorMatrix.cc:374
G4ErrorMatrix::trace
G4double trace() const
Definition: G4ErrorMatrix.cc:795
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
F32
#define F32
Definition: G4ErrorMatrix.cc:907
G4ErrorMatrix::apply
G4ErrorMatrix apply(G4double(*f)(G4double, G4int, G4int)) const
Definition: G4ErrorMatrix.cc:390
G4ErrorMatrix::dfinv_matrix
G4int dfinv_matrix(G4int *ir)
Definition: G4ErrorMatrix.cc:405
A53
#define A53
Definition: G4ErrorMatrix.cc:856
M23
#define M23
Definition: G4ErrorMatrix.cc:875
A04
#define A04
Definition: G4ErrorMatrix.cc:822
operator*
G4ErrorMatrix operator*(const G4ErrorMatrix &mat1, G4double t)
Definition: G4ErrorMatrix.cc:258
A40
#define A40
Definition: G4ErrorMatrix.cc:846
A05
#define A05
Definition: G4ErrorMatrix.cc:823
A51
#define A51
Definition: G4ErrorMatrix.cc:854
F02
#define F02
Definition: G4ErrorMatrix.cc:892
A41
#define A41
Definition: G4ErrorMatrix.cc:847
M22
#define M22
Definition: G4ErrorMatrix.cc:874
A34
#define A34
Definition: G4ErrorMatrix.cc:843
F01
#define F01
Definition: G4ErrorMatrix.cc:891
A11
#define A11
Definition: G4ErrorMatrix.cc:826
A33
#define A33
Definition: G4ErrorMatrix.cc:842
FatalException
Definition: G4ExceptionSeverity.hh:60
G4ErrorMatrix::invertHaywood4
virtual void invertHaywood4(G4int &ierr)
Definition: G4ErrorMatrix.cc:911
b1
static const G4double b1
Definition: G4BetheHeitlerModel.cc:76
G4ErrorMatrix::determinant
G4double determinant() const
Definition: G4ErrorMatrix.cc:776
A24
#define A24
Definition: G4ErrorMatrix.cc:836
G4endl
#define G4endl
Definition: G4ios.hh:61
G4ErrorMatrix::invertHaywood6
virtual void invertHaywood6(G4int &ierr)
Definition: G4ErrorMatrix.cc:1233
m
static const double m
Definition: G4SIunits.hh:128
F03
#define F03
Definition: G4ErrorMatrix.cc:893
F12
#define F12
Definition: G4ErrorMatrix.cc:897
G4ErrorMatrix::invertHaywood5
virtual void invertHaywood5(G4int &ierr)
Definition: G4ErrorMatrix.cc:1010
G4ErrorMatrix::operator-
G4ErrorMatrix operator-() const
Definition: G4ErrorMatrix.cc:215
A44
#define A44
Definition: G4ErrorMatrix.cc:850
G4double
double G4double
Definition: G4Types.hh:76
M10
#define M10
Definition: G4ErrorMatrix.cc:866
G4ErrorMatrix::operator+=
G4ErrorMatrix & operator+=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:306
M33
#define M33
Definition: G4ErrorMatrix.cc:881
G4ErrorMatrix::operator-=
G4ErrorMatrix & operator-=(const G4ErrorMatrix &m2)
Definition: G4ErrorMatrix.cc:313
M32
#define M32
Definition: G4ErrorMatrix.cc:880
A20
#define A20
Definition: G4ErrorMatrix.cc:832
M02
#define M02
Definition: G4ErrorMatrix.cc:862
F31
#define F31
Definition: G4ErrorMatrix.cc:906
F21
#define F21
Definition: G4ErrorMatrix.cc:901
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:95
operator-
G4ErrorMatrix operator-(const G4ErrorMatrix &mat1, const G4ErrorMatrix &mat2)
Definition: G4ErrorMatrix.cc:237
M11
#define M11
Definition: G4ErrorMatrix.cc:867
M20
#define M20
Definition: G4ErrorMatrix.cc:872