1 package org.lcsim.event.base;
2
3 import java.util.List;
4
5 import org.lcsim.event.CalorimeterHit;
6 import org.lcsim.event.Cluster;
7
8
9
10
11
12
13
14
15 public class TensorClusterPropertyCalculator extends AbstractClusterPropertyCalculator {
16
17 protected double[] mm_NE;
18 protected double[] mm_CE;
19 protected double[][] mm_PA;
20
21
22
23
24 public TensorClusterPropertyCalculator() {
25 }
26
27 protected void reset() {
28 super.reset();
29 mm_NE = new double[3];
30 mm_CE = new double[3];
31 mm_PA = new double[3][3];
32 }
33
34 public void calculateProperties(Cluster cluster) {
35 calculateProperties(cluster.getCalorimeterHits());
36 }
37
38 public void calculateProperties(List<CalorimeterHit> hits) {
39
40 reset();
41
42 for (int i = 0; i < 3; ++i) {
43 mm_NE[i] = 0.;
44 mm_CE[i] = 0.;
45 for (int j = 0; j < 3; ++j) {
46 mm_PA[i][j] = 0.;
47 }
48 }
49 double Etot = 0.0;
50 double Exx = 0.0;
51 double Eyy = 0.0;
52 double Ezz = 0.0;
53 double Exy = 0.0;
54 double Eyz = 0.0;
55 double Exz = 0.0;
56 double CEx = 0.0;
57 double CEy = 0.0;
58 double CEz = 0.0;
59 double CEr = 0.0;
60 double E1 = 0.0;
61 double E2 = 0.0;
62 double E3 = 0.0;
63 double NE1 = 0.0;
64 double NE2 = 0.0;
65 double NE3 = 0.0;
66 double Tr = 0.0;
67 double M = 0.0;
68 double Det = 0.0;
69 int nhits = hits.size();
70 for (int i = 0; i < hits.size(); i++) {
71 CalorimeterHit hit = hits.get(i);
72
73
74
75
76
77
78 double[] pos = hit.getPosition();
79 double E = hit.getCorrectedEnergy();
80 Etot += E;
81 CEx += E * pos[0];
82 CEy += E * pos[1];
83 CEz += E * pos[2];
84 Exx += E * (pos[1] * pos[1] + pos[2] * pos[2]);
85 Eyy += E * (pos[0] * pos[0] + pos[2] * pos[2]);
86 Ezz += E * (pos[1] * pos[1] + pos[0] * pos[0]);
87 Exy -= E * pos[0] * pos[1];
88 Eyz -= E * pos[1] * pos[2];
89 Exz -= E * pos[0] * pos[2];
90 }
91 CEx = CEx / Etot;
92 CEy = CEy / Etot;
93 CEz = CEz / Etot;
94 double CErSq = CEx * CEx + CEy * CEy + CEz * CEz;
95 CEr = Math.sqrt(CErSq);
96
97 if (nhits > 3) {
98 Exx = Exx - Etot * (CErSq - CEx * CEx);
99 Eyy = Eyy - Etot * (CErSq - CEy * CEy);
100 Ezz = Ezz - Etot * (CErSq - CEz * CEz);
101 Exy = Exy + Etot * CEx * CEy;
102 Eyz = Eyz + Etot * CEy * CEz;
103 Exz = Exz + Etot * CEz * CEx;
104
105
106 Tr = Exx + Eyy + Ezz;
107 double Dxx = Eyy * Ezz - Eyz * Eyz;
108 double Dyy = Ezz * Exx - Exz * Exz;
109 double Dzz = Exx * Eyy - Exy * Exy;
110 M = Dxx + Dyy + Dzz;
111 double Dxy = Exy * Ezz - Exz * Eyz;
112 double Dxz = Exy * Eyz - Exz * Eyy;
113 Det = Exx * Dxx - Exy * Dxy + Exz * Dxz;
114 double xt = Tr * Tr - 3 * M;
115 double sqrtxt = Math.sqrt(xt);
116
117
118
119 double mE1 = 0.;
120 double mE2 = 0.;
121 double mE3 = 0.;
122 double a = (3 * M - Tr * Tr) / 3.;
123 double b = (-2 * Tr * Tr * Tr + 9 * Tr * M - 27 * Det) / 27.;
124 double test = b * b / 4. + a * a * a / 27.;
125 if (test >= 0.01) {
126
127
128
129 } else {
130 double temp;
131 if (test >= 0.)
132 temp = 1.;
133 else
134 temp = Math.sqrt(b * b * 27. / (-a * a * a * 4.));
135 if (b > 0.)
136 temp = -temp;
137 double phi = Math.acos(temp);
138 double temp1 = 2. * Math.sqrt(-a / 3.);
139 mE1 = Tr / 3. + 2. * Math.sqrt(-a / 3.) * Math.cos(phi / 3.);
140 mE2 = Tr / 3. + 2. * Math.sqrt(-a / 3.) * Math.cos(phi / 3. + 2. * Math.PI / 3.);
141 mE3 = Tr / 3. + 2. * Math.sqrt(-a / 3.) * Math.cos(phi / 3. + 4. * Math.PI / 3.);
142 }
143 if (mE1 < mE2) {
144 if (mE1 < mE3) {
145 E1 = mE1;
146 if (mE2 < mE3) {
147 E2 = mE2;
148 E3 = mE3;
149 } else {
150 E2 = mE3;
151 E3 = mE2;
152 }
153 } else {
154 E1 = mE3;
155 E2 = mE1;
156 E3 = mE2;
157 }
158 } else {
159 if (mE2 < mE3) {
160 E1 = mE2;
161 if (mE1 < mE3) {
162 E2 = mE1;
163 E3 = mE3;
164 } else {
165 E2 = mE3;
166 E3 = mE1;
167 }
168 } else {
169 E1 = mE3;
170 E2 = mE2;
171 E3 = mE1;
172 }
173 }
174
175 NE1 = E1 / Etot;
176 NE2 = E2 / Etot;
177 NE3 = E3 / Etot;
178 double[] EV = new double[3];
179 EV[0] = E1;
180 EV[1] = E2;
181 EV[2] = E3;
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222 if (Exz == 0.0 && Eyz == 0.0) {
223
224 EV[0] = 0.5 * (Exx + Eyy) + 0.5 * Math.sqrt((Exx + Eyy) * (Exx + Eyy) + 4.0 * Exy * Exy);
225 EV[1] = 0.5 * (Exx + Eyy) - 0.5 * Math.sqrt((Exx + Eyy) * (Exx + Eyy) + 4.0 * Exy * Exy);
226 EV[2] = 0.0;
227 for (int i = 0; i < 2; i++) {
228 double nx_over_ny = Exy / (Exx - EV[i]);
229 double nx_unnormalized = nx_over_ny;
230 double ny_unnormalized = 1.0;
231 double norm = Math.sqrt(nx_unnormalized * nx_unnormalized + ny_unnormalized * ny_unnormalized);
232 mm_PA[i][0] = ny_unnormalized / norm;
233 mm_PA[i][1] = nx_unnormalized / norm;
234 mm_PA[i][2] = 0.0;
235 }
236
237 mm_PA[2][0] = 0.0;
238 mm_PA[2][1] = 0.0;
239 mm_PA[2][2] = 1.0;
240 } else if (Exz == 0.0 && Exy == 0.0) {
241
242 EV[0] = 0.5 * (Eyy + Ezz) + 0.5 * Math.sqrt((Eyy + Ezz) * (Eyy + Ezz) + 4.0 * Eyz * Eyz);
243 EV[1] = 0.5 * (Eyy + Ezz) - 0.5 * Math.sqrt((Eyy + Ezz) * (Eyy + Ezz) + 4.0 * Eyz * Eyz);
244 EV[2] = 0.0;
245 for (int i = 0; i < 2; i++) {
246 double ny_over_nz = Eyz / (Eyy - EV[i]);
247 double ny_unnormalized = ny_over_nz;
248 double nz_unnormalized = 1.0;
249 double norm = Math.sqrt(ny_unnormalized * ny_unnormalized + nz_unnormalized * nz_unnormalized);
250 mm_PA[i][0] = nz_unnormalized / norm;
251 mm_PA[i][1] = ny_unnormalized / norm;
252 mm_PA[i][2] = 0.0;
253 }
254 mm_PA[2][0] = 0.0;
255 mm_PA[2][1] = 0.0;
256 mm_PA[2][2] = 1.0;
257 } else {
258 for (int i = 0; i < 3; i++) {
259 double[] C = new double[3];
260 C[0] = 1.0;
261 C[2] = (Exy * Exy + (Eyy - EV[i]) * (EV[i] - Exx)) / ((Eyy - EV[i]) * Exz - Eyz * Exy);
262 C[1] = (EV[i] - Exx - Exz * C[2]) / Exy;
263 if (Exy == 0.0) {
264
265 if (Eyz != 0.0) {
266
267 C[1] = (Exz - (Ezz - EV[i]) * C[2]) / Eyz;
268 } else {
269
270 C[1] = -(Eyz * C[2]) / (Eyy - EV[i]);
271 }
272 }
273 double norm = Math.sqrt(C[0] * C[0] + C[1] * C[1] + C[2] * C[2]);
274 mm_PA[i][0] = C[0] / norm;
275 mm_PA[i][1] = C[1] / norm;
276 mm_PA[i][2] = C[2] / norm;
277 }
278 }
279 }
280 mm_NE[0] = NE1;
281 mm_NE[1] = NE2;
282 mm_NE[2] = NE3;
283 mm_CE[0] = CEx;
284 mm_CE[1] = CEy;
285 mm_CE[2] = CEz;
286 for (int i = 0; i < 6; i++) {
287 positionError[i] = 0.;
288 directionError[i] = 0.;
289 shapeParameters[i] = 0.;
290 }
291 for (int i = 0; i < 3; i++) {
292 position[i] = mm_CE[i];
293 shapeParameters[i] = mm_NE[i];
294 }
295 if (nhits > 3) {
296 double dr = Math.sqrt((position[0] + mm_PA[0][0]) * (position[0] + mm_PA[0][0]) + (position[1] + mm_PA[0][1]) * (position[1] + mm_PA[0][1]) + (position[2] + mm_PA[0][2]) * (position[2] + mm_PA[0][2])) - Math.sqrt((position[0]) * (position[0]) + (position[1]) * (position[1]) + (position[2]) * (position[2]));
297 double sign = 1.;
298 if (dr < 0.)
299 sign = -1.;
300 itheta = Math.acos(sign * mm_PA[0][2]);
301 iphi = Math.atan2(sign * mm_PA[0][1], sign * mm_PA[0][0]);
302 } else {
303 itheta = 999.;
304 iphi = 999.;
305 }
306 }
307
308 public double[] getPosition() {
309 return position;
310 }
311
312 public double[] getPositionError() {
313 return positionError;
314 }
315
316 public double getIPhi() {
317 return iphi;
318 }
319
320 public double getITheta() {
321 return itheta;
322 }
323
324 public double[] getDirectionError() {
325 return directionError;
326 }
327
328 public double[] getShapeParameters() {
329 return shapeParameters;
330 }
331
332 public double[] getNormalizedEigenvalues() {
333 return mm_NE;
334 }
335
336 public double[][] getPrincipleAxis() {
337 return mm_PA;
338 }
339
340 }