1 package org.lcsim.recon.tracking.trfbase;
2
3 import org.lcsim.recon.tracking.trfutil.TRFMath;
4 import Jama.Matrix;
5 /**
6 * Class TrackError is the corresponding symmetric 5x5 error matrix
7 * for TrackVector.
8 *<p>
9 * TrackError should support the following operations:
10 * <pre>
11 * construct: tv(), te(), td() (should fill with zeros)
12 * copy: tv2(tv1), te2(te1), td2(td1)
13 * assignment: tv2 = tv2, te2 = te1, td2 = td1
14 * indexing: tv(i), te(i,j), td(i,j)
15 * minimum element: tv.min(), te.min(), td.min()
16 * maximum element: tv.max(), te.max(), td.max()
17 * absolute minimum element: tv.amin(), te.amin(), td.amin()
18 * absolute maximum element: tv.amax(), te.amax(), td.amax()
19 * addition: tv2 += tv1, tv3 = tv1 + tv2,
20 * te2 += te1, te3 = tv1 + te2,
21 * td2 += td1, td3 = td1 + td2,
22 * subtraction: tv2 -= tv1, tv3 = tv1 - tv2,
23 * te2 -= te1, te3 = tv1 - te2,
24 * td2 -= td1, td3 = td1 - td2,
25 * inversion: int stat = invert( TrackError& te1 ) (0 for success)
26 * output stream: cout << tv << te << td
27 * equality: tv1==tv2, te1==te2, td1==td2
28 * inequality: tv!=tv2, te1!=te2, td1!=td2
29 * chisquare difference: chisq_diff(tv,te) = tv_T * te * tv
30 * (tv is the difference between two vectors and te is the inverse
31 * of the error matrix.)
32 * transpose: td.transpose()
33 * transform: te.xform(td)
34 * </pre>
35 *
36 * @author Norman A. Graf
37 * @version 1.0
38 */
39 public class TrackError
40 {
41 private Matrix _mat;
42 private int _size;
43
44 //
45
46 /**
47 * Default constructor creates 5x5 matrix of zeros.
48 *
49 */
50 public TrackError( )
51 {
52 _size = 5;
53 _mat = new Matrix(5,5);
54 }
55
56 //
57
58 /**
59 * Constructor from array.
60 *
61 * @param mat 5x5 array
62 */
63 public TrackError(double[][] mat)
64 {
65 _size = mat.length;
66 for (int i = 0; i < _size; i++)
67 {
68 if (mat[i].length != _size)
69 {
70 throw new IllegalArgumentException("All rows must have the same length.");
71 }
72 }
73 _mat = new Matrix(mat, _size, _size);
74
75 }
76
77 //
78
79 /**
80 * Constructor from Jama Matrix.
81 *
82 * @param mat 5x5 Jama Matrix of doubles representing the Track error matrix.
83 */
84 public TrackError( Matrix mat)
85 {
86 if(mat.getRowDimension()!=5 && mat.getColumnDimension() != 5)
87 {
88 throw new IllegalArgumentException("TrackError Matrix must be 5x5.");
89 }
90 _size = mat.getRowDimension();
91 _mat = new Matrix(mat.getArrayCopy());
92 }
93
94 //
95
96 /**
97 * Convenience method to set the error matrix to the identity matrix.
98 *
99 */
100 public void setIdentity()
101 {
102 _mat = Matrix.identity(_size, _size);
103 }
104
105 //
106
107 /**
108 * Copy Constructor
109 *
110 * @param te TrackError to copy.
111 */
112 public TrackError( TrackError te)
113 {
114 _size = 5;
115 _mat = te.getMatrix();
116 }
117
118 //
119
120 /**
121 * Return the underlying vector.
122 *
123 * @return 5x5 double array.
124 */
125 public double[][] matrix()
126 { //need to fix this!
127 return _mat.getArrayCopy();
128 }
129
130
131 /**
132 * String representation of TrackError.
133 *
134 * @return String representation of TrackError.
135 */
136 public String toString()
137 {
138 String className = getClass().getName();
139 int lastDot = className.lastIndexOf('.');
140 if(lastDot!=-1)className = className.substring(lastDot+1);
141
142 StringBuffer sb = new StringBuffer(className+"\n");
143 for(int i=0; i<_size; ++i)
144 {
145 for(int j = 0; j<=i; ++j)
146 {
147 sb.append(_mat.get(i,j)).append(" ");
148 }
149 sb.append("\n");
150 }
151 return sb.append("\n").toString();
152 }
153
154 /**
155 * Get the track error matrix.
156 *
157 * @return Jama Matrix.
158 */
159 public Matrix getMatrix()
160 {
161 return _mat.copy();
162 }
163
164
165 /**
166 * Set an element of the error matrix. This should be used
167 * with caution, since the error matrix should have qualities
168 * which may not be preserved under arbitrary element insertion.
169 * Inserting an element (i,j) will also set the element (j,i), since
170 * the matrix has to be at least symmetric.
171 * @param i row index.
172 * @param j column index.
173 * @param val value for element (i,j).
174 */
175 public void set(int i, int j, double val)
176 {
177 // This is a symmetric matrix, so always set the
178 // other element as well.
179 _mat.set(i,j,val);
180 if(i!=j) _mat.set(j,i,val);
181 }
182
183
184 /**
185 * Get an element of the error matrix.
186 *
187 * @param i row index.
188 * @param j column index.
189 * @return value of element (i,j).
190 */
191 public double get(int i, int j)
192 {
193 return _mat.get(i, j);
194 }
195
196 //cng what is the signature of this method?
197 // public void Xform( TrackDerivative deriv)
198 // {
199 // _mat = deriv.getMatrix().times(getMatrix().times(deriv.getMatrix().transpose()));
200 // }
201
202
203 /**
204 * TrackDerivative * TrackError * TrackDerivative(transpose).
205 *
206 * @param deriv TrackDerivative by which to transform the TrackError .
207 * @return TrackError transformed by Trackderivative.
208 */
209 public TrackError Xform( TrackDerivative deriv)
210 {
211 return new TrackError(deriv.getMatrix().times(getMatrix().times(deriv.getMatrix().transpose())));
212 }
213
214
215 /**
216 * Math function minus.
217 *
218 * @param te TrackError to be subtracted.
219 * @return <em>this</em> minus TrackError te.
220 */
221 public TrackError minus(TrackError te)
222 {
223 return new TrackError(getMatrix().minus(te.getMatrix()));
224 }
225
226
227 /**
228 *Math function plus.
229 *
230 * @param te TrackError to be added.
231 * @return <em>this</em> plus TrackError te.
232 */
233 public TrackError plus(TrackError te)
234 {
235 return new TrackError(getMatrix().plus(te.getMatrix()));
236 }
237
238
239 /**
240 * Math function multiply.
241 *
242 * @param te TrackError to be multiplied by.
243 * @return <em>this</em> times TrackError te.
244 */
245 public TrackError times(TrackError te)
246 {
247 return new TrackError( _mat.times(te._mat) );
248 }
249
250
251 /**
252 * Matrix inverse.
253 *
254 * @return Matrix inverse of <em>this</em>.
255 */
256 public TrackError inverse( )
257 {
258 return new TrackError(getMatrix().inverse());
259 }
260
261
262 /**
263 * Maximum TrackError element.
264 *
265 * @return value of maximum TrackError element.
266 */
267 public double max()
268 {
269 double[][] tmp = matrix();
270 double max = tmp[0][0];
271 for(int i = 0; i<5; ++i)
272 {
273 for(int j = 0; j<5; ++j)
274 {
275 if( tmp[i][j]>max) max = tmp[i][j];
276 }
277 }
278 return max;
279 }
280
281
282 /**
283 * Absolute maximum TrackError element.
284 *
285 * @return absolute value of absolute maximum TrackError element.
286 */
287 public double amax()
288 {
289 double[][] tmp = matrix();
290 double amax = Math.abs(tmp[0][0]);
291 for(int i = 0; i<5; ++i)
292 {
293 for(int j = 0; j<5; ++j)
294 {
295 if( Math.abs(tmp[i][j])>amax) amax = Math.abs(tmp[i][j]);
296 }
297 }
298 return amax;
299 }
300
301
302 /**
303 * Minimum TrackError element.
304 *
305 * @return value of minimum TrackError element.
306 */
307 public double min()
308 {
309 double[][] tmp = matrix();
310 double min = tmp[0][0];
311 for(int i = 0; i<5; ++i)
312 {
313 for(int j = 0; j<5; ++j)
314 {
315 if( tmp[i][j]<min) min = tmp[i][j];
316 }
317 }
318 return min;
319 }
320
321
322 /**
323 * Absolute minimum TrackError element.
324 *
325 * @return absolute value of absolute minimum TrackError element.
326 */
327 public double amin()
328 {
329 double[][] tmp = matrix();
330 double amin = Math.abs(tmp[0][0]);
331 for(int i = 0; i<5; ++i)
332 {
333 for(int j = 0; j<5; ++j)
334 {
335 if( Math.abs(tmp[i][j])<amin) amin = Math.abs(tmp[i][j]);
336 }
337 }
338 return amin;
339 }
340
341
342
343 /**
344 * Equality.
345 *
346 * @param te TrackError to check equality against.
347 * @return true if TrackErrors are equal.
348 */
349 public boolean equals( TrackError te)
350 {
351 if( _size != te._size ) return false;
352 for (int i = 0 ;i < _size ; ++i )
353 {
354 for (int j = 0; j < _size ; ++j )
355 {
356 if( _mat.get(i,j) != te._mat.get(i,j) ) return false;
357 }
358 }
359 return true;
360 }
361
362
363 /**
364 * Inequality convenience method.
365 *
366 * @param te TrackError to check equality against.
367 * @return true if TrackErrors are <b> not </b> equal.
368 */
369 public boolean notEquals(TrackError te)
370 {
371 return !equals(te);
372 }
373
374
375 /**
376 * Equality within tolerance.
377 *
378 * @param te TrackError to check equality against.
379 * @return true if TrackErrors are equal within tolerance.
380 */
381 public boolean isEqual( TrackError te)
382 {
383 if( _size != te._size ) return false;
384 for (int i = 0 ;i < _size ; ++i )
385 {
386 for (int j = 0; j < _size ; ++j )
387 {
388 if( !TRFMath.isEqual(_mat.get(i,j), te._mat.get(i,j)) ) return false;
389 }
390 }
391 return true;
392 }
393
394
395 /**
396 * Invert matrix in place.
397 *
398 * @return nonzero integer if error.
399 */
400 public int invert()
401 {
402 int stat = 0;
403 _mat = _mat.inverse();
404 return stat;
405 }
406
407 /**
408 * Normalize a matrix.
409 * Transformation which makes all diagonal elements unity.
410 *
411 */
412 public void normalize()
413 {
414
415 double[] vec = new double[_size];
416 for(int irow=0; irow<_size; irow++ )
417 {
418 vec[irow] = 1.0/Math.sqrt( _mat.get(irow,irow) );
419 double rowfac = Math.sqrt( vec[irow] );
420 for ( int icol=0; icol<=irow; icol++ )
421 {
422 double tmp = _mat.get(irow,icol);
423 tmp *= vec[irow]*vec[icol];
424 _mat.set(irow,icol, tmp);
425 _mat.set(icol,irow, tmp);
426 }
427 }
428 }
429 }