lcsim-contrib/src/main/java/org/lcsim/contrib/jeremym/hps
diff -N HPSEcalClusterer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HPSEcalClusterer.java 28 Apr 2011 00:34:30 -0000 1.1
@@ -0,0 +1,206 @@
+package org.lcsim.contrib.jeremym.hps;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+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.HPSEcal;
+import org.lcsim.geometry.subdetector.HPSEcal.XYSide;
+import org.lcsim.geometry.util.IDEncoder;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.util.Driver;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ * Creates clusters from CalorimeterHits in the HPSEcal detector.
+ *
+ * Clustering algorithm is from pages 83 and 84 of the HPS Proposal.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @version $Id: HPSEcalClusterer.java,v 1.1 2011/04/28 00:34:30 jeremy Exp $
+ */
+public class HPSEcalClusterer extends Driver
+{
+ HPSEcal ecal;
+
+ String ecalCollectionName;
+ String ecalName;
+ String clusterCollectionName = "EcalClusters";
+
+ // Minimum E for cluster seed.
+ double seedEMin = .05;
+
+ // Minimum E to add hit to cluster.
+ double addEMin = .03;
+
+ // Odd or even number of crystals in X.
+ boolean oddX;
+
+ public HPSEcalClusterer()
+ {}
+
+ public void setClusterCollectionName(String clusterCollectionName)
+ {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ public void setSeedEMin(double seedEMin)
+ {
+ this.seedEMin = seedEMin;
+ }
+
+ public void setAddEMin(double addEMin)
+ {
+ this.addEMin = addEMin;
+ }
+
+ 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 (ecalCollectionName == null)
+ throw new RuntimeException("The parameter ecalName was not set!");
+
+ System.out.println(this.getClass().getCanonicalName());
+ System.out.println(" seedEMin="+seedEMin);
+ System.out.println(" addEMin="+addEMin);
+ System.out.println();
+ }
+
+ public void detectorChanged(Detector detector)
+ {
+ ecal = (HPSEcal)detector.getSubdetector(ecalName);
+
+ 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());
+ }
+
+ public void process(EventHeader event)
+ {
+ //System.out.println("HPSEcalClusterer.process");
+
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ if (hits == null)
+ throw new RuntimeException("Event is missing ECal hits collection!");
+
+ IDDecoder dec = ecal.getIDDecoder();
+
+ // Hit map.
+ Map<Long,CalorimeterHit> hitMap = new HashMap<Long,CalorimeterHit>();
+
+ // Make map of x, y, and side to hit.
+ for (CalorimeterHit hit : hits)
+ {
+ dec.setID(hit.getCellID());
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ // Setup an encoder from the decoder.
+ IDEncoder enc = new IDEncoder(dec.getIDDescription());
+
+ // Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ // Loop over ECal hits.
+ for (CalorimeterHit hit : hits)
+ {
+ // Cut on min seed E.
+ if (hit.getRawEnergy() < seedEMin)
+ continue;
+
+ // Set decoder's ID.
+ dec.setID(hit.getCellID());
+
+ // Get ID values.
+ int ix = dec.getValue("ix");
+ int iy = dec.getValue("iy");
+ int side = dec.getValue("side");
+
+ // Get neighboring crystals to hit.
+ Set<XYSide> neighbors = ecal.getNeighbors(ix, iy, side);
+
+ //System.out.println("neighbors of " + ix + ", " + iy + ", " + side);
+
+ int[] buffer = new int[dec.getFieldCount()];
+ dec.getValues(buffer);
+
+ // Loop over neighbors.
+ List<CalorimeterHit> neighborHits = new ArrayList<CalorimeterHit>();
+ for (XYSide neighbor : neighbors)
+ {
+ // debug print
+ //System.out.println(" " + neighbor.x() + ", " + neighbor.y() + ", " + neighbor.side());
+
+ buffer[dec.getFieldIndex("ix")] = neighbor.x();
+ buffer[dec.getFieldIndex("iy")] = neighbor.y();
+ buffer[dec.getFieldIndex("side")] = neighbor.side();
+
+ // Create ID for this crystal.
+ // FIXME This is redundant and inefficient.
+ long neighborId = enc.setValues(buffer);
+
+ // Find neighbor hit if it exists.
+ CalorimeterHit neighborHit = hitMap.get(neighborId);
+
+ // Got a hit?
+ if (neighborHit != null)
+ {
+ // debug print
+ /*
+ dec.setID(neighborHit.getCellID());
+ System.out.println(" found neighbor hit w/ x, y, side ( "
+ + dec.getValue("ix") + ", "
+ + dec.getValue("iy") + ", "
+ + dec.getValue("side") + ")");
+ */
+
+ // Neighbor is hotter so skip this hit.
+ if (neighborHit.getRawEnergy() > hit.getRawEnergy())
+ {
+ //System.out.println("skipping hit because neighbor has more E");
+ continue;
+ }
+
+ // Add to cluster if above min E.
+ if (neighborHit.getRawEnergy() >= addEMin)
+ neighborHits.add(neighborHit);
+ }
+ }
+
+ // Make cluster.
+ BasicCluster cluster = new BasicCluster();
+ cluster.addHit(hit);
+ for (CalorimeterHit clusHit : neighborHits)
+ {
+ cluster.addHit(clusHit);
+ }
+ clusters.add(cluster);
+ //System.out.println("cluster nhits = " + cluster.getCalorimeterHits().size());
+ //System.out.println("cluster E = " + cluster.getEnergy());
+ }
+
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, clusters, Cluster.class, flag);
+ }
+}
\ No newline at end of file