Author: [log in to unmask]
Date: Mon Dec 15 15:54:53 2014
New Revision: 1738
Log:
Add new cluster package to ecal recon.
Added:
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/AbstractClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/CTPEcalClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterDriver.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterUtilities.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/Clusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClustererFactory.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/EcalClusterIC.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPEcalClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/LegacyClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleInnerCalClusterer.java
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/AbstractClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/AbstractClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/AbstractClusterer.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,21 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.List;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.geometry.subdetector.HPSEcal3.NeighborMap;
+
+public abstract class AbstractClusterer implements Clusterer {
+
+ HPSEcal3 ecal;
+ NeighborMap neighborMap;
+
+ public void setEcalSubdetector(HPSEcal3 ecal) {
+ this.ecal = ecal;
+ this.neighborMap = ecal.getNeighborMap();
+ }
+
+ public abstract List<Cluster> createClusters(List<CalorimeterHit> hits);
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/CTPEcalClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/CTPEcalClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/CTPEcalClusterer.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,383 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+import java.util.Collection;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.PriorityQueue;
+import java.util.Set;
+
+import org.lcsim.event.CalorimeterHit;
+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.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: CTPEcalClusterer.java,v 1.1 2013/02/25 22:39:24 meeg Exp $
+ */
+public class CTPEcalClusterer extends Driver {
+
+ Detector detector = null;
+ // The calorimeter object.
+ HPSEcal3 ecal;
+ IDDecoder dec;
+ // The identification name for getting the calorimeter object.
+ String ecalName;
+ // The collection name for the calorimeter hits.
+ String ecalCollectionName;
+ // The collection name in which to store clusters.
+ 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;
+ // The time period in which clusters may be formed. A negative value means that all hits
+ // will always be used in cluster finding, regardless of the time difference between them.
+ double clusterWindow = -1;
+ // The minimum energy needed for a hit to be considered.
+ double addEMin = 0;
+
+ public CTPEcalClusterer() {
+ }
+
+ public void setAddEMin(double addEMin) {
+ this.addEMin = addEMin;
+ }
+
+ public void setClusterWindow(double clusterWindow) {
+ this.clusterWindow = clusterWindow;
+ }
+
+ public void setClusterCollectionName(String clusterCollectionName) {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ public void setEcalCollectionName(String ecalCollectionName) {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setEcalName(String ecalName) {
+ this.ecalName = ecalName;
+ }
+
+ @Override
+ public void startOfData() {
+ // Make sure that there is a cluster collection name into which clusters may be placed.
+ if (ecalCollectionName == null) {
+ throw new RuntimeException("The parameter ecalCollectionName was not set!");
+ }
+
+ // Make sure that there is a calorimeter detector.
+ if (ecalName == null) {
+ throw new RuntimeException("The parameter ecalName was not set!");
+ }
+ }
+
+ @Override
+ public void detectorChanged(Detector detector) {
+
+ this.detector = detector;
+
+ // Get the Subdetector.
+ ecal = (HPSEcal3) detector.getSubdetector(ecalName);
+
+ // Get the decoder for the ECal IDs.
+ dec = ecal.getIDDecoder();
+
+ // Cache reference to the neighbor map.
+ neighborMap = ecal.getNeighborMap();
+
+ clusterCenters = new HashSet<Long>();
+ // Make set of valid cluster centers.
+ // Exclude edge crystals as good 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);
+ }
+ }
+ }
+
+ @Override
+ public void process(EventHeader event) {
+ // Make sure that this event has calorimeter hits.
+ if (event.hasCollection(CalorimeterHit.class, ecalCollectionName)) {
+ // Get the list of calorimeter hits from the event.
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+
+ // Define a list of clusters to be filled.
+ List<HPSEcalCluster> clusters;
+
+ // If there is a cluster window, run the cluster window code. A cluster window is a time
+ // period in nanoseconds during which hits can be applied to the same cluster.
+ if (clusterWindow >= 0) {
+ // Create priority queues. These are sorted by the time variable associated with each hit.
+ PriorityQueue<CalorimeterHit> futureHits = new PriorityQueue<CalorimeterHit>(10, new TimeComparator());
+ PriorityQueue<CalorimeterHit> pastHits = new PriorityQueue<CalorimeterHit>(10, new TimeComparator());
+
+ // Initialize the cluster list variable.
+ clusters = new ArrayList<HPSEcalCluster>();
+
+ // Populate the list of unprocessed hits with the calorimeter hits. These will then be sorted
+ // by time, from first to last, automatically by the priority queue.
+ for (CalorimeterHit hit : hits) {
+ if (hit.getRawEnergy() > addEMin) {
+ futureHits.add(hit);
+ }
+ }
+
+ // We process the unprocessed hits...
+ while (!futureHits.isEmpty()) {
+ // Move the first occurring hit from the unprocessed list to the processed list.
+ CalorimeterHit nextHit = futureHits.poll();
+ pastHits.add(nextHit);
+
+ // Add any hits that occurred at the same time as the hit we just added to the processed list.
+ while (!futureHits.isEmpty() && futureHits.peek().getTime() == nextHit.getTime()) {
+ pastHits.add(futureHits.poll());
+ }
+
+ // Remove hits that happened earlier than the cluster window period.
+ while (pastHits.peek().getTime() < nextHit.getTime() - clusterWindow) {
+ pastHits.poll();
+ }
+
+ // Calculate the cluster energy for each crystal. This should be the
+ // total energy for the 3x3 crystal collection sorrounding the center
+ // crystal.
+ sumHits(pastHits);
+
+ // Choose which crystal is the appropriate cluster crystal.
+ clusters.addAll(createClusters());
+ }
+ // If there is no cluster window, then all the hits in the event are visible simultaneously.
+ } else {
+ // Calculate the cluster energy of each crystal.
+ sumHits(hits);
+ // Generate the clusters.
+ clusters = createClusters();
+ }
+
+ // Put the cluster collection into the event.
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, clusters, HPSEcalCluster.class, flag);
+ }
+ }
+
+ public void sumHits(Collection<CalorimeterHit> hits) {
+ // Store the latest hit on each crystal in a map for later reference in
+ // the clustering algorithm.
+ hitMap = new HashMap<Long, CalorimeterHit>();
+ // Store the cluster energy for each crystal. Cluster energy represents
+ // the total energy of the 3x3 crystal set.
+ hitSums = new HashMap<Long, Double>();
+
+ // Loop over the active calorimeter hits to compute the cluster energies.
+ for (CalorimeterHit hit : hits) {
+ // Make a hit map for quick lookup by ID.
+ hitMap.put(hit.getCellID(), hit);
+
+ // Get the cell ID for the current crystal's neighbors.
+ Set<Long> neighbors = neighborMap.get(hit.getCellID());
+
+ // If there are no neighbors, something is rather wrong.
+ if (neighbors == null) {
+ throw new RuntimeException("Oops! Set of neighbors is null!");
+ }
+
+ // Store the energy of the current calorimeter hit.
+ Double hitSum;
+
+ // We are only interested in this crystal's cluster energy if it is
+ // a valid cluster crystal. Edge crystals are not allowed to be clusters,
+ // so these are ignored.
+ if (clusterCenters.contains(hit.getCellID())) {
+ // Check if an energy has been assigned to the crystal.
+ hitSum = hitSums.get(hit.getCellID());
+
+ // If not, then the crystal's cluster energy is equal to this hit's energy.
+ if (hitSum == null) {
+ hitSums.put(hit.getCellID(), hit.getRawEnergy());
+ }
+ // Otherwise, add the energy of this hit to the total crystal cluster energy.
+ else {
+ hitSums.put(hit.getCellID(), hitSum + hit.getRawEnergy());
+ }
+ }
+
+ // Loop over neighbors to add the current hit's energy to the neighbor's
+ // cluster energy.
+ for (Long neighborId : neighbors) {
+ // If the crystal is not an edge crystal, ignore its hit energy.
+ if (!clusterCenters.contains(neighborId)) {
+ continue;
+ }
+
+ // Get the cluster energy of the neighboring crystals.
+ hitSum = hitSums.get(neighborId);
+
+ // If the neighbor crystal has no cluster energy, then set the
+ // cluster energy to the current hit's energy.
+ if (hitSum == null) {
+ hitSums.put(neighborId, hit.getRawEnergy());
+ }
+ // Otherwise, add the hit's energy to the neighbor's cluster energy.
+ else {
+ hitSums.put(neighborId, hitSum + hit.getRawEnergy());
+ }
+ }
+ }
+ }
+
+ public List<HPSEcalCluster> createClusters() {
+ // Create a list of clusters to be added to the event,
+ List<HPSEcalCluster> clusters = new ArrayList<HPSEcalCluster>();
+
+ // We examine each crystal with a non-zero cluster energy.
+ for(Long possibleCluster : hitSums.keySet()) {
+ // Get the luster energy for the crystal this hit is assocaite with.
+ Double thisSum = hitSums.get(possibleCluster);
+
+ // Get neighboring crystals' IDs.
+ Set<Long> neighbors = neighborMap.get(possibleCluster);
+
+ // If there are no neighbors, throw an error.
+ if (neighbors == null) {
+ throw new RuntimeException("Oops! Set of neighbors is null!");
+ }
+
+ // Get the x/y position of the hit's associated crystal.
+ dec.setID(possibleCluster);
+ int x1 = dec.getValue("ix");
+ int y1 = dec.getValue("iy");
+
+ // Store whether it is a valid cluster or not.
+ boolean isCluster = true;
+
+ // Check to see if any of the crystal's neighbors preclude the current crystal
+ // from being a proper cluster. The cluster crystal should have the highest
+ // energy among its neighbors
+ for (Long neighborId : neighbors) {
+ // Get the x/y position of the neighbor's associated crystal.
+ dec.setID(neighborId);
+ int x2 = dec.getValue("ix");
+ int y2 = dec.getValue("iy");
+
+ // If the neighbor's energy value does not exist, we don't need to perform
+ // any additional checks for this neighbor. A crystal with no energy can
+ // not be the center of a cluster.
+ Double neighborSum = hitSums.get(neighborId);
+ if (neighborSum == null) {
+ continue;
+ }
+
+ // If the neighbor's energy value is greater than this crystal's value,
+ // then this crystal is not the cluster and we may terminate the check.
+ if (neighborSum > thisSum) {
+ isCluster = false;
+ break;
+ }
+ // If the crystals each have the same energy, we choose the crystal that
+ // is closest to the electron side of the detector. If both crystals are
+ // equally close, we choose the crystal closest to the beam gap. If the
+ // neighbor fits these parameters better, this is not a crystal and we
+ // may skip any further checks.
+ else if (neighborSum.equals(thisSum) && (x1 > x2 || (x1 == x2 && Math.abs(y1) < Math.abs(y2)))) {
+ isCluster = false;
+ break;
+ }
+ }
+
+ // If the crystal was not invalidated by the any of the neighboring crystals,
+ // then it is a cluster crystal and should be processed.
+ if (isCluster) {
+ // Make a list to store the hits that are part of this cluster.
+ List<CalorimeterHit> hits = new ArrayList<CalorimeterHit>();
+
+ // Store the time at which the cluster occurred.
+ double clusterTime = Double.NEGATIVE_INFINITY;
+
+ // Get the last hit on this crystal.
+ CalorimeterHit hit = hitMap.get(possibleCluster);
+
+ // If the hit exists, add it to the list of associated hits.
+ if (hit != null) {
+ hits.add(hit);
+
+ // If the latest hit's time is later than the current cluster time,
+ // set the cluster time to the latest hit's time.
+ if (hit.getTime() > clusterTime) {
+ clusterTime = hit.getTime();
+ }
+ }
+
+ // Add all of the neighboring crystals to the cluster, if they have a
+ // hit associated with them. Crystals with no hits are not actually part
+ // of a cluster.
+ for (Long neighborId : neighbors) {
+ hit = hitMap.get(neighborId);
+ if (hit != null) {
+ hits.add(hit);
+ if (hit.getTime() > clusterTime) {
+ clusterTime = hit.getTime();
+ }
+ }
+ }
+
+ // Generate a new cluster seed hit from the above results.
+ HPSCalorimeterHit seedHit = new HPSCalorimeterHit(0.0, clusterTime, possibleCluster, hits.get(0).getType());
+ seedHit.setMetaData(hits.get(0).getMetaData());
+
+ // Generate a new cluster from the seed hit.
+ HPSEcalCluster cluster = new HPSEcalCluster();
+ cluster.setSeedHit(seedHit);
+ // Populate the cluster with each of the chosen neighbors.
+ for (CalorimeterHit clusterHit : hits) {
+ cluster.addHit(clusterHit);
+ }
+ // Add the cluster to the cluster list.
+ clusters.add(cluster);
+ }
+ }
+
+ // Return the list of clusters.
+ return clusters;
+ }
+
+ static class TimeComparator implements Comparator<CalorimeterHit> {
+ // Compare by time with the earlier coming before the later.
+ public int compare(CalorimeterHit o1, CalorimeterHit o2) {
+ if (o1.getTime() == o2.getTime()) {
+ return 0;
+ } else {
+ return (o1.getTime() > o2.getTime()) ? 1 : -1;
+ }
+ }
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterDriver.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterDriver.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterDriver.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,128 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.List;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.Subdetector;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.lcio.LCIOConstants;
+import org.lcsim.util.Driver;
+
+/**
+ * This is a basic Driver that creates Cluster collections through the
+ * Clusterer interface.
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public class ClusterDriver extends Driver {
+
+ protected String ecalName = "Ecal";
+ protected HPSEcal3 ecal;
+ protected String outputClusterCollectionName = "EcalClusters";
+ protected String inputHitCollectionName = "EcalCalHits";
+ protected Clusterer clusterer;
+ protected boolean createEmptyClusterCollection = true;
+ protected boolean raiseErrorNoHitCollection = false;
+ protected boolean skipNoClusterEvents = false;
+ protected boolean writeClusterCollection = true;
+ protected boolean storeHits = true;
+
+ protected ClusterDriver() {
+ }
+
+ public void setEcalName(String ecalName) {
+ this.ecalName = ecalName;
+ }
+
+ public void setOutputClusterCollectionName(String outputClusterCollectionName) {
+ this.outputClusterCollectionName = outputClusterCollectionName;
+ }
+
+ public void setInputHitCollectionName(String inputHitCollectionName) {
+ this.inputHitCollectionName = inputHitCollectionName;
+ }
+
+ public void setSkipNoClusterEvents(boolean skipNoClusterEvents) {
+ this.skipNoClusterEvents = skipNoClusterEvents;
+ }
+
+ public void setWriteClusterCollection(boolean writeClusterCollection) {
+ this.writeClusterCollection = writeClusterCollection;
+ }
+
+ public void setRaiseErrorNoHitCollection(boolean raiseErrorNoHitCollection) {
+ this.raiseErrorNoHitCollection = raiseErrorNoHitCollection;
+ }
+
+ public void setStoreHits(boolean storeHits) {
+ this.storeHits = storeHits;
+ }
+
+ public void setClusterer(String name) {
+ clusterer = ClustererFactory.create(name);
+ }
+
+ public void setClusterer(Clusterer clusterer) {
+ this.clusterer = clusterer;
+ }
+
+ public void setCreateEmptyClusterCollection(boolean createEmptyClusterCollection) {
+ this.createEmptyClusterCollection = createEmptyClusterCollection;
+ }
+
+ public void detectorChanged(Detector detector) {
+ Subdetector subdetector = detector.getSubdetector(ecalName);
+ if (subdetector == null) {
+ throw new RuntimeException("There is no subdetector called " + ecalName + " in the detector.");
+ }
+ if (!(subdetector instanceof HPSEcal3)) {
+ throw new RuntimeException("Ther subdetector " + ecalName + " does not have the right type.");
+ }
+ ecal = (HPSEcal3) subdetector;
+ }
+
+ public void startOfData() {
+ if (this.clusterer == null) {
+ throw new RuntimeException("The clusterer was never initialized.");
+ }
+ if (this.clusterer instanceof AbstractClusterer) {
+ ((AbstractClusterer)clusterer).setEcalSubdetector(ecal);
+ }
+ }
+
+ /**
+ * This method implements the default clustering procedure based on input parameters.
+ */
+ public void process(EventHeader event) {
+ if (event.hasCollection(CalorimeterHit.class, inputHitCollectionName)) {
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, inputHitCollectionName);
+ List<Cluster> clusters = clusterer.createClusters(hits);
+ if (clusters == null) {
+ throw new RuntimeException("The clusterer returned a null pointer.");
+ }
+ if (clusters.isEmpty() && this.skipNoClusterEvents) {
+ throw new NextEventException();
+ }
+ if (event.hasCollection(Cluster.class, this.outputClusterCollectionName)) {
+ throw new RuntimeException("There is already a cluster collection called " + this.outputClusterCollectionName);
+ }
+ int flags = 0;
+ if (this.storeHits) {
+ flags = 1 << LCIOConstants.CLBIT_HITS;
+ }
+ if (!clusters.isEmpty() || this.createEmptyClusterCollection) {
+ event.put(outputClusterCollectionName, clusters, Cluster.class, flags);
+ if (!this.writeClusterCollection) {
+ event.getMetaData(clusters).setTransient(true);
+ }
+ }
+ } else {
+ this.getLogger().warning("The input hit collection " + this.inputHitCollectionName + " is missing from the event.");
+ if (this.raiseErrorNoHitCollection) {
+ throw new RuntimeException("The expected hit collection " + this.inputHitCollectionName + " is missing from the event.");
+ }
+ }
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterUtilities.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterUtilities.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterUtilities.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,41 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.HashSet;
+import java.util.LinkedHashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+
+public final class ClusterUtilities {
+
+ private ClusterUtilities() {
+ }
+
+ public static Map<Long, CalorimeterHit> createHitMap(List<CalorimeterHit> hits) {
+ Map<Long, CalorimeterHit> hitMap = new LinkedHashMap<Long, CalorimeterHit>();
+ for (CalorimeterHit hit : hits) {
+ hitMap.put(hit.getCellID(), hit);
+ }
+ return hitMap;
+ }
+
+ /**
+ * Given a hit, find its list of neighboring crystals that have hits and return their IDs.
+ * @param hit The input hit.
+ * @param hitMap The hit map with all the collection's hits.
+ * @return The set of neighboring hit IDs.
+ */
+ public static Set<Long> findNeighborHitIDs(HPSEcal3 ecal, CalorimeterHit hit, Map<Long, CalorimeterHit> hitMap) {
+ Set<Long> neigbhors = ecal.getNeighborMap().get(hit.getCellID());
+ Set<Long> neighborHitIDs = new HashSet<Long>();
+ for (long neighborID : neigbhors) {
+ if (hitMap.containsKey(neighborID)) {
+ neighborHitIDs.add(neighborID);
+ }
+ }
+ return neighborHitIDs;
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/Clusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/Clusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/Clusterer.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,12 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.List;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+
+public interface Clusterer {
+
+ List<Cluster> createClusters(List<CalorimeterHit> hits);
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClustererFactory.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClustererFactory.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClustererFactory.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,18 @@
+package org.hps.recon.ecal.cluster;
+
+public final class ClustererFactory {
+
+ private ClustererFactory() {
+ }
+
+ public static Clusterer create(String name) {
+ if (LegacyClusterer.class.getSimpleName().equals(name)) {
+ return new LegacyClusterer();
+ } if (SimpleInnerCalClusterer.class.getSimpleName().equals(name)) {
+ return new SimpleInnerCalClusterer();
+ } else {
+ throw new IllegalArgumentException("Unknown clusterer: " + name);
+ }
+ }
+
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/EcalClusterIC.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/EcalClusterIC.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/EcalClusterIC.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,868 @@
+package org.hps.recon.ecal.cluster;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.awt.Point;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import org.lcsim.detector.IGeometryInfo;
+import org.lcsim.detector.solids.Trd;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.geometry.subdetector.HPSEcal3.NeighborMap;
+import org.lcsim.lcio.LCIOConstants;
+import org.lcsim.util.Driver;
+
+
+/**
+ * This Driver creates clusters from the CalorimeterHits of an
+ * {@link org.lcsim.geometry.subdetectur.HPSEcal3} detector.
+ *
+ * This clustering logic is based on that from the CLAS-Note-2005-001.
+ * This sorts hits from highest to lowest energy and build clusters around
+ * each local maximum/seed hit. Common hits are distributed between clusters
+ * when minimum between two clusters. There is a threshold cut for minimum
+ * hit energy, minimum cluster energy, and minimum seed hit energy. There is
+ * also a timing threshold with respect to the seed hit. All of these parameters
+ * are tunable and should be refined with more analysis.
+ *
+ *
+ * @author Holly Szumila-Vance <[log in to unmask]>
+ * @author Kyle McCarty <[log in to unmask]>
+ *
+ */
+public class EcalClusterIC extends Driver {
+ // File writer to output cluster results.
+ FileWriter writeHits = null;
+ // LCIO collection name for calorimeter hits.
+ String ecalCollectionName;
+ // Name of the calorimeter detector object.
+ String ecalName = "Ecal";
+ // LCIO cluster collection name to which to write.
+ String clusterCollectionName = "EcalClustersIC";
+ // Collection name for rejected hits
+ String rejectedHitName = "RejectedHits";
+ // File path to which to write event display output.
+ String outfile = null;
+ // Map of crystals to their neighbors.
+ NeighborMap neighborMap = null;
+ // Minimum energy threshold for hits; lower energy hits will be
+ // excluded from clustering. Units in GeV.
+ double hitEnergyThreshold = 0.0075;
+ // Minimum energy threshold for seed hits; if seed hit is below
+ // cluster is excluded from output. Units in GeV.
+ double seedEnergyThreshold = 0.1;
+ // Minimum energy threshold for cluster hits; if total cluster
+ // energy is below, the cluster is excluded. Units in GeV.
+ double clusterEnergyThreshold = 0.3;
+ // A Comparator that sorts CalorimeterHit objects from highest to
+ // lowest energy.
+ private static final EnergyComparator ENERGY_COMP = new EnergyComparator();
+ // Track the event number for the purpose of outputting to event
+ // display format.
+ private int eventNum = -1;
+ // Apply time cut to hits
+ boolean timeCut = false;
+ // Minimum time cut window range. Units in ns.
+ double minTime = 0.0;
+ // Maximum time cut window range. Units in ns.
+ double timeWindow = 20.0;
+ // Variables for electron energy corrections
+ static final double ELECTRON_ENERGY_A = -0.0027;
+ static final double ELECTRON_ENERGY_B = -0.06;
+ static final double ELECTRON_ENERGY_C = 0.95;
+ // Variables for positron energy corrections
+ static final double POSITRON_ENERGY_A = -0.0096;
+ static final double POSITRON_ENERGY_B = -0.042;
+ static final double POSITRON_ENERGY_C = 0.94;
+ // Variables for photon energy corrections
+ static final double PHOTON_ENERGY_A = 0.0015;
+ static final double PHOTON_ENERGY_B = -0.047;
+ static final double PHOTON_ENERGY_C = 0.94;
+ // Variables for electron position corrections
+ static final double ELECTRON_POS_A = 0.0066;
+ static final double ELECTRON_POS_B = -0.03;
+ static final double ELECTRON_POS_C = 0.028;
+ static final double ELECTRON_POS_D = -0.45;
+ static final double ELECTRON_POS_E = 0.465;
+ // Variables for positron position corrections
+ static final double POSITRON_POS_A = 0.0072;
+ static final double POSITRON_POS_B = -0.031;
+ static final double POSITRON_POS_C = 0.007;
+ static final double POSITRON_POS_D = 0.342;
+ static final double POSITRON_POS_E = 0.108;
+ // Variables for photon position corrections
+ static final double PHOTON_POS_A = 0.005;
+ static final double PHOTON_POS_B = -0.032;
+ static final double PHOTON_POS_C = 0.011;
+ static final double PHOTON_POS_D = -0.037;
+ static final double PHOTON_POS_E = 0.294;
+
+
+
+ public void setClusterCollectionName(String clusterCollectionName) {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ public void setEcalCollectionName(String ecalCollectionName) {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setEcalName(String ecalName) {
+ this.ecalName = ecalName;
+ }
+
+ /**
+ * Output file name for event display output.
+ * @param outfile
+ */
+ public void setOutfile(String outfile) {
+ this.outfile = outfile;
+ }
+
+ public void setRejectedHitName(String rejectedHitName){
+ this.rejectedHitName = rejectedHitName;
+ }
+
+ /**
+ * Minimum energy for a hit to be used in a cluster. Default of 0.0075 GeV
+ *
+ * @param hitEnergyThreshold
+ */
+ public void sethitEnergyThreshold(double hitEnergyThreshold) {
+ this.hitEnergyThreshold = hitEnergyThreshold;
+ }
+
+ /**
+ * Minimum energy for a seed hit. Default of 0.1 GeV
+ *
+ * @param seedEnergyThreshold
+ */
+ public void setseedEnergyThreshold(double seedEnergyThreshold) {
+ this.seedEnergyThreshold = seedEnergyThreshold;
+ }
+
+ /**
+ * Minimum energy for a cluster. Default of 0.3 GeV
+ *
+ * @param clusterEnergyThreshold
+ */
+ public void setclusterEnergyThreshold(double clusterEnergyThreshold) {
+ this.clusterEnergyThreshold = clusterEnergyThreshold;
+ }
+
+ /**
+ * Apply time cuts to hits. Defaults to false.
+ *
+ * @param timeCut
+ */
+ public void setTimeCut(boolean timeCut) {
+ this.timeCut = timeCut;
+ }
+
+ /**
+ * Minimum hit time, if timeCut is true. Default of 0 ns.
+ *
+ * @param minTime
+ */
+ public void setMinTime(double minTime) {
+ this.minTime = minTime;
+ }
+
+ /**
+ * Width of time window, if timeCut is true. Default of 20 ns.
+ *
+ * @param timeWindow
+ */
+ public void setTimeWindow(double timeWindow) {
+ this.timeWindow = timeWindow;
+ }
+
+ // Make a map for quick calculation of the x-y position of crystal face
+ public Map<Point, double[]> correctedPositionMap = new HashMap<Point, double[]>();
+
+ public void startOfData() {
+ // Make sure that the calorimeter hit collection name is defined.
+ if (ecalCollectionName == null) {
+ throw new RuntimeException("The parameter ecalCollectionName was not set!");
+ }
+
+ // Make sure the name of calorimeter detector is defined.
+ if (ecalName == null) {
+ throw new RuntimeException("The parameter ecalName was not set!");
+ }
+
+ if (outfile!=null) {
+ // Create a file writer and clear the output file, if it exists.
+ try {
+ writeHits = new FileWriter(outfile);
+ writeHits.write("");
+ } catch (IOException e) {
+ }
+ }
+ }
+
+ public void detectorChanged(Detector detector) {
+ // Get the calorimeter.
+ HPSEcal3 ecal = (HPSEcal3) detector.getSubdetector(ecalName);
+
+ // Store the map of neighbor crystals for the current calorimeter set-up.
+ neighborMap = ecal.getNeighborMap();
+ }
+
+ public void process(EventHeader event) {
+ // Make sure the current event contains calorimeter hits.
+ if (event.hasCollection(CalorimeterHit.class, ecalCollectionName)) {
+
+ // Generate clusters from the calorimeter hits.
+ //List<HPSEcalClusterIC> clusterList = null;
+ try { createClusters(event); }
+ catch(IOException e) { }
+ }
+ }
+
+ public void createClusters(EventHeader event) throws IOException {
+
+ // Create a list to store the event hits in.
+ List<CalorimeterHit> hitList = new ArrayList<CalorimeterHit>();
+ List<CalorimeterHit> baseList = event.get(CalorimeterHit.class, ecalCollectionName);
+ for(CalorimeterHit r : baseList) {
+ hitList.add(r);
+ }
+
+ //Create a list to store the newly created clusters in.
+ ArrayList<HPSEcalClusterIC> clusterList = new ArrayList<HPSEcalClusterIC>();
+
+ //Create a list to store the rejected hits in.
+ ArrayList<CalorimeterHit> rejectedHitList = new ArrayList<CalorimeterHit>();
+
+ // Sort the list of hits by energy.
+ Collections.sort(hitList, ENERGY_COMP);
+
+ // Filter the hit list of any hits that fail to pass the
+ // designated threshold.
+ filterLoop:
+ for(int index = hitList.size() - 1; index >= 0; index--) {
+ // If the hit is below threshold or outside of time window, kill it.
+ if((hitList.get(index).getCorrectedEnergy() < hitEnergyThreshold)||
+ (timeCut && (hitList.get(index).getTime() < minTime || hitList.get(index).getTime() > (minTime + timeWindow)))) {
+ rejectedHitList.add(hitList.get(index));
+ hitList.remove(index);
+ }
+
+ // Since the hits are sorted by energy from highest to
+ // lowest, any hit that is above threshold means that all
+ // subsequent hits will also be above threshold. Continue through
+ // list to check in time window.
+ else { continue; }
+ }
+
+ //Create a map to connect the cell ID of a calorimeter crystal to the hit which occurred in that crystal.
+ HashMap<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
+ for (CalorimeterHit hit : hitList) { hitMap.put(hit.getCellID(), hit); }
+
+ //Map a crystal to a list of all clusters in which it is a member.
+ Map<CalorimeterHit, List<CalorimeterHit>> commonHits = new HashMap<CalorimeterHit, List<CalorimeterHit>>();
+
+ //Map a crystal to the seed of the cluster of which it is a member.
+ HashMap<CalorimeterHit, CalorimeterHit> hitSeedMap = new HashMap<CalorimeterHit, CalorimeterHit>();
+
+ // Loop through all calorimeter hits to locate seeds and perform
+ // first pass calculations for component and common hits.
+ for (int ii = 0; ii <= hitList.size() - 1; ii ++){
+ CalorimeterHit hit = hitList.get(ii);
+ // Get the set of all neighboring crystals to the current hit.
+ Set<Long> neighbors = neighborMap.get(hit.getCellID());
+
+ // Generate a list to store any neighboring hits in.
+ ArrayList<CalorimeterHit> neighborHits = new ArrayList<CalorimeterHit>();
+
+ // Sort through the set of neighbors and, if a hit exists
+ // which corresponds to a neighbor, add it to the list of
+ // neighboring hits.
+ for (Long neighbor : neighbors) {
+ // Get the neighboring hit.
+ CalorimeterHit neighborHit = hitMap.get(neighbor);
+
+ // If it exists, add it to the list.
+ if(neighborHit != null) { neighborHits.add(neighborHit); }
+ }
+
+ // Track whether the current hit is a seed hit or not.
+ boolean isSeed = true;
+
+ // Loops through all the neighboring hits to determine if
+ // the current hit is the local maximum within its set of
+ // neighboring hits.
+ seedHitLoop:
+ for(CalorimeterHit neighbor : neighborHits) {
+ if(!equalEnergies(hit, neighbor)) {
+ isSeed = false;
+ break seedHitLoop;
+ }
+ }
+ // If this hit is a seed hit, just map it to itself.
+ if (isSeed && hit.getCorrectedEnergy() >= seedEnergyThreshold) { hitSeedMap.put(hit, hit); }
+
+ // If this hit is a local maximum but does not pass seed threshold,
+ // remove from hit list and do not cluster.
+ else if (isSeed && hit.getCorrectedEnergy() < seedEnergyThreshold){
+ hitList.remove(ii);
+ rejectedHitList.add(hit);
+ ii --;
+ }
+
+ // If this hit is not a seed hit, see if it should be
+ // attached to any neighboring seed hits.
+ else {
+ // Sort through the list of neighboring hits.
+ for (CalorimeterHit neighborHit : neighborHits) {
+ // Check whether the neighboring hit is a seed.
+ if(hitSeedMap.get(neighborHit) == neighborHit) {
+ // If the neighboring hit is a seed hit and the
+ // current hit has been associated with a cluster,
+ // then it is a common hit between its previous
+ // seed and the neighboring seed.
+ if (hitSeedMap.containsKey(hit)) {
+ // Check and see if a list of common seeds
+ // for this hit already exists or not.
+ List<CalorimeterHit> commonHitList = commonHits.get(hit);
+
+ // If it does not, make a new one.
+ if(commonHitList == null) { commonHitList = new ArrayList<CalorimeterHit>(); }
+
+ // Add the neighbors to the seeds to set of
+ // common seeds.
+ commonHitList.add(neighborHit);
+ commonHitList.add(hitSeedMap.get(hit));
+
+
+ // Put the common seed list back into the set.
+ commonHits.put(hit, commonHitList);
+ }
+
+ // If the neighboring hit is a seed hit and the
+ // current hit has not been added to a cluster yet
+ // associate it with the neighboring seed and note
+ // that it has been clustered.
+ else {
+ hitSeedMap.put(hit, neighborHit);
+ }
+ }
+ }
+ }
+ } // End primary seed loop.
+
+ // Performs second pass calculations for component hits.
+ secondaryHitsLoop:
+ for (CalorimeterHit secondaryHit : hitList) {
+ // Look for hits that already have an associated seed/clustering.
+ if(!hitSeedMap.containsKey(secondaryHit)) { continue secondaryHitsLoop; }
+
+ // Get the secondary hit's neighboring crystals.
+ Set<Long> secondaryNeighbors = neighborMap.get(secondaryHit.getCellID());
+
+ // Make a list to store the hits associated with the
+ // neighboring crystals.
+ List<CalorimeterHit> secondaryNeighborHits = new ArrayList<CalorimeterHit>();
+
+ // Loop through the neighboring crystals.
+ for (Long secondaryNeighbor : secondaryNeighbors) {
+ // Get the hit associated with the neighboring crystal.
+ CalorimeterHit secondaryNeighborHit = hitMap.get(secondaryNeighbor);
+
+ // If the neighboring crystal exists and is not already
+ // in a cluster, add it to the list of neighboring hits.
+ if (secondaryNeighborHit != null && !hitSeedMap.containsKey(secondaryNeighborHit)) {
+ secondaryNeighborHits.add(secondaryNeighborHit);
+ }
+ }
+
+ // Loop over the secondary neighbor hits.
+ for (CalorimeterHit secondaryNeighborHit : secondaryNeighborHits) {
+ // If the neighboring hit is of lower energy than the
+ // current secondary hit, then associate the neighboring
+ // hit with the current secondary hit's seed.
+ if(!equalEnergies(secondaryNeighborHit, secondaryHit)){
+ hitSeedMap.put(secondaryNeighborHit, hitSeedMap.get(secondaryHit));}
+ else {continue;}
+ }
+ } // End component hits loop.
+
+
+ // Performs second pass calculations for common hits.
+ commonHitsLoop:
+ for (CalorimeterHit clusteredHit : hitSeedMap.keySet()) {
+
+ // Seed hits are never common hits and can be skipped.
+ if(hitSeedMap.get(clusteredHit) == clusteredHit) { continue commonHitsLoop; }
+
+ // Get the current clustered hit's neighboring crystals.
+ Set<Long> clusteredNeighbors = neighborMap.get(clusteredHit.getCellID());
+
+ // Store a list of all the clustered hits neighboring
+ // crystals which themselves contain hits.
+ List<CalorimeterHit> clusteredNeighborHits = new ArrayList<CalorimeterHit>();
+
+ // Loop through the neighbors and see if they have hits.
+ for (Long neighbor : clusteredNeighbors) {
+ // Get the hit associated with the neighbor.
+ CalorimeterHit clusteredNeighborHit = hitMap.get(neighbor);
+
+ // If it exists, add it to the neighboring hit list.
+
+ if (clusteredNeighborHit != null && hitSeedMap.get(clusteredNeighborHit) != null) {
+ clusteredNeighborHits.add(clusteredNeighborHit);
+ }
+ }
+
+ // Get the seed hit associated with this clustered hit.
+ CalorimeterHit clusteredHitSeed = hitSeedMap.get(clusteredHit);
+
+
+ // Loop over the clustered neighbor hits.
+ for (CalorimeterHit clusteredNeighborHit : clusteredNeighborHits) {
+ // Check to make sure that the clustered neighbor hit
+ // is not already associated with the current clustered
+ // hit's seed.
+
+ if ((hitSeedMap.get(clusteredNeighborHit) != clusteredHitSeed)){
+ // Check for lowest energy hit and that comparison hit is not already common.
+ // If already common, this boundary is already accounted for.
+ if(!equalEnergies(clusteredHit, clusteredNeighborHit)
+ && !commonHits.containsKey(clusteredNeighborHit)){
+
+ // Check and see if a list of common seeds
+ // for this hit already exists or not.
+ List<CalorimeterHit> commonHitList = commonHits.get(clusteredHit);
+
+ // If it does not, make a new one.
+ if(commonHitList == null) { commonHitList = new ArrayList<CalorimeterHit>();}
+
+ // Add the neighbors to the seeds to set of
+ // common seeds.
+ commonHitList.add(clusteredHitSeed);
+ commonHitList.add(hitSeedMap.get(clusteredNeighborHit));
+
+ // Put the common seed list back into the set.
+ commonHits.put(clusteredHit, commonHitList);
+
+ }
+ }
+ }
+ } // End common hits loop.
+
+
+ // Remove any common hits from the clustered hits collection.
+ for(CalorimeterHit commonHit : commonHits.keySet()) {
+ hitSeedMap.remove(commonHit);
+ }
+
+
+ /*
+ * All hits are sorted from above. The next part of the code is for calculating energies and positions.
+ */
+
+ //Create map to contain the total energy of each cluster
+ Map<CalorimeterHit, Double> seedEnergy = new HashMap<CalorimeterHit, Double>();
+
+ // Get energy of each cluster, excluding common hits
+ for (CalorimeterHit iSeed : hitList) {
+ if(hitSeedMap.get(iSeed) == iSeed) {
+ seedEnergy.put(iSeed, 0.0);
+ }
+ }
+
+ //Putting total cluster energies excluding common hit energies into map with seed keys
+ for (Map.Entry<CalorimeterHit, CalorimeterHit> entry : hitSeedMap.entrySet()) {
+ CalorimeterHit eSeed = entry.getValue();
+ double eEnergy = seedEnergy.get(eSeed);
+ eEnergy += entry.getKey().getCorrectedEnergy();
+ seedEnergy.put(eSeed, eEnergy);
+ }
+
+ // Create a map to contain final uncorrected cluster energies including common hit distributions.
+ Map<CalorimeterHit, Double> seedEnergyTot = seedEnergy;
+
+ //Distribute common hit energies with clusters
+ for (Map.Entry<CalorimeterHit, List<CalorimeterHit>> entry1 : commonHits.entrySet()) {
+ CalorimeterHit commonCell = entry1.getKey();
+ CalorimeterHit seedA = entry1.getValue().get(0);
+ CalorimeterHit seedB = entry1.getValue().get(1);
+ double eFractionA = (seedEnergy.get(seedA))/((seedEnergy.get(seedA)+seedEnergy.get(seedB)));
+ double eFractionB = (seedEnergy.get(seedB))/((seedEnergy.get(seedA)+seedEnergy.get(seedB)));
+ double currEnergyA = seedEnergyTot.get(seedA);
+ double currEnergyB = seedEnergyTot.get(seedB);
+ currEnergyA += eFractionA * commonCell.getCorrectedEnergy();
+ currEnergyB += eFractionB * commonCell.getCorrectedEnergy();
+
+ seedEnergyTot.put(seedA, currEnergyA);
+ seedEnergyTot.put(seedB, currEnergyB);
+ }
+
+
+ // Make mapping to contain clusters with corrected energy.
+ Map<CalorimeterHit, Double> seedEnergyCorr = new HashMap<CalorimeterHit, Double>();
+
+ // Energy Corrections as per HPS Note 2014-001
+ // Iterate through known clusters with energies and apply correction.
+ for (Map.Entry<CalorimeterHit, Double> entryC : seedEnergyTot.entrySet()) {
+ double rawEnergy = entryC.getValue();
+
+ // Energy correction for initial guess of electron:
+ int pdg = 11;
+ double corrEnergy = enCorrection(pdg, rawEnergy);
+
+ seedEnergyCorr.put(entryC.getKey(), corrEnergy);
+ }// end of energy corrections
+
+
+ // Cluster Position as per HPS Note 2014-001
+ // Create map with seed as key to position/centroid value
+ Map<CalorimeterHit, double[]> rawSeedPosition = new HashMap<CalorimeterHit, double[]>();
+ Map<CalorimeterHit, double[]> corrSeedPosition = new HashMap<CalorimeterHit, double[]>();
+
+ // top level iterates through seeds
+ for (Map.Entry<CalorimeterHit, Double> entryS : seedEnergyTot.entrySet()) {
+ //get the seed for this iteration
+ CalorimeterHit seedP = entryS.getKey();
+
+ double xCl = 0.0; // calculated cluster x position, prior to correction
+ double yCl = 0.0; // calculated cluster y position
+ double eNumX = 0.0;
+ double eNumY = 0.0;
+ double eDen = 0.0;
+ double w0 = 3.1;
+
+ // iterate through hits corresponding to seed
+ for (Map.Entry<CalorimeterHit, CalorimeterHit> entryP : hitSeedMap.entrySet()) {
+ if(entryP.getValue() == seedP){
+ ///////////////////////////////
+ // This block fills a map with crystal to center of face of crystal
+ // Get the hit indices as a Point.
+ int ix = entryP.getKey().getIdentifierFieldValue("ix");
+ int iy = entryP.getKey().getIdentifierFieldValue("iy");
+ Point hitIndex = new Point(ix, iy);
+
+ // Get the corrected position for this index pair.
+ double[] position = correctedPositionMap.get(hitIndex);
+
+ // If the result is null, it hasn't been calculated yet.
+ if(position == null) {
+ // Calculate the corrected position.
+ IGeometryInfo geom = entryP.getKey().getDetectorElement().getGeometry();
+ double[] pos = geom.transformLocalToGlobal(VecOp.add(geom.transformGlobalToLocal(geom.getPosition()),(Hep3Vector)new BasicHep3Vector(0,0,-1*((Trd)geom.getLogicalVolume().getSolid()).getZHalfLength()))).v();
+
+ // Convert the result to a Double[] array.
+ position = new double[3];
+ position[0] = pos[0];
+ position[1] = pos[1];
+ position[2] = pos[2];
+
+ // Store the result in the map.
+ correctedPositionMap.put(hitIndex, position);
+ }
+ ///////////////////////////////
+
+ //Use Method 3 weighting scheme described in note:
+ eNumX += Math.max(0.0,(w0+Math.log(entryP.getKey().getCorrectedEnergy()
+ /seedEnergyTot.get(seedP))))*(correctedPositionMap.get(hitIndex)[0]/10.0);
+ eNumY += Math.max(0.0,(w0+Math.log(entryP.getKey().getCorrectedEnergy()
+ /seedEnergyTot.get(seedP))))*(correctedPositionMap.get(hitIndex)[1]/10.0);
+ eDen += Math.max(0.0, (w0+Math.log(entryP.getKey().getCorrectedEnergy()/
+ seedEnergyTot.get(seedP))));
+ }
+
+ }
+
+ xCl = eNumX/eDen;
+ yCl = eNumY/eDen;
+
+ double[] rawPosition = new double[3];
+ rawPosition[0] = xCl*10.0;//mm
+ rawPosition[1] = yCl*10.0;//mm
+ int ix = seedP.getIdentifierFieldValue("ix");
+ int iy = seedP.getIdentifierFieldValue("iy");
+ Point hitIndex = new Point(ix, iy);
+ rawPosition[2] = correctedPositionMap.get(hitIndex)[2];
+
+
+
+ // Apply position correction factors:
+ // Position correction for electron:
+ int pdg = 11;
+ double xCorr = posCorrection(pdg, xCl*10.0, seedEnergyTot.get(seedP));
+
+ double[] corrPosition = new double[3];
+ corrPosition[0] = xCorr*10.0;//mm
+ corrPosition[1] = yCl*10.0;//mm
+ corrPosition[2] = correctedPositionMap.get(hitIndex)[2];
+
+ corrSeedPosition.put(seedP, corrPosition);
+ rawSeedPosition.put(seedP, rawPosition);
+
+
+ }// end of cluster position calculation
+
+
+ /*
+ * Outputs results to cluster collection.
+ */
+ // Only write output if something actually exists.
+ if (hitMap.size() != 0) {
+ // Loop over seeds
+ for (Map.Entry<CalorimeterHit, CalorimeterHit> entry2 : hitSeedMap.entrySet()) {
+ if (entry2.getKey() == entry2.getValue()){
+ if(seedEnergyCorr.get(entry2.getKey())<clusterEnergyThreshold)
+ {
+ //Not clustered for not passing cuts
+ rejectedHitList.add(entry2.getKey());
+ }
+
+ else{
+ // New cluster
+ HPSEcalClusterIC cluster = new HPSEcalClusterIC(entry2.getKey());
+ clusterList.add(cluster);
+ // Loop over hits belonging to seeds
+ for (Map.Entry<CalorimeterHit, CalorimeterHit> entry3 : hitSeedMap.entrySet()) {
+ if (entry3.getValue() == entry2.getValue()) {
+ if(rejectedHitList.contains(entry2.getValue())){
+ rejectedHitList.add(entry3.getKey());
+ }
+ else{
+ // Add hit to cluster
+ cluster.addHit(entry3.getKey());
+ }
+ }
+ }
+
+ for (Map.Entry<CalorimeterHit, List<CalorimeterHit>> entry4 : commonHits.entrySet()) {
+ if (entry4.getValue().contains(entry2.getKey())) {
+ // Add shared hits for energy distribution between clusters
+ cluster.addSharedHit(entry4.getKey());
+ }
+ }
+
+ //Input both raw and corrected cluster energies
+ if (seedEnergyCorr.values().size() > 0){
+ cluster.setEnergy(seedEnergyCorr.get(entry2.getKey()));
+ cluster.setRawEnergy(seedEnergyTot.get(entry2.getKey()));
+ }
+
+ //Input both uncorrected and corrected cluster positions.
+ cluster.setCorrPosition(corrSeedPosition.get(entry2.getKey()));
+ cluster.setRawPosition(rawSeedPosition.get(entry2.getKey()));
+
+
+ }// End checking thresholds and write out.
+ }
+ } //End cluster loop
+// System.out.println("Number of clusters: "+clusterList.size());
+
+
+ } //End event output loop.
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, clusterList, HPSEcalClusterIC.class, flag);
+ event.put(rejectedHitName, rejectedHitList, CalorimeterHit.class, flag);
+
+
+ }
+
+ public void endOfData() {
+ if (writeHits != null) {
+ // Close the event display output writer.
+ try {
+ writeHits.close();
+ } catch (IOException e) {
+ }
+ }
+ }
+
+
+ private static class EnergyComparator implements Comparator<CalorimeterHit> {
+ /**
+ * Compares the first hit with respect to the second. This
+ * method will compare hits first by energy, and then spatially.
+ * In the case of equal energy hits, the hit closest to the
+ * beam gap and closest to the positron side of the detector
+ * will be selected. If all of these conditions are true, the
+ * hit with the positive y-index will be selected. Hits with
+ * all four conditions matching are the same hit.
+ * @param hit1 The hit to compare.
+ * @param hit2 The hit with respect to which the first should
+ * be compared.
+ */
+ public int compare(CalorimeterHit hit1, CalorimeterHit hit2) {
+ // Hits are sorted on a hierarchy by three conditions. First,
+ // the hits with the highest energy come first. Next, they
+ // are ranked by vertical proximity to the beam gap, and
+ // lastly, they are sorted by horizontal proximity to the
+ // positron side of the detector.
+
+ // Get the hit energies.
+ double[] e = { hit1.getCorrectedEnergy(), hit2.getCorrectedEnergy() };
+
+ // Perform the energy comparison. The higher energy hit
+ // will be ordered first.
+ if(e[0] < e[1]) { return 1; }
+ else if(e[0] > e[1]) { return -1; }
+
+ // If the hits are the same energy, we must perform the
+ // spatial comparisons.
+ else {
+ // Get the position with respect to the beam gap.
+ int[] iy = { Math.abs(hit1.getIdentifierFieldValue("iy")), Math.abs(hit2.getIdentifierFieldValue("iy")) };
+
+ // The closest hit is first.
+ if(iy[0] > iy[1]) { return -1; }
+ else if(iy[0] < iy[1]) { return 1; }
+
+ // Hits that are identical in vertical distance from
+ // beam gap and energy are differentiated with distance
+ // horizontally from the positron side of the detector.
+ else {
+ // Get the position from the positron side.
+ int[] ix = { hit1.getIdentifierFieldValue("ix"), hit2.getIdentifierFieldValue("ix") };
+
+ // The closest hit is first.
+ if(ix[0] > ix[1]) { return 1; }
+ else if(ix[0] < ix[1]) { return -1; }
+
+ // If all of these checks are the same, compare
+ // the raw value for iy. If these are identical,
+ // then the two hits are the same. Otherwise, sort
+ // the numerical value of iy. (This removes the
+ // issue where hits (x, y) and (x, -y) can have
+ // the same energy and be otherwise seen as the
+ // same hit from the above checks.
+ else { return Integer.compare(hit1.getIdentifierFieldValue("iy"), hit2.getIdentifierFieldValue("iy")); }
+ }
+ }
+ }
+}
+
+
+
+ /**
+ * Handles pathological case where multiple neighboring crystals have EXACTLY the same energy.
+ * @param hit
+ * @param neighbor Neighbor to hit
+ * @return boolean value of if the hit is a seed
+ */
+ private boolean equalEnergies(CalorimeterHit hit, CalorimeterHit neighbor){
+ boolean isSeed = true;
+
+ int hix = hit.getIdentifierFieldValue("ix");
+ int hiy = Math.abs(hit.getIdentifierFieldValue("iy"));
+ int nix = neighbor.getIdentifierFieldValue("ix");
+ int niy = Math.abs(neighbor.getIdentifierFieldValue("iy"));
+ double hE = hit.getCorrectedEnergy();
+ double nE = neighbor.getCorrectedEnergy();
+ if(hE < nE) {
+ isSeed = false;
+ }
+ else if((hE == nE) && (hiy > niy)) {
+ isSeed = false;
+ }
+ else if((hE == nE) && (hiy == niy) && (hix > nix)) {
+ isSeed = false;
+ }
+ return isSeed;
+ }
+ /**
+ * Calculates energy correction based on cluster raw energy and particle type as per
+ *<a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param pdg Particle id as per PDG
+ * @param rawEnergy Raw Energy of the cluster (sum of hits with shared hit distribution)
+ * @return Corrected Energy
+ */
+ public double enCorrection(int pdg, double rawEnergy){
+ if (pdg == 11) { // Particle is electron
+ return energyCorrection(rawEnergy, ELECTRON_ENERGY_A, ELECTRON_ENERGY_B, ELECTRON_ENERGY_C);
+ }
+ else if (pdg == -11) { //Particle is positron
+ return energyCorrection(rawEnergy, POSITRON_ENERGY_A, POSITRON_ENERGY_B, POSITRON_ENERGY_C);
+ }
+ else if (pdg == 22) { //Particle is photon
+ return energyCorrection(rawEnergy, PHOTON_ENERGY_A, PHOTON_ENERGY_B, PHOTON_ENERGY_C);
+ }
+ else { //Unknown
+ double corrEnergy = rawEnergy;
+ return corrEnergy;}
+
+ }
+
+ /**
+ * Calculates the energy correction to a cluster given the variables from the fit as per
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param rawEnergy Raw energy of the cluster
+ * @param A,B,C from fitting in note
+ * @return Corrected Energy
+ */
+ public double energyCorrection(double rawEnergy, double varA, double varB, double varC){
+ double corrEnergy = rawEnergy / (varA * rawEnergy + varB / (Math.sqrt(rawEnergy)) + varC);
+ return corrEnergy;
+ }
+
+
+ /**
+ * Calculates position correction based on cluster raw energy, x calculated position,
+ * and particle type as per
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param pdg Particle id as per PDG
+ * @param xCl Calculated x centroid position of the cluster, uncorrected, at face
+ * @param rawEnergy Raw energy of the cluster (sum of hits with shared hit distribution)
+ * @return Corrected x position
+ */
+ public double posCorrection(int pdg, double xPos, double rawEnergy){
+ double xCl = xPos/10.0;//convert to mm
+ if (pdg == 11) { //Particle is electron
+ double xCorr = positionCorrection(xCl, rawEnergy, ELECTRON_POS_A, ELECTRON_POS_B, ELECTRON_POS_C, ELECTRON_POS_D, ELECTRON_POS_E);
+ return xCorr*10.0;
+ }
+ else if (pdg == -11) {// Particle is positron
+ double xCorr = positionCorrection(xCl, rawEnergy, POSITRON_POS_A, POSITRON_POS_B, POSITRON_POS_C, POSITRON_POS_D, POSITRON_POS_E);
+ return xCorr*10.0;
+ }
+ else if (pdg == 22) {// Particle is photon
+ double xCorr = positionCorrection(xCl, rawEnergy, PHOTON_POS_A, PHOTON_POS_B, PHOTON_POS_C, PHOTON_POS_D, PHOTON_POS_E);
+ return xCorr*10.0;
+ }
+ else { //Unknown
+ double xCorr = xCl;
+ return xCorr*10.0;}
+ }
+
+
+ /**
+ * Calculates the position correction in cm using the raw energy and variables associated with the fit
+ * of the particle as described in
+ * <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>
+ * @param xCl
+ * @param rawEnergy
+ * @param varA
+ * @param varB
+ * @param varC
+ * @param varD
+ * @param varE
+ * @return
+ */
+ public double positionCorrection(double xCl, double rawEnergy, double varA, double varB, double varC, double varD, double varE){
+ double xCorr = xCl-(varA/Math.sqrt(rawEnergy) + varB )*xCl-
+ (varC*rawEnergy + varD/Math.sqrt(rawEnergy) + varE);
+ return xCorr;
+ }
+
+
+ }
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPEcalClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPEcalClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/GTPEcalClusterer.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,382 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.LinkedList;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.geometry.subdetector.HPSEcal3.NeighborMap;
+import org.lcsim.lcio.LCIOConstants;
+import org.lcsim.util.Driver;
+
+/**
+ * Class <code>GTPCalorimeterClusterer</code> processes events and converts hits
+ * into clusters, where appropriate. It uses the modified 2014 clustering algorithm.<br/>
+ * <br/>
+ * For a hit to be a cluster center, it is required to have an energy above some
+ * tunable minimum threshold. Additionally, the hit must be a local maximum with
+ * respect to its neighbors and itself over a tunable (default 2) clock cycles.
+ * Hits that pass these checks are then required to additional have a total
+ * cluster energy that exceeds another tunable minimum threshold.<br/>
+ * <br/>
+ * A hit is added to a cluster as a component if it has a non-zero energy and
+ * within the aforementioned tunable time buffer used for clustering and is
+ * either at the same location as the seed hit or is a neighbor to the seed hit.
+ * @author Kyle McCarty
+ * @author Sho Uemura
+ */
+public class GTPEcalClusterer extends Driver {
+ Detector detector = null;
+
+
+ /**
+ * <b>calorimeter</b><br/><br/>
+ * <code>private HPSEcal3 <b>calorimeter</b></code><br/><br/>
+ * The sub-detector representing the calorimeter.
+ */
+ private HPSEcal3 calorimeter;
+ /**
+ * <b>ecalName</b><br/><br/>
+ * <code>private String <b>ecalName</b></code><br/><br/>
+ * The name of the calorimeter sub-detector stored in the
+ * <code>
+ * Detector</code> object for this run.
+ */
+ private String ecalName;
+ /**
+ * <b>clusterCollectionName</b><br/><br/>
+ * <code>private String <b>clusterCollectionName</b></code><br/><br/>
+ * The name of the LCIO collection name in which the clusters will be
+ * stored.
+ */
+ private String clusterCollectionName = "EcalClusters";
+ /**
+ * <b>clusterWindow</b><br/><br/>
+ * <code>private int <b>clusterWindow</b></code><br/><br/>
+ * Indicates the number of FADC clock cycles (each cycle is 4 ns) before and
+ * after a given cycle that should be considered when checking if a cluster
+ * is a local maximum in space-time.
+ */
+ private int clusterWindow = 2;
+ /**
+ * <b>hitBuffer</b><br/><br/>
+ * <code>private LinkedList<List<CalorimeterHit>> <b>hitBuffer</b></code><br/><br/>
+ * Stores a set of all the hits occurring in each clock cycle for the number
+ * of clock cycles that should be considered for clustering.
+ */
+ private LinkedList<Map<Long, CalorimeterHit>> hitBuffer;
+ /**
+ * <b>ecalCollectionName</b><br/><br/>
+ * <code>private String <b>ecalCollectionName</b></code><br/><br/>
+ * The name of LCIO collection containing the calorimeter hits that are to
+ * be used for clustering.
+ */
+ private String ecalCollectionName;
+ /**
+ * <b>neighborMap</b><br/><br/>
+ * <code>private NeighborMap <b>neighborMap</b></code><br/><br/>
+ * Maps the
+ * <code>long</code> crystal ID to the set of crystal IDs of the crystals
+ * which are adjacent to the key.
+ */
+ private NeighborMap neighborMap = null;
+ /**
+ * <b>seedEnergyThreshold</b><br/><br/>
+ * <code>private double <b>seedEnergyThreshold</b></code><br/><br/>
+ * The minimum energy required for a hit to be considered as a cluster
+ * center. Hits with energy less than this value will be ignored.
+ */
+ private double seedEnergyThreshold = 0.05;
+ /**
+ * <b>limitClusterRange</b><br/><br/>
+ * <code>private boolean <b>limitClusterRange</b></code><br/><br/>
+ * Whether an asymmetric or symmetric window should be used for
+ * adding hits to a cluster.
+ */
+ private boolean limitClusterRange = false;
+
+ /**
+ * Initializes detector-dependent parameters for clustering. Method is
+ * responsible for determining which crystals are valid cluster centers
+ * (stored in
+ * <code>validCrystals</code>), defining the detector object (stored in
+ * <code>calorimeter</code>), defining the detector ID decoder (stored in
+ * <code>decoder</code>), and defining the mapping of crystal IDs to the
+ * crystal's neighbors (defined in
+ * <code>neighborMap</code>).
+ *
+ * @param detector - The new detector to use.
+ */
+ @Override
+ public void detectorChanged(Detector detector) {
+
+ this.detector = detector;
+
+ // Get the calorimeter object.
+ calorimeter = (HPSEcal3) detector.getSubdetector(ecalName);
+
+ // Get a map to associate crystals with their neighbors.
+ neighborMap = calorimeter.getNeighborMap();
+ }
+
+ /**
+ * Generates a list of clusters from the current hit buffer. The "present"
+ * event is taken to be the list of hits occurring at index
+ * <code>clusterWindow</code>, which is the middle of the buffer.
+ *
+ * @return Returns a <code>List</code> of <code>HPSEcalCluster
+ * </code> objects generated from the current event.
+ */
+ public List<HPSEcalCluster> getClusters() {
+ // Generate a list for storing clusters.
+ List<HPSEcalCluster> clusters = new ArrayList<HPSEcalCluster>();
+
+ // Get the list of hits at the current time in the event buffer.
+ Map<Long, CalorimeterHit> currentHits = hitBuffer.get(clusterWindow);
+
+ // For a hit to be a cluster center, it must be a local maximum
+ // both with respect to its neighbors and itself both in the
+ // present time and at all times within the event buffer.
+ seedLoop:
+ for (Long currentID : currentHits.keySet()) {
+ // Get the actual hit object.
+ CalorimeterHit currentHit = currentHits.get(currentID);
+
+ // Store the energy of the current hit.
+ double currentEnergy = currentHit.getRawEnergy();
+
+ // If the hit energy is lower than the minimum threshold,
+ // then we immediately reject this hit as a possible cluster.
+ if (currentEnergy < seedEnergyThreshold) {
+ continue seedLoop;
+ }
+
+ // Store the crystals that are part of this potential cluster,
+ // starting with the cluster seed candidate.
+ HPSEcalCluster cluster = new HPSEcalCluster();
+ cluster.setSeedHit(currentHit);
+ cluster.addHit(currentHit);
+
+ // Get the set of neighbors for this hit.
+ Set<Long> neighbors = neighborMap.get(currentHit.getCellID());
+
+ // Sort through each event stored in the buffer.
+ int bufferIndex = 0;
+ for (Map<Long, CalorimeterHit> bufferHits : hitBuffer) {
+ // Get the hit energy at the current hit's position in
+ // the buffer, if it exists. Ignore the current seed candidate.
+ CalorimeterHit bufferHit = bufferHits.get(currentID);
+ if (bufferHit != null && bufferHit != currentHit) {
+ double bufferHitEnergy = bufferHit.getRawEnergy();
+
+ // Check to see if the hit at this point in the buffer
+ // is larger than then original hit. If it is, we may
+ // stop the comparison because this is not a cluster.
+ if (bufferHitEnergy > currentEnergy) {
+ continue seedLoop;
+ }
+
+ // If the buffer hit is smaller, then add its energy
+ // to the cluster total energy.
+ else {
+ if(limitClusterRange && bufferIndex <= clusterWindow + 1) { cluster.addHit(bufferHit); }
+ else if(!limitClusterRange) { cluster.addHit(bufferHit); }
+ }
+ }
+
+ // We must also make sure that the original hit is
+ // larger than all of the neighboring hits at this
+ // point in the buffer as well.
+ for (Long neighborID : neighbors) {
+ // Get the neighbor hit energy if it exists.
+ CalorimeterHit neighborHit = bufferHits.get(neighborID);
+ if (neighborHit != null) {
+ double neighborHitEnergy = neighborHit.getRawEnergy();
+
+ // Check to see if the neighbor hit at this point
+ // in the buffer is larger than then original hit.
+ // If it is, we may stop the comparison because this
+ // is not a cluster.
+ if (neighborHitEnergy > currentEnergy) {
+ continue seedLoop;
+ }
+
+ // If the buffer neighbor hit is smaller, then
+ // add its energy to the cluster total energy.
+ else {
+ if(limitClusterRange && bufferIndex <= clusterWindow + 1) { cluster.addHit(neighborHit); }
+ else if(!limitClusterRange) { cluster.addHit(neighborHit); }
+ }
+ }
+ }
+
+ // Increment the buffer index.
+ bufferIndex++;
+ }
+
+ // Add the cluster to the list of clusters.
+ clusters.add(cluster);
+ }
+
+ // Return the generated list of clusters.
+ return clusters;
+ }
+
+ /**
+ * Places hits from the current event into the event hit buffer and
+ * processes the buffer to extract clusters. Clusters are then stored in the
+ * event object.
+ *
+ * @param event - The event to process.
+ */
+ @Override
+ public void process(EventHeader event) {
+ // Only run the clusterer if this event has the calorimeter hit collection.
+ // This ensures this driver makes clusters in sync with the FADC readout cycle.
+ if (event.hasCollection(CalorimeterHit.class, ecalCollectionName)) {
+ // Get the list of calorimeter hits from the event.
+ List<CalorimeterHit> hitList = event.get(CalorimeterHit.class, ecalCollectionName);
+
+ // Store each hit in a set by its cell ID so that it may be
+ // easily acquired later.
+ HashMap<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
+ for (CalorimeterHit hit : hitList) {
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ // Remove the last event from the hit buffer and add the new one.
+ hitBuffer.removeLast();
+ hitBuffer.addFirst(hitMap);
+
+ // Run the clustering algorithm on the buffer.
+ List<HPSEcalCluster> clusterList = getClusters();
+
+ // Store the cluster list in the LCIO collection.
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, clusterList, HPSEcalCluster.class, flag);
+ }
+ }
+
+ /**
+ * Sets the name of the LCIO collection containing the calorimeter hits
+ * which should be used for clustering.
+ *
+ * @param ecalCollectionName - The appropriate collection name.
+ */
+ public void setEcalCollectionName(String ecalCollectionName) {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ /**
+ * Sets the name of the LCIO collection in which the clusters should be
+ * stored.
+ *
+ * @param clusterCollectionName - The desired collection name.
+ */
+ public void setClusterCollectionName(String clusterCollectionName) {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ /**
+ * Sets the name of the calorimeter sub-detector stored in the
+ * <code>Detector</code> object that is to be used for this run.
+ *
+ * @param ecalName - The calorimeter's name.
+ */
+ public void setEcalName(String ecalName) {
+ this.ecalName = ecalName;
+ }
+
+ /**
+ * Sets the number of clock cycles before and after a given cycle that will
+ * be used when checking whether a given hit is a local maximum in both time
+ * and space. Note that a value of
+ * <code>N
+ * </code> indicates that
+ * <code>N</code> clock cycles before and
+ * <code>N</code> clock cycles after will be considered. Thusly, a total of
+ * <code>2N + 1</code> clock cycles will be used total.
+ *
+ * @param clusterWindow - The number of additional clock cycles to include
+ * in the clustering checks. A negative value will be treated as zero.
+ */
+ public void setClusterWindow(int clusterWindow) {
+ // The cluster window of must always be at least zero.
+ if (clusterWindow < 0) {
+ this.clusterWindow = 0;
+ }
+
+ // If the cluster window is non-zero, then store it.
+ else {
+ this.clusterWindow = clusterWindow;
+ }
+ }
+
+ /**
+ * Sets whether hits should be added to a cluster from the entire
+ * cluster window or just the "future" hits, plus one clock-cycle
+ * of "past" hits as a safety buffer to account for time uncertainty.
+ *
+ * @param limitClusterRange - <code>true</code> indicates that
+ * the asymmetric clustering window should be used and <code>
+ * false</code> that the symmetric window should be used.
+ */
+ public void setLimitClusterRange(boolean limitClusterRange) {
+ this.limitClusterRange = limitClusterRange;
+ }
+
+ /**
+ * Sets the minimum energy threshold below which hits will not be considered
+ * as cluster centers.
+ *
+ * @param seedEnergyThreshold - The minimum energy for a cluster center.
+ */
+ public void setSeedEnergyThreshold(double seedEnergyThreshold) {
+ // A negative energy threshold is non-physical. All thresholds
+ // be at least zero.
+ if (seedEnergyThreshold < 0.0) {
+ this.seedEnergyThreshold = 0.0;
+ } // If the energy threshold is valid, then use it.
+ else {
+ this.seedEnergyThreshold = seedEnergyThreshold;
+ }
+ }
+
+ /**
+ * Initializes the clusterer. This ensures that the collection name for the
+ * calorimeter hits from which clusters are to be generated, along with the
+ * calorimeter name, have been defined. Method also initializes the event
+ * hit buffer.
+ */
+ @Override
+ public void startOfData() {
+ // Make sure that there is a cluster collection name into
+ // which clusters may be placed.
+ if (ecalCollectionName == null) {
+ throw new RuntimeException("The parameter ecalCollectionName was not set!");
+ }
+
+ // Make sure that there is a calorimeter detector.
+ if (ecalName == null) {
+ throw new RuntimeException("The parameter ecalName was not set!");
+ }
+
+ // Initiate the hit buffer.
+ hitBuffer = new LinkedList<Map<Long, CalorimeterHit>>();
+
+ // Populate the event buffer with (2 * clusterWindow + 1)
+ // empty events. These empty events represent the fact that
+ // the first few events will not have any events in the past
+ // portion of the buffer.
+ int bufferSize = (2 * clusterWindow) + 1;
+ for (int i = 0; i < bufferSize; i++) {
+ hitBuffer.add(new HashMap<Long, CalorimeterHit>(0));
+ }
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/LegacyClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/LegacyClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/LegacyClusterer.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,97 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import org.hps.recon.ecal.ECalUtils;
+import org.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+
+/**
+ * This Driver creates clusters from the CalorimeterHits of an
+ * {@link org.lcsim.geometry.subdetectur.HPSEcal3} detector.
+ *
+ * The clustering algorithm is from pages 83 and 84 of the HPS Proposal.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Tim Nelson <[log in to unmask]>
+ */
+public class LegacyClusterer extends AbstractClusterer {
+
+ // Minimum E for cluster seed.
+ double minimumClusterSeedEnergy = 0.05 * ECalUtils.GeV;
+
+ // Minimum E to add hit to cluster.
+ double minimumHitEnergy = 0.03 * ECalUtils.GeV;
+
+ void setMinimumClusterSeedEnergy(double minimumClusterSeedEnergy) {
+ this.minimumClusterSeedEnergy = minimumClusterSeedEnergy;
+ }
+
+ void setMinimumHitEnergy(double minimumHitEnergy) {
+ this.minimumHitEnergy = minimumHitEnergy;
+ if (minimumClusterSeedEnergy < minimumHitEnergy) {
+ minimumClusterSeedEnergy = minimumHitEnergy;
+ }
+ }
+
+ public List<Cluster> createClusters(List<CalorimeterHit> hits) {
+
+ Map<Long, CalorimeterHit> hitMap = ClusterUtilities.createHitMap(hits);
+
+ // New Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ // Loop over ECal hits to find cluster seeds.
+ for (CalorimeterHit hit : hitMap.values()) {
+ // Cut on min seed E.
+ if (hit.getRawEnergy() < minimumClusterSeedEnergy) {
+ continue;
+ }
+
+ // Get neighbor crystal IDs.
+ Set<Long> neighbors = neighborMap.get(hit.getCellID());
+
+ // List for neighboring hits.
+ List<CalorimeterHit> neighborHits = new ArrayList<CalorimeterHit>();
+
+ // Loop over neighbors to make hit list for cluster.
+ boolean isSeed = true;
+ for (Long neighborId : neighbors) {
+ // Find the neighbor hit in the event if it exists.
+ CalorimeterHit neighborHit = hitMap.get(neighborId);
+
+ // Was this neighbor cell hit?
+ if (neighborHit != null) {
+ // Check if neighbor cell has more energy.
+ if (neighborHit.getRawEnergy() > hit.getRawEnergy()) {
+ // Neighbor has more energy, so cell is not a seed.
+ isSeed = false;
+ break;
+ }
+
+ // Add to cluster if above min E.
+ if (neighborHit.getRawEnergy() >= minimumHitEnergy) {
+ neighborHits.add(neighborHit);
+ }
+ }
+ }
+
+ // Did we find a seed?
+ if (isSeed) {
+ // Make a cluster from the hit list.
+ HPSEcalCluster cluster = new HPSEcalCluster();
+ cluster.setSeedHit(hit);
+ cluster.addHit(hit);
+ for (CalorimeterHit clusHit : neighborHits) {
+ cluster.addHit(clusHit);
+ }
+ clusters.add(cluster);
+ }
+ }
+ return clusters;
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleInnerCalClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleInnerCalClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleInnerCalClusterer.java Mon Dec 15 15:54:53 2014
@@ -0,0 +1,151 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import org.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+
+/**
+ * This Driver creates clusters from the CalorimeterHits of an
+ * {@link org.lcsim.geometry.subdetector.HPSEcal3} detector.
+ *
+ * Uses basic IC clustering algorithm as given in CLAS note 2004-040: no common
+ * hits (hits are assigned to cluster with largest seed hit energy).
+ *
+ * Hit time information is not used, and multiple hits in the same crystal are
+ * not handled correctly (a warning is printed); optional time cut is applied at
+ * the beginning to discard hits too far from t0.
+ *
+ * @author Holly Szumila-Vance <[log in to unmask]>
+ * @author Sho Uemura <[log in to unmask]>
+ *
+ */
+public class SimpleInnerCalClusterer extends AbstractClusterer {
+
+ //Minimum energy that counts as hit
+ double Emin = 0.001;
+ boolean timeCut = false;
+ double minTime = 0.0;
+ double timeWindow = 20.0;
+
+ /**
+ * Minimum energy for a hit to be used in a cluster. Default of 0.001 GeV..
+ * @param Emin
+ */
+ public void setEmin(double Emin) {
+ this.Emin = Emin;
+ }
+
+ /**
+ * Apply time cuts to hits. Defaults to false.
+ * @param timeCut
+ */
+ public void setTimeCut(boolean timeCut) {
+ this.timeCut = timeCut;
+ }
+
+ /**
+ * Minimum hit time, if timeCut is true. Default of 0 ns.
+ * @param minTime
+ */
+ public void setMinTime(double minTime) {
+ this.minTime = minTime;
+ }
+
+ /**
+ * Width of time window, if timeCut is true. Default of 20 ns.
+ * @param timeWindow
+ */
+ public void setTimeWindow(double timeWindow) {
+ this.timeWindow = timeWindow;
+ }
+
+ public List<Cluster> createClusters(List<CalorimeterHit> allHits) {
+
+ // New Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ //Create a Calorimeter hit list in each event, then sort with highest energy first
+ ArrayList<CalorimeterHit> sortedHitList = new ArrayList<CalorimeterHit>(allHits.size());
+ for (CalorimeterHit h : allHits) {
+ //reject hits below the energy cut
+ if (h.getCorrectedEnergy() < Emin) {
+ continue;
+ }
+ //if time cut is being used, reject hits outside the time window
+ if (timeCut && (h.getTime() < minTime || h.getTime() > minTime + timeWindow)) {
+ continue;
+ }
+ sortedHitList.add(h);
+ }
+
+ //sort the list, highest energy first
+ Collections.sort(sortedHitList, Collections.reverseOrder(new EnergyComparator()));
+
+ //map from seed hit to cluster
+ Map<CalorimeterHit, HPSEcalCluster> seedToCluster = new HashMap<CalorimeterHit, HPSEcalCluster>();
+
+ //Quick Map to access hits from cell IDs
+ Map<Long, CalorimeterHit> idToHit = new HashMap<Long, CalorimeterHit>();
+
+ //map from each hit to its cluster seed
+ Map<CalorimeterHit, CalorimeterHit> hitToSeed = new HashMap<CalorimeterHit, CalorimeterHit>();
+
+ //Fill Map with cell ID and hit
+ for (CalorimeterHit hit : sortedHitList) {
+ //if (idToHit.containsKey(hit.getCellID())) {
+ // System.out.println(this.getName() + ": multiple CalorimeterHits in same crystal");
+ //}
+ idToHit.put(hit.getCellID(), hit);
+ }
+
+ for (CalorimeterHit hit : sortedHitList) {
+ CalorimeterHit biggestSeed = null;
+
+ for (Long neighbor : neighborMap.get(hit.getCellID())) {
+ CalorimeterHit neighborHit = idToHit.get(neighbor);
+ if (neighborHit == null) {
+ continue;
+ }
+
+ if (neighborHit.getCorrectedEnergy() > hit.getCorrectedEnergy()) {
+ CalorimeterHit neighborSeed = hitToSeed.get(neighborHit);
+ if (biggestSeed == null || neighborHit.getCorrectedEnergy() > biggestSeed.getCorrectedEnergy()) {
+ biggestSeed = neighborSeed;
+ }
+ }
+ }
+ if (biggestSeed == null) { //if no neighbors had more energy than this hit, this hit is a seed
+ hitToSeed.put(hit, hit);
+ HPSEcalCluster cluster = new HPSEcalCluster(hit.getCellID());
+ clusters.add(cluster);
+ seedToCluster.put(hit, cluster);
+ } else {
+ hitToSeed.put(hit, biggestSeed);
+ }
+ }
+
+ //add all hits to clusters
+ for (CalorimeterHit hit : sortedHitList) {
+ CalorimeterHit seed = hitToSeed.get(hit);
+ HPSEcalCluster cluster = seedToCluster.get(seed);
+ cluster.addHit(hit);
+ }
+
+ return clusters;
+ }
+
+ private class EnergyComparator implements Comparator<CalorimeterHit> {
+
+ @Override
+ public int compare(CalorimeterHit o1, CalorimeterHit o2) {
+ return Double.compare(o1.getCorrectedEnergy(), o2.getCorrectedEnergy());
+ }
+ }
+}
