View Javadoc

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   * Implementation of ClusterPropertyCalculator Interface using inertia Tensor calculation 
11   * to derive parameters.
12   * 
13   * @author cassell
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       * Calculate the cluster properties from the hits
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              // CalorimeterIDDecoder decoder = hit.getDecoder();
73              // decoder.setID(hit.getCellID());
74              // double[] pos = new double[3];
75              // pos[0] = decoder.getX();
76              // pos[1] = decoder.getY();
77              // pos[2] = decoder.getZ();
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          // now go to center of energy coords.
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             // eqn to solve for eigenvalues is x**3 - x**2*Tr + x*M - Det = 0
117             // crosses y axis at -Det and inflection points are
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                 // System.out.println("AbstractCluster: Only 1 real root!!!");
127                 // System.out.println("  nhits = " + nhits + "\n");
128                 // System.out.println(" a,b,test = " + a + " " + b + " " + test + "\n");
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             // Now calculate principal axes
183             // For eigenvalue EV, the axis is (nx, ny, nz) where:
184             // (Exx - EV)nx + (Exy)ny + (Exz)nz = 0
185             // (Eyx)nx + (Eyy - EV)ny + (Eyz)nz = 0
186             // (Ezx)nx + (Ezy)ny + (Ezz - EV)nz = 0
187             // Setting nx = 1, we have:
188             // (Exx - EV) + (Exy)ny + (Exz)nz = 0
189             // (Eyx) + (Eyy - EV)ny + (Eyz)nz = 0
190             // (Ezx) + (Ezy)ny + (Ezz - EV)nz = 0
191             // and so
192             // (Exy)ny = EV - Exx - (Exz)nz => ny = (EV - Exx - Exz*nz)/Exy
193             // What if Exy = 0? Then provided Eyz is non-zero we can write:
194             // (Ezx) + (Ezy)ny + (Ezz - EV)nz = 0
195             // ny = (Exz - (Ezz-EV)*nz)/Eyz
196             // What if Exy = 0 and Eyz = 0 but (Eyy - EV) is non-zero?
197             // (Eyy - EV)ny + (Eyz)nz = 0
198             // ny = -(Eyz*nz)/(Eyy-EV)
199 
200             // In the pathological case where Exz = Eyz = Ezz = 0:
201             // (Exx - EV)nx + (Exy)ny = 0 => ny/nx = -(Exx-EV)/Exy
202             // (Eyx)nx + (Eyy - EV)ny = 0 => ny/nx = -Eyx/(Eyy-EV)
203             // (EV)nz = 0
204             // so
205             // -ny/nx = (EV-Exx)/Exy = Eyx/(EV-Eyy)
206             // But watch out for order! Recalculate eigenvalues for this pathological case.
207             // (EV-Exx)(EV-Eyy) = Eyx*Exy
208             // EV^2 - EV(Exx+Eyy) + Exx*Eyy - Eyx*Exy = 0
209             //
210             // In another pathological case, Exz = Exy = 0:
211             // (Exx - EV)nx = 0
212             // (Eyy - EV)ny + (Eyz)nz = 0 => ny/nz = -(Eyz)/(Eyy-EV)
213             // (Ezy)ny + (Ezz - EV)nz = 0 => ny/nz = -(Ezz-EV)/(Ezy)
214             // so we cannot set nx = 1. Instead, write:
215             // -ny/nz = (Eyz)/(Eyy-EV) = (Ezz-EV)/(Ezy)
216             // Then
217             // (Eyz)(Ezy) = (Eyy-EV)(Ezz-EV)
218             // (Eyz)^2 = (Eyy)(Ezz) - (Eyy)(EV) - (Ezz)(EV) + (EV)^2
219             // EV^2 - EV(Eyy+Ezz) + Eyy*Ezz - Eyz*Eyz = 0
220 
221             // Handle pathological case
222             if (Exz == 0.0 && Eyz == 0.0) {
223                 // Recompute eigenvectors.
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                 // ... and now set third eigenvector to the z direction:
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                 // Another pathological case
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                         // Recompute
265                         if (Eyz != 0.0) {
266                             // ny = (Exz - (Ezz-EV)*nz)/Eyz
267                             C[1] = (Exz - (Ezz - EV[i]) * C[2]) / Eyz;
268                         } else {
269                             // ny = -(Eyz*nz)/(Eyy-EV)
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 }