1 /*
2 * Class DigitalReadoutChip
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 * Digital readout chip class. This class supports the minimal functions
26 * expected of a digital readout chip. The charge on a strip/pixel is digitized
27 * to give a specified ADC value if the charge is above the noise threshold.
28 * Noise is added to strips/pixels with charge, and random noise hits are generated
29 * as well. Methods are provided to decode the charge and time (although the
30 * current implementation always returns a time of 0).
31 *
32 * @author Richard Partridge
33 */
34 public class DigitalReadoutChip 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 DigitalChannel _channel = new DigitalChannel();
42
43 /** Creates a new instance of DigitalReadoutChip */
44 public DigitalReadoutChip() {
45 }
46
47 /**
48 * Set the noise intercept (i.e., the noise for 0 strip/pixel capacitance).
49 * Units are electrons of noise.
50 *
51 * @param noise_intercept noise for 0 capacitance
52 */
53 public void setNoiseIntercept(double noise_intercept) {
54 _channel.setNoiseIntercept(noise_intercept);
55 }
56
57 /**
58 * Set the noise slope (i.e., the proportionality between noise and capacitance).
59 * Units are electrons of noise per fF of capacitance.
60 *
61 * @param noise_slope noise slope per unit capacitance
62 */
63 public void setNoiseSlope(double noise_slope) {
64 _channel.setNoiseSlope(noise_slope);
65 }
66
67 /**
68 * Set the threshold for reading out a channel. Units are electrons.
69 *
70 * @param noise_threshold
71 */
72 public void setNoiseThreshold(double noise_threshold) {
73 _noise_threshold = noise_threshold;
74 _channel.setNoiseThreshold(noise_threshold);
75 }
76
77 /**
78 * Set the threshold for reading a channel if it's neighbor is
79 * above the noise threshold. Units are electrons.
80 *
81 * @param neighbor_threshold
82 */
83 public void setNeighborThreshold(double neighbor_threshold) {
84 _neighbor_threshold = neighbor_threshold;
85 }
86
87 /**
88 * Set the ADC output value for a hit.
89 *
90 * @param adc_for_hit
91 */
92 public void setConversionConstant(int adc_for_hit) {
93 _channel.setConversionConstant(adc_for_hit);
94 }
95
96 /**
97 * Return the DigitalChannel associated with a given channel number.
98 * For the digital readout, there is a single instance of DigitalChannel
99 * and thus the channel number is ignored.
100 *
101 * @param channel_number channel number
102 * @return associated DigitalChannel
103 */
104 public DigitalChannel getChannel(int channel_number) {
105 return _channel;
106 }
107
108 /**
109 * Given a collection of electrode data (i.e., charge on strips/pixels),
110 * return a map associating the channel and it's list of raw data.
111 *
112 * @param data electrode data from the charge distribution
113 * @param electrodes strip or pixel electrodes
114 * @return map containing the ADC counts for this sensor
115 */
116 public SortedMap<Integer, List<Integer>> readout(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
117
118 // If there is no electrode data for this readout chip, create an empty
119 // electrode data collection
120 if (data == null) {
121 data = new SiElectrodeDataCollection();
122 }
123
124 // Add noise hits to the electrode data collection
125 addNoise(data, electrodes);
126
127 // return the digitized charge data as a map that associates a hit
128 // channel with a list of raw data for the channel
129 return digitize(data, electrodes);
130 }
131
132 /**
133 * Decode the hit charge stored in the RawTrackerHit. Since this method
134 * isn't well defined for digital readout, return the noise threshold if
135 * there is a hit.
136 *
137 * @param hit raw hit
138 * @return hit pseudo-charge in units of electrons
139 */
140 public double decodeCharge(RawTrackerHit hit) {
141
142 // Get the ADC value
143 int adc = hit.getADCValues()[0];
144
145 // If there is a hit (i.e., adc > 0) return the charge as being the
146 // value of the noise threshold, otherwise return 0
147 if (adc > 0) return _noise_threshold;
148 else return 0.;
149 }
150
151 /**
152 * Decode the hit time. Currently, the digital readout chip ignores the
153 * hit time and returns 0.
154 *
155 * @param hit raw hit data
156 * @return hit time
157 */
158 public int decodeTime(RawTrackerHit hit) {
159 return 0;
160 }
161
162 /**
163 * Digitizes the hit channels in a SiElectrodeDataCollection.
164 *
165 * The SiElectrodeDataCollection is a map that associates a given channel with
166 * it's SiElectrodeData. The SiElectrodeData encapsulates the deposited charge
167 * on an strip/pixel and any associated SimTrackerHits.
168 *
169 * The output of this class is a map that associates a channel number with
170 * a list of raw data
171 *
172 * @param data electrode data collection
173 * @return map associating channels with a list of raw data
174 */
175 private SortedMap<Integer, List<Integer>> digitize(SiElectrodeDataCollection data,
176 SiSensorElectrodes electrodes) {
177
178 // Create the map that associates a given sensor channel with it's list of raw data
179 SortedMap<Integer, List<Integer>> chip_data = new TreeMap<Integer, List<Integer>>();
180
181 // Loop over the channels contained in the SiElectrodeDataCollection
182 for (Integer channel : data.keySet()) {
183
184 // Fetch the electrode data for this channel
185 SiElectrodeData eldata = data.get(channel);
186
187 // Get the charge in units of electrons
188 double charge = eldata.getCharge();
189
190 // If the charge is below the neighbor threshold, don't digitize it
191 if (charge < _neighbor_threshold) continue;
192
193 // If charge is between neighbor and noise thresholds, check it's neighbors
194 if (charge < _noise_threshold) {
195
196 // Loop over neighbors and look for a neighbor with charge above the noise
197 boolean nbrhit = false;
198 for (Integer nbr : electrodes.getNearestNeighborCells(channel)) {
199
200 // See if we have electrode data for this neighbor
201 SiElectrodeData nbrdata = data.get(nbr);
202 if (nbrdata == null) continue;
203
204 // See if we have found a neigbor above the noise threshold
205 if (nbrdata.getCharge() >= _noise_threshold) {
206 nbrhit = true;
207 break;
208 }
209 }
210
211 // If there were no neighbor channels above threshold, don't digitize it
212 if (!nbrhit) continue;
213 }
214
215 // Calculate the ADC value for this channel and make sure it is positive
216 int adc = getChannel(channel).computeAdcValue(charge);
217 if (adc <= 0) {
218 continue;
219 }
220
221 // Create a list containing the adc value - for the digital readout
222 // there is only 1 word of raw data
223 List<Integer> channel_data = new ArrayList<Integer>();
224 channel_data.add(adc);
225
226 // Save the list of raw data in the chip_data map
227 chip_data.put(channel, channel_data);
228 }
229
230 return chip_data;
231 }
232
233 /**
234 * Add noise hits for this readout chip
235 *
236 * @param data electrode data collection
237 * @param electrodes strip or pixel electrodes
238 */
239 private void addNoise(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
240
241 // First add noise to the strips/pixels in the SiElectrodeDataCollection
242 // Loop over the entries in the SiElectrodeDataCollection (which extends TreeMap)
243 for (Entry datum : data.entrySet()) {
244
245 // Get the channel number and electrode data for this entry
246 int channel = (Integer) datum.getKey();
247 SiElectrodeData eldata = (SiElectrodeData) datum.getValue();
248
249 // Get the RMS noise for this channel in units of electrons
250 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
251
252 // Add readout noise to the deposited charge
253 int noise_charge = (int) Math.round(_random.nextGaussian() * noise);
254 eldata.addCharge(noise_charge);
255 }
256
257 // Add random noise hits where the noise charge exceeds the noise threshold
258
259 // Find the number of pixels/strips that are not currently hit
260 int nelectrodes = electrodes.getNCells();
261 int nelectrodes_empty = nelectrodes - data.size();
262
263 // Get the noise threshold in units of the noise charge
264 double noiseRMS = _channel.computeNoise(electrodes.getCapacitance());
265 double normalized_integration_limit = _noise_threshold / noiseRMS;
266
267 // Calculate how many channels should get noise hits
268 double integral = Erf.phic(normalized_integration_limit);
269 int nchannels_throw = drawBinomial(nelectrodes_empty, integral);
270
271 // Now throw Gaussian randoms above the seed threshold and put signals on unoccupied channels
272 for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
273 // Throw to get a channel number
274 int channel = _random.nextInt(nelectrodes);
275 while (data.keySet().contains(channel)) {
276 channel = _random.nextInt(nelectrodes);
277 }
278
279 // Calculate the noise for this channel in units of electrons
280 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
281
282 // Throw Gaussian above threshold
283 int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
284 data.add(channel, new SiElectrodeData(charge));
285 }
286
287 // Now throw to lower threshold on channels that neighbor hits until we are exhausted
288 //-----------------------------------------------------------------------------------
289 nchannels_throw = 1;
290 while (nchannels_throw > 0) {
291 // System.out.println("\n"+"Throw nieghbors...");
292
293 // Get neighbor channels
294 Set<Integer> neighbors = new HashSet<Integer>();
295 for (int channel : data.keySet()) {
296 neighbors.addAll(electrodes.getNearestNeighborCells(channel));
297 }
298 neighbors.removeAll(data.keySet());
299
300 nelectrodes_empty = neighbors.size();
301
302 // Get the noise threshold in units of the noise charge
303 normalized_integration_limit = _neighbor_threshold / noiseRMS;
304
305 integral = Erf.phic(normalized_integration_limit);
306 nchannels_throw = drawBinomial(nelectrodes_empty, integral);
307
308 // Now throw Gaussian randoms above a threshold and put signals on unoccupied channels
309 for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
310 // Throw to get a channel number
311 List<Integer> neighbors_list = new ArrayList<Integer>(neighbors);
312
313 int channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
314
315 while (data.keySet().contains(channel)) {
316 channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
317 }
318
319 double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
320
321
322 // Throw Gaussian above threshold
323 int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
324 data.add(channel, new SiElectrodeData(charge));
325 }
326
327 }
328
329 }
330
331 public static int drawBinomial(int ntrials, double probability) {
332 _binomial.setNumberOfTrials(ntrials);
333 _binomial.setProbabilityOfSuccess(probability);
334
335 int nsuccess = 0;
336 try {
337 nsuccess = _binomial.inverseCumulativeProbability(_random.nextDouble());
338 } catch (MathException exception) {
339 throw new RuntimeException("Kpix failed to calculate inverse cumulative probability of binomial!");
340 }
341 return nsuccess;
342 }
343
344 /**
345 * Return a random variable following normal distribution, but beyond
346 * threshold provided during initialization.
347 *
348 * @param prob_above_threshold
349 * @return
350 */
351 public static double drawGaussianAboveThreshold(double prob_above_threshold) {
352 double draw, cumulative_probability;
353
354 draw = prob_above_threshold * _random.nextDouble();
355 cumulative_probability = 1.0 - prob_above_threshold + draw;
356
357 assert cumulative_probability < 1.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
358 assert cumulative_probability >= 0.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
359
360 double gaussian_random = 0;
361 try {
362 gaussian_random = _gaussian.inverseCumulativeProbability(cumulative_probability);
363 } catch (MathException e) {
364 System.out.println("MathException caught: " + e);
365 }
366
367 return gaussian_random;
368 }
369
370 /**
371 * DigitalChannel class representing a single channel's behavior
372 */
373 private class DigitalChannel implements ReadoutChannel {
374
375 private double _noise_intercept = 0.;
376 private double _noise_slope = 0.;
377 private double _noise_threshold = 100.;
378 private int _adc_for_hit = 100;
379
380 /**
381 * Set the conversion between ADC counts and charge in fC
382 *
383 * @param adc_per_fC conversion constant
384 */
385 private void setConversionConstant(int adc_for_hit) {
386 _adc_for_hit = adc_for_hit;
387 }
388
389 /**
390 * Return the conversion constant between ADC counts and charge in fC
391 *
392 * @return conversion constant
393 */
394 private double getConversionConstant() {
395 return _adc_for_hit;
396 }
397
398 /**
399 * Set the noise (in electrons) for 0 capacitance
400 *
401 * @param noise_intercept noise intercept
402 */
403 private void setNoiseIntercept(double noise_intercept) {
404 _noise_intercept = noise_intercept;
405 }
406
407 /**
408 * Set the capacitative noise slope (in electrons / pF)
409 *
410 * @param noise_slope noise slope
411 */
412 private void setNoiseSlope(double noise_slope) {
413 _noise_slope = noise_slope;
414 }
415
416 /**
417 * Set the noise threshold in units of electrons.
418 *
419 * @param noise_threshold
420 */
421 private void setNoiseThreshold(double noise_threshold) {
422 _noise_threshold = noise_threshold;
423 }
424
425 /**
426 * Return the noise in electrons for a given strip/pixel capacitance
427 *
428 * @param capacitance capacitance in pF
429 * @return noise in electrons
430 */
431 public double computeNoise(double capacitance) {
432 return _noise_intercept + _noise_slope * capacitance;
433 }
434
435 /**
436 * Calculate the ADC value associated with a pixel/strip charge deposit
437 *
438 * @param data electrode data
439 * @return charge
440 */
441 private int computeAdcValue(double charge) {
442
443 // Set to the specified value for a hit if we are above the
444 // noise threshold
445 int adc = 0;
446 if (charge > _noise_threshold) adc = _adc_for_hit;
447
448 // Don't allow negative adc values
449 if (adc < 0) {
450 adc = 0;
451 }
452
453 // Check for overflow - will be stored as a 16 bit integer
454 if (adc > 32767) {
455 adc = 32767;
456 }
457
458 return adc;
459 }
460 }
461 }