View Javadoc

1   package org.lcsim.mc.fast.tracking;
2   
3   import java.io.BufferedReader;
4   import java.io.IOException;
5   
6   import java.util.Arrays;
7   import java.util.StringTokenizer;
8   
9   public class LookupTable {
10      // cosine theta
11      private double[] m_key1;
12      // momentum
13      private double[] m_key2;
14      private double[][] m_matrix;
15      private int m_numBins1;
16      private int m_numBins2;
17  
18      // vector units
19      // 1 dr [ cm ] 10.
20      // 2 dphi [ ] 1.
21      // 3 domega [ cm-1 ] .1
22      // 4 dz [ cm ] 10.
23      // 5 dlambda [ ] 1.
24      //
25      double[] conversionFromCmToMm = { 10.0, 1.0, 0.1, 10.0, 1.0 };
26  
27      double[][] conversionFromCmToMmMatrix = { { 100.00, 10.00, 1.00, 100.00, 10.00 }, { 10.00, 1.00, 0.10, 10.00, 1.00 }, { 1.00, 0.10, 0.01, 1.00, 0.10 }, { 100.00, 10.00, 1.00, 100.00, 10.00 }, { 10.00, 1.00, 0.10, 10.00, 1.00 } };
28  
29      LookupTable(BufferedReader in, int iTerm, int jTerm) throws IOException {
30          // read in the number of cosine theta points
31          int m_numBins1 = Integer.parseInt(in.readLine());
32  
33          // read in the number of momentum points
34          int m_numBins2 = Integer.parseInt(in.readLine());
35  
36          m_matrix = new double[m_numBins1][m_numBins2];
37          m_key1 = new double[m_numBins1];
38          m_key2 = new double[m_numBins2];
39  
40          for (int i = 0; i < m_numBins1; i++) // i is # of cosine theta bin
41          {
42              m_key1[i] = Double.valueOf(in.readLine()).doubleValue(); // cosine theta
43              for (int j = 0; j < m_numBins2; j++) // j is # of momentum bin
44              {
45                  StringTokenizer t = new StringTokenizer(in.readLine());
46                  m_key2[j] = Double.valueOf(t.nextToken()).doubleValue();
47                  m_matrix[i][j] = Double.valueOf(t.nextToken()).doubleValue() * conversionFromCmToMmMatrix[iTerm][jTerm]; // momentum
48              }
49          }
50          if (!in.readLine().equals("end")) {
51              throw new IOException("Missing end in lookup table");
52          }
53      }
54  
55      public double interpolateVal(double val1, double val2) {
56          int index1 = binarySearch(m_key1, val1);
57          // implement cut-off
58          double t;
59          if (index1 < 0) {
60              t = m_key1[0];
61              index1 = 0;
62          } else if (index1 >= m_key1.length - 1) {
63              t = m_key1[m_key1.length - 1];
64              index1 = m_key1.length - 1;
65          } else
66              t = (val1 - m_key1[index1]) / (m_key1[index1 + 1] - m_key1[index1]);
67  
68          double u;
69          int index2 = binarySearch(m_key2, val2);
70          if (index2 < 0) {
71              u = m_key2[0];
72              index2 = 0;
73          } else if (index2 >= m_key2.length - 1) {
74              u = m_key2[m_key2.length - 1];
75              index2 = m_key2.length - 1;
76          } else
77              u = (val2 - m_key2[index2]) / (m_key2[index2 + 1] - m_key2[index2]);
78  
79          double y1 = m_matrix[index1][index2];
80          double y2 = m_matrix[index1 + 1][index2];
81          double y3 = m_matrix[index1 + 1][index2 + 1];
82          double y4 = m_matrix[index1][index2 + 1];
83  
84          return ((1 - t) * (1 - u) * y1) + (t * (1 - u) * y2) + (t * u * y3) + ((1 - t) * u * y4);
85      }
86  
87      private int binarySearch(double[] key, double value)
88      // {
89      // int result = binarySearchX(key,value);
90      // System.out.print("Looking for "+value+" in [");
91      // for (int i=0; i<key.length; i++) System.out.print(key[i]+",");
92      // System.out.print("] ");
93      // System.out.println("result="+result);
94      // return result;
95      // }
96      //
97      // private int binarySearchX(double[] key, double value)
98      {
99          if (value < key[0]) {
100             // throw new RuntimeException("Interpolation out of range: lower: "+value+" < "+key[0]);
101             return 0;
102         }
103 
104         int pos = Arrays.binarySearch(key, value);
105         if (pos > 0) {
106             return Math.min(pos, key.length - 2);
107         } else {
108             return Math.min(-pos - 2, key.length - 2);
109         }
110 
111         // // Ok, this isn't really a binary search, probably doesn't matter
112         // for (int i=1; i<key.length; i++) if (value<key[i]) return i-1;
113         //
114         //
115         // throw new LCDException("Interpolation out of range: upper: "
116         // +value+" >= "+key[key.length-1]);
117     }
118 }