1
2
3
4
5
6
7
8
9
10 package org.lcsim.recon.tracking.digitization.sisim;
11
12 import org.lcsim.detector.solids.GeomOp3D;
13 import org.lcsim.event.SimTrackerHit;
14 import org.lcsim.detector.tracker.silicon.SiSensor;
15 import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
16 import org.lcsim.detector.tracker.silicon.ChargeCarrier;
17 import org.lcsim.detector.ITransform3D;
18 import org.lcsim.detector.IRotation3D;
19 import org.lcsim.detector.tracker.silicon.ChargeDistribution;
20 import org.lcsim.detector.tracker.silicon.GaussianDistribution2D;
21
22 import hep.physics.vec.Hep3Vector;
23 import hep.physics.vec.BasicHep3Vector;
24 import hep.physics.matrix.BasicMatrix;
25 import hep.physics.vec.VecOp;
26
27 import java.util.ArrayList;
28 import java.util.List;
29 import java.util.Map;
30 import java.util.EnumMap;
31 import java.util.Map.Entry;
32 import java.util.SortedMap;
33
34 import org.lcsim.detector.Transform3D;
35 import org.lcsim.detector.Translation3D;
36 import org.lcsim.detector.solids.Line3D;
37 import org.lcsim.detector.solids.Point3D;
38
39
40
41
42
43 public class CDFSiSensorSim implements SiSensorSim
44 {
45
46
47 SiSensor _sensor = null;
48
49 Map<ChargeCarrier,Hep3Vector> _drift_direction = null;
50 Map<ChargeCarrier,SiElectrodeDataCollection> _sense_data = null;
51 Map<ChargeCarrier,SiElectrodeDataCollection> _readout_data = null;
52
53
54
55
56
57
58 double _trapping = 0.0;
59
60
61
62
63 private static double _DEPOSITION_GRANULARITY = 0.10;
64 private static final double DISTANCE_ERROR_THRESHOLD = 0.001;
65
66 private final boolean debug = false;
67
68
69
70
71 public CDFSiSensorSim()
72 {
73 _drift_direction = new EnumMap<ChargeCarrier,Hep3Vector>(ChargeCarrier.class);
74 _sense_data = new EnumMap<ChargeCarrier,SiElectrodeDataCollection>(ChargeCarrier.class);
75 _readout_data = new EnumMap<ChargeCarrier,SiElectrodeDataCollection>(ChargeCarrier.class);
76 for (ChargeCarrier carrier : ChargeCarrier.values())
77 {
78 _sense_data.put(carrier,new SiElectrodeDataCollection());
79 _readout_data.put(carrier,new SiElectrodeDataCollection());
80 }
81 }
82
83
84 public void setTrapping(double trapping)
85 {
86 _trapping = trapping;
87 }
88
89
90
91
92
93 public SiElectrodeDataCollection getReadoutData(ChargeCarrier carrier)
94 {
95 return _readout_data.get(carrier);
96 }
97
98
99
100
101
102 public void setSensor(SiSensor sensor)
103 {
104 _sensor = sensor;
105 }
106
107 public Map<ChargeCarrier,SiElectrodeDataCollection> computeElectrodeData()
108 {
109 if(debug) {
110 System.out.printf("%s: computeElectrodeData for sensor %s\n", this.getClass().getSimpleName(),this._sensor.getName());
111 System.out.printf("%s: # Sense strips: %d\n", this.getClass().getSimpleName(),_sensor.getSenseElectrodes(ChargeCarrier.HOLE).getNCells(0));
112 System.out.printf("%s: # Readout strips: %d\n", this.getClass().getSimpleName(),_sensor.getReadoutElectrodes(ChargeCarrier.HOLE).getNCells(0));
113 }
114 depositChargeOnSense();
115 transferChargeToReadout();
116
117 return _readout_data;
118 }
119
120
121 public void clearReadout()
122 {
123 for (ChargeCarrier carrier : ChargeCarrier.values())
124 {
125 _readout_data.get(carrier).clear();
126 }
127 }
128
129 public void lorentzCorrect(Hep3Vector position, ChargeCarrier carrier)
130 {
131 Line3D drift_line = new Line3D(new Point3D(position), driftDirection(carrier,new BasicHep3Vector(0,0,0)));
132
133 List<Point3D> intersection_points = new ArrayList<Point3D>();
134
135 for (ChargeCarrier bias_carrier : ChargeCarrier.values())
136 {
137 intersection_points.add(GeomOp3D.intersection(drift_line,_sensor.getBiasSurface(bias_carrier).getPlane()));
138 }
139
140 Hep3Vector corrected_position = VecOp.mult(0.5,VecOp.add(intersection_points.get(0),intersection_points.get(1)));
141
142 ITransform3D transform = new Transform3D(new Translation3D(VecOp.sub(corrected_position,position)));
143
144 transform.transform(position);
145 }
146
147
148 private void clearSense()
149 {
150 for (ChargeCarrier carrier : ChargeCarrier.values())
151 {
152 _sense_data.get(carrier).clear();
153 }
154 }
155
156 private void depositChargeOnSense()
157 {
158
159 if(debug)
160 System.out.printf("%s: depositChargeOnSense for sensor %s\n", this.getClass().getSimpleName(),this._sensor.getName());
161
162
163 for (ChargeCarrier carrier : ChargeCarrier.values())
164 {
165 if (_sensor.hasElectrodesOnSide(carrier))
166 {
167 _drift_direction.put( carrier, driftDirection(carrier,new BasicHep3Vector(0.0,0.0,0.0)) );
168 }
169 }
170
171 ITransform3D global_to_sensor = _sensor.getGeometry().getGlobalToLocal();
172 List<SimTrackerHit> hits = _sensor.getReadout().getHits(SimTrackerHit.class);
173
174 for (SimTrackerHit hit : hits)
175 {
176 if(debug)
177 System.out.printf("%s: depositChargeOnSense for sim tracker hit at %s \n", this.getClass().getSimpleName(),hit.getPositionVec().toString());
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193 TrackSegment track = new TrackSegment(hit);
194 track.transform(global_to_sensor);
195
196 int nsegments = 0;
197
198
199 for (ChargeCarrier carrier : ChargeCarrier.values())
200 {
201 if (_sensor.hasElectrodesOnSide(carrier))
202 {
203
204 nsegments = Math.max(nsegments,nSegments(track,carrier, _DEPOSITION_GRANULARITY));
205 }
206 }
207
208
209
210
211
212
213
214
215
216
217
218
219
220 double segment_length = track.getLength()/nsegments;
221 double segment_charge = track.getEloss()/nsegments/_sensor.getBulk().ENERGY_EHPAIR;
222
223
224
225
226 Hep3Vector segment_step = VecOp.mult(segment_length,track.getDirection());
227 Hep3Vector segment_center = VecOp.add( track.getP1(),VecOp.mult(0.5,segment_step) );
228
229
230
231
232 for (int iseg = 0; iseg < nsegments; iseg++)
233 {
234
235
236
237
238
239 for (ChargeCarrier carrier : ChargeCarrier.values())
240 {
241 if (_sensor.hasElectrodesOnSide(carrier))
242 {
243 SiSensorElectrodes electrodes = _sensor.getSenseElectrodes(carrier);
244
245
246 double collection_efficiency = 1.0 - 10*_trapping*
247 driftVector(segment_center,carrier).magnitude();
248 collection_efficiency = Math.max(0.0,Math.min(1.0,collection_efficiency));
249 segment_charge *= (collection_efficiency);
250
251 ChargeDistribution charge_distribution = diffusionDistribution(segment_charge,segment_center,carrier);
252 charge_distribution.transform(electrodes.getParentToLocal());
253
254 SortedMap<Integer,Integer> sense_charge = electrodes.computeElectrodeData(charge_distribution);
255
256
257 _sense_data.get(carrier).add(new SiElectrodeDataCollection(sense_charge,hit));
258 }
259 }
260
261
262 segment_center = VecOp.add(segment_center, segment_step);
263 }
264
265 }
266
267 if(debug) {
268 System.out.printf("%s: Final sense charge map:\n", this.getClass().getSimpleName());
269 for(Map.Entry<Integer, Integer> entry : _sense_data.get(ChargeCarrier.HOLE).getChargeMap().entrySet()) {
270 System.out.printf("%s: cell %d -> %d \n", this.getClass().getSimpleName(),entry.getKey(),entry.getValue());
271 }
272 }
273
274 }
275
276 private void transferChargeToReadout()
277 {
278 for (ChargeCarrier carrier : ChargeCarrier.values())
279 {
280 if (_sensor.hasElectrodesOnSide(carrier))
281 {
282 if (_sensor.isACCoupled(carrier))
283 {
284
285
286 SiSensorElectrodes sense_electrodes = _sensor.getSenseElectrodes(carrier);
287 SiSensorElectrodes readout_electrodes = _sensor.getReadoutElectrodes(carrier);
288 BasicMatrix transfer_efficiencies = _sensor.getTransferEfficiencies(carrier);
289
290 SiElectrodeDataCollection sense_data = _sense_data.get(carrier);
291 SiElectrodeDataCollection readout_data = _readout_data.get(carrier);
292
293 for (Integer sense_cell : sense_data.keySet())
294 {
295
296
297
298 SiElectrodeData sense_cell_data = sense_data.get(sense_cell);
299
300
301
302 int sense_row = sense_electrodes.getRowNumber(sense_cell);
303 int sense_col = sense_electrodes.getColumnNumber(sense_cell);
304 int row_steps = transfer_efficiencies.getNRows()-1;
305 int col_steps = transfer_efficiencies.getNColumns()-1;
306
307
308
309
310
311
312
313
314 for (int irow = sense_row - row_steps; irow <= sense_row + row_steps; irow++)
315 {
316
317
318
319 for (int icol = sense_col - col_steps; icol <= sense_col + col_steps; icol++)
320 {
321
322
323
324 int sense_id = sense_electrodes.getCellID(irow,icol);
325 if (sense_id < 0) continue;
326 Hep3Vector sense_position = sense_electrodes.getCellPosition(sense_id);
327 int readout_cell = readout_electrodes.getCellID(sense_position);
328
329
330
331
332
333
334
335
336
337
338 if ( readout_electrodes.isValidCell(readout_cell) &&
339 readout_electrodes.getPositionInCell(sense_position).x() == 0.0 &&
340 readout_electrodes.getPositionInCell(sense_position).y() == 0.0 )
341 {
342
343 double transfer_efficiency = transfer_efficiencies.e(Math.abs(irow-sense_row),Math.abs(icol-sense_col));
344
345
346
347 SiElectrodeData readout_datum = new SiElectrodeData((int)Math.round(transfer_efficiency*sense_cell_data.getCharge()),
348 sense_cell_data.getSimulatedHits());
349 readout_data.add(readout_cell,readout_datum);
350
351
352 }
353
354 }
355
356 }
357
358 }
359
360 }
361 else
362 {
363 _readout_data.put(carrier,_sense_data.get(carrier));
364 }
365
366
367 }
368 }
369
370 clearSense();
371
372 }
373
374 private int nSegments(TrackSegment track, ChargeCarrier carrier, double deposition_granularity)
375 {
376
377 int nsegments = 0;
378 if (!_sensor.hasElectrodesOnSide(carrier)) return nsegments;
379 SiSensorElectrodes electrodes = _sensor.getSenseElectrodes(carrier);
380
381
382
383
384
385
386 nsegments = (int)Math.ceil(Math.abs(VecOp.dot(track.getVector(),electrodes.getGeometry().getNormal()))/(_sensor.getThickness()*deposition_granularity));
387
388
389
390 Hep3Vector deposition_line = VecOp.sub( driftDestination(track.getP2(),carrier),
391 driftDestination(track.getP1(),carrier) );
392
393 int naxes = electrodes.getNAxes();
394 for (int iaxis = 0; iaxis < naxes; iaxis++)
395 {
396 IRotation3D electrode_to_sensor = electrodes.getParentToLocal().getRotation().inverse();
397 Hep3Vector measured_coordinate = electrode_to_sensor.rotated(electrodes.getMeasuredCoordinate(iaxis));
398
399 double projected_deposition_length = Math.abs(VecOp.dot(deposition_line,measured_coordinate));
400
401
402
403
404
405 int required_segments = (int)Math.ceil(projected_deposition_length/
406 (deposition_granularity*electrodes.getPitch(iaxis)));
407 nsegments = Math.max(nsegments,required_segments);
408 }
409
410
411
412
413
414
415
416
417
418
419
420 return nsegments;
421 }
422
423
424 private Hep3Vector driftDestination(Hep3Vector origin, ChargeCarrier carrier)
425 {
426
427 return VecOp.add(origin,driftVector(origin, carrier));
428 }
429
430 private Hep3Vector driftVector(Hep3Vector origin, ChargeCarrier carrier)
431 {
432
433
434 double drift_vector_scale = _sensor.distanceFromSide(origin,carrier)/VecOp.dot(_drift_direction.get(carrier),_sensor.getBiasSurface(carrier).getNormal());
435
436
437 return VecOp.mult(drift_vector_scale,_drift_direction.get(carrier));
438 }
439
440 private Hep3Vector driftDirection(ChargeCarrier carrier, Hep3Vector local_position)
441 {
442
443
444
445
446
447
448
449 Hep3Vector b_field = _sensor.getBField(local_position);
450
451
452
453 Hep3Vector e_field = VecOp.mult(1000., _sensor.electricField(local_position));
454
455
456
457 double mobility = 1.0e-4 * _sensor.getBulk().mobility(carrier);
458
459
460 double qmu = carrier.charge() * mobility;
461
462
463
464 Hep3Vector vpar = VecOp.mult(qmu, e_field);
465
466
467
468 Hep3Vector vperp = VecOp.mult(qmu*qmu, VecOp.cross(e_field, b_field));
469
470
471
472 Hep3Vector drift_direction = VecOp.unit(VecOp.add(vpar, vperp));
473
474
475
476 return drift_direction;
477
478 }
479
480 private ChargeDistribution diffusionDistribution(double segment_charge, Hep3Vector origin, ChargeCarrier carrier)
481 {
482
483
484
485
486 double distance = _sensor.distanceFromSide(origin,carrier);
487 double thickness = _sensor.getThickness();
488
489
490
491
492 if (distance < -DISTANCE_ERROR_THRESHOLD || distance > thickness + DISTANCE_ERROR_THRESHOLD){
493 throw new RuntimeException("Distance is outside of sensor by more than "+DISTANCE_ERROR_THRESHOLD+". Distance = "+distance+
494 ". If this is an isolated event, then perhaps DISTANCE_ERROR_THRESHOLD must be increased in CDFSiSensorSim");
495 }
496 else if (distance < 0) distance = 0.;
497 else if (distance > thickness) distance = thickness;
498
499
500
501 double bias_voltage = _sensor.getBiasVoltage();
502 double depletion_voltage = _sensor.getDepletionVoltage();
503
504
505 double difference_V = bias_voltage - depletion_voltage;
506 double sum_V = bias_voltage + depletion_voltage;
507 double common_factor = 2.0 * distance * depletion_voltage / thickness;
508
509
510
511
512
513 double sigmasq = _sensor.getBulk().K_BOLTZMANN * _sensor.getBulk().getTemperature() *
514 _sensor.getThickness()*_sensor.getThickness() / _sensor.getDepletionVoltage();
515 if (_sensor.getBulk().isNtype() == (carrier==ChargeCarrier.HOLE))
516 {
517 sigmasq *= Math.log( sum_V / (sum_V - common_factor));
518 }
519 else
520 {
521 sigmasq *= Math.log( (difference_V + common_factor) / difference_V );
522 }
523
524 double sigma = Math.sqrt(sigmasq);
525
526
527
528
529
530
531 Hep3Vector drift_direction = _drift_direction.get(carrier);
532 Hep3Vector bias_surface_normal = _sensor.getBiasSurface(carrier).getNormal();
533
534 Hep3Vector major_axis;
535 Hep3Vector minor_axis;
536
537 double major_axis_length;
538 double minor_axis_length;
539
540 if (VecOp.cross(drift_direction,bias_surface_normal).magnitude() < GeomOp3D.ANGULAR_TOLERANCE)
541 {
542
543
544
545
546
547
548 major_axis = VecOp.cross(bias_surface_normal,_sensor.getBiasSurface(carrier).getEdges().get(0).getDirection());
549 minor_axis = VecOp.cross(bias_surface_normal,major_axis);
550
551
552
553
554 major_axis_length = minor_axis_length = sigma;
555
556 }
557 else
558 {
559
560
561 major_axis = VecOp.unit(VecOp.cross(bias_surface_normal,VecOp.cross(_drift_direction.get(carrier),bias_surface_normal)));
562 minor_axis = VecOp.cross(bias_surface_normal,major_axis);
563
564
565
566
567 double cos_theta_lorentz = VecOp.dot(_drift_direction.get(carrier),_sensor.getBiasSurface(carrier).getNormal());
568
569
570
571 minor_axis_length = sigma*(1.0/cos_theta_lorentz);
572 major_axis_length = minor_axis_length*(1.0/cos_theta_lorentz);
573 }
574
575
576 major_axis = VecOp.mult(major_axis_length,major_axis);
577 minor_axis = VecOp.mult(minor_axis_length,minor_axis);
578
579
580
581
582
583 Hep3Vector drift_destination = driftDestination(origin,carrier);
584
585
586
587 ChargeDistribution distribution = new GaussianDistribution2D(segment_charge, drift_destination,major_axis,minor_axis);
588
589
590
591
592
593
594
595
596 return distribution;
597
598 }
599
600 }
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618