diff -N HPSEcalCTPClusterer.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ HPSEcalCTPClusterer.java	15 Mar 2012 23:37:28 -0000	1.1
@@ -0,0 +1,237 @@

+package org.lcsim.hps.recon.ecal;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.geometry.subdetector.HPSEcal3.NeighborMap;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.util.Driver;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ * Creates clusters from CalorimeterHits in the HPSEcal detector.
+ *
+ * The clustering algorithm is from JLab Hall B 6 GeV DVCS Trigger Design doc.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Tim Nelson <[log in to unmask]>
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: HPSEcalCTPClusterer.java,v 1.1 2012/03/15 23:37:28 meeg Exp $
+ */
+public class HPSEcalCTPClusterer extends Driver {
+
+	HPSEcal3 ecal;
+	IDDecoder dec;
+	String ecalName;
+	String ecalCollectionName;
+	String clusterCollectionName = "EcalClusters";
+	Set<Long> clusterCenters = null;
+	Map<Long, Double> hitSums = null;
+	Map<Long, CalorimeterHit> hitMap = null;
+	// Map of crystals to their neighbors.
+	NeighborMap neighborMap = null;
+
+	public HPSEcalCTPClusterer() {
+	}
+
+	public void setClusterCollectionName(String clusterCollectionName) {
+		this.clusterCollectionName = clusterCollectionName;
+	}
+
+	public void setEcalCollectionName(String ecalCollectionName) {
+		this.ecalCollectionName = ecalCollectionName;
+	}
+
+	public void setEcalName(String ecalName) {
+		this.ecalName = ecalName;
+	}
+
+	public void startOfData() {
+		if (ecalCollectionName == null) {
+			throw new RuntimeException("The parameter ecalCollectionName was not set!");
+		}
+
+		if (ecalName == null) {
+			throw new RuntimeException("The parameter ecalName was not set!");
+		}
+	}
+
+	public void detectorChanged(Detector detector) {
+		// Get the Subdetector.
+		ecal = (HPSEcal3) detector.getSubdetector(ecalName);
+
+		// Get the decoder for the ECal IDs.
+		dec = ecal.getIDDecoder();
+
+		// Cache ref to neighbor map.
+		neighborMap = ecal.getNeighborMap();
+
+		clusterCenters = new HashSet<Long>();
+		//Make set of valid cluster centers.
+		for (Long cellID : neighborMap.keySet()) {
+			boolean isValidCenter = true;
+			Set<Long> neighbors = neighborMap.get(cellID);
+			for (Long neighborID : neighbors) {
+				Set<Long> neighborneighbors = new HashSet<Long>();
+				neighborneighbors.addAll(neighborMap.get(neighborID));
+				neighborneighbors.add(neighborID);
+
+				if (neighborneighbors.containsAll(neighbors)) {
+					isValidCenter = false;
+					break;
+				}
+			}
+			if (isValidCenter) {
+				clusterCenters.add(cellID);
+			}
+		}
+
+		//System.out.println(ecal.getName());
+		//System.out.println("  nx="+ecal.nx());
+		//System.out.println("  ny="+ecal.ny());
+		//System.out.println("  beamgap="+ecal.beamGap());
+		//System.out.println("  dface="+ecal.distanceToFace());
+
+		//System.out.println(neighborMap.toString());
+	}
+
+	public void process(EventHeader event) {
+		//System.out.println(this.getClass().getCanonicalName() + " - process");
+
+		// Get the list of raw ECal hits.
+		List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+		if (hits == null) {
+			throw new RuntimeException("Event is missing ECal raw hits collection!");
+		}
+
+		// Get the list of raw ECal hits.
+
+		sumHits(hits);
+
+		// Put Cluster collection into event.
+		int flag = 1 << LCIOConstants.CLBIT_HITS;
+		event.put(clusterCollectionName, createClusters(), Cluster.class, flag);
+	}
+
+	public void sumHits(List<CalorimeterHit> hits) {
+		// Hit map.
+		hitMap = new HashMap<Long, CalorimeterHit>();
+
+		hitSums = new HashMap<Long, Double>();
+		// Loop over ECal hits to compute energy sums
+		for (CalorimeterHit hit : hits) {
+			// Make a hit map for quick lookup by ID.
+			hitMap.put(hit.getCellID(), hit);
+
+			// Get neighbor crystal IDs. 
+			Set<Long> neighbors = neighborMap.get(hit.getCellID());
+
+			if (neighbors == null) {
+				throw new RuntimeException("Oops!  Set of neighbors is null!");
+			}
+
+			Double hitSum;
+
+			if (clusterCenters.contains(hit.getCellID())) {
+				hitSum = hitSums.get(hit.getCellID());
+				if (hitSum == null) {
+					hitSums.put(hit.getCellID(), hit.getRawEnergy());
+				} else {
+					hitSums.put(hit.getCellID(), hitSum + hit.getRawEnergy());
+				}
+			}
+
+			// Loop over neighbors to make hit list for cluster.
+			for (Long neighborId : neighbors) {
+				if (!clusterCenters.contains(neighborId)) {
+					continue;
+				}
+				hitSum = hitSums.get(neighborId);
+				if (hitSum == null) {
+					hitSums.put(neighborId, hit.getRawEnergy());
+				} else {
+					hitSums.put(neighborId, hitSum + hit.getRawEnergy());
+				}
+			}
+		}
+	}
+
+	public List<Cluster> createClusters() {
+//		boolean printClusters;
+		// New Cluster list to be added to event.
+		List<Cluster> clusters = new ArrayList<Cluster>();
+		//System.out.println("New event");
+		//for each crystal with a nonzero hit count, test for cluster
+		for (Long possibleCluster : hitSums.keySet()) {
+			Double thisSum = hitSums.get(possibleCluster);
+
+			// Get neighbor crystal IDs.
+			Set<Long> neighbors = neighborMap.get(possibleCluster);
+
+			if (neighbors == null) {
+				throw new RuntimeException("Oops!  Set of neighbors is null!");
+			}
+
+			//Apply peak detector scheme.
+			// Set the ID.
+//			dec.setID(possibleCluster);
+//			int x1 = dec.getValue("ix");
+//			int y1 = dec.getValue("iy");
+//			System.out.printf("\nThis cluster: E= %f, ID=%d, x=%d, y=%d, neighbors=%d\n", thisSum, possibleCluster, x1, y1, neighbors.size());
+			boolean isCluster = true;
+			for (Long neighborId : neighbors) {
+				// Set the ID.
+//				dec.setID(neighborId);
+//				int x2 = dec.getValue("ix");
+//				int y2 = dec.getValue("iy");
+
+				Double neighborSum = hitSums.get(neighborId);
+				if (neighborSum == null) {
+					continue;
+				}
+
+//				System.out.printf("Neighbor cluster: E= %f, ID=%d, x=%d, y=%d, neighbors=%d\n", neighborSum, neighborId, x2, y2, neighborMap.get(neighborId).size());
+
+				if (neighborSum > thisSum) {
+//					System.out.println("Reject cluster: sum cut");
+					isCluster = false;
+					break;
+				} else if (neighborSum.equals(thisSum) && neighborId > possibleCluster) {
+//					System.out.println("Reject cluster: tie-breaker cut");
+					isCluster = false;
+					break;
+				}
+			}
+
+			if (isCluster) {
+//				System.out.println("Accept this cluster");
+				HPSEcalCluster cluster = new HPSEcalCluster(possibleCluster);
+				CalorimeterHit hit = hitMap.get(possibleCluster);
+				if (hit != null /*&& hit.getRawEnergy() > hitEMin && hit.getRawEnergy() < hitEMax*/) {
+//					System.out.printf("Adding hit, E = %f\n", hit.getRawEnergy());
+					cluster.addHit(hit);
+				}
+
+				for (Long neighborId : neighbors) {
+					// Find the neighbor hit in the event if it exists.
+					hit = hitMap.get(neighborId);
+					if (hit != null /*&& hit.getRawEnergy() > hitEMin && hit.getRawEnergy() < hitEMax*/) {
+//						System.out.printf("Adding hit, E = %f\n", hit.getRawEnergy());
+						cluster.addHit(hit);
+					}
+				}
+				clusters.add(cluster);
+			}
+		}
+		return clusters;
+	}
+}