1 /*
2 * StripHit2DMaker.java
3 *
4 * Created on December 20, 2007, 11:39 PM
5 *
6 * To change this template, choose Tools | Template Manager
7 * and open the template in the editor.
8 */
9
10 package org.lcsim.recon.tracking.digitization.sisim;
11
12 import hep.physics.matrix.SymmetricMatrix;
13
14 import java.util.ArrayList;
15 import java.util.HashMap;
16 import java.util.List;
17 import java.util.Map;
18
19 import org.lcsim.detector.IDetectorElement;
20 import org.lcsim.detector.IReadout;
21 import org.lcsim.detector.solids.GeomOp3D;
22 import org.lcsim.detector.tracker.silicon.SiSensor;
23 import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
24 import org.lcsim.detector.tracker.silicon.SiTrackerModule;
25 import org.lcsim.event.RawTrackerHit;
26 import org.lcsim.event.SimTrackerHit;
27
28 /**
29 *
30 * @author tknelson
31 */
32 public class StripHit2DMaker implements StripHitCombiner
33 {
34
35 private static String _NAME = "RawTrackerHitMaker";
36
37 /**
38 * Creates a new instance of StripHit2DMaker
39 */
40 public StripHit2DMaker()
41 {
42 }
43
44 public String getName()
45 {
46 return _NAME;
47 }
48
49 public List<SiTrackerHitStrip2D> makeHits(IDetectorElement detector)
50 {
51
52 List<SiTrackerHitStrip2D> hits_2D = new ArrayList<SiTrackerHitStrip2D>();
53
54 List<SiTrackerModule> modules = detector.findDescendants(SiTrackerModule.class);
55
56 // Loop over all modules
57 for (SiTrackerModule module : modules)
58 {
59 List<SiSensor> sensors = module.findDescendants(SiSensor.class);
60
61 boolean doublesided_sensors = false;
62 boolean pixel_sensors = false;
63
64 // Process double-sided sensors first
65 for (SiSensor sensor : sensors)
66 {
67 if (sensor.isDoubleSided())
68 {
69 doublesided_sensors = true;
70 // if (sensor.hasPixels())
71 // {
72 // pixel_sensors = true;
73 // }
74 // {
75 hits_2D.addAll(this.makeHits(sensor));
76 // }
77 }
78 }
79
80 // If have two sensors and neither is double-sided then make 2d hits from single-sided sensor hits
81 if (sensors.size() == 2 && !doublesided_sensors && !pixel_sensors)
82 {
83 hits_2D.addAll(this.makeHits(module));
84 }
85 }
86 return hits_2D;
87 }
88
89 public List<SiTrackerHitStrip2D> makeHits(SiSensor sensor)
90 {
91 List<SiTrackerHitStrip2D> hits_2D = new ArrayList<SiTrackerHitStrip2D>();
92
93 if (sensor.isDoubleSided() && !sensor.hasPixels())
94 {
95 // Get hits for this sensor
96 IReadout ro = sensor.getReadout();
97 List<SiTrackerHitStrip1D> hits = ro.getHits(SiTrackerHitStrip1D.class);
98
99 // Make map from electrodes to lists of hits
100 Map<SiSensorElectrodes,List<SiTrackerHitStrip1D>> side_hitlists = new HashMap<SiSensorElectrodes,List<SiTrackerHitStrip1D>>();
101 for (SiTrackerHitStrip1D hit : hits)
102 {
103 SiSensorElectrodes electrodes = hit.getReadoutElectrodes();
104 if (side_hitlists.get(electrodes) == null)
105 {
106 side_hitlists.put(electrodes,new ArrayList<SiTrackerHitStrip1D>());
107 }
108
109 side_hitlists.get(electrodes).add(hit);
110 }
111
112 // If have hits on both sides, make 2d strip hits
113 if (side_hitlists.keySet().size() == 2)
114 {
115
116 // System.out.println("Found sensor with hits on both sides!");
117
118 List<List<SiTrackerHitStrip1D>> hitlists = new ArrayList<List<SiTrackerHitStrip1D>>(side_hitlists.values());
119
120 for (SiTrackerHitStrip1D hit1: hitlists.get(0))
121 {
122 for (SiTrackerHitStrip1D hit2: hitlists.get(1))
123 {
124 List<SiTrackerHitStrip1D> hitpair = new ArrayList<SiTrackerHitStrip1D>();
125 hitpair.add(hit1);
126 hitpair.add(hit2);
127 hits_2D.add(makeTrackerHit2D(hitpair));
128 }
129 }
130 }
131 }
132
133 return hits_2D;
134
135 }
136
137
138 public List<SiTrackerHitStrip2D> makeHits(SiTrackerModule module)
139 {
140 List<SiTrackerHitStrip2D> hits_2D = new ArrayList<SiTrackerHitStrip2D>();
141 List<SiSensor> sensors = module.findDescendants(SiSensor.class);
142
143 if (sensors.size() == 2)
144 {
145
146 // Map to lists of hits for both sensors
147 Map<SiSensor,List<SiTrackerHitStrip1D>> sensor_hitlists = new HashMap<SiSensor,List<SiTrackerHitStrip1D>>();
148
149 // Loop over sensors
150 for (SiSensor sensor : sensors)
151 {
152 // bail out if we find a double-sided sensor
153 if (sensor.isDoubleSided() || sensor.hasPixels())
154 {
155 return hits_2D;
156 }
157
158 // Get hits for this sensor
159 IReadout ro = sensor.getReadout();
160 List<SiTrackerHitStrip1D> hits = ro.getHits(SiTrackerHitStrip1D.class);
161
162 // Add hits to map
163 if (hits.size() != 0)
164 {
165 sensor_hitlists.put(sensor,new ArrayList<SiTrackerHitStrip1D>());
166 sensor_hitlists.get(sensor).addAll(hits);
167 }
168
169 }
170
171 // If we have hits on both sensors, make 2D hits
172 if (sensor_hitlists.keySet().size() == 2)
173 {
174 // System.out.println("Found module with hits on both sides. Module: "+module.getName());
175
176 List<List<SiTrackerHitStrip1D>> hitlists = new ArrayList<List<SiTrackerHitStrip1D>>(sensor_hitlists.values());
177
178 for (SiTrackerHitStrip1D hit1: hitlists.get(0))
179 {
180 for (SiTrackerHitStrip1D hit2: hitlists.get(1))
181 {
182 List<SiTrackerHitStrip1D> hitpair = new ArrayList<SiTrackerHitStrip1D>();
183 hitpair.add(hit1);
184 hitpair.add(hit2);
185 hits_2D.add(makeTrackerHit2D(hitpair));
186 }
187 }
188 }
189
190 }
191
192 return hits_2D;
193 }
194
195
196 private SiTrackerHitStrip2D makeTrackerHit2D(List<SiTrackerHitStrip1D> hits_1D)
197 {
198 double energy = 0;
199 double time = 0;
200 List<RawTrackerHit> raw_hits = new ArrayList<RawTrackerHit>();
201
202 for (SiTrackerHitStrip1D hit_1D : hits_1D)
203 {
204 energy += hit_1D.getdEdx();
205 time += hit_1D.getdEdx()*hit_1D.getTime();
206 raw_hits.addAll(hit_1D.getRawHits());
207 }
208 time /= energy;
209
210 TrackerHitType decoded_type = new TrackerHitType(TrackerHitType.CoordinateSystem.GLOBAL,
211 TrackerHitType.MeasurementType.STRIP_2D);
212
213 return new SiTrackerHitStrip2D(calculateGlobalCovariance(hits_1D), energy, time,
214 raw_hits, decoded_type, hits_1D);
215
216 }
217
218 SymmetricMatrix calculateGlobalCovariance(List<SiTrackerHitStrip1D> hits_1D)
219 {
220 double[] covariance_inverted_sum = new double[6];
221
222 for (SiTrackerHitStrip1D hit_1D : hits_1D)
223 {
224
225 //System.out.println("proc hit: " + ((SimTrackerHit)hit_1D.getSimHits().toArray()[0]).getDetectorElement().getName());
226
227 // SymmetricMatrix covariance = hit_1D.getTransformedHit(TrackerHitType.CoordinateSystem.GLOBAL).getCovarianceAsMatrix();
228 // SymmetricMatrix covariance_inverted = new SymmetricMatrix(3);
229 //
230 // for (int i = 0; i<3; i++)
231 // {
232 // System.out.println("covariance.diagonal(" +i+ "): "+covariance.diagonal(i));
233 //
234 // if (covariance.diagonal(i) == 0)
235 // {
236 // System.out.println("Found zero element for i: "+i);
237 // covariance.setElement(i,i,1);
238 // }
239 // }
240 //
241 // System.out.println("Covariance matrix: "+covariance);
242 //
243 // inverse(covariance,covariance_inverted);
244 //
245 // System.out.println("Covariance matrix inverted: "+covariance_inverted);
246
247
248 // Get covariance and invert to get matrices that add linearly in contrbuting to chisquared
249 // FIXME for now we cheat and attribute a negligible measurement error to the third coordinate
250 // to avoid an indeteminate (non-invertable) matrix.
251 // Ideally, we should treat only the 2X2 submatrix
252
253 SymmetricMatrix covariance_inverted = hit_1D.getTransformedHit(TrackerHitType.CoordinateSystem.GLOBAL).getCovarianceAsMatrix();
254 for (int i = 0; i<covariance_inverted.getNRows(); i++)
255 {
256 if (covariance_inverted.diagonal(i) < GeomOp3D.DISTANCE_TOLERANCE*GeomOp3D.DISTANCE_TOLERANCE)
257 {
258 // System.out.println("Found zero element for i: "+i);
259 covariance_inverted.setElement(i,i,GeomOp3D.DISTANCE_TOLERANCE*GeomOp3D.DISTANCE_TOLERANCE);
260 }
261 }
262
263 //System.out.println("Covariance matrix: "+covariance_inverted);
264
265 covariance_inverted.invert();
266
267 // System.out.println("Covariance matrix inverted: "+covariance_inverted);
268
269
270
271
272 // Sum the matrices
273 for (int i = 0; i < covariance_inverted_sum.length; i++)
274 {
275 covariance_inverted_sum[i] += covariance_inverted.asPackedArray(true)[i];
276 }
277
278 }
279
280 SymmetricMatrix covariance = new SymmetricMatrix(3,covariance_inverted_sum,true);
281 covariance.invert();
282
283 //System.out.println("Covariance matrix after: "+covariance);
284
285 return covariance;
286 }
287
288
289 // public static void inverse(Matrix mIn, MutableMatrix mOut) throws InvalidMatrixException
290 // {
291 // int order = mIn.getNRows();
292 // if (order != mIn.getNColumns()) throw new InvalidMatrixException("Matrix.inverse only supports square matrices");
293 // if (order != mOut.getNColumns() && order != mOut.getNRows()) throw new InvalidMatrixException("mOut must be same size as mIn");
294 //
295 // int[] ik = new int[order];
296 // int[] jk = new int[order];
297 // double[][] array = new double[order][order];
298 // for (int i=0;i<order;i++)
299 // {
300 // System.out.println(" i: "+i);
301 // for (int j=0; j<order; j++)
302 // {
303 // System.out.println(" j: "+j);
304 // array[i][j] = mIn.e(i,j);
305 // System.out.println("array(i,j): "+array[i][j]);
306 // }
307 // }
308 //
309 // for (int k=0; k<order; k++)
310 // {
311 // // Find largest element in rest of matrix
312 // double amax = 0;
313 // for (int i=k; i<order; i++)
314 // {
315 // for (int j=k; j<order; j++)
316 // {
317 //
318 // double test = array[i][j];
319 //
320 // for (int i2=0;i2<order;i2++)
321 // {
322 // System.out.println(" i2: "+i2);
323 // for (int j2=0; j2<order; j2++)
324 // {
325 // System.out.println(" j2: "+j2);
326 // System.out.println("array(i2,j2): "+array[i2][j2]);
327 // }
328 // }
329 //
330 //
331 // System.out.println("k: "+k+" i: "+i+" j:"+j);
332 // System.out.println("test "+test);
333 // System.out.println("array[i][j]: "+array[i][j]);
334 // System.out.println("Math.abs(array[i][j]): "+Math.abs(array[i][j]));
335 // System.out.println("Math.abs(amax): "+Math.abs(amax));
336 //
337 // if (Math.abs(array[i][j]) > Math.abs(amax))
338 // {
339 // amax = array[i][j];
340 // ik[k] = i;
341 // jk[k] = j;
342 // }
343 // }
344 // }
345 //
346 // // Interchange rows and columns to put max in array[k][k]
347 //
348 // System.out.println("FINAL amax: "+amax);
349 //
350 // if (amax == 0)
351 // {
352 // throw new IndeterminateMatrixException();
353 // }
354 //
355 // {
356 // int i = ik[k];
357 // assert(k <= i);
358 // if (i > k)
359 // {
360 // for (int j=0; j<order; j++)
361 // {
362 // double save = array[k][j];
363 // array[k][j] = array[i][j];
364 // array[i][j] = -save;
365 // }
366 // }
367 // }
368 // {
369 // int j = jk[k];
370 // assert(k <= j);
371 // if (j > k)
372 // {
373 // for (int i=0; i<order; i++)
374 // {
375 // double save = array[i][k];
376 // array[i][k] = array[i][j];
377 // array[i][j] = -save;
378 // }
379 // }
380 // }
381 //
382 // // Accumulate elements of inverse matrix
383 //
384 // for (int i=0; i<order; i++)
385 // {
386 // if (i == k) continue;
387 // array[i][k] = -array[i][k]/amax;
388 // }
389 // for (int i=0; i<order; i++)
390 // {
391 // if (i == k) continue;
392 // for (int j=0; j<order; j++)
393 // {
394 // if (j == k) continue;
395 // array[i][j] += array[i][k]*array[k][j];
396 // }
397 // }
398 // for (int j=0; j<order; j++)
399 // {
400 // if (j == k) continue;
401 // array[k][j] = array[k][j]/amax;
402 // }
403 // array[k][k] = 1/amax;
404 // }
405 //
406 // // restore ordering of matrix
407 //
408 // for (int l=0; l<order; l++)
409 // {
410 // int k = order - l - 1;
411 // {
412 // int j = ik[k];
413 // if (j>k)
414 // {
415 // for (int i=0; i<order; i++)
416 // {
417 // double save = array[i][k];
418 // array[i][k] = -array[i][j];
419 // array[i][j] = save;
420 // }
421 // }
422 // }
423 // {
424 // int i = jk[k];
425 // if (i>k)
426 // {
427 // for (int j=0; j<order; j++)
428 // {
429 // double save = array[k][j];
430 // array[k][j] = -array[i][j];
431 // array[i][j] = save;
432 // }
433 // }
434 // }
435 // }
436 // for (int i=0;i<order;i++)
437 // {
438 // for (int j=0; j<order; j++)
439 // {
440 // mOut.setElement(i,j,array[i][j]);
441 // }
442 // }
443 // }
444
445
446 }
447
448
449 // // Get hits for this sensor
450 // IReadout ro = sensor.getReadout();
451 // List<SiTrackerHitStrip1D> hits = ro.getHits(SiTrackerHitStrip1D.class);
452 //
453 // // Make map from electrodes to lists of hits
454 // Map<SiSensorElectrodes,List<SiTrackerHitStrip1D>> side_hitlists = new HashMap<SiSensorElectrodes,List<SiTrackerHitStrip1D>>();
455 // for (SiTrackerHitStrip1D hit : hits)
456 // {
457 // SiSensorElectrodes electrodes = hit.getReadoutElectrodes();
458 // if (side_hitlists.get(electrodes) == null)
459 // {
460 // side_hitlists.put(electrodes,new ArrayList<SiTrackerHitStrip1D>());
461 // }
462 //
463 // side_hitlists.get(electrodes).add(hit);
464 // }
465 //
466 // // If have hits on both sides, make 2d strip hits
467 // if (side_hitlists.keySet().size() == 2)
468 // {
469 // List<List<SiTrackerHitStrip1D>> hitlists = new ArrayList<List<SiTrackerHitStrip1D>>(side_hitlists.values());
470 //
471 // for (SiTrackerHitStrip1D hit1: hitlists.get(0))
472 // {
473 // for (SiTrackerHitStrip1D hit2: hitlists.get(1))
474 // {
475 // List<SiTrackerHitStrip1D> hitpair = new ArrayList<SiTrackerHitStrip1D>();
476 // hitpair.add(hit1);
477 // hitpair.add(hit2);
478 // hits_2D.add(makeTrackerHit2D(hitpair));
479 // }
480 // }
481 // }
482
483
484 // Map<SiSensor,List<SiTrackerHitStrip1D>> sensor_hitlists = new HashMap<SiSensor,List<SiTrackerHitStrip1D>>();
485 //
486 // for (SiSensor sensor : sensors)
487 // {
488 // // Get hits for this sensor
489 // IReadout ro = sensor.getReadout();
490 // List<SiTrackerHitStrip1D> hits = ro.getHits(SiTrackerHitStrip1D.class);
491 //
492 // // If sensor is single-sided strip sensor, add hits to map
493 // if (hits.size() != 0)
494 // {
495 // sensor_hitlists.put(sensor,new ArrayList<SiTrackerHitStrip1D>());
496 // sensor_hitlists.get(sensor).addAll(hits);
497 // }
498 //
499 // }
500 //
501 // if (sensor_hitlists.keySet().size() == 2)
502 // {
503 // List<List<SiTrackerHitStrip1D>> hitlists = new ArrayList<List<SiTrackerHitStrip1D>>(sensor_hitlists.values());
504 //
505 // for (SiTrackerHitStrip1D hit1: hitlists.get(0))
506 // {
507 // for (SiTrackerHitStrip1D hit2: hitlists.get(1))
508 // {
509 // List<SiTrackerHitStrip1D> hitpair = new ArrayList<SiTrackerHitStrip1D>();
510 // hitpair.add(hit1);
511 // hitpair.add(hit2);
512 // hits_2D.add(makeTrackerHit2D(hitpair));
513 // }
514 // }
515 // }