// @(#)root/base:$Id$ // Authors: Fons Rademakers, Eddy Offermann Jan 2004 /************************************************************************* * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. * * All rights reserved. * * * * For the licensing terms see $ROOTSYS/LICENSE. * * For the list of contributors see $ROOTSYS/README/CREDITS. * *************************************************************************/ ////////////////////////////////////////////////////////////////////////// // // // TMatrixTCramerInv // // // // Encapsulate templates of Cramer Inversion routines. // // // // The 4x4, 5x5 and 6x6 are adapted from routines written by // // Mark Fischler and Steven Haywood as part of the CLHEP package // // // // Although for sizes <= 6x6 the Cramer Inversion has a gain in speed // // compared to factorization schemes (like LU) , one pays a price in // // accuracy . // // // // For Example: // // H * H^-1 = U, where H is a 5x5 Hilbert matrix // // U is a 5x5 Unity matrix // // // // LU : |U_jk| < 10e-13 for j!=k // // Cramer: |U_jk| < 10e-7 for j!=k // // // // however Cramer algorithm is about 10 (!) times faster // ////////////////////////////////////////////////////////////////////////// #include "TMatrixTCramerInv.h" #if !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD) NamespaceImp(TMatrixTCramerInv); #endif //////////////////////////////////////////////////////////////////////////////// template Bool_t TMatrixTCramerInv::Inv2x2(TMatrixT &m,Double_t *determ) { if (m.GetNrows() != 2 || m.GetNcols() != 2 || m.GetRowLwb() != m.GetColLwb()) { Error("Inv2x2","matrix should be square 2x2"); return kFALSE; } Element *pM = m.GetMatrixArray(); const Double_t det = pM[0] * pM[3] - pM[2] * pM[1]; if (determ) *determ = det; const Double_t s = 1./det; if ( det == 0 ) { Error("Inv2x2","matrix is singular"); return kFALSE; } const Double_t tmp = s*pM[3]; pM[1] *= -s; pM[2] *= -s; pM[3] = s*pM[0]; pM[0] = tmp; return kTRUE; } //////////////////////////////////////////////////////////////////////////////// template Bool_t TMatrixTCramerInv::Inv3x3(TMatrixT &m,Double_t *determ) { if (m.GetNrows() != 3 || m.GetNcols() != 3 || m.GetRowLwb() != m.GetColLwb()) { Error("Inv3x3","matrix should be square 3x3"); return kFALSE; } Element *pM = m.GetMatrixArray(); const Double_t c00 = pM[4] * pM[8] - pM[5] * pM[7]; const Double_t c01 = pM[5] * pM[6] - pM[3] * pM[8]; const Double_t c02 = pM[3] * pM[7] - pM[4] * pM[6]; const Double_t c10 = pM[7] * pM[2] - pM[8] * pM[1]; const Double_t c11 = pM[8] * pM[0] - pM[6] * pM[2]; const Double_t c12 = pM[6] * pM[1] - pM[7] * pM[0]; const Double_t c20 = pM[1] * pM[5] - pM[2] * pM[4]; const Double_t c21 = pM[2] * pM[3] - pM[0] * pM[5]; const Double_t c22 = pM[0] * pM[4] - pM[1] * pM[3]; const Double_t t0 = TMath::Abs(pM[0]); const Double_t t1 = TMath::Abs(pM[3]); const Double_t t2 = TMath::Abs(pM[6]); Double_t det; Double_t tmp; if (t0 >= t1) { if (t2 >= t0) { tmp = pM[6]; det = c12*c01-c11*c02; } else { tmp = pM[0]; det = c11*c22-c12*c21; } } else if (t2 >= t1) { tmp = pM[6]; det = c12*c01-c11*c02; } else { tmp = pM[3]; det = c02*c21-c01*c22; } if ( det == 0 || tmp == 0) { Error("Inv3x3","matrix is singular"); return kFALSE; } const Double_t s = tmp/det; if (determ) *determ = 1./s; pM[0] = s*c00; pM[1] = s*c10; pM[2] = s*c20; pM[3] = s*c01; pM[4] = s*c11; pM[5] = s*c21; pM[6] = s*c02; pM[7] = s*c12; pM[8] = s*c22; return kTRUE; } // GFij are indices for a 4x4 matrix. #define GF00 0 #define GF01 1 #define GF02 2 #define GF03 3 #define GF10 4 #define GF11 5 #define GF12 6 #define GF13 7 #define GF20 8 #define GF21 9 #define GF22 10 #define GF23 11 #define GF30 12 #define GF31 13 #define GF32 14 #define GF33 15 //////////////////////////////////////////////////////////////////////////////// template Bool_t TMatrixTCramerInv::Inv4x4(TMatrixT &m,Double_t *determ) { if (m.GetNrows() != 4 || m.GetNcols() != 4 || m.GetRowLwb() != m.GetColLwb()) { Error("Inv4x4","matrix should be square 4x4"); return kFALSE; } Element *pM = m.GetMatrixArray(); // Find all NECESSARY 2x2 dets: (18 of them) const Double_t det2_12_01 = pM[GF10]*pM[GF21] - pM[GF11]*pM[GF20]; const Double_t det2_12_02 = pM[GF10]*pM[GF22] - pM[GF12]*pM[GF20]; const Double_t det2_12_03 = pM[GF10]*pM[GF23] - pM[GF13]*pM[GF20]; const Double_t det2_12_13 = pM[GF11]*pM[GF23] - pM[GF13]*pM[GF21]; const Double_t det2_12_23 = pM[GF12]*pM[GF23] - pM[GF13]*pM[GF22]; const Double_t det2_12_12 = pM[GF11]*pM[GF22] - pM[GF12]*pM[GF21]; const Double_t det2_13_01 = pM[GF10]*pM[GF31] - pM[GF11]*pM[GF30]; const Double_t det2_13_02 = pM[GF10]*pM[GF32] - pM[GF12]*pM[GF30]; const Double_t det2_13_03 = pM[GF10]*pM[GF33] - pM[GF13]*pM[GF30]; const Double_t det2_13_12 = pM[GF11]*pM[GF32] - pM[GF12]*pM[GF31]; const Double_t det2_13_13 = pM[GF11]*pM[GF33] - pM[GF13]*pM[GF31]; const Double_t det2_13_23 = pM[GF12]*pM[GF33] - pM[GF13]*pM[GF32]; const Double_t det2_23_01 = pM[GF20]*pM[GF31] - pM[GF21]*pM[GF30]; const Double_t det2_23_02 = pM[GF20]*pM[GF32] - pM[GF22]*pM[GF30]; const Double_t det2_23_03 = pM[GF20]*pM[GF33] - pM[GF23]*pM[GF30]; const Double_t det2_23_12 = pM[GF21]*pM[GF32] - pM[GF22]*pM[GF31]; const Double_t det2_23_13 = pM[GF21]*pM[GF33] - pM[GF23]*pM[GF31]; const Double_t det2_23_23 = pM[GF22]*pM[GF33] - pM[GF23]*pM[GF32]; // Find all NECESSARY 3x3 dets: (16 of them) const Double_t det3_012_012 = pM[GF00]*det2_12_12 - pM[GF01]*det2_12_02 + pM[GF02]*det2_12_01; const Double_t det3_012_013 = pM[GF00]*det2_12_13 - pM[GF01]*det2_12_03 + pM[GF03]*det2_12_01; const Double_t det3_012_023 = pM[GF00]*det2_12_23 - pM[GF02]*det2_12_03 + pM[GF03]*det2_12_02; const Double_t det3_012_123 = pM[GF01]*det2_12_23 - pM[GF02]*det2_12_13 + pM[GF03]*det2_12_12; const Double_t det3_013_012 = pM[GF00]*det2_13_12 - pM[GF01]*det2_13_02 + pM[GF02]*det2_13_01; const Double_t det3_013_013 = pM[GF00]*det2_13_13 - pM[GF01]*det2_13_03 + pM[GF03]*det2_13_01; const Double_t det3_013_023 = pM[GF00]*det2_13_23 - pM[GF02]*det2_13_03 + pM[GF03]*det2_13_02; const Double_t det3_013_123 = pM[GF01]*det2_13_23 - pM[GF02]*det2_13_13 + pM[GF03]*det2_13_12; const Double_t det3_023_012 = pM[GF00]*det2_23_12 - pM[GF01]*det2_23_02 + pM[GF02]*det2_23_01; const Double_t det3_023_013 = pM[GF00]*det2_23_13 - pM[GF01]*det2_23_03 + pM[GF03]*det2_23_01; const Double_t det3_023_023 = pM[GF00]*det2_23_23 - pM[GF02]*det2_23_03 + pM[GF03]*det2_23_02; const Double_t det3_023_123 = pM[GF01]*det2_23_23 - pM[GF02]*det2_23_13 + pM[GF03]*det2_23_12; const Double_t det3_123_012 = pM[GF10]*det2_23_12 - pM[GF11]*det2_23_02 + pM[GF12]*det2_23_01; const Double_t det3_123_013 = pM[GF10]*det2_23_13 - pM[GF11]*det2_23_03 + pM[GF13]*det2_23_01; const Double_t det3_123_023 = pM[GF10]*det2_23_23 - pM[GF12]*det2_23_03 + pM[GF13]*det2_23_02; const Double_t det3_123_123 = pM[GF11]*det2_23_23 - pM[GF12]*det2_23_13 + pM[GF13]*det2_23_12; // Find the 4x4 det: const Double_t det = pM[GF00]*det3_123_123 - pM[GF01]*det3_123_023 + pM[GF02]*det3_123_013 - pM[GF03]*det3_123_012; if (determ) *determ = det; if ( det == 0 ) { Error("Inv4x4","matrix is singular"); return kFALSE; } const Double_t oneOverDet = 1.0/det; const Double_t mn1OverDet = - oneOverDet; pM[GF00] = det3_123_123 * oneOverDet; pM[GF01] = det3_023_123 * mn1OverDet; pM[GF02] = det3_013_123 * oneOverDet; pM[GF03] = det3_012_123 * mn1OverDet; pM[GF10] = det3_123_023 * mn1OverDet; pM[GF11] = det3_023_023 * oneOverDet; pM[GF12] = det3_013_023 * mn1OverDet; pM[GF13] = det3_012_023 * oneOverDet; pM[GF20] = det3_123_013 * oneOverDet; pM[GF21] = det3_023_013 * mn1OverDet; pM[GF22] = det3_013_013 * oneOverDet; pM[GF23] = det3_012_013 * mn1OverDet; pM[GF30] = det3_123_012 * mn1OverDet; pM[GF31] = det3_023_012 * oneOverDet; pM[GF32] = det3_013_012 * mn1OverDet; pM[GF33] = det3_012_012 * oneOverDet; return kTRUE; } // GMij are indices for a 5x5 matrix. #define GM00 0 #define GM01 1 #define GM02 2 #define GM03 3 #define GM04 4 #define GM10 5 #define GM11 6 #define GM12 7 #define GM13 8 #define GM14 9 #define GM20 10 #define GM21 11 #define GM22 12 #define GM23 13 #define GM24 14 #define GM30 15 #define GM31 16 #define GM32 17 #define GM33 18 #define GM34 19 #define GM40 20 #define GM41 21 #define GM42 22 #define GM43 23 #define GM44 24 //////////////////////////////////////////////////////////////////////////////// template Bool_t TMatrixTCramerInv::Inv5x5(TMatrixT &m,Double_t *determ) { if (m.GetNrows() != 5 || m.GetNcols() != 5 || m.GetRowLwb() != m.GetColLwb()) { Error("Inv5x5","matrix should be square 5x5"); return kFALSE; } Element *pM = m.GetMatrixArray(); // Find all NECESSARY 2x2 dets: (30 of them) const Double_t det2_23_01 = pM[GM20]*pM[GM31] - pM[GM21]*pM[GM30]; const Double_t det2_23_02 = pM[GM20]*pM[GM32] - pM[GM22]*pM[GM30]; const Double_t det2_23_03 = pM[GM20]*pM[GM33] - pM[GM23]*pM[GM30]; const Double_t det2_23_04 = pM[GM20]*pM[GM34] - pM[GM24]*pM[GM30]; const Double_t det2_23_12 = pM[GM21]*pM[GM32] - pM[GM22]*pM[GM31]; const Double_t det2_23_13 = pM[GM21]*pM[GM33] - pM[GM23]*pM[GM31]; const Double_t det2_23_14 = pM[GM21]*pM[GM34] - pM[GM24]*pM[GM31]; const Double_t det2_23_23 = pM[GM22]*pM[GM33] - pM[GM23]*pM[GM32]; const Double_t det2_23_24 = pM[GM22]*pM[GM34] - pM[GM24]*pM[GM32]; const Double_t det2_23_34 = pM[GM23]*pM[GM34] - pM[GM24]*pM[GM33]; const Double_t det2_24_01 = pM[GM20]*pM[GM41] - pM[GM21]*pM[GM40]; const Double_t det2_24_02 = pM[GM20]*pM[GM42] - pM[GM22]*pM[GM40]; const Double_t det2_24_03 = pM[GM20]*pM[GM43] - pM[GM23]*pM[GM40]; const Double_t det2_24_04 = pM[GM20]*pM[GM44] - pM[GM24]*pM[GM40]; const Double_t det2_24_12 = pM[GM21]*pM[GM42] - pM[GM22]*pM[GM41]; const Double_t det2_24_13 = pM[GM21]*pM[GM43] - pM[GM23]*pM[GM41]; const Double_t det2_24_14 = pM[GM21]*pM[GM44] - pM[GM24]*pM[GM41]; const Double_t det2_24_23 = pM[GM22]*pM[GM43] - pM[GM23]*pM[GM42]; const Double_t det2_24_24 = pM[GM22]*pM[GM44] - pM[GM24]*pM[GM42]; const Double_t det2_24_34 = pM[GM23]*pM[GM44] - pM[GM24]*pM[GM43]; const Double_t det2_34_01 = pM[GM30]*pM[GM41] - pM[GM31]*pM[GM40]; const Double_t det2_34_02 = pM[GM30]*pM[GM42] - pM[GM32]*pM[GM40]; const Double_t det2_34_03 = pM[GM30]*pM[GM43] - pM[GM33]*pM[GM40]; const Double_t det2_34_04 = pM[GM30]*pM[GM44] - pM[GM34]*pM[GM40]; const Double_t det2_34_12 = pM[GM31]*pM[GM42] - pM[GM32]*pM[GM41]; const Double_t det2_34_13 = pM[GM31]*pM[GM43] - pM[GM33]*pM[GM41]; const Double_t det2_34_14 = pM[GM31]*pM[GM44] - pM[GM34]*pM[GM41]; const Double_t det2_34_23 = pM[GM32]*pM[GM43] - pM[GM33]*pM[GM42]; const Double_t det2_34_24 = pM[GM32]*pM[GM44] - pM[GM34]*pM[GM42]; const Double_t det2_34_34 = pM[GM33]*pM[GM44] - pM[GM34]*pM[GM43]; // Find all NECESSARY 3x3 dets: (40 of them) const Double_t det3_123_012 = pM[GM10]*det2_23_12 - pM[GM11]*det2_23_02 + pM[GM12]*det2_23_01; const Double_t det3_123_013 = pM[GM10]*det2_23_13 - pM[GM11]*det2_23_03 + pM[GM13]*det2_23_01; const Double_t det3_123_014 = pM[GM10]*det2_23_14 - pM[GM11]*det2_23_04 + pM[GM14]*det2_23_01; const Double_t det3_123_023 = pM[GM10]*det2_23_23 - pM[GM12]*det2_23_03 + pM[GM13]*det2_23_02; const Double_t det3_123_024 = pM[GM10]*det2_23_24 - pM[GM12]*det2_23_04 + pM[GM14]*det2_23_02; const Double_t det3_123_034 = pM[GM10]*det2_23_34 - pM[GM13]*det2_23_04 + pM[GM14]*det2_23_03; const Double_t det3_123_123 = pM[GM11]*det2_23_23 - pM[GM12]*det2_23_13 + pM[GM13]*det2_23_12; const Double_t det3_123_124 = pM[GM11]*det2_23_24 - pM[GM12]*det2_23_14 + pM[GM14]*det2_23_12; const Double_t det3_123_134 = pM[GM11]*det2_23_34 - pM[GM13]*det2_23_14 + pM[GM14]*det2_23_13; const Double_t det3_123_234 = pM[GM12]*det2_23_34 - pM[GM13]*det2_23_24 + pM[GM14]*det2_23_23; const Double_t det3_124_012 = pM[GM10]*det2_24_12 - pM[GM11]*det2_24_02 + pM[GM12]*det2_24_01; const Double_t det3_124_013 = pM[GM10]*det2_24_13 - pM[GM11]*det2_24_03 + pM[GM13]*det2_24_01; const Double_t det3_124_014 = pM[GM10]*det2_24_14 - pM[GM11]*det2_24_04 + pM[GM14]*det2_24_01; const Double_t det3_124_023 = pM[GM10]*det2_24_23 - pM[GM12]*det2_24_03 + pM[GM13]*det2_24_02; const Double_t det3_124_024 = pM[GM10]*det2_24_24 - pM[GM12]*det2_24_04 + pM[GM14]*det2_24_02; const Double_t det3_124_034 = pM[GM10]*det2_24_34 - pM[GM13]*det2_24_04 + pM[GM14]*det2_24_03; const Double_t det3_124_123 = pM[GM11]*det2_24_23 - pM[GM12]*det2_24_13 + pM[GM13]*det2_24_12; const Double_t det3_124_124 = pM[GM11]*det2_24_24 - pM[GM12]*det2_24_14 + pM[GM14]*det2_24_12; const Double_t det3_124_134 = pM[GM11]*det2_24_34 - pM[GM13]*det2_24_14 + pM[GM14]*det2_24_13; const Double_t det3_124_234 = pM[GM12]*det2_24_34 - pM[GM13]*det2_24_24 + pM[GM14]*det2_24_23; const Double_t det3_134_012 = pM[GM10]*det2_34_12 - pM[GM11]*det2_34_02 + pM[GM12]*det2_34_01; const Double_t det3_134_013 = pM[GM10]*det2_34_13 - pM[GM11]*det2_34_03 + pM[GM13]*det2_34_01; const Double_t det3_134_014 = pM[GM10]*det2_34_14 - pM[GM11]*det2_34_04 + pM[GM14]*det2_34_01; const Double_t det3_134_023 = pM[GM10]*det2_34_23 - pM[GM12]*det2_34_03 + pM[GM13]*det2_34_02; const Double_t det3_134_024 = pM[GM10]*det2_34_24 - pM[GM12]*det2_34_04 + pM[GM14]*det2_34_02; const Double_t det3_134_034 = pM[GM10]*det2_34_34 - pM[GM13]*det2_34_04 + pM[GM14]*det2_34_03; const Double_t det3_134_123 = pM[GM11]*det2_34_23 - pM[GM12]*det2_34_13 + pM[GM13]*det2_34_12; const Double_t det3_134_124 = pM[GM11]*det2_34_24 - pM[GM12]*det2_34_14 + pM[GM14]*det2_34_12; const Double_t det3_134_134 = pM[GM11]*det2_34_34 - pM[GM13]*det2_34_14 + pM[GM14]*det2_34_13; const Double_t det3_134_234 = pM[GM12]*det2_34_34 - pM[GM13]*det2_34_24 + pM[GM14]*det2_34_23; const Double_t det3_234_012 = pM[GM20]*det2_34_12 - pM[GM21]*det2_34_02 + pM[GM22]*det2_34_01; const Double_t det3_234_013 = pM[GM20]*det2_34_13 - pM[GM21]*det2_34_03 + pM[GM23]*det2_34_01; const Double_t det3_234_014 = pM[GM20]*det2_34_14 - pM[GM21]*det2_34_04 + pM[GM24]*det2_34_01; const Double_t det3_234_023 = pM[GM20]*det2_34_23 - pM[GM22]*det2_34_03 + pM[GM23]*det2_34_02; const Double_t det3_234_024 = pM[GM20]*det2_34_24 - pM[GM22]*det2_34_04 + pM[GM24]*det2_34_02; const Double_t det3_234_034 = pM[GM20]*det2_34_34 - pM[GM23]*det2_34_04 + pM[GM24]*det2_34_03; const Double_t det3_234_123 = pM[GM21]*det2_34_23 - pM[GM22]*det2_34_13 + pM[GM23]*det2_34_12; const Double_t det3_234_124 = pM[GM21]*det2_34_24 - pM[GM22]*det2_34_14 + pM[GM24]*det2_34_12; const Double_t det3_234_134 = pM[GM21]*det2_34_34 - pM[GM23]*det2_34_14 + pM[GM24]*det2_34_13; const Double_t det3_234_234 = pM[GM22]*det2_34_34 - pM[GM23]*det2_34_24 + pM[GM24]*det2_34_23; // Find all NECESSARY 4x4 dets: (25 of them) const Double_t det4_0123_0123 = pM[GM00]*det3_123_123 - pM[GM01]*det3_123_023 + pM[GM02]*det3_123_013 - pM[GM03]*det3_123_012; const Double_t det4_0123_0124 = pM[GM00]*det3_123_124 - pM[GM01]*det3_123_024 + pM[GM02]*det3_123_014 - pM[GM04]*det3_123_012; const Double_t det4_0123_0134 = pM[GM00]*det3_123_134 - pM[GM01]*det3_123_034 + pM[GM03]*det3_123_014 - pM[GM04]*det3_123_013; const Double_t det4_0123_0234 = pM[GM00]*det3_123_234 - pM[GM02]*det3_123_034 + pM[GM03]*det3_123_024 - pM[GM04]*det3_123_023; const Double_t det4_0123_1234 = pM[GM01]*det3_123_234 - pM[GM02]*det3_123_134 + pM[GM03]*det3_123_124 - pM[GM04]*det3_123_123; const Double_t det4_0124_0123 = pM[GM00]*det3_124_123 - pM[GM01]*det3_124_023 + pM[GM02]*det3_124_013 - pM[GM03]*det3_124_012; const Double_t det4_0124_0124 = pM[GM00]*det3_124_124 - pM[GM01]*det3_124_024 + pM[GM02]*det3_124_014 - pM[GM04]*det3_124_012; const Double_t det4_0124_0134 = pM[GM00]*det3_124_134 - pM[GM01]*det3_124_034 + pM[GM03]*det3_124_014 - pM[GM04]*det3_124_013; const Double_t det4_0124_0234 = pM[GM00]*det3_124_234 - pM[GM02]*det3_124_034 + pM[GM03]*det3_124_024 - pM[GM04]*det3_124_023; const Double_t det4_0124_1234 = pM[GM01]*det3_124_234 - pM[GM02]*det3_124_134 + pM[GM03]*det3_124_124 - pM[GM04]*det3_124_123; const Double_t det4_0134_0123 = pM[GM00]*det3_134_123 - pM[GM01]*det3_134_023 + pM[GM02]*det3_134_013 - pM[GM03]*det3_134_012; const Double_t det4_0134_0124 = pM[GM00]*det3_134_124 - pM[GM01]*det3_134_024 + pM[GM02]*det3_134_014 - pM[GM04]*det3_134_012; const Double_t det4_0134_0134 = pM[GM00]*det3_134_134 - pM[GM01]*det3_134_034 + pM[GM03]*det3_134_014 - pM[GM04]*det3_134_013; const Double_t det4_0134_0234 = pM[GM00]*det3_134_234 - pM[GM02]*det3_134_034 + pM[GM03]*det3_134_024 - pM[GM04]*det3_134_023; const Double_t det4_0134_1234 = pM[GM01]*det3_134_234 - pM[GM02]*det3_134_134 + pM[GM03]*det3_134_124 - pM[GM04]*det3_134_123; const Double_t det4_0234_0123 = pM[GM00]*det3_234_123 - pM[GM01]*det3_234_023 + pM[GM02]*det3_234_013 - pM[GM03]*det3_234_012; const Double_t det4_0234_0124 = pM[GM00]*det3_234_124 - pM[GM01]*det3_234_024 + pM[GM02]*det3_234_014 - pM[GM04]*det3_234_012; const Double_t det4_0234_0134 = pM[GM00]*det3_234_134 - pM[GM01]*det3_234_034 + pM[GM03]*det3_234_014 - pM[GM04]*det3_234_013; const Double_t det4_0234_0234 = pM[GM00]*det3_234_234 - pM[GM02]*det3_234_034 + pM[GM03]*det3_234_024 - pM[GM04]*det3_234_023; const Double_t det4_0234_1234 = pM[GM01]*det3_234_234 - pM[GM02]*det3_234_134 + pM[GM03]*det3_234_124 - pM[GM04]*det3_234_123; const Double_t det4_1234_0123 = pM[GM10]*det3_234_123 - pM[GM11]*det3_234_023 + pM[GM12]*det3_234_013 - pM[GM13]*det3_234_012; const Double_t det4_1234_0124 = pM[GM10]*det3_234_124 - pM[GM11]*det3_234_024 + pM[GM12]*det3_234_014 - pM[GM14]*det3_234_012; const Double_t det4_1234_0134 = pM[GM10]*det3_234_134 - pM[GM11]*det3_234_034 + pM[GM13]*det3_234_014 - pM[GM14]*det3_234_013; const Double_t det4_1234_0234 = pM[GM10]*det3_234_234 - pM[GM12]*det3_234_034 + pM[GM13]*det3_234_024 - pM[GM14]*det3_234_023; const Double_t det4_1234_1234 = pM[GM11]*det3_234_234 - pM[GM12]*det3_234_134 + pM[GM13]*det3_234_124 - pM[GM14]*det3_234_123; // Find the 5x5 det: const Double_t det = pM[GM00]*det4_1234_1234 - pM[GM01]*det4_1234_0234 + pM[GM02]*det4_1234_0134 - pM[GM03]*det4_1234_0124 + pM[GM04]*det4_1234_0123; if (determ) *determ = det; if ( det == 0 ) { Error("Inv5x5","matrix is singular"); return kFALSE; } const Double_t oneOverDet = 1.0/det; const Double_t mn1OverDet = - oneOverDet; pM[GM00] = det4_1234_1234 * oneOverDet; pM[GM01] = det4_0234_1234 * mn1OverDet; pM[GM02] = det4_0134_1234 * oneOverDet; pM[GM03] = det4_0124_1234 * mn1OverDet; pM[GM04] = det4_0123_1234 * oneOverDet; pM[GM10] = det4_1234_0234 * mn1OverDet; pM[GM11] = det4_0234_0234 * oneOverDet; pM[GM12] = det4_0134_0234 * mn1OverDet; pM[GM13] = det4_0124_0234 * oneOverDet; pM[GM14] = det4_0123_0234 * mn1OverDet; pM[GM20] = det4_1234_0134 * oneOverDet; pM[GM21] = det4_0234_0134 * mn1OverDet; pM[GM22] = det4_0134_0134 * oneOverDet; pM[GM23] = det4_0124_0134 * mn1OverDet; pM[GM24] = det4_0123_0134 * oneOverDet; pM[GM30] = det4_1234_0124 * mn1OverDet; pM[GM31] = det4_0234_0124 * oneOverDet; pM[GM32] = det4_0134_0124 * mn1OverDet; pM[GM33] = det4_0124_0124 * oneOverDet; pM[GM34] = det4_0123_0124 * mn1OverDet; pM[GM40] = det4_1234_0123 * oneOverDet; pM[GM41] = det4_0234_0123 * mn1OverDet; pM[GM42] = det4_0134_0123 * oneOverDet; pM[GM43] = det4_0124_0123 * mn1OverDet; pM[GM44] = det4_0123_0123 * oneOverDet; return kTRUE; } // Aij are indices for a 6x6 matrix. #define GA00 0 #define GA01 1 #define GA02 2 #define GA03 3 #define GA04 4 #define GA05 5 #define GA10 6 #define GA11 7 #define GA12 8 #define GA13 9 #define GA14 10 #define GA15 11 #define GA20 12 #define GA21 13 #define GA22 14 #define GA23 15 #define GA24 16 #define GA25 17 #define GA30 18 #define GA31 19 #define GA32 20 #define GA33 21 #define GA34 22 #define GA35 23 #define GA40 24 #define GA41 25 #define GA42 26 #define GA43 27 #define GA44 28 #define GA45 29 #define GA50 30 #define GA51 31 #define GA52 32 #define GA53 33 #define GA54 34 #define GA55 35 //////////////////////////////////////////////////////////////////////////////// template Bool_t TMatrixTCramerInv::Inv6x6(TMatrixT &m,Double_t *determ) { if (m.GetNrows() != 6 || m.GetNcols() != 6 || m.GetRowLwb() != m.GetColLwb()) { Error("Inv6x6","matrix should be square 6x6"); return kFALSE; } Element *pM = m.GetMatrixArray(); // Find all NECESSGARY 2x2 dets: (45 of them) const Double_t det2_34_01 = pM[GA30]*pM[GA41] - pM[GA31]*pM[GA40]; const Double_t det2_34_02 = pM[GA30]*pM[GA42] - pM[GA32]*pM[GA40]; const Double_t det2_34_03 = pM[GA30]*pM[GA43] - pM[GA33]*pM[GA40]; const Double_t det2_34_04 = pM[GA30]*pM[GA44] - pM[GA34]*pM[GA40]; const Double_t det2_34_05 = pM[GA30]*pM[GA45] - pM[GA35]*pM[GA40]; const Double_t det2_34_12 = pM[GA31]*pM[GA42] - pM[GA32]*pM[GA41]; const Double_t det2_34_13 = pM[GA31]*pM[GA43] - pM[GA33]*pM[GA41]; const Double_t det2_34_14 = pM[GA31]*pM[GA44] - pM[GA34]*pM[GA41]; const Double_t det2_34_15 = pM[GA31]*pM[GA45] - pM[GA35]*pM[GA41]; const Double_t det2_34_23 = pM[GA32]*pM[GA43] - pM[GA33]*pM[GA42]; const Double_t det2_34_24 = pM[GA32]*pM[GA44] - pM[GA34]*pM[GA42]; const Double_t det2_34_25 = pM[GA32]*pM[GA45] - pM[GA35]*pM[GA42]; const Double_t det2_34_34 = pM[GA33]*pM[GA44] - pM[GA34]*pM[GA43]; const Double_t det2_34_35 = pM[GA33]*pM[GA45] - pM[GA35]*pM[GA43]; const Double_t det2_34_45 = pM[GA34]*pM[GA45] - pM[GA35]*pM[GA44]; const Double_t det2_35_01 = pM[GA30]*pM[GA51] - pM[GA31]*pM[GA50]; const Double_t det2_35_02 = pM[GA30]*pM[GA52] - pM[GA32]*pM[GA50]; const Double_t det2_35_03 = pM[GA30]*pM[GA53] - pM[GA33]*pM[GA50]; const Double_t det2_35_04 = pM[GA30]*pM[GA54] - pM[GA34]*pM[GA50]; const Double_t det2_35_05 = pM[GA30]*pM[GA55] - pM[GA35]*pM[GA50]; const Double_t det2_35_12 = pM[GA31]*pM[GA52] - pM[GA32]*pM[GA51]; const Double_t det2_35_13 = pM[GA31]*pM[GA53] - pM[GA33]*pM[GA51]; const Double_t det2_35_14 = pM[GA31]*pM[GA54] - pM[GA34]*pM[GA51]; const Double_t det2_35_15 = pM[GA31]*pM[GA55] - pM[GA35]*pM[GA51]; const Double_t det2_35_23 = pM[GA32]*pM[GA53] - pM[GA33]*pM[GA52]; const Double_t det2_35_24 = pM[GA32]*pM[GA54] - pM[GA34]*pM[GA52]; const Double_t det2_35_25 = pM[GA32]*pM[GA55] - pM[GA35]*pM[GA52]; const Double_t det2_35_34 = pM[GA33]*pM[GA54] - pM[GA34]*pM[GA53]; const Double_t det2_35_35 = pM[GA33]*pM[GA55] - pM[GA35]*pM[GA53]; const Double_t det2_35_45 = pM[GA34]*pM[GA55] - pM[GA35]*pM[GA54]; const Double_t det2_45_01 = pM[GA40]*pM[GA51] - pM[GA41]*pM[GA50]; const Double_t det2_45_02 = pM[GA40]*pM[GA52] - pM[GA42]*pM[GA50]; const Double_t det2_45_03 = pM[GA40]*pM[GA53] - pM[GA43]*pM[GA50]; const Double_t det2_45_04 = pM[GA40]*pM[GA54] - pM[GA44]*pM[GA50]; const Double_t det2_45_05 = pM[GA40]*pM[GA55] - pM[GA45]*pM[GA50]; const Double_t det2_45_12 = pM[GA41]*pM[GA52] - pM[GA42]*pM[GA51]; const Double_t det2_45_13 = pM[GA41]*pM[GA53] - pM[GA43]*pM[GA51]; const Double_t det2_45_14 = pM[GA41]*pM[GA54] - pM[GA44]*pM[GA51]; const Double_t det2_45_15 = pM[GA41]*pM[GA55] - pM[GA45]*pM[GA51]; const Double_t det2_45_23 = pM[GA42]*pM[GA53] - pM[GA43]*pM[GA52]; const Double_t det2_45_24 = pM[GA42]*pM[GA54] - pM[GA44]*pM[GA52]; const Double_t det2_45_25 = pM[GA42]*pM[GA55] - pM[GA45]*pM[GA52]; const Double_t det2_45_34 = pM[GA43]*pM[GA54] - pM[GA44]*pM[GA53]; const Double_t det2_45_35 = pM[GA43]*pM[GA55] - pM[GA45]*pM[GA53]; const Double_t det2_45_45 = pM[GA44]*pM[GA55] - pM[GA45]*pM[GA54]; // Find all NECESSGARY 3x3 dets: (80 of them) const Double_t det3_234_012 = pM[GA20]*det2_34_12 - pM[GA21]*det2_34_02 + pM[GA22]*det2_34_01; const Double_t det3_234_013 = pM[GA20]*det2_34_13 - pM[GA21]*det2_34_03 + pM[GA23]*det2_34_01; const Double_t det3_234_014 = pM[GA20]*det2_34_14 - pM[GA21]*det2_34_04 + pM[GA24]*det2_34_01; const Double_t det3_234_015 = pM[GA20]*det2_34_15 - pM[GA21]*det2_34_05 + pM[GA25]*det2_34_01; const Double_t det3_234_023 = pM[GA20]*det2_34_23 - pM[GA22]*det2_34_03 + pM[GA23]*det2_34_02; const Double_t det3_234_024 = pM[GA20]*det2_34_24 - pM[GA22]*det2_34_04 + pM[GA24]*det2_34_02; const Double_t det3_234_025 = pM[GA20]*det2_34_25 - pM[GA22]*det2_34_05 + pM[GA25]*det2_34_02; const Double_t det3_234_034 = pM[GA20]*det2_34_34 - pM[GA23]*det2_34_04 + pM[GA24]*det2_34_03; const Double_t det3_234_035 = pM[GA20]*det2_34_35 - pM[GA23]*det2_34_05 + pM[GA25]*det2_34_03; const Double_t det3_234_045 = pM[GA20]*det2_34_45 - pM[GA24]*det2_34_05 + pM[GA25]*det2_34_04; const Double_t det3_234_123 = pM[GA21]*det2_34_23 - pM[GA22]*det2_34_13 + pM[GA23]*det2_34_12; const Double_t det3_234_124 = pM[GA21]*det2_34_24 - pM[GA22]*det2_34_14 + pM[GA24]*det2_34_12; const Double_t det3_234_125 = pM[GA21]*det2_34_25 - pM[GA22]*det2_34_15 + pM[GA25]*det2_34_12; const Double_t det3_234_134 = pM[GA21]*det2_34_34 - pM[GA23]*det2_34_14 + pM[GA24]*det2_34_13; const Double_t det3_234_135 = pM[GA21]*det2_34_35 - pM[GA23]*det2_34_15 + pM[GA25]*det2_34_13; const Double_t det3_234_145 = pM[GA21]*det2_34_45 - pM[GA24]*det2_34_15 + pM[GA25]*det2_34_14; const Double_t det3_234_234 = pM[GA22]*det2_34_34 - pM[GA23]*det2_34_24 + pM[GA24]*det2_34_23; const Double_t det3_234_235 = pM[GA22]*det2_34_35 - pM[GA23]*det2_34_25 + pM[GA25]*det2_34_23; const Double_t det3_234_245 = pM[GA22]*det2_34_45 - pM[GA24]*det2_34_25 + pM[GA25]*det2_34_24; const Double_t det3_234_345 = pM[GA23]*det2_34_45 - pM[GA24]*det2_34_35 + pM[GA25]*det2_34_34; const Double_t det3_235_012 = pM[GA20]*det2_35_12 - pM[GA21]*det2_35_02 + pM[GA22]*det2_35_01; const Double_t det3_235_013 = pM[GA20]*det2_35_13 - pM[GA21]*det2_35_03 + pM[GA23]*det2_35_01; const Double_t det3_235_014 = pM[GA20]*det2_35_14 - pM[GA21]*det2_35_04 + pM[GA24]*det2_35_01; const Double_t det3_235_015 = pM[GA20]*det2_35_15 - pM[GA21]*det2_35_05 + pM[GA25]*det2_35_01; const Double_t det3_235_023 = pM[GA20]*det2_35_23 - pM[GA22]*det2_35_03 + pM[GA23]*det2_35_02; const Double_t det3_235_024 = pM[GA20]*det2_35_24 - pM[GA22]*det2_35_04 + pM[GA24]*det2_35_02; const Double_t det3_235_025 = pM[GA20]*det2_35_25 - pM[GA22]*det2_35_05 + pM[GA25]*det2_35_02; const Double_t det3_235_034 = pM[GA20]*det2_35_34 - pM[GA23]*det2_35_04 + pM[GA24]*det2_35_03; const Double_t det3_235_035 = pM[GA20]*det2_35_35 - pM[GA23]*det2_35_05 + pM[GA25]*det2_35_03; const Double_t det3_235_045 = pM[GA20]*det2_35_45 - pM[GA24]*det2_35_05 + pM[GA25]*det2_35_04; const Double_t det3_235_123 = pM[GA21]*det2_35_23 - pM[GA22]*det2_35_13 + pM[GA23]*det2_35_12; const Double_t det3_235_124 = pM[GA21]*det2_35_24 - pM[GA22]*det2_35_14 + pM[GA24]*det2_35_12; const Double_t det3_235_125 = pM[GA21]*det2_35_25 - pM[GA22]*det2_35_15 + pM[GA25]*det2_35_12; const Double_t det3_235_134 = pM[GA21]*det2_35_34 - pM[GA23]*det2_35_14 + pM[GA24]*det2_35_13; const Double_t det3_235_135 = pM[GA21]*det2_35_35 - pM[GA23]*det2_35_15 + pM[GA25]*det2_35_13; const Double_t det3_235_145 = pM[GA21]*det2_35_45 - pM[GA24]*det2_35_15 + pM[GA25]*det2_35_14; const Double_t det3_235_234 = pM[GA22]*det2_35_34 - pM[GA23]*det2_35_24 + pM[GA24]*det2_35_23; const Double_t det3_235_235 = pM[GA22]*det2_35_35 - pM[GA23]*det2_35_25 + pM[GA25]*det2_35_23; const Double_t det3_235_245 = pM[GA22]*det2_35_45 - pM[GA24]*det2_35_25 + pM[GA25]*det2_35_24; const Double_t det3_235_345 = pM[GA23]*det2_35_45 - pM[GA24]*det2_35_35 + pM[GA25]*det2_35_34; const Double_t det3_245_012 = pM[GA20]*det2_45_12 - pM[GA21]*det2_45_02 + pM[GA22]*det2_45_01; const Double_t det3_245_013 = pM[GA20]*det2_45_13 - pM[GA21]*det2_45_03 + pM[GA23]*det2_45_01; const Double_t det3_245_014 = pM[GA20]*det2_45_14 - pM[GA21]*det2_45_04 + pM[GA24]*det2_45_01; const Double_t det3_245_015 = pM[GA20]*det2_45_15 - pM[GA21]*det2_45_05 + pM[GA25]*det2_45_01; const Double_t det3_245_023 = pM[GA20]*det2_45_23 - pM[GA22]*det2_45_03 + pM[GA23]*det2_45_02; const Double_t det3_245_024 = pM[GA20]*det2_45_24 - pM[GA22]*det2_45_04 + pM[GA24]*det2_45_02; const Double_t det3_245_025 = pM[GA20]*det2_45_25 - pM[GA22]*det2_45_05 + pM[GA25]*det2_45_02; const Double_t det3_245_034 = pM[GA20]*det2_45_34 - pM[GA23]*det2_45_04 + pM[GA24]*det2_45_03; const Double_t det3_245_035 = pM[GA20]*det2_45_35 - pM[GA23]*det2_45_05 + pM[GA25]*det2_45_03; const Double_t det3_245_045 = pM[GA20]*det2_45_45 - pM[GA24]*det2_45_05 + pM[GA25]*det2_45_04; const Double_t det3_245_123 = pM[GA21]*det2_45_23 - pM[GA22]*det2_45_13 + pM[GA23]*det2_45_12; const Double_t det3_245_124 = pM[GA21]*det2_45_24 - pM[GA22]*det2_45_14 + pM[GA24]*det2_45_12; const Double_t det3_245_125 = pM[GA21]*det2_45_25 - pM[GA22]*det2_45_15 + pM[GA25]*det2_45_12; const Double_t det3_245_134 = pM[GA21]*det2_45_34 - pM[GA23]*det2_45_14 + pM[GA24]*det2_45_13; const Double_t det3_245_135 = pM[GA21]*det2_45_35 - pM[GA23]*det2_45_15 + pM[GA25]*det2_45_13; const Double_t det3_245_145 = pM[GA21]*det2_45_45 - pM[GA24]*det2_45_15 + pM[GA25]*det2_45_14; const Double_t det3_245_234 = pM[GA22]*det2_45_34 - pM[GA23]*det2_45_24 + pM[GA24]*det2_45_23; const Double_t det3_245_235 = pM[GA22]*det2_45_35 - pM[GA23]*det2_45_25 + pM[GA25]*det2_45_23; const Double_t det3_245_245 = pM[GA22]*det2_45_45 - pM[GA24]*det2_45_25 + pM[GA25]*det2_45_24; const Double_t det3_245_345 = pM[GA23]*det2_45_45 - pM[GA24]*det2_45_35 + pM[GA25]*det2_45_34; const Double_t det3_345_012 = pM[GA30]*det2_45_12 - pM[GA31]*det2_45_02 + pM[GA32]*det2_45_01; const Double_t det3_345_013 = pM[GA30]*det2_45_13 - pM[GA31]*det2_45_03 + pM[GA33]*det2_45_01; const Double_t det3_345_014 = pM[GA30]*det2_45_14 - pM[GA31]*det2_45_04 + pM[GA34]*det2_45_01; const Double_t det3_345_015 = pM[GA30]*det2_45_15 - pM[GA31]*det2_45_05 + pM[GA35]*det2_45_01; const Double_t det3_345_023 = pM[GA30]*det2_45_23 - pM[GA32]*det2_45_03 + pM[GA33]*det2_45_02; const Double_t det3_345_024 = pM[GA30]*det2_45_24 - pM[GA32]*det2_45_04 + pM[GA34]*det2_45_02; const Double_t det3_345_025 = pM[GA30]*det2_45_25 - pM[GA32]*det2_45_05 + pM[GA35]*det2_45_02; const Double_t det3_345_034 = pM[GA30]*det2_45_34 - pM[GA33]*det2_45_04 + pM[GA34]*det2_45_03; const Double_t det3_345_035 = pM[GA30]*det2_45_35 - pM[GA33]*det2_45_05 + pM[GA35]*det2_45_03; const Double_t det3_345_045 = pM[GA30]*det2_45_45 - pM[GA34]*det2_45_05 + pM[GA35]*det2_45_04; const Double_t det3_345_123 = pM[GA31]*det2_45_23 - pM[GA32]*det2_45_13 + pM[GA33]*det2_45_12; const Double_t det3_345_124 = pM[GA31]*det2_45_24 - pM[GA32]*det2_45_14 + pM[GA34]*det2_45_12; const Double_t det3_345_125 = pM[GA31]*det2_45_25 - pM[GA32]*det2_45_15 + pM[GA35]*det2_45_12; const Double_t det3_345_134 = pM[GA31]*det2_45_34 - pM[GA33]*det2_45_14 + pM[GA34]*det2_45_13; const Double_t det3_345_135 = pM[GA31]*det2_45_35 - pM[GA33]*det2_45_15 + pM[GA35]*det2_45_13; const Double_t det3_345_145 = pM[GA31]*det2_45_45 - pM[GA34]*det2_45_15 + pM[GA35]*det2_45_14; const Double_t det3_345_234 = pM[GA32]*det2_45_34 - pM[GA33]*det2_45_24 + pM[GA34]*det2_45_23; const Double_t det3_345_235 = pM[GA32]*det2_45_35 - pM[GA33]*det2_45_25 + pM[GA35]*det2_45_23; const Double_t det3_345_245 = pM[GA32]*det2_45_45 - pM[GA34]*det2_45_25 + pM[GA35]*det2_45_24; const Double_t det3_345_345 = pM[GA33]*det2_45_45 - pM[GA34]*det2_45_35 + pM[GA35]*det2_45_34; // Find all NECESSGARY 4x4 dets: (75 of them) const Double_t det4_1234_0123 = pM[GA10]*det3_234_123 - pM[GA11]*det3_234_023 + pM[GA12]*det3_234_013 - pM[GA13]*det3_234_012; const Double_t det4_1234_0124 = pM[GA10]*det3_234_124 - pM[GA11]*det3_234_024 + pM[GA12]*det3_234_014 - pM[GA14]*det3_234_012; const Double_t det4_1234_0125 = pM[GA10]*det3_234_125 - pM[GA11]*det3_234_025 + pM[GA12]*det3_234_015 - pM[GA15]*det3_234_012; const Double_t det4_1234_0134 = pM[GA10]*det3_234_134 - pM[GA11]*det3_234_034 + pM[GA13]*det3_234_014 - pM[GA14]*det3_234_013; const Double_t det4_1234_0135 = pM[GA10]*det3_234_135 - pM[GA11]*det3_234_035 + pM[GA13]*det3_234_015 - pM[GA15]*det3_234_013; const Double_t det4_1234_0145 = pM[GA10]*det3_234_145 - pM[GA11]*det3_234_045 + pM[GA14]*det3_234_015 - pM[GA15]*det3_234_014; const Double_t det4_1234_0234 = pM[GA10]*det3_234_234 - pM[GA12]*det3_234_034 + pM[GA13]*det3_234_024 - pM[GA14]*det3_234_023; const Double_t det4_1234_0235 = pM[GA10]*det3_234_235 - pM[GA12]*det3_234_035 + pM[GA13]*det3_234_025 - pM[GA15]*det3_234_023; const Double_t det4_1234_0245 = pM[GA10]*det3_234_245 - pM[GA12]*det3_234_045 + pM[GA14]*det3_234_025 - pM[GA15]*det3_234_024; const Double_t det4_1234_0345 = pM[GA10]*det3_234_345 - pM[GA13]*det3_234_045 + pM[GA14]*det3_234_035 - pM[GA15]*det3_234_034; const Double_t det4_1234_1234 = pM[GA11]*det3_234_234 - pM[GA12]*det3_234_134 + pM[GA13]*det3_234_124 - pM[GA14]*det3_234_123; const Double_t det4_1234_1235 = pM[GA11]*det3_234_235 - pM[GA12]*det3_234_135 + pM[GA13]*det3_234_125 - pM[GA15]*det3_234_123; const Double_t det4_1234_1245 = pM[GA11]*det3_234_245 - pM[GA12]*det3_234_145 + pM[GA14]*det3_234_125 - pM[GA15]*det3_234_124; const Double_t det4_1234_1345 = pM[GA11]*det3_234_345 - pM[GA13]*det3_234_145 + pM[GA14]*det3_234_135 - pM[GA15]*det3_234_134; const Double_t det4_1234_2345 = pM[GA12]*det3_234_345 - pM[GA13]*det3_234_245 + pM[GA14]*det3_234_235 - pM[GA15]*det3_234_234; const Double_t det4_1235_0123 = pM[GA10]*det3_235_123 - pM[GA11]*det3_235_023 + pM[GA12]*det3_235_013 - pM[GA13]*det3_235_012; const Double_t det4_1235_0124 = pM[GA10]*det3_235_124 - pM[GA11]*det3_235_024 + pM[GA12]*det3_235_014 - pM[GA14]*det3_235_012; const Double_t det4_1235_0125 = pM[GA10]*det3_235_125 - pM[GA11]*det3_235_025 + pM[GA12]*det3_235_015 - pM[GA15]*det3_235_012; const Double_t det4_1235_0134 = pM[GA10]*det3_235_134 - pM[GA11]*det3_235_034 + pM[GA13]*det3_235_014 - pM[GA14]*det3_235_013; const Double_t det4_1235_0135 = pM[GA10]*det3_235_135 - pM[GA11]*det3_235_035 + pM[GA13]*det3_235_015 - pM[GA15]*det3_235_013; const Double_t det4_1235_0145 = pM[GA10]*det3_235_145 - pM[GA11]*det3_235_045 + pM[GA14]*det3_235_015 - pM[GA15]*det3_235_014; const Double_t det4_1235_0234 = pM[GA10]*det3_235_234 - pM[GA12]*det3_235_034 + pM[GA13]*det3_235_024 - pM[GA14]*det3_235_023; const Double_t det4_1235_0235 = pM[GA10]*det3_235_235 - pM[GA12]*det3_235_035 + pM[GA13]*det3_235_025 - pM[GA15]*det3_235_023; const Double_t det4_1235_0245 = pM[GA10]*det3_235_245 - pM[GA12]*det3_235_045 + pM[GA14]*det3_235_025 - pM[GA15]*det3_235_024; const Double_t det4_1235_0345 = pM[GA10]*det3_235_345 - pM[GA13]*det3_235_045 + pM[GA14]*det3_235_035 - pM[GA15]*det3_235_034; const Double_t det4_1235_1234 = pM[GA11]*det3_235_234 - pM[GA12]*det3_235_134 + pM[GA13]*det3_235_124 - pM[GA14]*det3_235_123; const Double_t det4_1235_1235 = pM[GA11]*det3_235_235 - pM[GA12]*det3_235_135 + pM[GA13]*det3_235_125 - pM[GA15]*det3_235_123; const Double_t det4_1235_1245 = pM[GA11]*det3_235_245 - pM[GA12]*det3_235_145 + pM[GA14]*det3_235_125 - pM[GA15]*det3_235_124; const Double_t det4_1235_1345 = pM[GA11]*det3_235_345 - pM[GA13]*det3_235_145 + pM[GA14]*det3_235_135 - pM[GA15]*det3_235_134; const Double_t det4_1235_2345 = pM[GA12]*det3_235_345 - pM[GA13]*det3_235_245 + pM[GA14]*det3_235_235 - pM[GA15]*det3_235_234; const Double_t det4_1245_0123 = pM[GA10]*det3_245_123 - pM[GA11]*det3_245_023 + pM[GA12]*det3_245_013 - pM[GA13]*det3_245_012; const Double_t det4_1245_0124 = pM[GA10]*det3_245_124 - pM[GA11]*det3_245_024 + pM[GA12]*det3_245_014 - pM[GA14]*det3_245_012; const Double_t det4_1245_0125 = pM[GA10]*det3_245_125 - pM[GA11]*det3_245_025 + pM[GA12]*det3_245_015 - pM[GA15]*det3_245_012; const Double_t det4_1245_0134 = pM[GA10]*det3_245_134 - pM[GA11]*det3_245_034 + pM[GA13]*det3_245_014 - pM[GA14]*det3_245_013; const Double_t det4_1245_0135 = pM[GA10]*det3_245_135 - pM[GA11]*det3_245_035 + pM[GA13]*det3_245_015 - pM[GA15]*det3_245_013; const Double_t det4_1245_0145 = pM[GA10]*det3_245_145 - pM[GA11]*det3_245_045 + pM[GA14]*det3_245_015 - pM[GA15]*det3_245_014; const Double_t det4_1245_0234 = pM[GA10]*det3_245_234 - pM[GA12]*det3_245_034 + pM[GA13]*det3_245_024 - pM[GA14]*det3_245_023; const Double_t det4_1245_0235 = pM[GA10]*det3_245_235 - pM[GA12]*det3_245_035 + pM[GA13]*det3_245_025 - pM[GA15]*det3_245_023; const Double_t det4_1245_0245 = pM[GA10]*det3_245_245 - pM[GA12]*det3_245_045 + pM[GA14]*det3_245_025 - pM[GA15]*det3_245_024; const Double_t det4_1245_0345 = pM[GA10]*det3_245_345 - pM[GA13]*det3_245_045 + pM[GA14]*det3_245_035 - pM[GA15]*det3_245_034; const Double_t det4_1245_1234 = pM[GA11]*det3_245_234 - pM[GA12]*det3_245_134 + pM[GA13]*det3_245_124 - pM[GA14]*det3_245_123; const Double_t det4_1245_1235 = pM[GA11]*det3_245_235 - pM[GA12]*det3_245_135 + pM[GA13]*det3_245_125 - pM[GA15]*det3_245_123; const Double_t det4_1245_1245 = pM[GA11]*det3_245_245 - pM[GA12]*det3_245_145 + pM[GA14]*det3_245_125 - pM[GA15]*det3_245_124; const Double_t det4_1245_1345 = pM[GA11]*det3_245_345 - pM[GA13]*det3_245_145 + pM[GA14]*det3_245_135 - pM[GA15]*det3_245_134; const Double_t det4_1245_2345 = pM[GA12]*det3_245_345 - pM[GA13]*det3_245_245 + pM[GA14]*det3_245_235 - pM[GA15]*det3_245_234; const Double_t det4_1345_0123 = pM[GA10]*det3_345_123 - pM[GA11]*det3_345_023 + pM[GA12]*det3_345_013 - pM[GA13]*det3_345_012; const Double_t det4_1345_0124 = pM[GA10]*det3_345_124 - pM[GA11]*det3_345_024 + pM[GA12]*det3_345_014 - pM[GA14]*det3_345_012; const Double_t det4_1345_0125 = pM[GA10]*det3_345_125 - pM[GA11]*det3_345_025 + pM[GA12]*det3_345_015 - pM[GA15]*det3_345_012; const Double_t det4_1345_0134 = pM[GA10]*det3_345_134 - pM[GA11]*det3_345_034 + pM[GA13]*det3_345_014 - pM[GA14]*det3_345_013; const Double_t det4_1345_0135 = pM[GA10]*det3_345_135 - pM[GA11]*det3_345_035 + pM[GA13]*det3_345_015 - pM[GA15]*det3_345_013; const Double_t det4_1345_0145 = pM[GA10]*det3_345_145 - pM[GA11]*det3_345_045 + pM[GA14]*det3_345_015 - pM[GA15]*det3_345_014; const Double_t det4_1345_0234 = pM[GA10]*det3_345_234 - pM[GA12]*det3_345_034 + pM[GA13]*det3_345_024 - pM[GA14]*det3_345_023; const Double_t det4_1345_0235 = pM[GA10]*det3_345_235 - pM[GA12]*det3_345_035 + pM[GA13]*det3_345_025 - pM[GA15]*det3_345_023; const Double_t det4_1345_0245 = pM[GA10]*det3_345_245 - pM[GA12]*det3_345_045 + pM[GA14]*det3_345_025 - pM[GA15]*det3_345_024; const Double_t det4_1345_0345 = pM[GA10]*det3_345_345 - pM[GA13]*det3_345_045 + pM[GA14]*det3_345_035 - pM[GA15]*det3_345_034; const Double_t det4_1345_1234 = pM[GA11]*det3_345_234 - pM[GA12]*det3_345_134 + pM[GA13]*det3_345_124 - pM[GA14]*det3_345_123; const Double_t det4_1345_1235 = pM[GA11]*det3_345_235 - pM[GA12]*det3_345_135 + pM[GA13]*det3_345_125 - pM[GA15]*det3_345_123; const Double_t det4_1345_1245 = pM[GA11]*det3_345_245 - pM[GA12]*det3_345_145 + pM[GA14]*det3_345_125 - pM[GA15]*det3_345_124; const Double_t det4_1345_1345 = pM[GA11]*det3_345_345 - pM[GA13]*det3_345_145 + pM[GA14]*det3_345_135 - pM[GA15]*det3_345_134; const Double_t det4_1345_2345 = pM[GA12]*det3_345_345 - pM[GA13]*det3_345_245 + pM[GA14]*det3_345_235 - pM[GA15]*det3_345_234; const Double_t det4_2345_0123 = pM[GA20]*det3_345_123 - pM[GA21]*det3_345_023 + pM[GA22]*det3_345_013 - pM[GA23]*det3_345_012; const Double_t det4_2345_0124 = pM[GA20]*det3_345_124 - pM[GA21]*det3_345_024 + pM[GA22]*det3_345_014 - pM[GA24]*det3_345_012; const Double_t det4_2345_0125 = pM[GA20]*det3_345_125 - pM[GA21]*det3_345_025 + pM[GA22]*det3_345_015 - pM[GA25]*det3_345_012; const Double_t det4_2345_0134 = pM[GA20]*det3_345_134 - pM[GA21]*det3_345_034 + pM[GA23]*det3_345_014 - pM[GA24]*det3_345_013; const Double_t det4_2345_0135 = pM[GA20]*det3_345_135 - pM[GA21]*det3_345_035 + pM[GA23]*det3_345_015 - pM[GA25]*det3_345_013; const Double_t det4_2345_0145 = pM[GA20]*det3_345_145 - pM[GA21]*det3_345_045 + pM[GA24]*det3_345_015 - pM[GA25]*det3_345_014; const Double_t det4_2345_0234 = pM[GA20]*det3_345_234 - pM[GA22]*det3_345_034 + pM[GA23]*det3_345_024 - pM[GA24]*det3_345_023; const Double_t det4_2345_0235 = pM[GA20]*det3_345_235 - pM[GA22]*det3_345_035 + pM[GA23]*det3_345_025 - pM[GA25]*det3_345_023; const Double_t det4_2345_0245 = pM[GA20]*det3_345_245 - pM[GA22]*det3_345_045 + pM[GA24]*det3_345_025 - pM[GA25]*det3_345_024; const Double_t det4_2345_0345 = pM[GA20]*det3_345_345 - pM[GA23]*det3_345_045 + pM[GA24]*det3_345_035 - pM[GA25]*det3_345_034; const Double_t det4_2345_1234 = pM[GA21]*det3_345_234 - pM[GA22]*det3_345_134 + pM[GA23]*det3_345_124 - pM[GA24]*det3_345_123; const Double_t det4_2345_1235 = pM[GA21]*det3_345_235 - pM[GA22]*det3_345_135 + pM[GA23]*det3_345_125 - pM[GA25]*det3_345_123; const Double_t det4_2345_1245 = pM[GA21]*det3_345_245 - pM[GA22]*det3_345_145 + pM[GA24]*det3_345_125 - pM[GA25]*det3_345_124; const Double_t det4_2345_1345 = pM[GA21]*det3_345_345 - pM[GA23]*det3_345_145 + pM[GA24]*det3_345_135 - pM[GA25]*det3_345_134; const Double_t det4_2345_2345 = pM[GA22]*det3_345_345 - pM[GA23]*det3_345_245 + pM[GA24]*det3_345_235 - pM[GA25]*det3_345_234; // Find all NECESSGARY 5x5 dets: (36 of them) const Double_t det5_01234_01234 = pM[GA00]*det4_1234_1234 - pM[GA01]*det4_1234_0234 + pM[GA02]*det4_1234_0134 - pM[GA03]*det4_1234_0124 + pM[GA04]*det4_1234_0123; const Double_t det5_01234_01235 = pM[GA00]*det4_1234_1235 - pM[GA01]*det4_1234_0235 + pM[GA02]*det4_1234_0135 - pM[GA03]*det4_1234_0125 + pM[GA05]*det4_1234_0123; const Double_t det5_01234_01245 = pM[GA00]*det4_1234_1245 - pM[GA01]*det4_1234_0245 + pM[GA02]*det4_1234_0145 - pM[GA04]*det4_1234_0125 + pM[GA05]*det4_1234_0124; const Double_t det5_01234_01345 = pM[GA00]*det4_1234_1345 - pM[GA01]*det4_1234_0345 + pM[GA03]*det4_1234_0145 - pM[GA04]*det4_1234_0135 + pM[GA05]*det4_1234_0134; const Double_t det5_01234_02345 = pM[GA00]*det4_1234_2345 - pM[GA02]*det4_1234_0345 + pM[GA03]*det4_1234_0245 - pM[GA04]*det4_1234_0235 + pM[GA05]*det4_1234_0234; const Double_t det5_01234_12345 = pM[GA01]*det4_1234_2345 - pM[GA02]*det4_1234_1345 + pM[GA03]*det4_1234_1245 - pM[GA04]*det4_1234_1235 + pM[GA05]*det4_1234_1234; const Double_t det5_01235_01234 = pM[GA00]*det4_1235_1234 - pM[GA01]*det4_1235_0234 + pM[GA02]*det4_1235_0134 - pM[GA03]*det4_1235_0124 + pM[GA04]*det4_1235_0123; const Double_t det5_01235_01235 = pM[GA00]*det4_1235_1235 - pM[GA01]*det4_1235_0235 + pM[GA02]*det4_1235_0135 - pM[GA03]*det4_1235_0125 + pM[GA05]*det4_1235_0123; const Double_t det5_01235_01245 = pM[GA00]*det4_1235_1245 - pM[GA01]*det4_1235_0245 + pM[GA02]*det4_1235_0145 - pM[GA04]*det4_1235_0125 + pM[GA05]*det4_1235_0124; const Double_t det5_01235_01345 = pM[GA00]*det4_1235_1345 - pM[GA01]*det4_1235_0345 + pM[GA03]*det4_1235_0145 - pM[GA04]*det4_1235_0135 + pM[GA05]*det4_1235_0134; const Double_t det5_01235_02345 = pM[GA00]*det4_1235_2345 - pM[GA02]*det4_1235_0345 + pM[GA03]*det4_1235_0245 - pM[GA04]*det4_1235_0235 + pM[GA05]*det4_1235_0234; const Double_t det5_01235_12345 = pM[GA01]*det4_1235_2345 - pM[GA02]*det4_1235_1345 + pM[GA03]*det4_1235_1245 - pM[GA04]*det4_1235_1235 + pM[GA05]*det4_1235_1234; const Double_t det5_01245_01234 = pM[GA00]*det4_1245_1234 - pM[GA01]*det4_1245_0234 + pM[GA02]*det4_1245_0134 - pM[GA03]*det4_1245_0124 + pM[GA04]*det4_1245_0123; const Double_t det5_01245_01235 = pM[GA00]*det4_1245_1235 - pM[GA01]*det4_1245_0235 + pM[GA02]*det4_1245_0135 - pM[GA03]*det4_1245_0125 + pM[GA05]*det4_1245_0123; const Double_t det5_01245_01245 = pM[GA00]*det4_1245_1245 - pM[GA01]*det4_1245_0245 + pM[GA02]*det4_1245_0145 - pM[GA04]*det4_1245_0125 + pM[GA05]*det4_1245_0124; const Double_t det5_01245_01345 = pM[GA00]*det4_1245_1345 - pM[GA01]*det4_1245_0345 + pM[GA03]*det4_1245_0145 - pM[GA04]*det4_1245_0135 + pM[GA05]*det4_1245_0134; const Double_t det5_01245_02345 = pM[GA00]*det4_1245_2345 - pM[GA02]*det4_1245_0345 + pM[GA03]*det4_1245_0245 - pM[GA04]*det4_1245_0235 + pM[GA05]*det4_1245_0234; const Double_t det5_01245_12345 = pM[GA01]*det4_1245_2345 - pM[GA02]*det4_1245_1345 + pM[GA03]*det4_1245_1245 - pM[GA04]*det4_1245_1235 + pM[GA05]*det4_1245_1234; const Double_t det5_01345_01234 = pM[GA00]*det4_1345_1234 - pM[GA01]*det4_1345_0234 + pM[GA02]*det4_1345_0134 - pM[GA03]*det4_1345_0124 + pM[GA04]*det4_1345_0123; const Double_t det5_01345_01235 = pM[GA00]*det4_1345_1235 - pM[GA01]*det4_1345_0235 + pM[GA02]*det4_1345_0135 - pM[GA03]*det4_1345_0125 + pM[GA05]*det4_1345_0123; const Double_t det5_01345_01245 = pM[GA00]*det4_1345_1245 - pM[GA01]*det4_1345_0245 + pM[GA02]*det4_1345_0145 - pM[GA04]*det4_1345_0125 + pM[GA05]*det4_1345_0124; const Double_t det5_01345_01345 = pM[GA00]*det4_1345_1345 - pM[GA01]*det4_1345_0345 + pM[GA03]*det4_1345_0145 - pM[GA04]*det4_1345_0135 + pM[GA05]*det4_1345_0134; const Double_t det5_01345_02345 = pM[GA00]*det4_1345_2345 - pM[GA02]*det4_1345_0345 + pM[GA03]*det4_1345_0245 - pM[GA04]*det4_1345_0235 + pM[GA05]*det4_1345_0234; const Double_t det5_01345_12345 = pM[GA01]*det4_1345_2345 - pM[GA02]*det4_1345_1345 + pM[GA03]*det4_1345_1245 - pM[GA04]*det4_1345_1235 + pM[GA05]*det4_1345_1234; const Double_t det5_02345_01234 = pM[GA00]*det4_2345_1234 - pM[GA01]*det4_2345_0234 + pM[GA02]*det4_2345_0134 - pM[GA03]*det4_2345_0124 + pM[GA04]*det4_2345_0123; const Double_t det5_02345_01235 = pM[GA00]*det4_2345_1235 - pM[GA01]*det4_2345_0235 + pM[GA02]*det4_2345_0135 - pM[GA03]*det4_2345_0125 + pM[GA05]*det4_2345_0123; const Double_t det5_02345_01245 = pM[GA00]*det4_2345_1245 - pM[GA01]*det4_2345_0245 + pM[GA02]*det4_2345_0145 - pM[GA04]*det4_2345_0125 + pM[GA05]*det4_2345_0124; const Double_t det5_02345_01345 = pM[GA00]*det4_2345_1345 - pM[GA01]*det4_2345_0345 + pM[GA03]*det4_2345_0145 - pM[GA04]*det4_2345_0135 + pM[GA05]*det4_2345_0134; const Double_t det5_02345_02345 = pM[GA00]*det4_2345_2345 - pM[GA02]*det4_2345_0345 + pM[GA03]*det4_2345_0245 - pM[GA04]*det4_2345_0235 + pM[GA05]*det4_2345_0234; const Double_t det5_02345_12345 = pM[GA01]*det4_2345_2345 - pM[GA02]*det4_2345_1345 + pM[GA03]*det4_2345_1245 - pM[GA04]*det4_2345_1235 + pM[GA05]*det4_2345_1234; const Double_t det5_12345_01234 = pM[GA10]*det4_2345_1234 - pM[GA11]*det4_2345_0234 + pM[GA12]*det4_2345_0134 - pM[GA13]*det4_2345_0124 + pM[GA14]*det4_2345_0123; const Double_t det5_12345_01235 = pM[GA10]*det4_2345_1235 - pM[GA11]*det4_2345_0235 + pM[GA12]*det4_2345_0135 - pM[GA13]*det4_2345_0125 + pM[GA15]*det4_2345_0123; const Double_t det5_12345_01245 = pM[GA10]*det4_2345_1245 - pM[GA11]*det4_2345_0245 + pM[GA12]*det4_2345_0145 - pM[GA14]*det4_2345_0125 + pM[GA15]*det4_2345_0124; const Double_t det5_12345_01345 = pM[GA10]*det4_2345_1345 - pM[GA11]*det4_2345_0345 + pM[GA13]*det4_2345_0145 - pM[GA14]*det4_2345_0135 + pM[GA15]*det4_2345_0134; const Double_t det5_12345_02345 = pM[GA10]*det4_2345_2345 - pM[GA12]*det4_2345_0345 + pM[GA13]*det4_2345_0245 - pM[GA14]*det4_2345_0235 + pM[GA15]*det4_2345_0234; const Double_t det5_12345_12345 = pM[GA11]*det4_2345_2345 - pM[GA12]*det4_2345_1345 + pM[GA13]*det4_2345_1245 - pM[GA14]*det4_2345_1235 + pM[GA15]*det4_2345_1234; // Find the determinant const Double_t det = pM[GA00]*det5_12345_12345 - pM[GA01]*det5_12345_02345 + pM[GA02]*det5_12345_01345 - pM[GA03]*det5_12345_01245 + pM[GA04]*det5_12345_01235 - pM[GA05]*det5_12345_01234; if (determ) *determ = det; if ( det == 0 ) { Error("Inv6x6","matrix is singular"); return kFALSE; } const Double_t oneOverDet = 1.0/det; const Double_t mn1OverDet = - oneOverDet; pM[GA00] = det5_12345_12345*oneOverDet; pM[GA01] = det5_02345_12345*mn1OverDet; pM[GA02] = det5_01345_12345*oneOverDet; pM[GA03] = det5_01245_12345*mn1OverDet; pM[GA04] = det5_01235_12345*oneOverDet; pM[GA05] = det5_01234_12345*mn1OverDet; pM[GA10] = det5_12345_02345*mn1OverDet; pM[GA11] = det5_02345_02345*oneOverDet; pM[GA12] = det5_01345_02345*mn1OverDet; pM[GA13] = det5_01245_02345*oneOverDet; pM[GA14] = det5_01235_02345*mn1OverDet; pM[GA15] = det5_01234_02345*oneOverDet; pM[GA20] = det5_12345_01345*oneOverDet; pM[GA21] = det5_02345_01345*mn1OverDet; pM[GA22] = det5_01345_01345*oneOverDet; pM[GA23] = det5_01245_01345*mn1OverDet; pM[GA24] = det5_01235_01345*oneOverDet; pM[GA25] = det5_01234_01345*mn1OverDet; pM[GA30] = det5_12345_01245*mn1OverDet; pM[GA31] = det5_02345_01245*oneOverDet; pM[GA32] = det5_01345_01245*mn1OverDet; pM[GA33] = det5_01245_01245*oneOverDet; pM[GA34] = det5_01235_01245*mn1OverDet; pM[GA35] = det5_01234_01245*oneOverDet; pM[GA40] = det5_12345_01235*oneOverDet; pM[GA41] = det5_02345_01235*mn1OverDet; pM[GA42] = det5_01345_01235*oneOverDet; pM[GA43] = det5_01245_01235*mn1OverDet; pM[GA44] = det5_01235_01235*oneOverDet; pM[GA45] = det5_01234_01235*mn1OverDet; pM[GA50] = det5_12345_01234*mn1OverDet; pM[GA51] = det5_02345_01234*oneOverDet; pM[GA52] = det5_01345_01234*mn1OverDet; pM[GA53] = det5_01245_01234*oneOverDet; pM[GA54] = det5_01235_01234*mn1OverDet; pM[GA55] = det5_01234_01234*oneOverDet; return kTRUE; } #ifndef ROOT_TMatrixFfwd #include "TMatrixFfwd.h" #endif template Bool_t TMatrixTCramerInv::Inv2x2(TMatrixF&,Double_t*); template Bool_t TMatrixTCramerInv::Inv3x3(TMatrixF&,Double_t*); template Bool_t TMatrixTCramerInv::Inv4x4(TMatrixF&,Double_t*); template Bool_t TMatrixTCramerInv::Inv5x5(TMatrixF&,Double_t*); template Bool_t TMatrixTCramerInv::Inv6x6(TMatrixF&,Double_t*); #ifndef ROOT_TMatrixDfwd #include "TMatrixDfwd.h" #endif template Bool_t TMatrixTCramerInv::Inv2x2(TMatrixD&,Double_t*); template Bool_t TMatrixTCramerInv::Inv3x3(TMatrixD&,Double_t*); template Bool_t TMatrixTCramerInv::Inv4x4(TMatrixD&,Double_t*); template Bool_t TMatrixTCramerInv::Inv5x5(TMatrixD&,Double_t*); template Bool_t TMatrixTCramerInv::Inv6x6(TMatrixD&,Double_t*);