lcsim/src/org/lcsim/recon/tracking/digitization/sistripsim
diff -u -r1.1 -r1.2
--- GenericReadoutChip.java 22 Apr 2009 23:10:38 -0000 1.1
+++ GenericReadoutChip.java 23 Apr 2009 23:05:23 -0000 1.2
@@ -4,16 +4,21 @@
package org.lcsim.recon.tracking.digitization.sistripsim;
import java.util.ArrayList;
+import java.util.HashSet;
import java.util.List;
+import java.util.Map.Entry;
import java.util.Random;
+import java.util.Set;
import java.util.SortedMap;
import java.util.TreeMap;
+import org.apache.commons.math.MathException;
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.math.probability.Erf;
import org.lcsim.recon.tracking.digitization.sistripsim.ReadoutChip.ReadoutChannel;
/**
@@ -111,7 +116,9 @@
// If there is no electrode data for this readout chip, create an empty
// electrode data collection
- if (data == null) data = new SiElectrodeDataCollection();
+ if (data == null) {
+ data = new SiElectrodeDataCollection();
+ }
// Add noise hits to the electrode data collection
// if (_strip_clusterer != null) {
@@ -175,7 +182,9 @@
// Calculate the ADC value for this channel and make sure it is positive
int adc = getChannel(channel).computeAdcValue(eldata);
- if (adc <= 0) continue;
+ if (adc <= 0) {
+ continue;
+ }
// Create a list containing the adc value - for the generic readout
// there is only 1 word of raw data
@@ -189,158 +198,136 @@
return chip_data;
}
- /*
+ /**
+ * Add noise hits for this readout chip
+ *
+ * @param data electrode data collection
+ * @param electrodes strip or pixel electrodes
+ */
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));
- }
+ // First add noise to the strips/pixels in the SiElectrodeDataCollection
+ // Loop over the entries in the SiElectrodeDataCollection (which extends TreeMap)
+ for (Entry datum : data.entrySet()) {
+
+ // Get the channel number and electrode data for this entry
+ int channel = (Integer) datum.getKey();
+ SiElectrodeData eldata = (SiElectrodeData) datum.getValue();
+
+ // Get the RMS noise for this channel in units of electrons
+ double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
+
+ // Add readout noise to the deposited charge
+ int noise_charge = (int) Math.round(_random.nextGaussian() * noise);
+ eldata.addCharge(noise_charge);
+ }
- double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
+ // Add random noise hits where the noise charge exceeds the cluster seed threshold
- // System.out.println(" noise: "+noise);
- // System.out.println(" Gaussian above threshold: "+drawGaussianAboveThreshold(integral));
+ // Find the number of pixels/strips that are not currently hit
+ int nelectrodes = electrodes.getNCells();
+ int nelectrodes_empty = nelectrodes - data.size();
+
+ // Get the seed threshold from the clusterer
+ double normalized_integration_limit = _clusterer.getSeedThreshold();
+
+ // Calculate how many channels should get noise hits
+ double integral = Erf.phic(normalized_integration_limit);
+ int nchannels_throw = drawBinomial(nelectrodes_empty, integral);
+
+ // Now throw Gaussian randoms above the seed 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);
+ }
+
+ // Calculate the noise for this channel in units of electrons
+ double noise = getChannel(channel).computeNoise(electrodes.getCapacitance(channel));
+
+ // Throw Gaussian above threshold
+ int charge = (int) Math.round(drawGaussianAboveThreshold(integral) * noise);
+ data.add(channel, new SiElectrodeData(charge));
+ }
- // 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 = _clusterer.getNeighborThreshold(); // We should get this from same place as clustering code
+
+ integral = Erf.phic(normalized_integration_limit);
+ nchannels_throw = drawBinomial(nelectrodes_empty, integral);
+
+ // 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));
+
+
+ // 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);
+ _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;
+ 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;
+ 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;
+ 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;
+ 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);
- }
+ double gaussian_random = 0;
+ try {
+ gaussian_random = _gaussian.inverseCumulativeProbability(cumulative_probability);
+ } catch (MathException e) {
+ System.out.println("MathException caught: " + e);
+ }
- return gaussian_random;
+ return gaussian_random;
}
- */
/**
* GenericChannel class representing a single channel's behavior
@@ -412,10 +399,14 @@
int adc = (int) Math.floor(charge * 1.6e-4 * _adc_per_fC);
// Don't allow negative adc values
- if (adc < 0) adc = 0;
+ if (adc < 0) {
+ adc = 0;
+ }
// Check for overflow - will be stored as a 16 bit integer
- if (adc > 32767) adc = 32767;
+ if (adc > 32767) {
+ adc = 32767;
+ }
return adc;
}