1
2
3
4
5
6
7
8
9 package org.lcsim.recon.emid.hmatrix;
10
11 import hep.physics.matrix.SymmetricMatrix;
12 import org.lcsim.math.interpolation.BilinearInterpolator;
13
14
15
16
17
18
19 public class InterpolatedHMatrix
20 {
21
22 private BilinearInterpolator[] _valInterpolators;
23
24
25 private BilinearInterpolator[] _covInterpolators;
26
27 private int _dim;
28 private int _covDim;
29
30
31
32 public InterpolatedHMatrix(int dim, BilinearInterpolator[] valInterpolators, BilinearInterpolator[] covInterpolators)
33 {
34 _dim = dim;
35 _covDim = (_dim*(_dim+1))/2;
36 _valInterpolators = valInterpolators;
37 _covInterpolators = covInterpolators;
38 }
39
40
41
42
43
44
45
46
47
48
49
50 public double zeta(double theta, double energy, double[] dat)
51 {
52 double[] means = new double[_dim];
53 double[] cov = new double[_covDim];
54 double[] tmp = new double[_dim];
55 double[] tmp2 = new double[_dim];
56
57 for(int i=0; i<_dim; ++i)
58 {
59 means[i] = _valInterpolators[i].interpolateValueAt(theta, energy);
60 }
61
62 for(int i=0; i<_covDim; ++i)
63 {
64 cov[i] = _covInterpolators[i].interpolateValueAt(theta, energy);
65 }
66
67
68 for(int i=0; i<_dim; ++i)
69 {
70 tmp[i] = dat[i] - means[i];
71 }
72
73 double zeta = 0.;
74
75
76 double[][] invcov = new double[_dim][_dim];
77 int count = 0;
78 for(int i=0; i<_dim; ++i)
79 {
80 for(int j=0; j<i+1; ++j)
81 {
82 invcov[i][j] = cov[count];
83 invcov[j][i] = cov[count];
84 count++;
85 }
86 }
87
88 for(int i=0; i<_dim; ++i)
89 {
90 tmp2[i] = 0.;
91 for(int j=0; j<_dim; ++j)
92 {
93 tmp2[i]+=invcov[j][i]*tmp[j];
94 }
95 zeta += tmp[i]*tmp2[i];
96 }
97 return zeta;
98 }
99 }