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

Last change on this file since 1258 was 1228, checked in by garnier, 14 years ago

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