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
11 private double[] m_key1;
12
13 private double[] m_key2;
14 private double[][] m_matrix;
15 private int m_numBins1;
16 private int m_numBins2;
17
18
19
20
21
22
23
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
31 int m_numBins1 = Integer.parseInt(in.readLine());
32
33
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++)
41 {
42 m_key1[i] = Double.valueOf(in.readLine()).doubleValue();
43 for (int j = 0; j < m_numBins2; j++)
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];
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
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
90
91
92
93
94
95
96
97
98 {
99 if (value < key[0]) {
100
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
112
113
114
115
116
117 }
118 }