1 /*
2 * Class BasicReadoutChip
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 * Basic readout chip class. This class supports the minimal functions expected of
26 * a readout chip. The charge on a strip/pixel is digitized as an integer number
27 * with a simple ADC with programmable resolution and dynamic range. A chip with
28 * 1-bit ADC resolution (binary readout) is treated as a special case.
29 *
30 * Noise is added to strips with charge and random noise hits are generated as well.
31 * Methods are provided to decode the charge and time (although the current
32 * implementation always returns a time of 0).
33 *
34 * This implementation has thresholds that are settable in units of RMS noise of
35 * each channel to enable simluation of highly optimized readout chains. If
36 * absolute thresholds are desired, GenericReadoutChip should be used instead.
37 *
38 * @author Tim Nelson
39 */
40 public class BasicReadoutChip implements ReadoutChip
41 {
42
43 private static Random _random = new Random();
44 private static NormalDistribution _gaussian = new NormalDistributionImpl(0.0, 1.0);
45 private static BinomialDistribution _binomial = new BinomialDistributionImpl(1, 1);
46 private double _noise_threshold = 4;
47 private double _neighbor_threshold = 4;
48 private BasicChannel _channel = new BasicChannel();
49 private ADC _adc = new ADC();
50
51 /** Creates a new instance of BasicReadoutChip */
52 public BasicReadoutChip()
53 {
54 }
55
56 /**
57 * Set the noise intercept (i.e., the noise for 0 strip/pixel capacitance).
58 * Units are electrons of noise.
59 *
60 * @param noise_intercept noise for 0 capacitance
61 */
62 public void setNoiseIntercept(double noise_intercept)
63 {
64 _channel.setNoiseIntercept(noise_intercept);
65 }
66
67 /**
68 * Set the noise slope (i.e., the proportionality between noise and capacitance).
69 * Units are electrons of noise per fF of capacitance.
70 *
71 * @param noise_slope noise slope per unit capacitance
72 */
73 public void setNoiseSlope(double noise_slope)
74 {
75 _channel.setNoiseSlope(noise_slope);
76 }
77
78 /**
79 * Set the threshold for reading out a channel. Units are multiples of RMS noise.
80 *
81 * @param noise_threshold
82 */
83 public void setNoiseThreshold(double noise_threshold)
84 {
85 _noise_threshold = noise_threshold;
86 }
87
88 /**
89 * Set the threshold for reading a channel if its neighbor is
90 * above the noise threshold. Units are multiples of RMS noise.
91 *
92 * @param neighbor_threshold
93 */
94 public void setNeighborThreshold(double neighbor_threshold)
95 {
96 _neighbor_threshold = neighbor_threshold;
97 }
98
99 /**
100 * Set the number of bits of ADC resolution
101 *
102 * @param nbits
103 */
104 public void setNbits(int nbits)
105 {
106 getADC().setNbits(nbits);
107 }
108
109 /**
110 * Set the dynamic range of the ADC
111 *
112 * @param dynamic_range in fC
113 */
114 public void setDynamicRange(double dynamic_range)
115 {
116 getADC().setDynamicRange(dynamic_range);
117 }
118
119 /**
120 * Return the BasicChannel associated with a given channel number.
121 * For the basic readout, there is a single instance of BasicChannel
122 * and thus the channel number is ignored.
123 *
124 * @param channel_number channel number
125 * @return associated BasicReadoutChannel
126 */
127 public BasicChannel getChannel(int channel_number)
128 {
129 return _channel;
130 }
131
132 private ADC getADC()
133 {
134 return _adc;
135 }
136
137 /**
138 * Given a collection of electrode data (i.e., charge on strips/pixels),
139 * return a map associating the channel and it's list of raw data.
140 *
141 * @param data electrode data from the charge distribution
142 * @param electrodes strip or pixel electrodes
143 * @return map containing the ADC counts for this sensor
144 */
145 public SortedMap<Integer, List<Integer>> readout(SiElectrodeDataCollection data, SiSensorElectrodes electrodes)
146 {
147
148 // If there is no electrode data for this readout chip, create an empty
149 // electrode data collection
150 if (data == null)
151 {
152 data = new SiElectrodeDataCollection();
153 }
154
155 // Add noise hits to the electrode data collection
156 addNoise(data, electrodes);
157
158 // return the digitized charge data as a map that associates a hit
159 // channel with a list of raw data for the channel
160 return digitize(data, electrodes);
161 }
162
163 /**
164 * Decode the hit charge stored in the RawTrackerHit
165 *
166 * @param hit raw hit
167 * @return hit charge in units of electrons
168 */
169 public double decodeCharge(RawTrackerHit hit)
170 {
171 return getADC().decodeCharge(hit.getADCValues()[0]);
172 }
173
174 /**
175 * Decode the hit time. Currently, the basic readout chip ignores the
176 * hit time and returns 0.
177 *
178 * @param hit raw hit data
179 * @return hit time
180 */
181 public int decodeTime(RawTrackerHit hit)
182 {
183 return 0;
184 }
185
186 /**
187 * Digitizes the hit channels in a SiElectrodeDataCollection.
188 *
189 * The SiElectrodeDataCollection is a map that associates a given channel with
190 * it's SiElectrodeData. The SiElectrodeData encapsulates the deposited charge
191 * on an strip/pixel and any associated SimTrackerHits.
192 *
193 * The output of this class is a map that associates a channel number with
194 * a list of raw data
195 *
196 * @param data electrode data collection
197 * @return map associating channels with a list of raw data
198 */
199 private SortedMap<Integer, List<Integer>> digitize(SiElectrodeDataCollection data,
200 SiSensorElectrodes electrodes)
201 {
202
203 // Create the map that associates a given sensor channel with it's list of raw data
204 SortedMap<Integer, List<Integer>> chip_data = new TreeMap<Integer, List<Integer>>();
205
206 // Loop over the channels contained in the SiElectrodeDataCollection
207 for (Integer channel : data.keySet())
208 {
209
210 // Fetch the electrode data for this channel
211 SiElectrodeData eldata = data.get(channel);
212
213 // Get the charge in units of electrons
214 double charge = eldata.getCharge();
215
216 // Get the noise RMS in units of electrons
217 double noiseRMS = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
218
219 // If the charge is below the neighbor threshold, don't digitize it
220 if (charge < _neighbor_threshold*noiseRMS)
221 {
222 continue;
223 }
224
225 // If charge is between neighbor and noise thresholds, check it's neighbors
226 if (charge < _noise_threshold*noiseRMS)
227 {
228
229 // Loop over neighbors and look for a neighbor with charge above the noise
230 boolean nbrhit = false;
231 for (Integer nbr : electrodes.getNearestNeighborCells(channel))
232 {
233
234 // See if we have electrode data for this neighbor
235 SiElectrodeData nbrdata = data.get(nbr);
236 if (nbrdata == null)
237 {
238 continue;
239 }
240
241 // See if we have found a neigbor above the noise threshold
242 if (nbrdata.getCharge() >= _noise_threshold)
243 {
244 nbrhit = true;
245 break;
246 }
247 }
248
249 // If there were no neighbor channels above threshold, don't digitize it
250 if (!nbrhit)
251 {
252 continue;
253 }
254 }
255
256 // Calculate the ADC value for this channel and make sure it is positive
257 int adc = getADC().convert(charge);
258 if (adc <= 0)
259 {
260 continue;
261 }
262
263 // Create a list containing the adc value - for the basic readout
264 // there is only 1 word of raw data
265 List<Integer> channel_data = new ArrayList<Integer>();
266 channel_data.add(adc);
267
268 // Save the list of raw data in the chip_data map
269 chip_data.put(channel, channel_data);
270 }
271
272 return chip_data;
273 }
274
275 /**
276 * Add noise hits for this readout chip
277 *
278 * @param data electrode data collection
279 * @param electrodes strip or pixel electrodes
280 */
281 private void addNoise(SiElectrodeDataCollection data, SiSensorElectrodes electrodes)
282 {
283
284 // First add noise to the strips/pixels in the SiElectrodeDataCollection
285 // Loop over the entries in the SiElectrodeDataCollection (which extends TreeMap)
286 for (Entry datum : data.entrySet())
287 {
288
289 // Get the channel number and electrode data for this entry
290 int channel = (Integer) datum.getKey();
291 SiElectrodeData eldata = (SiElectrodeData) datum.getValue();
292
293 // Get the RMS noise for this channel in units of electrons
294 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
295
296 // Add readout noise to the deposited charge
297 int noise_charge = (int) Math.round(_random.nextGaussian() * noise);
298 eldata.addCharge(noise_charge);
299 }
300
301 // Add random noise hits where the noise charge exceeds the noise threshold
302
303 // Find the number of pixels/strips that are not currently hit
304 int nelectrodes = electrodes.getNCells();
305 int nelectrodes_empty = nelectrodes - data.size();
306
307 // Get the noise threshold in units of the noise charge
308 double normalized_integration_limit = _noise_threshold;
309
310 // Calculate how many channels should get noise hits
311 double integral = Erf.phic(normalized_integration_limit);
312 int nchannels_throw = drawBinomial(nelectrodes_empty, integral);
313
314 // Now throw Gaussian randoms above the seed threshold and put signals on unoccupied channels
315 for (int ithrow = 0; ithrow < nchannels_throw; ithrow++)
316 {
317 // Throw to get a channel number
318 int channel = _random.nextInt(nelectrodes);
319 while (data.keySet().contains(channel))
320 {
321 channel = _random.nextInt(nelectrodes);
322 }
323
324 // Calculate the noise for this channel in units of electrons
325 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
326
327 // Throw Gaussian above threshold
328 int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
329 data.add(channel, new SiElectrodeData(charge));
330 }
331
332 // Now throw to lower threshold on channels that neighbor hits until we are exhausted
333 //-----------------------------------------------------------------------------------
334 nchannels_throw = 1;
335 while (nchannels_throw > 0)
336 {
337 // System.out.println("\n"+"Throw nieghbors...");
338
339 // Get neighbor channels
340 Set<Integer> neighbors = new HashSet<Integer>();
341 for (int channel : data.keySet())
342 {
343 neighbors.addAll(electrodes.getNearestNeighborCells(channel));
344 }
345 neighbors.removeAll(data.keySet());
346
347 nelectrodes_empty = neighbors.size();
348
349 // Get the noise threshold in units of the noise charge
350 normalized_integration_limit = _neighbor_threshold;
351
352 integral = Erf.phic(normalized_integration_limit);
353 nchannels_throw = drawBinomial(nelectrodes_empty, integral);
354
355 // Now throw Gaussian randoms above a threshold and put signals on unoccupied channels
356 for (int ithrow = 0; ithrow < nchannels_throw; ithrow++)
357 {
358 // Throw to get a channel number
359 List<Integer> neighbors_list = new ArrayList<Integer>(neighbors);
360
361 int channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
362
363 while (data.keySet().contains(channel))
364 {
365 channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
366 }
367
368 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
369
370
371 // Throw Gaussian above threshold
372 int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
373 data.add(channel, new SiElectrodeData(charge));
374 }
375
376 }
377
378 }
379
380 public static int drawBinomial(int ntrials, double probability)
381 {
382 _binomial.setNumberOfTrials(ntrials);
383 _binomial.setProbabilityOfSuccess(probability);
384
385 int nsuccess = 0;
386 try
387 {
388 nsuccess = _binomial.inverseCumulativeProbability(_random.nextDouble());
389 } catch (MathException exception)
390 {
391 throw new RuntimeException("BasicReadoutChip failed to calculate inverse cumulative probability of binomial!");
392 }
393 return nsuccess;
394 }
395
396 /**
397 * Return a random variable following normal distribution, but beyond
398 * threshold provided during initialization.
399 */
400 public static double drawGaussianAboveThreshold(double prob_above_threshold)
401 {
402 double draw, cumulative_probability;
403
404 draw = prob_above_threshold * _random.nextDouble();
405 cumulative_probability = 1.0 - prob_above_threshold + draw;
406
407 assert cumulative_probability < 1.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
408 assert cumulative_probability >= 0.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
409
410 double gaussian_random = 0;
411 try
412 {
413 gaussian_random = _gaussian.inverseCumulativeProbability(cumulative_probability);
414 } catch (MathException e)
415 {
416 System.out.println("MathException caught: " + e);
417 }
418
419 return gaussian_random;
420 }
421
422 /**
423 * BasicChannel class representing a single channel's behavior
424 *
425 * Note that binary readout is a special case. Anything positive value
426 * passed to a binary ADC for digitization is assumed to have crossed t
427 * hreshold and is assigned a value of 1. Decoding binary readout results
428 * in either 0 or dynamic_range.
429 */
430 private class BasicChannel implements ReadoutChannel
431 {
432
433 private double _noise_intercept = 0.;
434 private double _noise_slope = 0.;
435
436 /**
437 * Set the noise (in electrons) for 0 capacitance
438 *
439 * @param noise_intercept noise intercept
440 */
441 private void setNoiseIntercept(double noise_intercept)
442 {
443 _noise_intercept = noise_intercept;
444 }
445
446 /**
447 * Set the capacitative noise slope (in electrons / pF)
448 *
449 * @param noise_slope noise slope
450 */
451 private void setNoiseSlope(double noise_slope)
452 {
453 _noise_slope = noise_slope;
454 }
455
456 /**
457 * Return the noise in electrons for a given strip/pixel capacitance
458 *
459 * @param capacitance capacitance in pF
460 * @return noise in electrons
461 */
462 public double computeNoise(double capacitance)
463 {
464 return _noise_intercept + _noise_slope * capacitance;
465 }
466 }
467
468 /**
469 * ADC class representing analog to digital converter.
470 */
471 private class ADC
472 {
473
474 private int _nbits = 8;
475 private double _dynamic_range = 20.;
476
477 /**
478 * Set the ADC resolution in number of bits.
479 *
480 * @param nbits number of bits
481 */
482 private void setNbits(int nbits)
483 {
484 _nbits = nbits;
485 }
486
487 /**
488 * Set the dynamic range in fC
489 *
490 * @param dynamic range
491 */
492 private void setDynamicRange(double dynamic_range)
493 {
494 _dynamic_range = dynamic_range;
495 }
496
497 /**
498 * Compute the maximum ADC value
499 *
500 * @return largest possible ADC value according to # of bits
501 */
502 private int maxADCValue()
503 {
504 return (int) Math.pow(2, _nbits) - 1;
505 }
506
507 /**
508 * Compute the conversion constant in ADC/fC
509 *
510 * @return conversion constant for ADC
511 */
512 private double conversionConstant()
513 {
514 return maxADCValue() / _dynamic_range;
515 }
516
517 /**
518 * Perform analog to digital conversion
519 *
520 * @return digital ADC output between 0 and maxADCValue
521 */
522 public int convert(double charge)
523 {
524 if (_nbits != 1)
525 {
526 return Math.max(0, Math.min(maxADCValue(), (int) Math.floor(charge * 1.602e-4 * conversionConstant())));
527 }
528 else
529 {
530 if (charge <= 0.0)
531 {
532 return 0;
533 }
534 else
535 {
536 return 1;
537 }
538 }
539 }
540
541 /**
542 * Decode charge from ADC value
543 *
544 * @return charge specified by a given ADC value
545 */
546 public double decodeCharge(int adc_value)
547 {
548 if (_nbits != 1)
549 {
550 return (adc_value + 0.5) / (1.602e-4 * conversionConstant());
551 }
552 else
553 {
554 return adc_value*_dynamic_range;
555 }
556
557 }
558 }
559 }