1 /*
2 * NearestNeighborRMS.java
3 */
4 package org.lcsim.recon.tracking.digitization.sisim;
5
6 import java.util.ArrayList;
7 import java.util.HashMap;
8 import java.util.LinkedList;
9 import java.util.List;
10 import java.util.Map;
11 import java.util.Set;
12 import org.lcsim.detector.identifier.IIdentifier;
13 import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
14 import org.lcsim.detector.tracker.silicon.SiTrackerIdentifierHelper;
15 import org.lcsim.event.RawTrackerHit;
16
17 /**
18 * This class uses a nearest neighbor algorithm to find clusters of hits on a
19 * set of silicon sensor electrodes (i.e., an instance of SiStrips or SiPixels).
20 *
21 * The algorithm first finds possible cluster seeds that are above the seed
22 * threshold. Starting from a seed, the neighbor strips/pixels are searched
23 * to see if we have hits above the neighbor threshold. Neighbor channels
24 * are added until we have a cluster where every neighboring channel for the
25 * cluster has charge lying below the neighbor threshold.
26 *
27 * Thresholds are specified in units of RMS noise of each electrode.
28 *
29 * HashMaps are used to speed up the algorithm when there are large numbers
30 * of hits, and the algorithm is expected to scale linearly with the number
31 * of RawTrackerHits on the electrodes.
32 *
33 * @author Tim Nelson (adapted from NearestNeighbor code of Rich Partridge)
34 */
35 public class NearestNeighborRMS implements ClusteringAlgorithm {
36
37 private static String _NAME = "NearestNeighborRMS";
38 private double _seed_threshold;
39 private double _neighbor_threshold;
40 private double _cluster_threshold;
41
42 /**
43 * Instantiate NearestNeighborRMS with specified thresholds.
44 * Seed threshold is the minimum charge to initiate a cluster. Neighbor
45 * threshold is the minimum charge to add a neighboring cell to a cluster.
46 * Cluster threshold is minimum charge of the entire cluster.
47 * All thresholds are in units of RMS noise of the channel(s).
48 *
49 * @param seed_threshold seed threhold
50 * @param neighbor_threshold neighbor threshold
51 * @param cluster_threshold cluster threshold
52 */
53 public NearestNeighborRMS(double seed_threshold, double neighbor_threshold, double cluster_threshold) {
54 _seed_threshold = seed_threshold;
55 _neighbor_threshold = neighbor_threshold;
56 _cluster_threshold = cluster_threshold;
57 }
58
59 /**
60 * Instantiate NearestNeighborRMS with default thresholds:
61 *
62 * seed_threshold = 4*RMS noise
63 * neighbor_threhold = 3*RMS noise
64 * cluster_threshold = 4*RMS noise
65 */
66 public NearestNeighborRMS() {
67 this(4.0, 3.0, 4.0);
68 }
69
70 /**
71 * Set the seed threshold. Units are RMS noise.
72 *
73 * @param seed_threshold seed threshold
74 */
75 public void setSeedThreshold(double seed_threshold) {
76 _seed_threshold = seed_threshold;
77 }
78
79 /**
80 * Set the neighbor threshold. Units are RMS noise.
81 *
82 * @param neighbor_threshold neighbor threshold
83 */
84 public void setNeighborThreshold(double neighbor_threshold) {
85 _neighbor_threshold = neighbor_threshold;
86 }
87
88 /**
89 * Set the cluster threshold. Units are RMS noise.
90 *
91 * @param cluster_threshold cluster threshold
92 */
93 public void setClusterThreshold(double cluster_threshold) {
94 _cluster_threshold = cluster_threshold;
95 }
96
97 /**
98 * Return the seed threshold. Units are RMS noise.
99 *
100 * @return seed threshold
101 */
102 public double getSeedThreshold() {
103 return _seed_threshold;
104 }
105
106 /**
107 * Return the neighbor threshold. Units are RMS noise.
108 *
109 * @return neighbor threshold
110 */
111 public double getNeighborThreshold() {
112 return _neighbor_threshold;
113 }
114
115 /**
116 * Return the cluster threshold. Units are RMS noise.
117 *
118 * @return cluster threshold
119 */
120 public double getClusterThreshold() {
121 return _cluster_threshold;
122 }
123
124 /**
125 * Return the name of the clustering algorithm.
126 *
127 * @return _name clusting algorithm name
128 */
129 public String getName() {
130 return _NAME;
131 }
132
133 /**
134 * Find clusters using the nearest neighbor algorithm.
135 *
136 * @param electrodes electrodes we are clustering
137 * @param readout_chip readout chip for these electrodes
138 * @param raw_hits List of RawTrackerHits to be clustered
139 * @return list of clusters, with a cluster being a list of RawTrackerHits
140 */
141 public List<List<RawTrackerHit>> findClusters(SiSensorElectrodes electrodes,
142 ReadoutChip readout_chip, List<RawTrackerHit> raw_hits) {
143
144 // Check that the seed threshold is at least as large as the neighbor threshold
145 if (_seed_threshold < _neighbor_threshold) {
146 throw new RuntimeException("Tracker hit clustering error: seed threshold below neighbor threshold");
147 }
148
149 // Create maps that show the channel status and relate the channel number to the raw hit and vice versa
150 int mapsize = 2 * raw_hits.size();
151 Map<Integer, Boolean> clusterable = new HashMap<Integer, Boolean>(mapsize);
152 Map<RawTrackerHit, Integer> hit_to_channel = new HashMap<RawTrackerHit, Integer>(mapsize);
153 Map<Integer, RawTrackerHit> channel_to_hit = new HashMap<Integer, RawTrackerHit>(mapsize);
154
155 // Create list of channel numbers to be used as cluster seeds
156 List<Integer> cluster_seeds = new ArrayList<Integer>();
157
158 // Loop over the raw hits and construct the maps used to relate cells and hits, initialize the
159 // clustering status map, and create a list of possible cluster seeds
160 for (RawTrackerHit raw_hit : raw_hits) {
161
162 // get the channel number for this hit
163 SiTrackerIdentifierHelper sid_helper = (SiTrackerIdentifierHelper) raw_hit.getIdentifierHelper();
164 IIdentifier id = raw_hit.getIdentifier();
165 int channel_number = sid_helper.getElectrodeValue(id);
166
167 // Check for duplicate RawTrackerHit
168 if (hit_to_channel.containsKey(raw_hit)) {
169 throw new RuntimeException("Duplicate hit: "+id.toString());
170 }
171
172 // Check for duplicate RawTrackerHits or channel numbers
173 if (channel_to_hit.containsKey(channel_number)) {
174 throw new RuntimeException("Duplicate channel number: "+channel_number);
175 }
176
177 // Add this hit to the maps that relate channels and hits
178 hit_to_channel.put(raw_hit, channel_number);
179 channel_to_hit.put(channel_number, raw_hit);
180
181 // Get the signal from the readout chip
182 double signal = readout_chip.decodeCharge(raw_hit);
183 double noiseRMS = readout_chip.getChannel(channel_number).computeNoise(electrodes.getCapacitance(channel_number));
184
185 // Mark this hit as available for clustering if it is above the neighbor threshold
186 clusterable.put(channel_number, signal/noiseRMS >= _neighbor_threshold);
187
188 // Add this hit to the list of seeds if it is above the seed threshold
189 if (signal/noiseRMS >= _seed_threshold) {
190 cluster_seeds.add(channel_number);
191 }
192 }
193
194 // Create a list of clusters
195 List<List<RawTrackerHit>> cluster_list = new ArrayList<List<RawTrackerHit>>();
196
197 // Now loop over the cluster seeds to form clusters
198 for (int seed_channel : cluster_seeds) {
199
200 // First check if this hit is still available for clustering
201 if (!clusterable.get(seed_channel)) continue;
202
203 // Create a new cluster
204 List<RawTrackerHit> cluster = new ArrayList<RawTrackerHit>();
205 double cluster_signal = 0.;
206 double cluster_noise_squared = 0.;
207
208 // Create a queue to hold channels whose neighbors need to be checked for inclusion
209 LinkedList<Integer> unchecked = new LinkedList<Integer>();
210
211 // Add the seed channel to the unchecked list and mark it as unavailable for clustering
212 unchecked.addLast(seed_channel);
213 clusterable.put(seed_channel, false);
214
215 // Check the neighbors of channels added to the cluster
216 while (unchecked.size() > 0) {
217
218 // Pull the next channel off the queue and add it's hit to the cluster
219 int clustered_cell = unchecked.removeFirst();
220 cluster.add(channel_to_hit.get(clustered_cell));
221 cluster_signal += readout_chip.decodeCharge(channel_to_hit.get(clustered_cell));
222 cluster_noise_squared += Math.pow(readout_chip.getChannel(clustered_cell).computeNoise(electrodes.getCapacitance(clustered_cell)),2);
223
224 // Get the neigbor channels
225 Set<Integer> neighbor_channels = electrodes.getNearestNeighborCells(clustered_cell);
226
227 // Now loop over the neighbors and see if we can add them to the cluster
228 for (int channel : neighbor_channels) {
229
230 // Get the status of this channel
231 Boolean addhit = clusterable.get(channel);
232
233 // If the map entry is null, there is no raw hit for this channel
234 if (addhit == null) continue;
235
236 // Check if this neighbor channel is still available for clustering
237 if (!addhit) continue;
238
239 // Add channel to the list of unchecked clustered channels
240 // and mark it unavailable for clustering
241 unchecked.addLast(channel);
242 clusterable.put(channel, false);
243
244 } // end of loop over neighbor cells
245 } // end of loop over unchecked cells
246
247 // Finished with this cluster, check cluster threshold and add it to the list of clusters
248 if (cluster.size() > 0 &&
249 cluster_signal/Math.sqrt(cluster_noise_squared) > _cluster_threshold)
250 {
251 cluster_list.add(cluster);
252 }
253
254 } // End of loop over seeds
255
256 // Finished finding clusters
257 return cluster_list;
258 }
259 }