hps-java/src/main/java/org/lcsim/hps/recon/ecal
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;
+ }
+}