View Javadoc

1   /*
2    * Class GenericReadoutChip
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   * Generic 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 programmable conversion constant.  Noise is added to strips with charge,
28   * and random noise hits are generated as well.  Methods are provided to decode
29   * the charge and time (although the current implementation always returns a time
30   * of 0).
31   *
32   * @author Richard Partridge
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      /** Creates a new instance of GenericReadoutChip */
44      public GenericReadoutChip() {
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      }
75  
76      /**
77       * Set the threshold for reading a channel if it's neighbor is
78       * above the noise threshold.  Units are electrons.
79       *
80       * @param neighbor_threshold
81       */
82      public void setNeighborThreshold(double neighbor_threshold) {
83          _neighbor_threshold = neighbor_threshold;
84      }
85  
86      /**
87       * Set the conversion between charge and ADC counts.  Units are ADC counts
88       * per fC of charge.
89       *
90       * @param adc_per_fC
91       */
92      public void setConversionConstant(double adc_per_fC) {
93          _channel.setConversionConstant(adc_per_fC);
94      }
95  
96      /**
97       * Return the GenericChannel associated with a given channel number.
98       * For the generic readout, there is a single instance of GenericChannel
99       * and thus the channel number is ignored.
100      *
101      * @param channel_number channel number
102      * @return associated GenericReadoutChannel
103      */
104     public GenericChannel 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
134      *
135      * @param hit raw hit
136      * @return hit charge in units of electrons
137      */
138     public double decodeCharge(RawTrackerHit hit) {
139 
140         //  Get the ADC value
141         int adc = hit.getADCValues()[0];
142 
143         //  Return the charge in electrons associated with the ADC value
144         return adc / (1.6e-4 * _channel.getConversionConstant());
145     }
146 
147     /**
148      * Decode the hit time.  Currently, the generic readout chip ignores the
149      * hit time and returns 0.
150      *
151      * @param hit raw hit data
152      * @return hit time
153      */
154     public int decodeTime(RawTrackerHit hit) {
155         return 0;
156     }
157 
158     /**
159      * Digitizes the hit channels in a SiElectrodeDataCollection.
160      *
161      * The SiElectrodeDataCollection is a map that associates a given channel with
162      * it's SiElectrodeData.  The SiElectrodeData encapsulates the deposited charge
163      * on an strip/pixel and any associated SimTrackerHits.
164      *
165      * The output of this class is a map that associates a channel number with
166      * a list of raw data
167      *
168      * @param data electrode data collection
169      * @return map associating channels with a list of raw data
170      */
171     private SortedMap<Integer, List<Integer>> digitize(SiElectrodeDataCollection data,
172             SiSensorElectrodes electrodes) {
173 
174         //  Create the map that associates a given sensor channel with it's list of raw data
175         SortedMap<Integer, List<Integer>> chip_data = new TreeMap<Integer, List<Integer>>();
176 
177         //  Loop over the channels contained in the SiElectrodeDataCollection
178         for (Integer channel : data.keySet()) {
179 
180             //  Fetch the electrode data for this channel
181             SiElectrodeData eldata = data.get(channel);
182 
183             //  Get the charge in units of electrons
184             double charge = eldata.getCharge();
185 
186             //  If the charge is below the neighbor threshold, don't digitize it
187             if (charge < _neighbor_threshold) continue;
188 
189             //  If charge is between neighbor and noise thresholds, check it's neighbors
190             if (charge < _noise_threshold) {
191 
192                 //  Loop over neighbors and look for a neighbor with charge above the noise
193                 boolean nbrhit = false;
194                 for (Integer nbr : electrodes.getNearestNeighborCells(channel)) {
195 
196                     //  See if we have electrode data for this neighbor
197                     SiElectrodeData nbrdata = data.get(nbr);
198                     if (nbrdata == null) continue;
199 
200                     //  See if we have found a neigbor above the noise threshold
201                     if (nbrdata.getCharge() >= _noise_threshold) {
202                         nbrhit = true;
203                         break;
204                     }
205                 }
206                 
207                 //  If there were no neighbor channels above threshold, don't digitize it
208                 if (!nbrhit) continue;
209             }
210 
211             //  Calculate the ADC value for this channel and make sure it is positive
212             int adc = getChannel(channel).computeAdcValue(charge);
213             if (adc <= 0) {
214                 continue;
215             }
216 
217             //  Create a list containing the adc value - for the generic readout
218             //  there is only 1 word of raw data
219             List<Integer> channel_data = new ArrayList<Integer>();
220             channel_data.add(adc);
221 
222             //  Save the list of raw data in the chip_data map
223             chip_data.put(channel, channel_data);
224         }
225 
226         return chip_data;
227     }
228 
229     /**
230      * Add noise hits for this readout chip
231      *
232      * @param data electrode data collection
233      * @param electrodes strip or pixel electrodes
234      */
235     private void addNoise(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
236 
237         //  First add noise to the strips/pixels in the SiElectrodeDataCollection
238         //  Loop over the entries in the SiElectrodeDataCollection (which extends TreeMap)
239         for (Entry datum : data.entrySet()) {
240 
241             //  Get the channel number and electrode data for this entry
242             int channel = (Integer) datum.getKey();
243             SiElectrodeData eldata = (SiElectrodeData) datum.getValue();
244 
245             //  Get the RMS noise for this channel in units of electrons
246             double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
247 
248             //  Add readout noise to the deposited charge
249             int noise_charge = (int) Math.round(_random.nextGaussian() * noise);
250             eldata.addCharge(noise_charge);
251         }
252 
253         //  Add random noise hits where the noise charge exceeds the noise threshold
254 
255         //  Find the number of pixels/strips that are not currently hit
256         int nelectrodes = electrodes.getNCells();
257         int nelectrodes_empty = nelectrodes - data.size();
258 
259         //  Get the noise threshold in units of the noise charge
260         double noiseRMS = _channel.computeNoise(electrodes.getCapacitance());
261         double normalized_integration_limit = _noise_threshold / noiseRMS;
262 
263         //  Calculate how many channels should get noise hits
264         double integral = Erf.phic(normalized_integration_limit);
265         int nchannels_throw = drawBinomial(nelectrodes_empty, integral);
266 
267         // Now throw Gaussian randoms above the seed threshold and put signals on unoccupied channels
268         for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
269             // Throw to get a channel number
270             int channel = _random.nextInt(nelectrodes);
271             while (data.keySet().contains(channel)) {
272                 channel = _random.nextInt(nelectrodes);
273             }
274 
275             //  Calculate the noise for this channel in units of electrons
276             double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
277 
278             // Throw Gaussian above threshold
279             int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
280             data.add(channel, new SiElectrodeData(charge));
281         }
282 
283         // Now throw to lower threshold on channels that neighbor hits until we are exhausted
284         //-----------------------------------------------------------------------------------
285         nchannels_throw = 1;
286         while (nchannels_throw > 0) {
287             //            System.out.println("\n"+"Throw nieghbors...");
288 
289             // Get neighbor channels
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             //  Get the noise threshold in units of the noise charge
299             normalized_integration_limit = _neighbor_threshold / noiseRMS;
300 
301             integral = Erf.phic(normalized_integration_limit);
302             nchannels_throw = drawBinomial(nelectrodes_empty, integral);
303 
304             // Now throw Gaussian randoms above a threshold and put signals on unoccupied channels
305             for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
306                 // Throw to get a channel number
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                 // Throw Gaussian above threshold
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      * Return a random variable following normal distribution, but beyond
342      * threshold provided during initialization.
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      * GenericChannel class representing a single channel's behavior
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          * Set the conversion between ADC counts and charge in fC
376          *
377          * @param adc_per_fC conversion constant
378          */
379         private void setConversionConstant(double adc_per_fC) {
380             _adc_per_fC = adc_per_fC;
381         }
382 
383         /**
384          * Return the conversion constant between ADC counts and charge in fC
385          *
386          * @return conversion constant
387          */
388         private double getConversionConstant() {
389             return _adc_per_fC;
390         }
391 
392         /**
393          * Set the noise (in electrons) for 0 capacitance
394          *
395          * @param noise_intercept noise intercept
396          */
397         private void setNoiseIntercept(double noise_intercept) {
398             _noise_intercept = noise_intercept;
399         }
400 
401         /**
402          * Set the capacitative noise slope (in electrons / pF)
403          *
404          * @param noise_slope noise slope
405          */
406         private void setNoiseSlope(double noise_slope) {
407             _noise_slope = noise_slope;
408         }
409 
410         /**
411          * Return the noise in electrons for a given strip/pixel capacitance
412          *
413          * @param capacitance capacitance in pF
414          * @return noise in electrons
415          */
416         public double computeNoise(double capacitance) {
417             return _noise_intercept + _noise_slope * capacitance;
418         }
419 
420         /**
421          * Calculate the ADC value associated with a pixel/strip charge deposit
422          * 
423          * @param data electrode data
424          * @return charge
425          */
426         private int computeAdcValue(double charge) {
427 
428             //  Convert from electrons to ADC counts (1 electron = 1.6 x 10^-4 fC)
429             int adc = (int) Math.floor(charge * 1.6e-4 * _adc_per_fC);
430 
431             //  Don't allow negative adc values
432             if (adc < 0) {
433                 adc = 0;
434             }
435 
436             //  Check for overflow - will be stored as a 16 bit integer
437             if (adc > 32767) {
438                 adc = 32767;
439             }
440 
441             return adc;
442         }
443     }
444 }