TMatrixTSymCramerInv.cxx

Go to the documentation of this file.
00001 // @(#)root/base:$Id: TMatrixTSymCramerInv.cxx 20882 2007-11-19 11:31:26Z rdm $
00002 // Authors: Fons Rademakers, Eddy Offermann  Oct 2004
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 //////////////////////////////////////////////////////////////////////////
00013 //                                                                      //
00014 // TMatrixTSymCramerInv                                                 //
00015 //                                                                      //
00016 // Encapsulate templates of Cramer Inversion routines.                  //
00017 //                                                                      //
00018 // The 4x4, 5x5 and 6x6 are adapted from routines written by            //
00019 // Mark Fischler and Steven Haywood as part of the CLHEP package        //
00020 //                                                                      //
00021 // Although for sizes <= 6x6 the Cramer Inversion has a gain in speed   //
00022 // compared to factorization schemes (like LU) , one pays a price in    //
00023 // accuracy  .                                                          //
00024 //                                                                      //
00025 // For Example:                                                         //
00026 //  H * H^-1 = U, where H is a 5x5 Hilbert matrix                       //
00027 //                      U is a 5x5 Unity matrix                         //
00028 //                                                                      //
00029 // LU    : |U_jk| < 10e-13 for  j!=k                                    //
00030 // Cramer: |U_jk| < 10e-7  for  j!=k                                    //
00031 //                                                                      //
00032 //  however Cramer algorithm is about 10 (!) times faster               //
00033 //////////////////////////////////////////////////////////////////////////
00034 
00035 #include "TMatrixTSymCramerInv.h"
00036 
00037 #if !defined(R__ALPHA) && !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
00038 NamespaceImp(TMatrixTSymCramerInv);
00039 #endif
00040 
00041 //______________________________________________________________________________
00042 template<class Element> 
00043 Bool_t TMatrixTSymCramerInv::Inv2x2(TMatrixTSym<Element> &m,Double_t *determ)
00044 {
00045    if (m.GetNrows() != 2) {
00046       Error("Inv2x2","matrix should be square 2x2");
00047       return kFALSE;
00048    }
00049 
00050    Element *pM = m.GetMatrixArray();
00051 
00052    const Double_t det = pM[0] * pM[3] - pM[1] * pM[1];
00053 
00054    if (determ)
00055       *determ = det;
00056 
00057    if ( det == 0 ) {
00058       Error("Inv2x2","matrix is singular");
00059       return kFALSE;
00060    }
00061 
00062    const Double_t tmp1 =   pM[3] / det;
00063    pM[3] = pM[0] / det;
00064    pM[2] = pM[1] = -pM[1] / det;
00065    pM[0] = tmp1;
00066 
00067    return kTRUE;
00068 }
00069 
00070 //______________________________________________________________________________
00071 template<class Element> 
00072 Bool_t TMatrixTSymCramerInv::Inv3x3(TMatrixTSym<Element> &m,Double_t *determ)
00073 {
00074    if (m.GetNrows() != 3) {
00075       Error("Inv3x3","matrix should be square 3x3");
00076       return kFALSE;
00077    }
00078 
00079    Element *pM = m.GetMatrixArray();
00080 
00081    const Double_t c00 = pM[4] * pM[8] - pM[5] * pM[5];
00082    const Double_t c01 = pM[5] * pM[2] - pM[1] * pM[8];
00083    const Double_t c02 = pM[1] * pM[5] - pM[4] * pM[2];
00084    const Double_t c11 = pM[8] * pM[0] - pM[2] * pM[2];
00085    const Double_t c12 = pM[2] * pM[1] - pM[5] * pM[0];
00086    const Double_t c22 = pM[0] * pM[4] - pM[1] * pM[1];
00087 
00088    const Double_t t0  = TMath::Abs(pM[0]);
00089    const Double_t t1  = TMath::Abs(pM[1]);
00090    const Double_t t2  = TMath::Abs(pM[2]);
00091 
00092    Double_t det;
00093    Double_t tmp;
00094 
00095    if (t0 >= t1) {
00096       if (t2 >= t0) {
00097          tmp = pM[2];
00098          det = c12*c01-c11*c02;
00099       } else {
00100          tmp = pM[0];
00101          det = c11*c22-c12*c12;
00102       }
00103    } else if (t2 >= t1) {
00104       tmp = pM[2];
00105       det = c12*c01-c11*c02;
00106    } else {
00107       tmp = pM[1];
00108       det = c02*c12-c01*c22;
00109    }
00110 
00111    if ( det == 0 || tmp == 0) {
00112       Error("Inv3x3","matrix is singular");
00113       return kFALSE;
00114    }
00115 
00116    Double_t s = tmp/det;
00117    if (determ)
00118       *determ = 1./s;
00119 
00120    pM[0] = s*c00;
00121    pM[1] = s*c01;
00122    pM[2] = s*c02;
00123    pM[3] = pM[1];
00124    pM[4] = s*c11;
00125    pM[5] = s*c12;
00126    pM[6] = pM[2];
00127    pM[7] = pM[5];
00128    pM[8] = s*c22;
00129   
00130    return kTRUE;
00131 }
00132 
00133 // SFij are indices for a 4x4 symmetric matrix.
00134 
00135 #define SF00 0
00136 #define SF01 1
00137 #define SF02 2
00138 #define SF03 3
00139 
00140 #define SF10 1
00141 #define SF11 5
00142 #define SF12 6
00143 #define SF13 7
00144 
00145 #define SF20 2
00146 #define SF21 6
00147 #define SF22 10
00148 #define SF23 11
00149 
00150 #define SF30 3
00151 #define SF31 7
00152 #define SF32 11
00153 #define SF33 15
00154 
00155 //______________________________________________________________________________
00156 template<class Element> 
00157 Bool_t TMatrixTSymCramerInv::Inv4x4(TMatrixTSym<Element> &m,Double_t *determ)
00158 {
00159    if (m.GetNrows() != 4) {
00160       Error("Inv4x4","matrix should be square 4x4");
00161       return kFALSE;
00162    }
00163 
00164    Element *pM = m.GetMatrixArray();
00165 
00166    // Find all NECESSARY 2x2 dets:  (14 of them)
00167 
00168    const Double_t mDet2_12_01 = pM[SF10]*pM[SF21] - pM[SF11]*pM[SF20];
00169    const Double_t mDet2_12_02 = pM[SF10]*pM[SF22] - pM[SF12]*pM[SF20];
00170    const Double_t mDet2_12_12 = pM[SF11]*pM[SF22] - pM[SF12]*pM[SF21];
00171    const Double_t mDet2_13_01 = pM[SF10]*pM[SF31] - pM[SF11]*pM[SF30];
00172    const Double_t mDet2_13_02 = pM[SF10]*pM[SF32] - pM[SF12]*pM[SF30];
00173    const Double_t mDet2_13_03 = pM[SF10]*pM[SF33] - pM[SF13]*pM[SF30];
00174    const Double_t mDet2_13_12 = pM[SF11]*pM[SF32] - pM[SF12]*pM[SF31];
00175    const Double_t mDet2_13_13 = pM[SF11]*pM[SF33] - pM[SF13]*pM[SF31];
00176    const Double_t mDet2_23_01 = pM[SF20]*pM[SF31] - pM[SF21]*pM[SF30];
00177    const Double_t mDet2_23_02 = pM[SF20]*pM[SF32] - pM[SF22]*pM[SF30];
00178    const Double_t mDet2_23_03 = pM[SF20]*pM[SF33] - pM[SF23]*pM[SF30];
00179    const Double_t mDet2_23_12 = pM[SF21]*pM[SF32] - pM[SF22]*pM[SF31];
00180    const Double_t mDet2_23_13 = pM[SF21]*pM[SF33] - pM[SF23]*pM[SF31];
00181    const Double_t mDet2_23_23 = pM[SF22]*pM[SF33] - pM[SF23]*pM[SF32];
00182 
00183   // SFind all NECESSSFRY 3x3 dets:   (10 of them)
00184   
00185    const Double_t mDet3_012_012 = pM[SF00]*mDet2_12_12 - pM[SF01]*mDet2_12_02
00186                                 + pM[SF02]*mDet2_12_01;
00187    const Double_t mDet3_013_012 = pM[SF00]*mDet2_13_12 - pM[SF01]*mDet2_13_02
00188                                 + pM[SF02]*mDet2_13_01;
00189    const Double_t mDet3_013_013 = pM[SF00]*mDet2_13_13 - pM[SF01]*mDet2_13_03
00190                                 + pM[SF03]*mDet2_13_01;
00191    const Double_t mDet3_023_012 = pM[SF00]*mDet2_23_12 - pM[SF01]*mDet2_23_02
00192                                 + pM[SF02]*mDet2_23_01;
00193    const Double_t mDet3_023_013 = pM[SF00]*mDet2_23_13 - pM[SF01]*mDet2_23_03
00194                                 + pM[SF03]*mDet2_23_01;
00195    const Double_t mDet3_023_023 = pM[SF00]*mDet2_23_23 - pM[SF02]*mDet2_23_03
00196                                 + pM[SF03]*mDet2_23_02;
00197    const Double_t mDet3_123_012 = pM[SF10]*mDet2_23_12 - pM[SF11]*mDet2_23_02
00198                                 + pM[SF12]*mDet2_23_01;
00199    const Double_t mDet3_123_013 = pM[SF10]*mDet2_23_13 - pM[SF11]*mDet2_23_03
00200                                 + pM[SF13]*mDet2_23_01;
00201    const Double_t mDet3_123_023 = pM[SF10]*mDet2_23_23 - pM[SF12]*mDet2_23_03
00202                                 + pM[SF13]*mDet2_23_02;
00203    const Double_t mDet3_123_123 = pM[SF11]*mDet2_23_23 - pM[SF12]*mDet2_23_13
00204                                 + pM[SF13]*mDet2_23_12;
00205 
00206    // Find the 4x4 det:
00207 
00208    const Double_t det = pM[SF00]*mDet3_123_123 - pM[SF01]*mDet3_123_023
00209                       + pM[SF02]*mDet3_123_013 - pM[SF03]*mDet3_123_012;
00210 
00211    if (determ)
00212       *determ = det;
00213 
00214    if ( det == 0 ) {
00215       Error("Inv4x4","matrix is singular");
00216       return kFALSE;
00217    }
00218 
00219    const Double_t oneOverDet = 1.0/det;
00220    const Double_t mn1OverDet = - oneOverDet;
00221 
00222    pM[SF00] =  mDet3_123_123 * oneOverDet;
00223    pM[SF01] =  mDet3_123_023 * mn1OverDet;
00224    pM[SF02] =  mDet3_123_013 * oneOverDet;
00225    pM[SF03] =  mDet3_123_012 * mn1OverDet;
00226 
00227    pM[SF11] =  mDet3_023_023 * oneOverDet;
00228    pM[SF12] =  mDet3_023_013 * mn1OverDet;
00229    pM[SF13] =  mDet3_023_012 * oneOverDet;
00230 
00231    pM[SF22] =  mDet3_013_013 * oneOverDet;
00232    pM[SF23] =  mDet3_013_012 * mn1OverDet;
00233 
00234    pM[SF33] =  mDet3_012_012 * oneOverDet;
00235 
00236    for (Int_t irow = 0; irow < 4; irow++) {
00237       const Int_t rowOff1 = irow*4; 
00238       for (Int_t icol = 0; icol < irow; icol++) {
00239          const Int_t rowOff2 = icol*4; 
00240          pM[rowOff1+icol] = pM[rowOff2+irow];
00241       }
00242    } 
00243 
00244    return kTRUE;
00245 }
00246 
00247 // Mij are indices for a 5x5 matrix.
00248     
00249 #define SM00 0
00250 #define SM01 1
00251 #define SM02 2
00252 #define SM03 3
00253 #define SM04 4
00254   
00255 #define SM10 1
00256 #define SM11 6
00257 #define SM12 7
00258 #define SM13 8
00259 #define SM14 9
00260 
00261 #define SM20 2
00262 #define SM21 7
00263 #define SM22 12
00264 #define SM23 13
00265 #define SM24 14
00266 
00267 #define SM30 3
00268 #define SM31 8
00269 #define SM32 13
00270 #define SM33 18
00271 #define SM34 19
00272 
00273 #define SM40 4
00274 #define SM41 9
00275 #define SM42 14
00276 #define SM43 19
00277 #define SM44 24
00278 
00279 //______________________________________________________________________________
00280 template<class Element> 
00281 Bool_t TMatrixTSymCramerInv::Inv5x5(TMatrixTSym<Element> &m,Double_t *determ)
00282 {
00283    if (m.GetNrows() != 5) {
00284       Error("Inv5x5","matrix should be square 5x5");
00285       return kFALSE;
00286    }
00287 
00288    Element *pM = m.GetMatrixArray();
00289 
00290    // Find all NECESSARY 2x2 dets:  (25 of them)
00291 
00292    const Double_t mDet2_23_01 = pM[SM20]*pM[SM31] - pM[SM21]*pM[SM30];
00293    const Double_t mDet2_23_02 = pM[SM20]*pM[SM32] - pM[SM22]*pM[SM30];
00294    const Double_t mDet2_23_03 = pM[SM20]*pM[SM33] - pM[SM23]*pM[SM30];
00295    const Double_t mDet2_23_12 = pM[SM21]*pM[SM32] - pM[SM22]*pM[SM31];
00296    const Double_t mDet2_23_13 = pM[SM21]*pM[SM33] - pM[SM23]*pM[SM31];
00297    const Double_t mDet2_23_23 = pM[SM22]*pM[SM33] - pM[SM23]*pM[SM32];
00298    const Double_t mDet2_24_01 = pM[SM20]*pM[SM41] - pM[SM21]*pM[SM40];
00299    const Double_t mDet2_24_02 = pM[SM20]*pM[SM42] - pM[SM22]*pM[SM40];
00300    const Double_t mDet2_24_03 = pM[SM20]*pM[SM43] - pM[SM23]*pM[SM40];
00301    const Double_t mDet2_24_04 = pM[SM20]*pM[SM44] - pM[SM24]*pM[SM40];
00302    const Double_t mDet2_24_12 = pM[SM21]*pM[SM42] - pM[SM22]*pM[SM41];
00303    const Double_t mDet2_24_13 = pM[SM21]*pM[SM43] - pM[SM23]*pM[SM41];
00304    const Double_t mDet2_24_14 = pM[SM21]*pM[SM44] - pM[SM24]*pM[SM41];
00305    const Double_t mDet2_24_23 = pM[SM22]*pM[SM43] - pM[SM23]*pM[SM42];
00306    const Double_t mDet2_24_24 = pM[SM22]*pM[SM44] - pM[SM24]*pM[SM42];
00307    const Double_t mDet2_34_01 = pM[SM30]*pM[SM41] - pM[SM31]*pM[SM40];
00308    const Double_t mDet2_34_02 = pM[SM30]*pM[SM42] - pM[SM32]*pM[SM40];
00309    const Double_t mDet2_34_03 = pM[SM30]*pM[SM43] - pM[SM33]*pM[SM40];
00310    const Double_t mDet2_34_04 = pM[SM30]*pM[SM44] - pM[SM34]*pM[SM40];
00311    const Double_t mDet2_34_12 = pM[SM31]*pM[SM42] - pM[SM32]*pM[SM41];
00312    const Double_t mDet2_34_13 = pM[SM31]*pM[SM43] - pM[SM33]*pM[SM41];
00313    const Double_t mDet2_34_14 = pM[SM31]*pM[SM44] - pM[SM34]*pM[SM41];
00314    const Double_t mDet2_34_23 = pM[SM32]*pM[SM43] - pM[SM33]*pM[SM42];
00315    const Double_t mDet2_34_24 = pM[SM32]*pM[SM44] - pM[SM34]*pM[SM42];
00316    const Double_t mDet2_34_34 = pM[SM33]*pM[SM44] - pM[SM34]*pM[SM43];
00317 
00318    // Find all NECESSARY 3x3 dets:   (30 of them)
00319 
00320    const Double_t mDet3_123_012 = pM[SM10]*mDet2_23_12 - pM[SM11]*mDet2_23_02 + pM[SM12]*mDet2_23_01;
00321    const Double_t mDet3_123_013 = pM[SM10]*mDet2_23_13 - pM[SM11]*mDet2_23_03 + pM[SM13]*mDet2_23_01;
00322    const Double_t mDet3_123_023 = pM[SM10]*mDet2_23_23 - pM[SM12]*mDet2_23_03 + pM[SM13]*mDet2_23_02;
00323    const Double_t mDet3_123_123 = pM[SM11]*mDet2_23_23 - pM[SM12]*mDet2_23_13 + pM[SM13]*mDet2_23_12;
00324    const Double_t mDet3_124_012 = pM[SM10]*mDet2_24_12 - pM[SM11]*mDet2_24_02 + pM[SM12]*mDet2_24_01;
00325    const Double_t mDet3_124_013 = pM[SM10]*mDet2_24_13 - pM[SM11]*mDet2_24_03 + pM[SM13]*mDet2_24_01;
00326    const Double_t mDet3_124_014 = pM[SM10]*mDet2_24_14 - pM[SM11]*mDet2_24_04 + pM[SM14]*mDet2_24_01;
00327    const Double_t mDet3_124_023 = pM[SM10]*mDet2_24_23 - pM[SM12]*mDet2_24_03 + pM[SM13]*mDet2_24_02;
00328    const Double_t mDet3_124_024 = pM[SM10]*mDet2_24_24 - pM[SM12]*mDet2_24_04 + pM[SM14]*mDet2_24_02;
00329    const Double_t mDet3_124_123 = pM[SM11]*mDet2_24_23 - pM[SM12]*mDet2_24_13 + pM[SM13]*mDet2_24_12;
00330    const Double_t mDet3_124_124 = pM[SM11]*mDet2_24_24 - pM[SM12]*mDet2_24_14 + pM[SM14]*mDet2_24_12;
00331    const Double_t mDet3_134_012 = pM[SM10]*mDet2_34_12 - pM[SM11]*mDet2_34_02 + pM[SM12]*mDet2_34_01;
00332    const Double_t mDet3_134_013 = pM[SM10]*mDet2_34_13 - pM[SM11]*mDet2_34_03 + pM[SM13]*mDet2_34_01;
00333    const Double_t mDet3_134_014 = pM[SM10]*mDet2_34_14 - pM[SM11]*mDet2_34_04 + pM[SM14]*mDet2_34_01;
00334    const Double_t mDet3_134_023 = pM[SM10]*mDet2_34_23 - pM[SM12]*mDet2_34_03 + pM[SM13]*mDet2_34_02;
00335    const Double_t mDet3_134_024 = pM[SM10]*mDet2_34_24 - pM[SM12]*mDet2_34_04 + pM[SM14]*mDet2_34_02;
00336    const Double_t mDet3_134_034 = pM[SM10]*mDet2_34_34 - pM[SM13]*mDet2_34_04 + pM[SM14]*mDet2_34_03;
00337    const Double_t mDet3_134_123 = pM[SM11]*mDet2_34_23 - pM[SM12]*mDet2_34_13 + pM[SM13]*mDet2_34_12;
00338    const Double_t mDet3_134_124 = pM[SM11]*mDet2_34_24 - pM[SM12]*mDet2_34_14 + pM[SM14]*mDet2_34_12;
00339    const Double_t mDet3_134_134 = pM[SM11]*mDet2_34_34 - pM[SM13]*mDet2_34_14 + pM[SM14]*mDet2_34_13;
00340    const Double_t mDet3_234_012 = pM[SM20]*mDet2_34_12 - pM[SM21]*mDet2_34_02 + pM[SM22]*mDet2_34_01;
00341    const Double_t mDet3_234_013 = pM[SM20]*mDet2_34_13 - pM[SM21]*mDet2_34_03 + pM[SM23]*mDet2_34_01;
00342    const Double_t mDet3_234_014 = pM[SM20]*mDet2_34_14 - pM[SM21]*mDet2_34_04 + pM[SM24]*mDet2_34_01;
00343    const Double_t mDet3_234_023 = pM[SM20]*mDet2_34_23 - pM[SM22]*mDet2_34_03 + pM[SM23]*mDet2_34_02;
00344    const Double_t mDet3_234_024 = pM[SM20]*mDet2_34_24 - pM[SM22]*mDet2_34_04 + pM[SM24]*mDet2_34_02;
00345    const Double_t mDet3_234_034 = pM[SM20]*mDet2_34_34 - pM[SM23]*mDet2_34_04 + pM[SM24]*mDet2_34_03;
00346    const Double_t mDet3_234_123 = pM[SM21]*mDet2_34_23 - pM[SM22]*mDet2_34_13 + pM[SM23]*mDet2_34_12;
00347    const Double_t mDet3_234_124 = pM[SM21]*mDet2_34_24 - pM[SM22]*mDet2_34_14 + pM[SM24]*mDet2_34_12;
00348    const Double_t mDet3_234_134 = pM[SM21]*mDet2_34_34 - pM[SM23]*mDet2_34_14 + pM[SM24]*mDet2_34_13;
00349    const Double_t mDet3_234_234 = pM[SM22]*mDet2_34_34 - pM[SM23]*mDet2_34_24 + pM[SM24]*mDet2_34_23;
00350 
00351    // Find all NECESSARY 4x4 dets:   (15 of them)
00352 
00353    const Double_t mDet4_0123_0123 = pM[SM00]*mDet3_123_123 - pM[SM01]*mDet3_123_023
00354                                   + pM[SM02]*mDet3_123_013 - pM[SM03]*mDet3_123_012;
00355    const Double_t mDet4_0124_0123 = pM[SM00]*mDet3_124_123 - pM[SM01]*mDet3_124_023
00356                                   + pM[SM02]*mDet3_124_013 - pM[SM03]*mDet3_124_012;
00357    const Double_t mDet4_0124_0124 = pM[SM00]*mDet3_124_124 - pM[SM01]*mDet3_124_024
00358                                   + pM[SM02]*mDet3_124_014 - pM[SM04]*mDet3_124_012;
00359    const Double_t mDet4_0134_0123 = pM[SM00]*mDet3_134_123 - pM[SM01]*mDet3_134_023
00360                                   + pM[SM02]*mDet3_134_013 - pM[SM03]*mDet3_134_012;
00361    const Double_t mDet4_0134_0124 = pM[SM00]*mDet3_134_124 - pM[SM01]*mDet3_134_024
00362                                   + pM[SM02]*mDet3_134_014 - pM[SM04]*mDet3_134_012;
00363    const Double_t mDet4_0134_0134 = pM[SM00]*mDet3_134_134 - pM[SM01]*mDet3_134_034
00364                                   + pM[SM03]*mDet3_134_014 - pM[SM04]*mDet3_134_013;
00365    const Double_t mDet4_0234_0123 = pM[SM00]*mDet3_234_123 - pM[SM01]*mDet3_234_023
00366                                   + pM[SM02]*mDet3_234_013 - pM[SM03]*mDet3_234_012;
00367    const Double_t mDet4_0234_0124 = pM[SM00]*mDet3_234_124 - pM[SM01]*mDet3_234_024
00368                                   + pM[SM02]*mDet3_234_014 - pM[SM04]*mDet3_234_012;
00369    const Double_t mDet4_0234_0134 = pM[SM00]*mDet3_234_134 - pM[SM01]*mDet3_234_034
00370                                   + pM[SM03]*mDet3_234_014 - pM[SM04]*mDet3_234_013;
00371    const Double_t mDet4_0234_0234 = pM[SM00]*mDet3_234_234 - pM[SM02]*mDet3_234_034
00372                                   + pM[SM03]*mDet3_234_024 - pM[SM04]*mDet3_234_023;
00373    const Double_t mDet4_1234_0123 = pM[SM10]*mDet3_234_123 - pM[SM11]*mDet3_234_023
00374                                   + pM[SM12]*mDet3_234_013 - pM[SM13]*mDet3_234_012;
00375    const Double_t mDet4_1234_0124 = pM[SM10]*mDet3_234_124 - pM[SM11]*mDet3_234_024
00376                                   + pM[SM12]*mDet3_234_014 - pM[SM14]*mDet3_234_012;
00377    const Double_t mDet4_1234_0134 = pM[SM10]*mDet3_234_134 - pM[SM11]*mDet3_234_034
00378                                   + pM[SM13]*mDet3_234_014 - pM[SM14]*mDet3_234_013;
00379    const Double_t mDet4_1234_0234 = pM[SM10]*mDet3_234_234 - pM[SM12]*mDet3_234_034
00380                                   + pM[SM13]*mDet3_234_024 - pM[SM14]*mDet3_234_023;
00381    const Double_t mDet4_1234_1234 = pM[SM11]*mDet3_234_234 - pM[SM12]*mDet3_234_134
00382                                   + pM[SM13]*mDet3_234_124 - pM[SM14]*mDet3_234_123;
00383 
00384    // Find the 5x5 det:
00385 
00386    const Double_t det = pM[SM00]*mDet4_1234_1234 - pM[SM01]*mDet4_1234_0234 + pM[SM02]*mDet4_1234_0134
00387                       - pM[SM03]*mDet4_1234_0124 + pM[SM04]*mDet4_1234_0123;
00388    if (determ)
00389       *determ = det;
00390 
00391    if ( det == 0 ) {
00392       Error("Inv5x5","matrix is singular");
00393       return kFALSE;
00394    }
00395 
00396    const Double_t oneOverDet = 1.0/det;
00397    const Double_t mn1OverDet = - oneOverDet;
00398 
00399    pM[SM00] = mDet4_1234_1234 * oneOverDet;
00400    pM[SM01] = mDet4_1234_0234 * mn1OverDet;
00401    pM[SM02] = mDet4_1234_0134 * oneOverDet;
00402    pM[SM03] = mDet4_1234_0124 * mn1OverDet;
00403    pM[SM04] = mDet4_1234_0123 * oneOverDet;
00404 
00405    pM[SM11] = mDet4_0234_0234 * oneOverDet;
00406    pM[SM12] = mDet4_0234_0134 * mn1OverDet;
00407    pM[SM13] = mDet4_0234_0124 * oneOverDet;
00408    pM[SM14] = mDet4_0234_0123 * mn1OverDet;
00409 
00410    pM[SM22] = mDet4_0134_0134 * oneOverDet;
00411    pM[SM23] = mDet4_0134_0124 * mn1OverDet;
00412    pM[SM24] = mDet4_0134_0123 * oneOverDet;
00413 
00414    pM[SM33] = mDet4_0124_0124 * oneOverDet;
00415    pM[SM34] = mDet4_0124_0123 * mn1OverDet;
00416 
00417    pM[SM44] = mDet4_0123_0123 * oneOverDet;
00418 
00419    for (Int_t irow = 0; irow < 5; irow++) {
00420       const Int_t rowOff1 = irow*5; 
00421       for (Int_t icol = 0; icol < irow; icol++) {
00422          const Int_t rowOff2 = icol*5; 
00423          pM[rowOff1+icol] = pM[rowOff2+irow];
00424       }
00425    } 
00426 
00427    return kTRUE;
00428 }
00429 
00430 // Aij are indices for a 6x6 symmetric matrix.
00431 
00432 #define SA00 0
00433 #define SA01 1
00434 #define SA02 2
00435 #define SA03 3
00436 #define SA04 4
00437 #define SA05 5
00438 
00439 #define SA10 1
00440 #define SA11 7
00441 #define SA12 8
00442 #define SA13 9
00443 #define SA14 10
00444 #define SA15 11
00445 
00446 #define SA20 2
00447 #define SA21 8
00448 #define SA22 14
00449 #define SA23 15
00450 #define SA24 16
00451 #define SA25 17
00452 
00453 #define SA30 3
00454 #define SA31 9
00455 #define SA32 15
00456 #define SA33 21
00457 #define SA34 22
00458 #define SA35 23
00459 
00460 #define SA40 4
00461 #define SA41 10
00462 #define SA42 16
00463 #define SA43 22
00464 #define SA44 28
00465 #define SA45 29
00466 
00467 #define SA50 5
00468 #define SA51 11
00469 #define SA52 17
00470 #define SA53 23
00471 #define SA54 29
00472 #define SA55 35
00473 
00474 //______________________________________________________________________________
00475 template<class Element> 
00476 Bool_t TMatrixTSymCramerInv::Inv6x6(TMatrixTSym<Element> &m,Double_t *determ)
00477 {
00478    if (m.GetNrows() != 6 || m.GetNcols() != 6 || m.GetRowLwb() != m.GetColLwb()) {
00479       Error("Inv6x6","matrix should be square 6x6");
00480       return kFALSE;
00481    }
00482 
00483    Element *pM = m.GetMatrixArray();
00484 
00485    // Find all NECESSSARY 2x2 dets:  (39 of them)
00486 
00487    const Double_t mDet2_34_01 = pM[SA30]*pM[SA41] - pM[SA31]*pM[SA40];
00488    const Double_t mDet2_34_02 = pM[SA30]*pM[SA42] - pM[SA32]*pM[SA40];
00489    const Double_t mDet2_34_03 = pM[SA30]*pM[SA43] - pM[SA33]*pM[SA40];
00490    const Double_t mDet2_34_04 = pM[SA30]*pM[SA44] - pM[SA34]*pM[SA40];
00491    const Double_t mDet2_34_12 = pM[SA31]*pM[SA42] - pM[SA32]*pM[SA41];
00492    const Double_t mDet2_34_13 = pM[SA31]*pM[SA43] - pM[SA33]*pM[SA41];
00493    const Double_t mDet2_34_14 = pM[SA31]*pM[SA44] - pM[SA34]*pM[SA41];
00494    const Double_t mDet2_34_23 = pM[SA32]*pM[SA43] - pM[SA33]*pM[SA42];
00495    const Double_t mDet2_34_24 = pM[SA32]*pM[SA44] - pM[SA34]*pM[SA42];
00496    const Double_t mDet2_34_34 = pM[SA33]*pM[SA44] - pM[SA34]*pM[SA43];
00497    const Double_t mDet2_35_01 = pM[SA30]*pM[SA51] - pM[SA31]*pM[SA50];
00498    const Double_t mDet2_35_02 = pM[SA30]*pM[SA52] - pM[SA32]*pM[SA50];
00499    const Double_t mDet2_35_03 = pM[SA30]*pM[SA53] - pM[SA33]*pM[SA50];
00500    const Double_t mDet2_35_04 = pM[SA30]*pM[SA54] - pM[SA34]*pM[SA50];
00501    const Double_t mDet2_35_05 = pM[SA30]*pM[SA55] - pM[SA35]*pM[SA50];
00502    const Double_t mDet2_35_12 = pM[SA31]*pM[SA52] - pM[SA32]*pM[SA51];
00503    const Double_t mDet2_35_13 = pM[SA31]*pM[SA53] - pM[SA33]*pM[SA51];
00504    const Double_t mDet2_35_14 = pM[SA31]*pM[SA54] - pM[SA34]*pM[SA51];
00505    const Double_t mDet2_35_15 = pM[SA31]*pM[SA55] - pM[SA35]*pM[SA51];
00506    const Double_t mDet2_35_23 = pM[SA32]*pM[SA53] - pM[SA33]*pM[SA52];
00507    const Double_t mDet2_35_24 = pM[SA32]*pM[SA54] - pM[SA34]*pM[SA52];
00508    const Double_t mDet2_35_25 = pM[SA32]*pM[SA55] - pM[SA35]*pM[SA52];
00509    const Double_t mDet2_35_34 = pM[SA33]*pM[SA54] - pM[SA34]*pM[SA53];
00510    const Double_t mDet2_35_35 = pM[SA33]*pM[SA55] - pM[SA35]*pM[SA53];
00511    const Double_t mDet2_45_01 = pM[SA40]*pM[SA51] - pM[SA41]*pM[SA50];
00512    const Double_t mDet2_45_02 = pM[SA40]*pM[SA52] - pM[SA42]*pM[SA50];
00513    const Double_t mDet2_45_03 = pM[SA40]*pM[SA53] - pM[SA43]*pM[SA50];
00514    const Double_t mDet2_45_04 = pM[SA40]*pM[SA54] - pM[SA44]*pM[SA50];
00515    const Double_t mDet2_45_05 = pM[SA40]*pM[SA55] - pM[SA45]*pM[SA50];
00516    const Double_t mDet2_45_12 = pM[SA41]*pM[SA52] - pM[SA42]*pM[SA51];
00517    const Double_t mDet2_45_13 = pM[SA41]*pM[SA53] - pM[SA43]*pM[SA51];
00518    const Double_t mDet2_45_14 = pM[SA41]*pM[SA54] - pM[SA44]*pM[SA51];
00519    const Double_t mDet2_45_15 = pM[SA41]*pM[SA55] - pM[SA45]*pM[SA51];
00520    const Double_t mDet2_45_23 = pM[SA42]*pM[SA53] - pM[SA43]*pM[SA52];
00521    const Double_t mDet2_45_24 = pM[SA42]*pM[SA54] - pM[SA44]*pM[SA52];
00522    const Double_t mDet2_45_25 = pM[SA42]*pM[SA55] - pM[SA45]*pM[SA52];
00523    const Double_t mDet2_45_34 = pM[SA43]*pM[SA54] - pM[SA44]*pM[SA53];
00524    const Double_t mDet2_45_35 = pM[SA43]*pM[SA55] - pM[SA45]*pM[SA53];
00525    const Double_t mDet2_45_45 = pM[SA44]*pM[SA55] - pM[SA45]*pM[SA54];
00526 
00527    // Find all NECESSSARY 3x3 dets:  (65 of them)
00528 
00529    const Double_t mDet3_234_012 = pM[SA20]*mDet2_34_12 - pM[SA21]*mDet2_34_02 + pM[SA22]*mDet2_34_01;
00530    const Double_t mDet3_234_013 = pM[SA20]*mDet2_34_13 - pM[SA21]*mDet2_34_03 + pM[SA23]*mDet2_34_01;
00531    const Double_t mDet3_234_014 = pM[SA20]*mDet2_34_14 - pM[SA21]*mDet2_34_04 + pM[SA24]*mDet2_34_01;
00532    const Double_t mDet3_234_023 = pM[SA20]*mDet2_34_23 - pM[SA22]*mDet2_34_03 + pM[SA23]*mDet2_34_02;
00533    const Double_t mDet3_234_024 = pM[SA20]*mDet2_34_24 - pM[SA22]*mDet2_34_04 + pM[SA24]*mDet2_34_02;
00534    const Double_t mDet3_234_034 = pM[SA20]*mDet2_34_34 - pM[SA23]*mDet2_34_04 + pM[SA24]*mDet2_34_03;
00535    const Double_t mDet3_234_123 = pM[SA21]*mDet2_34_23 - pM[SA22]*mDet2_34_13 + pM[SA23]*mDet2_34_12;
00536    const Double_t mDet3_234_124 = pM[SA21]*mDet2_34_24 - pM[SA22]*mDet2_34_14 + pM[SA24]*mDet2_34_12;
00537    const Double_t mDet3_234_134 = pM[SA21]*mDet2_34_34 - pM[SA23]*mDet2_34_14 + pM[SA24]*mDet2_34_13;
00538    const Double_t mDet3_234_234 = pM[SA22]*mDet2_34_34 - pM[SA23]*mDet2_34_24 + pM[SA24]*mDet2_34_23;
00539    const Double_t mDet3_235_012 = pM[SA20]*mDet2_35_12 - pM[SA21]*mDet2_35_02 + pM[SA22]*mDet2_35_01;
00540    const Double_t mDet3_235_013 = pM[SA20]*mDet2_35_13 - pM[SA21]*mDet2_35_03 + pM[SA23]*mDet2_35_01;
00541    const Double_t mDet3_235_014 = pM[SA20]*mDet2_35_14 - pM[SA21]*mDet2_35_04 + pM[SA24]*mDet2_35_01;
00542    const Double_t mDet3_235_015 = pM[SA20]*mDet2_35_15 - pM[SA21]*mDet2_35_05 + pM[SA25]*mDet2_35_01;
00543    const Double_t mDet3_235_023 = pM[SA20]*mDet2_35_23 - pM[SA22]*mDet2_35_03 + pM[SA23]*mDet2_35_02;
00544    const Double_t mDet3_235_024 = pM[SA20]*mDet2_35_24 - pM[SA22]*mDet2_35_04 + pM[SA24]*mDet2_35_02;
00545    const Double_t mDet3_235_025 = pM[SA20]*mDet2_35_25 - pM[SA22]*mDet2_35_05 + pM[SA25]*mDet2_35_02;
00546    const Double_t mDet3_235_034 = pM[SA20]*mDet2_35_34 - pM[SA23]*mDet2_35_04 + pM[SA24]*mDet2_35_03;
00547    const Double_t mDet3_235_035 = pM[SA20]*mDet2_35_35 - pM[SA23]*mDet2_35_05 + pM[SA25]*mDet2_35_03;
00548    const Double_t mDet3_235_123 = pM[SA21]*mDet2_35_23 - pM[SA22]*mDet2_35_13 + pM[SA23]*mDet2_35_12;
00549    const Double_t mDet3_235_124 = pM[SA21]*mDet2_35_24 - pM[SA22]*mDet2_35_14 + pM[SA24]*mDet2_35_12;
00550    const Double_t mDet3_235_125 = pM[SA21]*mDet2_35_25 - pM[SA22]*mDet2_35_15 + pM[SA25]*mDet2_35_12;
00551    const Double_t mDet3_235_134 = pM[SA21]*mDet2_35_34 - pM[SA23]*mDet2_35_14 + pM[SA24]*mDet2_35_13;
00552    const Double_t mDet3_235_135 = pM[SA21]*mDet2_35_35 - pM[SA23]*mDet2_35_15 + pM[SA25]*mDet2_35_13;
00553    const Double_t mDet3_235_234 = pM[SA22]*mDet2_35_34 - pM[SA23]*mDet2_35_24 + pM[SA24]*mDet2_35_23;
00554    const Double_t mDet3_235_235 = pM[SA22]*mDet2_35_35 - pM[SA23]*mDet2_35_25 + pM[SA25]*mDet2_35_23;
00555    const Double_t mDet3_245_012 = pM[SA20]*mDet2_45_12 - pM[SA21]*mDet2_45_02 + pM[SA22]*mDet2_45_01;
00556    const Double_t mDet3_245_013 = pM[SA20]*mDet2_45_13 - pM[SA21]*mDet2_45_03 + pM[SA23]*mDet2_45_01;
00557    const Double_t mDet3_245_014 = pM[SA20]*mDet2_45_14 - pM[SA21]*mDet2_45_04 + pM[SA24]*mDet2_45_01;
00558    const Double_t mDet3_245_015 = pM[SA20]*mDet2_45_15 - pM[SA21]*mDet2_45_05 + pM[SA25]*mDet2_45_01;
00559    const Double_t mDet3_245_023 = pM[SA20]*mDet2_45_23 - pM[SA22]*mDet2_45_03 + pM[SA23]*mDet2_45_02;
00560    const Double_t mDet3_245_024 = pM[SA20]*mDet2_45_24 - pM[SA22]*mDet2_45_04 + pM[SA24]*mDet2_45_02;
00561    const Double_t mDet3_245_025 = pM[SA20]*mDet2_45_25 - pM[SA22]*mDet2_45_05 + pM[SA25]*mDet2_45_02;
00562    const Double_t mDet3_245_034 = pM[SA20]*mDet2_45_34 - pM[SA23]*mDet2_45_04 + pM[SA24]*mDet2_45_03;
00563    const Double_t mDet3_245_035 = pM[SA20]*mDet2_45_35 - pM[SA23]*mDet2_45_05 + pM[SA25]*mDet2_45_03;
00564    const Double_t mDet3_245_045 = pM[SA20]*mDet2_45_45 - pM[SA24]*mDet2_45_05 + pM[SA25]*mDet2_45_04;
00565    const Double_t mDet3_245_123 = pM[SA21]*mDet2_45_23 - pM[SA22]*mDet2_45_13 + pM[SA23]*mDet2_45_12;
00566    const Double_t mDet3_245_124 = pM[SA21]*mDet2_45_24 - pM[SA22]*mDet2_45_14 + pM[SA24]*mDet2_45_12;
00567    const Double_t mDet3_245_125 = pM[SA21]*mDet2_45_25 - pM[SA22]*mDet2_45_15 + pM[SA25]*mDet2_45_12;
00568    const Double_t mDet3_245_134 = pM[SA21]*mDet2_45_34 - pM[SA23]*mDet2_45_14 + pM[SA24]*mDet2_45_13;
00569    const Double_t mDet3_245_135 = pM[SA21]*mDet2_45_35 - pM[SA23]*mDet2_45_15 + pM[SA25]*mDet2_45_13;
00570    const Double_t mDet3_245_145 = pM[SA21]*mDet2_45_45 - pM[SA24]*mDet2_45_15 + pM[SA25]*mDet2_45_14;
00571    const Double_t mDet3_245_234 = pM[SA22]*mDet2_45_34 - pM[SA23]*mDet2_45_24 + pM[SA24]*mDet2_45_23;
00572    const Double_t mDet3_245_235 = pM[SA22]*mDet2_45_35 - pM[SA23]*mDet2_45_25 + pM[SA25]*mDet2_45_23;
00573    const Double_t mDet3_245_245 = pM[SA22]*mDet2_45_45 - pM[SA24]*mDet2_45_25 + pM[SA25]*mDet2_45_24;
00574    const Double_t mDet3_345_012 = pM[SA30]*mDet2_45_12 - pM[SA31]*mDet2_45_02 + pM[SA32]*mDet2_45_01;
00575    const Double_t mDet3_345_013 = pM[SA30]*mDet2_45_13 - pM[SA31]*mDet2_45_03 + pM[SA33]*mDet2_45_01;
00576    const Double_t mDet3_345_014 = pM[SA30]*mDet2_45_14 - pM[SA31]*mDet2_45_04 + pM[SA34]*mDet2_45_01;
00577    const Double_t mDet3_345_015 = pM[SA30]*mDet2_45_15 - pM[SA31]*mDet2_45_05 + pM[SA35]*mDet2_45_01;
00578    const Double_t mDet3_345_023 = pM[SA30]*mDet2_45_23 - pM[SA32]*mDet2_45_03 + pM[SA33]*mDet2_45_02;
00579    const Double_t mDet3_345_024 = pM[SA30]*mDet2_45_24 - pM[SA32]*mDet2_45_04 + pM[SA34]*mDet2_45_02;
00580    const Double_t mDet3_345_025 = pM[SA30]*mDet2_45_25 - pM[SA32]*mDet2_45_05 + pM[SA35]*mDet2_45_02;
00581    const Double_t mDet3_345_034 = pM[SA30]*mDet2_45_34 - pM[SA33]*mDet2_45_04 + pM[SA34]*mDet2_45_03;
00582    const Double_t mDet3_345_035 = pM[SA30]*mDet2_45_35 - pM[SA33]*mDet2_45_05 + pM[SA35]*mDet2_45_03;
00583    const Double_t mDet3_345_045 = pM[SA30]*mDet2_45_45 - pM[SA34]*mDet2_45_05 + pM[SA35]*mDet2_45_04;
00584    const Double_t mDet3_345_123 = pM[SA31]*mDet2_45_23 - pM[SA32]*mDet2_45_13 + pM[SA33]*mDet2_45_12;
00585    const Double_t mDet3_345_124 = pM[SA31]*mDet2_45_24 - pM[SA32]*mDet2_45_14 + pM[SA34]*mDet2_45_12;
00586    const Double_t mDet3_345_125 = pM[SA31]*mDet2_45_25 - pM[SA32]*mDet2_45_15 + pM[SA35]*mDet2_45_12;
00587    const Double_t mDet3_345_134 = pM[SA31]*mDet2_45_34 - pM[SA33]*mDet2_45_14 + pM[SA34]*mDet2_45_13;
00588    const Double_t mDet3_345_135 = pM[SA31]*mDet2_45_35 - pM[SA33]*mDet2_45_15 + pM[SA35]*mDet2_45_13;
00589    const Double_t mDet3_345_145 = pM[SA31]*mDet2_45_45 - pM[SA34]*mDet2_45_15 + pM[SA35]*mDet2_45_14;
00590    const Double_t mDet3_345_234 = pM[SA32]*mDet2_45_34 - pM[SA33]*mDet2_45_24 + pM[SA34]*mDet2_45_23;
00591    const Double_t mDet3_345_235 = pM[SA32]*mDet2_45_35 - pM[SA33]*mDet2_45_25 + pM[SA35]*mDet2_45_23;
00592    const Double_t mDet3_345_245 = pM[SA32]*mDet2_45_45 - pM[SA34]*mDet2_45_25 + pM[SA35]*mDet2_45_24;
00593    const Double_t mDet3_345_345 = pM[SA33]*mDet2_45_45 - pM[SA34]*mDet2_45_35 + pM[SA35]*mDet2_45_34;
00594 
00595    // Find all NECESSSARY 4x4 dets:  (55 of them)
00596 
00597    const Double_t mDet4_1234_0123 = pM[SA10]*mDet3_234_123 - pM[SA11]*mDet3_234_023
00598                                   + pM[SA12]*mDet3_234_013 - pM[SA13]*mDet3_234_012;
00599    const Double_t mDet4_1234_0124 = pM[SA10]*mDet3_234_124 - pM[SA11]*mDet3_234_024
00600                                   + pM[SA12]*mDet3_234_014 - pM[SA14]*mDet3_234_012;
00601    const Double_t mDet4_1234_0134 = pM[SA10]*mDet3_234_134 - pM[SA11]*mDet3_234_034
00602                                   + pM[SA13]*mDet3_234_014 - pM[SA14]*mDet3_234_013;
00603    const Double_t mDet4_1234_0234 = pM[SA10]*mDet3_234_234 - pM[SA12]*mDet3_234_034
00604                                   + pM[SA13]*mDet3_234_024 - pM[SA14]*mDet3_234_023;
00605    const Double_t mDet4_1234_1234 = pM[SA11]*mDet3_234_234 - pM[SA12]*mDet3_234_134
00606                                   + pM[SA13]*mDet3_234_124 - pM[SA14]*mDet3_234_123;
00607    const Double_t mDet4_1235_0123 = pM[SA10]*mDet3_235_123 - pM[SA11]*mDet3_235_023
00608                                   + pM[SA12]*mDet3_235_013 - pM[SA13]*mDet3_235_012;
00609    const Double_t mDet4_1235_0124 = pM[SA10]*mDet3_235_124 - pM[SA11]*mDet3_235_024
00610                                   + pM[SA12]*mDet3_235_014 - pM[SA14]*mDet3_235_012;
00611    const Double_t mDet4_1235_0125 = pM[SA10]*mDet3_235_125 - pM[SA11]*mDet3_235_025
00612                                   + pM[SA12]*mDet3_235_015 - pM[SA15]*mDet3_235_012;
00613    const Double_t mDet4_1235_0134 = pM[SA10]*mDet3_235_134 - pM[SA11]*mDet3_235_034
00614                                   + pM[SA13]*mDet3_235_014 - pM[SA14]*mDet3_235_013;
00615    const Double_t mDet4_1235_0135 = pM[SA10]*mDet3_235_135 - pM[SA11]*mDet3_235_035
00616                                   + pM[SA13]*mDet3_235_015 - pM[SA15]*mDet3_235_013;
00617    const Double_t mDet4_1235_0234 = pM[SA10]*mDet3_235_234 - pM[SA12]*mDet3_235_034
00618                                   + pM[SA13]*mDet3_235_024 - pM[SA14]*mDet3_235_023;
00619    const Double_t mDet4_1235_0235 = pM[SA10]*mDet3_235_235 - pM[SA12]*mDet3_235_035
00620                                   + pM[SA13]*mDet3_235_025 - pM[SA15]*mDet3_235_023;
00621    const Double_t mDet4_1235_1234 = pM[SA11]*mDet3_235_234 - pM[SA12]*mDet3_235_134
00622                                   + pM[SA13]*mDet3_235_124 - pM[SA14]*mDet3_235_123;
00623    const Double_t mDet4_1235_1235 = pM[SA11]*mDet3_235_235 - pM[SA12]*mDet3_235_135
00624                                   + pM[SA13]*mDet3_235_125 - pM[SA15]*mDet3_235_123;
00625    const Double_t mDet4_1245_0123 = pM[SA10]*mDet3_245_123 - pM[SA11]*mDet3_245_023
00626                                   + pM[SA12]*mDet3_245_013 - pM[SA13]*mDet3_245_012;
00627    const Double_t mDet4_1245_0124 = pM[SA10]*mDet3_245_124 - pM[SA11]*mDet3_245_024
00628                                   + pM[SA12]*mDet3_245_014 - pM[SA14]*mDet3_245_012;
00629    const Double_t mDet4_1245_0125 = pM[SA10]*mDet3_245_125 - pM[SA11]*mDet3_245_025
00630                                   + pM[SA12]*mDet3_245_015 - pM[SA15]*mDet3_245_012;
00631    const Double_t mDet4_1245_0134 = pM[SA10]*mDet3_245_134 - pM[SA11]*mDet3_245_034
00632                                   + pM[SA13]*mDet3_245_014 - pM[SA14]*mDet3_245_013;
00633    const Double_t mDet4_1245_0135 = pM[SA10]*mDet3_245_135 - pM[SA11]*mDet3_245_035
00634                                   + pM[SA13]*mDet3_245_015 - pM[SA15]*mDet3_245_013;
00635    const Double_t mDet4_1245_0145 = pM[SA10]*mDet3_245_145 - pM[SA11]*mDet3_245_045
00636                                   + pM[SA14]*mDet3_245_015 - pM[SA15]*mDet3_245_014;
00637    const Double_t mDet4_1245_0234 = pM[SA10]*mDet3_245_234 - pM[SA12]*mDet3_245_034
00638                                   + pM[SA13]*mDet3_245_024 - pM[SA14]*mDet3_245_023;
00639    const Double_t mDet4_1245_0235 = pM[SA10]*mDet3_245_235 - pM[SA12]*mDet3_245_035
00640                                   + pM[SA13]*mDet3_245_025 - pM[SA15]*mDet3_245_023;
00641    const Double_t mDet4_1245_0245 = pM[SA10]*mDet3_245_245 - pM[SA12]*mDet3_245_045
00642                                   + pM[SA14]*mDet3_245_025 - pM[SA15]*mDet3_245_024;
00643    const Double_t mDet4_1245_1234 = pM[SA11]*mDet3_245_234 - pM[SA12]*mDet3_245_134
00644                                   + pM[SA13]*mDet3_245_124 - pM[SA14]*mDet3_245_123;
00645    const Double_t mDet4_1245_1235 = pM[SA11]*mDet3_245_235 - pM[SA12]*mDet3_245_135
00646                                   + pM[SA13]*mDet3_245_125 - pM[SA15]*mDet3_245_123;
00647    const Double_t mDet4_1245_1245 = pM[SA11]*mDet3_245_245 - pM[SA12]*mDet3_245_145
00648                                   + pM[SA14]*mDet3_245_125 - pM[SA15]*mDet3_245_124;
00649    const Double_t mDet4_1345_0123 = pM[SA10]*mDet3_345_123 - pM[SA11]*mDet3_345_023
00650                                   + pM[SA12]*mDet3_345_013 - pM[SA13]*mDet3_345_012;
00651    const Double_t mDet4_1345_0124 = pM[SA10]*mDet3_345_124 - pM[SA11]*mDet3_345_024
00652                                   + pM[SA12]*mDet3_345_014 - pM[SA14]*mDet3_345_012;
00653    const Double_t mDet4_1345_0125 = pM[SA10]*mDet3_345_125 - pM[SA11]*mDet3_345_025
00654                                   + pM[SA12]*mDet3_345_015 - pM[SA15]*mDet3_345_012;
00655    const Double_t mDet4_1345_0134 = pM[SA10]*mDet3_345_134 - pM[SA11]*mDet3_345_034
00656                                   + pM[SA13]*mDet3_345_014 - pM[SA14]*mDet3_345_013;
00657    const Double_t mDet4_1345_0135 = pM[SA10]*mDet3_345_135 - pM[SA11]*mDet3_345_035
00658                                   + pM[SA13]*mDet3_345_015 - pM[SA15]*mDet3_345_013;
00659    const Double_t mDet4_1345_0145 = pM[SA10]*mDet3_345_145 - pM[SA11]*mDet3_345_045
00660                                   + pM[SA14]*mDet3_345_015 - pM[SA15]*mDet3_345_014;
00661    const Double_t mDet4_1345_0234 = pM[SA10]*mDet3_345_234 - pM[SA12]*mDet3_345_034
00662                                   + pM[SA13]*mDet3_345_024 - pM[SA14]*mDet3_345_023;
00663    const Double_t mDet4_1345_0235 = pM[SA10]*mDet3_345_235 - pM[SA12]*mDet3_345_035
00664                                   + pM[SA13]*mDet3_345_025 - pM[SA15]*mDet3_345_023;
00665    const Double_t mDet4_1345_0245 = pM[SA10]*mDet3_345_245 - pM[SA12]*mDet3_345_045
00666                                   + pM[SA14]*mDet3_345_025 - pM[SA15]*mDet3_345_024;
00667    const Double_t mDet4_1345_0345 = pM[SA10]*mDet3_345_345 - pM[SA13]*mDet3_345_045
00668                                   + pM[SA14]*mDet3_345_035 - pM[SA15]*mDet3_345_034;
00669    const Double_t mDet4_1345_1234 = pM[SA11]*mDet3_345_234 - pM[SA12]*mDet3_345_134
00670                                   + pM[SA13]*mDet3_345_124 - pM[SA14]*mDet3_345_123;
00671    const Double_t mDet4_1345_1235 = pM[SA11]*mDet3_345_235 - pM[SA12]*mDet3_345_135
00672                                   + pM[SA13]*mDet3_345_125 - pM[SA15]*mDet3_345_123;
00673    const Double_t mDet4_1345_1245 = pM[SA11]*mDet3_345_245 - pM[SA12]*mDet3_345_145
00674                                   + pM[SA14]*mDet3_345_125 - pM[SA15]*mDet3_345_124;
00675    const Double_t mDet4_1345_1345 = pM[SA11]*mDet3_345_345 - pM[SA13]*mDet3_345_145
00676                                   + pM[SA14]*mDet3_345_135 - pM[SA15]*mDet3_345_134;
00677    const Double_t mDet4_2345_0123 = pM[SA20]*mDet3_345_123 - pM[SA21]*mDet3_345_023
00678                                   + pM[SA22]*mDet3_345_013 - pM[SA23]*mDet3_345_012;
00679    const Double_t mDet4_2345_0124 = pM[SA20]*mDet3_345_124 - pM[SA21]*mDet3_345_024
00680                                   + pM[SA22]*mDet3_345_014 - pM[SA24]*mDet3_345_012;
00681    const Double_t mDet4_2345_0125 = pM[SA20]*mDet3_345_125 - pM[SA21]*mDet3_345_025
00682                                   + pM[SA22]*mDet3_345_015 - pM[SA25]*mDet3_345_012;
00683    const Double_t mDet4_2345_0134 = pM[SA20]*mDet3_345_134 - pM[SA21]*mDet3_345_034
00684                                   + pM[SA23]*mDet3_345_014 - pM[SA24]*mDet3_345_013;
00685    const Double_t mDet4_2345_0135 = pM[SA20]*mDet3_345_135 - pM[SA21]*mDet3_345_035
00686                                   + pM[SA23]*mDet3_345_015 - pM[SA25]*mDet3_345_013;
00687    const Double_t mDet4_2345_0145 = pM[SA20]*mDet3_345_145 - pM[SA21]*mDet3_345_045
00688                                   + pM[SA24]*mDet3_345_015 - pM[SA25]*mDet3_345_014;
00689    const Double_t mDet4_2345_0234 = pM[SA20]*mDet3_345_234 - pM[SA22]*mDet3_345_034
00690                                   + pM[SA23]*mDet3_345_024 - pM[SA24]*mDet3_345_023;
00691    const Double_t mDet4_2345_0235 = pM[SA20]*mDet3_345_235 - pM[SA22]*mDet3_345_035
00692                                   + pM[SA23]*mDet3_345_025 - pM[SA25]*mDet3_345_023;
00693    const Double_t mDet4_2345_0245 = pM[SA20]*mDet3_345_245 - pM[SA22]*mDet3_345_045
00694                                   + pM[SA24]*mDet3_345_025 - pM[SA25]*mDet3_345_024;
00695    const Double_t mDet4_2345_0345 = pM[SA20]*mDet3_345_345 - pM[SA23]*mDet3_345_045
00696                                   + pM[SA24]*mDet3_345_035 - pM[SA25]*mDet3_345_034;
00697    const Double_t mDet4_2345_1234 = pM[SA21]*mDet3_345_234 - pM[SA22]*mDet3_345_134
00698                                   + pM[SA23]*mDet3_345_124 - pM[SA24]*mDet3_345_123;
00699    const Double_t mDet4_2345_1235 = pM[SA21]*mDet3_345_235 - pM[SA22]*mDet3_345_135
00700                                   + pM[SA23]*mDet3_345_125 - pM[SA25]*mDet3_345_123;
00701    const Double_t mDet4_2345_1245 = pM[SA21]*mDet3_345_245 - pM[SA22]*mDet3_345_145
00702                                   + pM[SA24]*mDet3_345_125 - pM[SA25]*mDet3_345_124;
00703    const Double_t mDet4_2345_1345 = pM[SA21]*mDet3_345_345 - pM[SA23]*mDet3_345_145
00704                                   + pM[SA24]*mDet3_345_135 - pM[SA25]*mDet3_345_134;
00705    const Double_t mDet4_2345_2345 = pM[SA22]*mDet3_345_345 - pM[SA23]*mDet3_345_245
00706                                   + pM[SA24]*mDet3_345_235 - pM[SA25]*mDet3_345_234;
00707 
00708    // Find all NECESSSARY 5x5 dets:  (19 of them)
00709 
00710    const Double_t mDet5_01234_01234 = pM[SA00]*mDet4_1234_1234 - pM[SA01]*mDet4_1234_0234
00711                                     + pM[SA02]*mDet4_1234_0134 - pM[SA03]*mDet4_1234_0124 + pM[SA04]*mDet4_1234_0123;
00712    const Double_t mDet5_01235_01234 = pM[SA00]*mDet4_1235_1234 - pM[SA01]*mDet4_1235_0234
00713                                     + pM[SA02]*mDet4_1235_0134 - pM[SA03]*mDet4_1235_0124 + pM[SA04]*mDet4_1235_0123;
00714    const Double_t mDet5_01235_01235 = pM[SA00]*mDet4_1235_1235 - pM[SA01]*mDet4_1235_0235
00715                                     + pM[SA02]*mDet4_1235_0135 - pM[SA03]*mDet4_1235_0125 + pM[SA05]*mDet4_1235_0123;
00716    const Double_t mDet5_01245_01234 = pM[SA00]*mDet4_1245_1234 - pM[SA01]*mDet4_1245_0234
00717                                     + pM[SA02]*mDet4_1245_0134 - pM[SA03]*mDet4_1245_0124 + pM[SA04]*mDet4_1245_0123;
00718    const Double_t mDet5_01245_01235 = pM[SA00]*mDet4_1245_1235 - pM[SA01]*mDet4_1245_0235
00719                                     + pM[SA02]*mDet4_1245_0135 - pM[SA03]*mDet4_1245_0125 + pM[SA05]*mDet4_1245_0123;
00720    const Double_t mDet5_01245_01245 = pM[SA00]*mDet4_1245_1245 - pM[SA01]*mDet4_1245_0245
00721                                     + pM[SA02]*mDet4_1245_0145 - pM[SA04]*mDet4_1245_0125 + pM[SA05]*mDet4_1245_0124;
00722    const Double_t mDet5_01345_01234 = pM[SA00]*mDet4_1345_1234 - pM[SA01]*mDet4_1345_0234
00723                                     + pM[SA02]*mDet4_1345_0134 - pM[SA03]*mDet4_1345_0124 + pM[SA04]*mDet4_1345_0123;
00724    const Double_t mDet5_01345_01235 = pM[SA00]*mDet4_1345_1235 - pM[SA01]*mDet4_1345_0235
00725                                     + pM[SA02]*mDet4_1345_0135 - pM[SA03]*mDet4_1345_0125 + pM[SA05]*mDet4_1345_0123;
00726    const Double_t mDet5_01345_01245 = pM[SA00]*mDet4_1345_1245 - pM[SA01]*mDet4_1345_0245
00727                                     + pM[SA02]*mDet4_1345_0145 - pM[SA04]*mDet4_1345_0125 + pM[SA05]*mDet4_1345_0124;
00728    const Double_t mDet5_01345_01345 = pM[SA00]*mDet4_1345_1345 - pM[SA01]*mDet4_1345_0345
00729                                     + pM[SA03]*mDet4_1345_0145 - pM[SA04]*mDet4_1345_0135 + pM[SA05]*mDet4_1345_0134;
00730    const Double_t mDet5_02345_01234 = pM[SA00]*mDet4_2345_1234 - pM[SA01]*mDet4_2345_0234
00731                                     + pM[SA02]*mDet4_2345_0134 - pM[SA03]*mDet4_2345_0124 + pM[SA04]*mDet4_2345_0123;
00732    const Double_t mDet5_02345_01235 = pM[SA00]*mDet4_2345_1235 - pM[SA01]*mDet4_2345_0235
00733                                     + pM[SA02]*mDet4_2345_0135 - pM[SA03]*mDet4_2345_0125 + pM[SA05]*mDet4_2345_0123;
00734    const Double_t mDet5_02345_01245 = pM[SA00]*mDet4_2345_1245 - pM[SA01]*mDet4_2345_0245
00735                                     + pM[SA02]*mDet4_2345_0145 - pM[SA04]*mDet4_2345_0125 + pM[SA05]*mDet4_2345_0124;
00736    const Double_t mDet5_02345_01345 = pM[SA00]*mDet4_2345_1345 - pM[SA01]*mDet4_2345_0345
00737                                     + pM[SA03]*mDet4_2345_0145 - pM[SA04]*mDet4_2345_0135 + pM[SA05]*mDet4_2345_0134;
00738    const Double_t mDet5_02345_02345 = pM[SA00]*mDet4_2345_2345 - pM[SA02]*mDet4_2345_0345
00739                                     + pM[SA03]*mDet4_2345_0245 - pM[SA04]*mDet4_2345_0235 + pM[SA05]*mDet4_2345_0234;
00740    const Double_t mDet5_12345_01234 = pM[SA10]*mDet4_2345_1234 - pM[SA11]*mDet4_2345_0234
00741                                     + pM[SA12]*mDet4_2345_0134 - pM[SA13]*mDet4_2345_0124 + pM[SA14]*mDet4_2345_0123;
00742    const Double_t mDet5_12345_01235 = pM[SA10]*mDet4_2345_1235 - pM[SA11]*mDet4_2345_0235
00743                                     + pM[SA12]*mDet4_2345_0135 - pM[SA13]*mDet4_2345_0125 + pM[SA15]*mDet4_2345_0123;
00744    const Double_t mDet5_12345_01245 = pM[SA10]*mDet4_2345_1245 - pM[SA11]*mDet4_2345_0245
00745                                     + pM[SA12]*mDet4_2345_0145 - pM[SA14]*mDet4_2345_0125 + pM[SA15]*mDet4_2345_0124;
00746    const Double_t mDet5_12345_01345 = pM[SA10]*mDet4_2345_1345 - pM[SA11]*mDet4_2345_0345
00747                                     + pM[SA13]*mDet4_2345_0145 - pM[SA14]*mDet4_2345_0135 + pM[SA15]*mDet4_2345_0134;
00748    const Double_t mDet5_12345_02345 = pM[SA10]*mDet4_2345_2345 - pM[SA12]*mDet4_2345_0345
00749                                     + pM[SA13]*mDet4_2345_0245 - pM[SA14]*mDet4_2345_0235 + pM[SA15]*mDet4_2345_0234;
00750    const Double_t mDet5_12345_12345 = pM[SA11]*mDet4_2345_2345 - pM[SA12]*mDet4_2345_1345
00751                                     + pM[SA13]*mDet4_2345_1245 - pM[SA14]*mDet4_2345_1235 + pM[SA15]*mDet4_2345_1234;
00752 
00753    // Find the determinant 
00754 
00755    const Double_t det = pM[SA00]*mDet5_12345_12345 - pM[SA01]*mDet5_12345_02345 + pM[SA02]*mDet5_12345_01345
00756                       - pM[SA03]*mDet5_12345_01245 + pM[SA04]*mDet5_12345_01235 - pM[SA05]*mDet5_12345_01234;
00757 
00758    if (determ)
00759       *determ = det;
00760 
00761    if ( det == 0 ) {
00762       Error("Inv6x6","matrix is singular");
00763       return kFALSE;
00764    }
00765 
00766    const Double_t oneOverDet = 1.0/det;
00767    const Double_t mn1OverDet = - oneOverDet;
00768 
00769    pM[SA00] =  mDet5_12345_12345*oneOverDet;
00770    pM[SA01] =  mDet5_12345_02345*mn1OverDet;
00771    pM[SA02] =  mDet5_12345_01345*oneOverDet;
00772    pM[SA03] =  mDet5_12345_01245*mn1OverDet;
00773    pM[SA04] =  mDet5_12345_01235*oneOverDet;
00774    pM[SA05] =  mDet5_12345_01234*mn1OverDet;
00775 
00776    pM[SA11] =  mDet5_02345_02345*oneOverDet;
00777    pM[SA12] =  mDet5_02345_01345*mn1OverDet;
00778    pM[SA13] =  mDet5_02345_01245*oneOverDet;
00779    pM[SA14] =  mDet5_02345_01235*mn1OverDet;
00780    pM[SA15] =  mDet5_02345_01234*oneOverDet;
00781 
00782    pM[SA22] =  mDet5_01345_01345*oneOverDet;
00783    pM[SA23] =  mDet5_01345_01245*mn1OverDet;
00784    pM[SA24] =  mDet5_01345_01235*oneOverDet;
00785    pM[SA25] =  mDet5_01345_01234*mn1OverDet;
00786 
00787    pM[SA33] =  mDet5_01245_01245*oneOverDet;
00788    pM[SA34] =  mDet5_01245_01235*mn1OverDet;
00789    pM[SA35] =  mDet5_01245_01234*oneOverDet;
00790 
00791    pM[SA44] =  mDet5_01235_01235*oneOverDet;
00792    pM[SA45] =  mDet5_01235_01234*mn1OverDet;
00793 
00794    pM[SA55] =  mDet5_01234_01234*oneOverDet;
00795 
00796    for (Int_t irow = 0; irow < 6; irow++) {
00797       const Int_t rowOff1 = irow*6; 
00798       for (Int_t icol = 0; icol < irow; icol++) {
00799          const Int_t rowOff2 = icol*6; 
00800          pM[rowOff1+icol] = pM[rowOff2+irow];
00801       }
00802    } 
00803 
00804    return kTRUE;
00805 }
00806 
00807 #ifndef ROOT_TMatrixFSymfwd
00808 #include "TMatrixFSymfwd.h"
00809 #endif
00810 
00811 template Bool_t TMatrixTSymCramerInv::Inv2x2<Float_t>(TMatrixFSym&,Double_t*);
00812 template Bool_t TMatrixTSymCramerInv::Inv3x3<Float_t>(TMatrixFSym&,Double_t*);
00813 template Bool_t TMatrixTSymCramerInv::Inv4x4<Float_t>(TMatrixFSym&,Double_t*);
00814 template Bool_t TMatrixTSymCramerInv::Inv5x5<Float_t>(TMatrixFSym&,Double_t*);
00815 template Bool_t TMatrixTSymCramerInv::Inv6x6<Float_t>(TMatrixFSym&,Double_t*);
00816 
00817 #ifndef ROOT_TMatrixDSymfwd
00818 #include "TMatrixDSymfwd.h"
00819 #endif
00820 
00821 template Bool_t TMatrixTSymCramerInv::Inv2x2<Double_t>(TMatrixDSym&,Double_t*);
00822 template Bool_t TMatrixTSymCramerInv::Inv3x3<Double_t>(TMatrixDSym&,Double_t*);
00823 template Bool_t TMatrixTSymCramerInv::Inv4x4<Double_t>(TMatrixDSym&,Double_t*);
00824 template Bool_t TMatrixTSymCramerInv::Inv5x5<Double_t>(TMatrixDSym&,Double_t*);
00825 template Bool_t TMatrixTSymCramerInv::Inv6x6<Double_t>(TMatrixDSym&,Double_t*);

Generated on Tue Jul 5 14:36:36 2011 for ROOT_528-00b_version by  doxygen 1.5.1