1
2
3
4 package org.lcsim.recon.tracking.digitization.sisim;
5
6 import java.util.ArrayList;
7 import java.util.HashSet;
8 import java.util.List;
9 import java.util.Map.Entry;
10 import java.util.Random;
11 import java.util.Set;
12 import java.util.SortedMap;
13 import java.util.TreeMap;
14 import org.apache.commons.math.MathException;
15 import org.apache.commons.math.distribution.BinomialDistribution;
16 import org.apache.commons.math.distribution.BinomialDistributionImpl;
17 import org.apache.commons.math.distribution.NormalDistribution;
18 import org.apache.commons.math.distribution.NormalDistributionImpl;
19 import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
20 import org.lcsim.event.RawTrackerHit;
21 import org.lcsim.math.probability.Erf;
22 import org.lcsim.recon.tracking.digitization.sisim.ReadoutChip.ReadoutChannel;
23
24
25
26
27
28
29
30
31
32
33
34 public class GenericReadoutChip implements ReadoutChip {
35
36 private static Random _random = new Random();
37 private static NormalDistribution _gaussian = new NormalDistributionImpl(0.0, 1.0);
38 private static BinomialDistribution _binomial = new BinomialDistributionImpl(1, 1);
39 private double _noise_threshold;
40 private double _neighbor_threshold;
41 private GenericChannel _channel = new GenericChannel();
42
43
44 public GenericReadoutChip() {
45 }
46
47
48
49
50
51
52
53 public void setNoiseIntercept(double noise_intercept) {
54 _channel.setNoiseIntercept(noise_intercept);
55 }
56
57
58
59
60
61
62
63 public void setNoiseSlope(double noise_slope) {
64 _channel.setNoiseSlope(noise_slope);
65 }
66
67
68
69
70
71
72 public void setNoiseThreshold(double noise_threshold) {
73 _noise_threshold = noise_threshold;
74 }
75
76
77
78
79
80
81
82 public void setNeighborThreshold(double neighbor_threshold) {
83 _neighbor_threshold = neighbor_threshold;
84 }
85
86
87
88
89
90
91
92 public void setConversionConstant(double adc_per_fC) {
93 _channel.setConversionConstant(adc_per_fC);
94 }
95
96
97
98
99
100
101
102
103
104 public GenericChannel getChannel(int channel_number) {
105 return _channel;
106 }
107
108
109
110
111
112
113
114
115
116 public SortedMap<Integer, List<Integer>> readout(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
117
118
119
120 if (data == null) {
121 data = new SiElectrodeDataCollection();
122 }
123
124
125 addNoise(data, electrodes);
126
127
128
129 return digitize(data, electrodes);
130 }
131
132
133
134
135
136
137
138 public double decodeCharge(RawTrackerHit hit) {
139
140
141 int adc = hit.getADCValues()[0];
142
143
144 return adc / (1.6e-4 * _channel.getConversionConstant());
145 }
146
147
148
149
150
151
152
153
154 public int decodeTime(RawTrackerHit hit) {
155 return 0;
156 }
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171 private SortedMap<Integer, List<Integer>> digitize(SiElectrodeDataCollection data,
172 SiSensorElectrodes electrodes) {
173
174
175 SortedMap<Integer, List<Integer>> chip_data = new TreeMap<Integer, List<Integer>>();
176
177
178 for (Integer channel : data.keySet()) {
179
180
181 SiElectrodeData eldata = data.get(channel);
182
183
184 double charge = eldata.getCharge();
185
186
187 if (charge < _neighbor_threshold) continue;
188
189
190 if (charge < _noise_threshold) {
191
192
193 boolean nbrhit = false;
194 for (Integer nbr : electrodes.getNearestNeighborCells(channel)) {
195
196
197 SiElectrodeData nbrdata = data.get(nbr);
198 if (nbrdata == null) continue;
199
200
201 if (nbrdata.getCharge() >= _noise_threshold) {
202 nbrhit = true;
203 break;
204 }
205 }
206
207
208 if (!nbrhit) continue;
209 }
210
211
212 int adc = getChannel(channel).computeAdcValue(charge);
213 if (adc <= 0) {
214 continue;
215 }
216
217
218
219 List<Integer> channel_data = new ArrayList<Integer>();
220 channel_data.add(adc);
221
222
223 chip_data.put(channel, channel_data);
224 }
225
226 return chip_data;
227 }
228
229
230
231
232
233
234
235 private void addNoise(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
236
237
238
239 for (Entry datum : data.entrySet()) {
240
241
242 int channel = (Integer) datum.getKey();
243 SiElectrodeData eldata = (SiElectrodeData) datum.getValue();
244
245
246 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
247
248
249 int noise_charge = (int) Math.round(_random.nextGaussian() * noise);
250 eldata.addCharge(noise_charge);
251 }
252
253
254
255
256 int nelectrodes = electrodes.getNCells();
257 int nelectrodes_empty = nelectrodes - data.size();
258
259
260 double noiseRMS = _channel.computeNoise(electrodes.getCapacitance());
261 double normalized_integration_limit = _noise_threshold / noiseRMS;
262
263
264 double integral = Erf.phic(normalized_integration_limit);
265 int nchannels_throw = drawBinomial(nelectrodes_empty, integral);
266
267
268 for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
269
270 int channel = _random.nextInt(nelectrodes);
271 while (data.keySet().contains(channel)) {
272 channel = _random.nextInt(nelectrodes);
273 }
274
275
276 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
277
278
279 int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
280 data.add(channel, new SiElectrodeData(charge));
281 }
282
283
284
285 nchannels_throw = 1;
286 while (nchannels_throw > 0) {
287
288
289
290 Set<Integer> neighbors = new HashSet<Integer>();
291 for (int channel : data.keySet()) {
292 neighbors.addAll(electrodes.getNearestNeighborCells(channel));
293 }
294 neighbors.removeAll(data.keySet());
295
296 nelectrodes_empty = neighbors.size();
297
298
299 normalized_integration_limit = _neighbor_threshold / noiseRMS;
300
301 integral = Erf.phic(normalized_integration_limit);
302 nchannels_throw = drawBinomial(nelectrodes_empty, integral);
303
304
305 for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
306
307 List<Integer> neighbors_list = new ArrayList<Integer>(neighbors);
308
309 int channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
310
311 while (data.keySet().contains(channel)) {
312 channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
313 }
314
315 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
316
317
318
319 int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
320 data.add(channel, new SiElectrodeData(charge));
321 }
322
323 }
324
325 }
326
327 public static int drawBinomial(int ntrials, double probability) {
328 _binomial.setNumberOfTrials(ntrials);
329 _binomial.setProbabilityOfSuccess(probability);
330
331 int nsuccess = 0;
332 try {
333 nsuccess = _binomial.inverseCumulativeProbability(_random.nextDouble());
334 } catch (MathException exception) {
335 throw new RuntimeException("Kpix failed to calculate inverse cumulative probability of binomial!");
336 }
337 return nsuccess;
338 }
339
340
341
342
343
344 public static double drawGaussianAboveThreshold(double prob_above_threshold) {
345 double draw, cumulative_probability;
346
347 draw = prob_above_threshold * _random.nextDouble();
348 cumulative_probability = 1.0 - prob_above_threshold + draw;
349
350 assert cumulative_probability < 1.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
351 assert cumulative_probability >= 0.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
352
353 double gaussian_random = 0;
354 try {
355 gaussian_random = _gaussian.inverseCumulativeProbability(cumulative_probability);
356 } catch (MathException e) {
357 System.out.println("MathException caught: " + e);
358 }
359
360 return gaussian_random;
361 }
362
363
364
365
366 private class GenericChannel implements ReadoutChannel {
367
368 private double _noise_intercept = 0.;
369 private double _noise_slope = 0.;
370 private double _noise_threshold = 1.;
371 private double _neighbor_threshold = 1.;
372 private double _adc_per_fC = 100.;
373
374
375
376
377
378
379 private void setConversionConstant(double adc_per_fC) {
380 _adc_per_fC = adc_per_fC;
381 }
382
383
384
385
386
387
388 private double getConversionConstant() {
389 return _adc_per_fC;
390 }
391
392
393
394
395
396
397 private void setNoiseIntercept(double noise_intercept) {
398 _noise_intercept = noise_intercept;
399 }
400
401
402
403
404
405
406 private void setNoiseSlope(double noise_slope) {
407 _noise_slope = noise_slope;
408 }
409
410
411
412
413
414
415
416 public double computeNoise(double capacitance) {
417 return _noise_intercept + _noise_slope * capacitance;
418 }
419
420
421
422
423
424
425
426 private int computeAdcValue(double charge) {
427
428
429 int adc = (int) Math.floor(charge * 1.6e-4 * _adc_per_fC);
430
431
432 if (adc < 0) {
433 adc = 0;
434 }
435
436
437 if (adc > 32767) {
438 adc = 32767;
439 }
440
441 return adc;
442 }
443 }
444 }