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