lcsim/src/org/lcsim/recon/tracking/digitization/sistripsim
diff -N GenericReadoutChip.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ GenericReadoutChip.java 22 Apr 2009 23:10:38 -0000 1.1
@@ -0,0 +1,423 @@
+/*
+ * Class GenericReadoutChip
+ */
+package org.lcsim.recon.tracking.digitization.sistripsim;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Random;
+import java.util.SortedMap;
+import java.util.TreeMap;
+import org.apache.commons.math.distribution.BinomialDistribution;
+import org.apache.commons.math.distribution.BinomialDistributionImpl;
+import org.apache.commons.math.distribution.NormalDistribution;
+import org.apache.commons.math.distribution.NormalDistributionImpl;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.recon.tracking.digitization.sistripsim.ReadoutChip.ReadoutChannel;
+
+/**
+ * Generic readout chip class. This class supports the minimal functions expected of
+ * a readout chip. The charge on a strip/pixel is digitized as an integer number
+ * with a programmable conversion constant. Noise is added to strips with charge,
+ * and random noise hits are generated as well. Methods are provided to decode
+ * the charge and time (although the current implementation always returns a time
+ * of 0).
+ *
+ * @author Richard Partridge
+ */
+public class GenericReadoutChip implements ReadoutChip {
+
+ private static Random _random = new Random();
+ private static NormalDistribution _gaussian = new NormalDistributionImpl(0.0, 1.0);
+ private static BinomialDistribution _binomial = new BinomialDistributionImpl(1, 1);
+ private StripClusterer _clusterer;
+ private GenericChannel _channel = new GenericChannel();
+ private double _capacitance;
+
+ /** Creates a new instance of GenericReadoutChip */
+ public GenericReadoutChip() {
+ }
+
+ /**
+ * Set the noise intercept (i.e., the noise for 0 strip/pixel capacitance).
+ * Units are electrons of noise.
+ *
+ * @param noise_intercept noise for 0 capacitance
+ */
+ public void setNoiseIntercept(double noise_intercept) {
+ _channel.setNoiseIntercept(noise_intercept);
+ }
+
+ /**
+ * Set the noise slope (i.e., the proportionality between noise and capacitance).
+ * Units are electrons of noise per fF of capacitance.
+ *
+ * @param noise_slope noise slope per unit capacitance
+ */
+ public void setNoiseSlope(double noise_slope) {
+ _channel.setNoiseSlope(noise_slope);
+ }
+
+ /**
+ * Set the capacitance of the strip/pixels. Units are fF.
+ *
+ * @param capacitance
+ */
+ public void setCapacitance(double capacitance) {
+ _capacitance = capacitance;
+ }
+
+ /**
+ * Set the conversion between charge and ADC counts. Units are ADC counts
+ * per fC of charge.
+ *
+ * @param adc_per_fC
+ */
+ public void setConversionConstant(double adc_per_fC) {
+ _channel.setConversionConstant(adc_per_fC);
+ }
+
+ /**
+ * Set the clusterer to be used for this readout chip.
+ *
+ * @param _clusterer
+ */
+ public void setStripClusterer(StripClusterer clusterer) {
+ _clusterer = clusterer;
+ }
+
+ /**
+ * Return the GenericChannel associated with a given channel number.
+ * For the generic readout, there is a single instance of GenericChannel
+ * and thus the channel number is ignored.
+ *
+ * @param channel_number channel number
+ * @return associated GenericReadoutChannel
+ */
+ public GenericChannel getChannel(int channel_number) {
+ return _channel;
+ }
+
+ /**
+ * Given a collection of electrode data (i.e., charge on strips/pixels),
+ * return a map associating the channel and it's list of raw data.
+ *
+ * @param data electrode data from the charge distribution
+ * @param electrodes strip or pixel electrodes
+ * @return map containing the ADC counts for this sensor
+ */
+ public SortedMap<Integer, List<Integer>> readout(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
+
+ // If there is no electrode data for this readout chip, create an empty
+ // electrode data collection
+ if (data == null) data = new SiElectrodeDataCollection();
+
+ // Add noise hits to the electrode data collection
+// if (_strip_clusterer != null) {
+// addNoise(data, electrodes);
+// }
+
+ // return the digitized charge data as a map that associates a hit
+ // channel with a list of raw data for the channel
+ return digitize(data);
+ }
+
+ /**
+ * Decode the hit charge stored in the RawTrackerHit
+ *
+ * @param hit raw hit
+ * @return hit charge
+ */
+ public double decodeCharge(RawTrackerHit hit) {
+
+ // Get the ADC value
+ int adc = hit.getADCValues()[0];
+
+ // Return the charge associated with the ADC value
+ return adc / _channel.getConversionConstant();
+ }
+
+ /**
+ * Decode the hit time. Currently, the generic readout chip ignores the
+ * hit time and returns 0.
+ *
+ * @param hit raw hit data
+ * @return hit time
+ */
+ public int decodeTime(RawTrackerHit hit) {
+ return 0;
+ }
+
+ /**
+ * Digitizes the hit channels in a SiElectrodeDataCollection.
+ *
+ * The SiElectrodeDataCollection is a map that associates a given channel with
+ * it's SiElectrodeData. The SiElectrodeData encapsulates the deposited charge
+ * on an strip/pixel and any associated SimTrackerHits.
+ *
+ * The output of this class is a map that associates a channel number with
+ * a list of raw data
+ *
+ * @param data electrode data collection
+ * @return map associating channels with a list of raw data
+ */
+ private SortedMap<Integer, List<Integer>> digitize(SiElectrodeDataCollection data) {
+
+ // Create the map that associates a given sensor channel with it's list of raw data
+ SortedMap<Integer, List<Integer>> chip_data = new TreeMap<Integer, List<Integer>>();
+
+ // Loop over the channels contained in the SiElectrodeDataCollection
+ for (Integer channel : data.keySet()) {
+
+ // Fetch the electrode data for this channel
+ SiElectrodeData eldata = data.get(channel);
+
+ // Calculate the ADC value for this channel and make sure it is positive
+ int adc = getChannel(channel).computeAdcValue(eldata);
+ if (adc <= 0) continue;
+
+ // Create a list containing the adc value - for the generic readout
+ // there is only 1 word of raw data
+ List<Integer> channel_data = new ArrayList<Integer>();
+ channel_data.add(adc);
+
+ // Save the list of raw data in the chip_data map
+ chip_data.put(channel, channel_data);
+ }
+
+ return chip_data;
+ }
+
+ /*
+ private void addNoise(SiElectrodeDataCollection data, SiSensorElectrodes electrodes) {
+
+ // System.out.println("\n"+"Adding noise...");
+
+ // Add full noise distribution to any cells with charge deposition
+ //----------------------------------------------------------------
+ for (Entry datum : data.entrySet()) {
+ int channel = (Integer) datum.getKey();
+ double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
+ int origCharge = ((SiElectrodeData) datum.getValue()).getCharge();
+ int addedNoise = (int) Math.round(_random.nextGaussian() * noise);
+ // System.out.println("Kpix::addNoise channel " + channel + " charge = " + origCharge + " noise = " + addedNoise);
+ if (addedNoise + origCharge < 0) {
+ // System.out.println("Kpix::addNoise preventing charge from going negative");
+ addedNoise = -origCharge;
+ }
+ ((SiElectrodeData) datum.getValue()).addCharge(addedNoise);
+ // System.out.println("Kpix::addNoise new charge = " + ((SiElectrodeData) datum.getValue()).getCharge());
+ }
+
+ // Throw cluster seeds on all channels
+ //------------------------------------
+ // System.out.println("\n"+"Throw flyers...");
+
+ int nelectrodes = electrodes.getNCells();
+ int nelectrodes_empty = nelectrodes - data.size();
+ double normalized_integration_limit = _strip_clusterer.getSeedThreshold(); // We should get this from same place as clustering code
+
+ double integral = normalCDF(normalized_integration_limit);
+ int nchannels_throw = drawBinomial(nelectrodes_empty, integral);
+
+ // System.out.println(" # Empty channels: "+nelectrodes_empty);
+ // System.out.println(" "+normalized_integration_limit+"-sigma integral: "+integral);
+ // System.out.println(" Mean # channels: "+nelectrodes_empty*integral);
+ // System.out.println(" Binomial draw: "+nchannels_throw);
+
+ // Now throw Gaussian randoms above a threshold and put signals on unoccupied channels
+ for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
+ // Throw to get a channel number
+ int channel = _random.nextInt(nelectrodes);
+ while (data.keySet().contains(channel)) {
+ channel = _random.nextInt(nelectrodes);
+ }
+
+ double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
+
+ // System.out.println(" noise: "+noise);
+ // System.out.println(" Gaussian above threshold: "+drawGaussianAboveThreshold(integral));
+
+ // Throw Gaussian above threshold
+ int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
+ data.add(channel, new SiElectrodeData(charge));
+ }
+
+ // Now throw to lower threshold on channels that neighbor hits until we are exhausted
+ //-----------------------------------------------------------------------------------
+ nchannels_throw = 1;
+ while (nchannels_throw > 0) {
+ // System.out.println("\n"+"Throw nieghbors...");
+
+ // Get neighbor channels
+ Set<Integer> neighbors = new HashSet<Integer>();
+ for (int channel : data.keySet()) {
+ neighbors.addAll(electrodes.getNearestNeighborCells(channel));
+ }
+ neighbors.removeAll(data.keySet());
+
+ nelectrodes_empty = neighbors.size();
+ normalized_integration_limit = _strip_clusterer.getNeighborThreshold(); // We should get this from same place as clustering code
+
+ integral = normalCDF(normalized_integration_limit);
+ nchannels_throw = drawBinomial(nelectrodes_empty, integral);
+
+ // System.out.println(" # Empty channels: "+nelectrodes_empty);
+ // System.out.println(" "+normalized_integration_limit+"-sigma integral: "+integral);
+ // System.out.println(" Mean # channels: "+nelectrodes_empty*integral);
+ // System.out.println(" Binomial draw: "+nchannels_throw);
+
+ // Now throw Gaussian randoms above a threshold and put signals on unoccupied channels
+ for (int ithrow = 0; ithrow < nchannels_throw; ithrow++) {
+ // Throw to get a channel number
+ List<Integer> neighbors_list = new ArrayList<Integer>(neighbors);
+
+ int channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
+
+ while (data.keySet().contains(channel)) {
+ channel = neighbors_list.get(_random.nextInt(nelectrodes_empty));
+ }
+
+ double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
+
+ // System.out.println(" noise: "+noise);
+ // System.out.println(" Gaussian above threshold: "+drawGaussianAboveThreshold(integral));
+
+ // Throw Gaussian above threshold
+ int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
+ data.add(channel, new SiElectrodeData(charge));
+ }
+
+ }
+
+ }
+ */
+ /*
+ public static double normalCDF(double normalized_integration_limit) {
+ double integral = 0;
+ try {
+ integral = (1.0 - Erf.erf(normalized_integration_limit / Math.sqrt(2.0))) / 2.0;
+ } catch (MathException no_convergence) {
+ System.out.println("Warning: erf fails to converge!! ");
+ System.out.println(" normalized integration limit: " + normalized_integration_limit);
+ }
+ return integral;
+ }
+
+ public static int drawBinomial(int ntrials, double probability) {
+ _binomial.setNumberOfTrials(ntrials);
+ _binomial.setProbabilityOfSuccess(probability);
+
+ int nsuccess = 0;
+ try {
+ nsuccess = _binomial.inverseCumulativeProbability(_random.nextDouble());
+ } catch (MathException exception) {
+ throw new RuntimeException("Kpix failed to calculate inverse cumulative probability of binomial!");
+ }
+ return nsuccess;
+ }
+
+ /**
+ * Return a random variable following normal distribution, but beyond
+ * threshold provided during initialization.
+ */
+ /* public static double drawGaussianAboveThreshold(double prob_above_threshold) {
+ double draw, cumulative_probability;
+
+ draw = prob_above_threshold * _random.nextDouble();
+ cumulative_probability = 1.0 - prob_above_threshold + draw;
+
+ assert cumulative_probability < 1.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
+ assert cumulative_probability >= 0.0 : "cumulProb=" + cumulative_probability + ", draw=" + draw + ", probAboveThreshold=" + prob_above_threshold;
+
+ double gaussian_random = 0;
+ try {
+ gaussian_random = _gaussian.inverseCumulativeProbability(cumulative_probability);
+ } catch (MathException e) {
+ System.out.println("MathException caught: " + e);
+ }
+
+ return gaussian_random;
+ }
+ */
+
+ /**
+ * GenericChannel class representing a single channel's behavior
+ */
+ private class GenericChannel implements ReadoutChannel {
+
+ private double _noise_intercept = 0.;
+ private double _noise_slope = 0.;
+ private double _adc_per_fC = 100.;
+
+ /**
+ * Set the conversion between ADC counts and charge in fC
+ *
+ * @param adc_per_fC conversion constant
+ */
+ private void setConversionConstant(double adc_per_fC) {
+ _adc_per_fC = adc_per_fC;
+ }
+
+ /**
+ * Return the conversion constant between ADC counts and charge in fC
+ *
+ * @return conversion constant
+ */
+ private double getConversionConstant() {
+ return _adc_per_fC;
+ }
+
+ /**
+ * Set the noise (in electrons) for 0 capacitance
+ *
+ * @param noise_intercept noise intercept
+ */
+ private void setNoiseIntercept(double noise_intercept) {
+ _noise_intercept = noise_intercept;
+ }
+
+ /**
+ * Set the capacitative noise slope (in electrons / fF)
+ *
+ * @param noise_slope noise slope
+ */
+ private void setNoiseSlope(double noise_slope) {
+ _noise_slope = noise_slope;
+ }
+
+ /**
+ * Return the noise for a given strip/pixel capacitance
+ *
+ * @param capacitance capacitance in fF
+ * @return noise in electrons
+ */
+ public double computeNoise(double capacitance) {
+ return _noise_intercept + _noise_slope * capacitance;
+ }
+
+ /**
+ * Calculate the ADC value associated with a pixel/strip charge deposit
+ *
+ * @param data electrode data
+ * @return charge
+ */
+ private int computeAdcValue(SiElectrodeData data) {
+
+ // Get the charge in units of electrons
+ double charge = data.getCharge();
+
+ // Convert from electrons to ADC counts (1 electron = 1.6 x 10^-4 fC)
+ int adc = (int) Math.floor(charge * 1.6e-4 * _adc_per_fC);
+
+ // Don't allow negative adc values
+ if (adc < 0) adc = 0;
+
+ // Check for overflow - will be stored as a 16 bit integer
+ if (adc > 32767) adc = 32767;
+
+ return adc;
+ }
+ }
+}
