View Javadoc

1   /*
2    * InterpolatedHMatrix.java
3    *
4    * Created on June 5, 2008, 10:34 AM
5    *
6    * $Id: InterpolatedHMatrix.java,v 1.1 2008/06/06 07:27:32 ngraf Exp $
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   * This class interpolates the average values and inverse covariance matrix elements
16   * for photon showers based on derivations at fixed energy and angle.
17   * @author Norman Graf
18   */
19  public class InterpolatedHMatrix
20  {
21      //Interpolation for the vector of measurement averages
22      private BilinearInterpolator[] _valInterpolators;
23      //Interpolation for the inverse covariance matrix elements represented as a
24      //packed lower-diagonal matrix
25      private BilinearInterpolator[] _covInterpolators;
26      
27      private int _dim;
28      private int _covDim;
29      
30      /** Creates a new instance of InterpolatedHMatrix */
31      // TODO check robustness, use deep copy if necessary.
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       * Returns the closeness of the vector dat to the expectations based on an 
42       * interpolation between discrete values of angle and energy at which showers
43       * were simulated.
44       *
45       * @param theta  the polar angle
46       * @param energy the energy of the electromagnetic shower
47       * @param dat the vector of fractional energies by layer
48       * @return the distance of this shower to that expected, chi-squared for normal distributions of layer energies.
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          // vector of measured-predicted values
68          for(int i=0; i<_dim; ++i)
69          {
70              tmp[i] = dat[i] - means[i];
71          }
72  
73          double zeta = 0.;
74          // reexpand lower-diagonal into full double array
75          // could this be optimized?
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          // covariance matrix times difference vector
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  }