Author: [log in to unmask]
Date: Fri Dec 19 16:22:04 2014
New Revision: 1796
Log:
Omnibus update to new clustering package. Add javadoc. Improve interfaces and implementations. Add cosmic clusterers from analysis module.
Added:
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterDriver.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NearestNeighborClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCuts.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCutsImpl.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleCosmicClusterer.java
Modified:
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/ClasInnerCalClusterDriver.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterer.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/LegacyClusterer.java
java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleClasInnerCalClusterer.java
Modified: 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 (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/AbstractClusterer.java Fri Dec 19 16:22:04 2014
@@ -4,6 +4,7 @@
import org.hps.conditions.database.DatabaseConditionsManager;
import org.lcsim.conditions.ConditionsEvent;
+import org.lcsim.detector.identifier.IIdentifierHelper;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.Cluster;
import org.lcsim.event.EventHeader;
@@ -14,7 +15,10 @@
* This is an abstract class that {@link Clusterer} classes should implement
* to perform a clustering algorithm on a <code>CalorimeterHit</code> collection.
* The sub-class should implement {@link #createClusters(List)} which is
- * the method that should perform the clustering.
+ * the method that should perform the clustering algorithm.
+ *
+ * @see Clusterer
+ * @see org.lcsim.event.Cluster
*
* @author Jeremy McCormick <[log in to unmask]>
*/
@@ -22,9 +26,42 @@
protected HPSEcal3 ecal;
protected NeighborMap neighborMap;
- protected double[] cuts;
- protected double[] defaultCuts;
- protected String[] cutNames;
+ protected NumericalCuts cuts;
+
+ /**
+ * Default constructor which takes names of cuts and their default values.
+ * These arguments cannot be null. (Instead use the no-arg constructor.)
+ * @param cutNames The names of the cuts for this clustering algorithm.
+ * @param defaultCuts The default cut values for the algorithm matching the cutNames ordering.
+ * @throw IllegalArgumentException if the arguments are null or the arrays are different lengths.
+ */
+ protected AbstractClusterer(String cutNames[], double[] defaultCuts) {
+ if (cutNames == null) {
+ throw new IllegalArgumentException("The cutNames is set to null.");
+ }
+ if (defaultCuts == null) {
+ throw new IllegalArgumentException("The defaultCuts is set to null.");
+ }
+ if (cutNames.length != defaultCuts.length) {
+ throw new IllegalArgumentException("The cutNames and defaultCuts are not the same length.");
+ }
+ cuts = new NumericalCutsImpl(cutNames, defaultCuts);
+ }
+
+ /**
+ * Constructor with cuts set.
+ * @param cuts The numerical cuts.
+ */
+ protected AbstractClusterer(NumericalCuts cuts) {
+ this.cuts = cuts;
+ }
+
+ /**
+ * Default no-arg constructor.
+ */
+ protected AbstractClusterer() {
+ // No cuts are set.
+ }
/**
* This is the primary method for sub-classes to implement their clustering algorithm.
@@ -50,83 +87,20 @@
*/
public void initialize() {
}
-
+
/**
- * Default constructor which takes names of cuts and their default values.
- * Even if there are no cuts, these should be arrays of length 0 instead of null.
- * @param cutNames The names of the cuts for this clustering algorithm.
- * @param defaultCuts The default cut values for the algorithm matching the cutNames ordering.
- * @throw IllegalArgumentException if the arguments are null or the arrays are different lengths.
+ * Get the numerical cut settings.
+ * @return The numerical cuts.
*/
- protected AbstractClusterer(String cutNames[], double[] defaultCuts) {
- if (cutNames == null) {
- throw new IllegalArgumentException("The cutNames is set to null.");
- }
- if (defaultCuts == null) {
- throw new IllegalArgumentException("The defaultCuts is set to null.");
- }
- if (cutNames.length != defaultCuts.length) {
- throw new IllegalArgumentException("The cutNames and defaultCuts are not the same length.");
- }
- this.cutNames = cutNames;
- this.defaultCuts = defaultCuts;
- this.cuts = defaultCuts;
- }
-
- public void setCuts(double[] cuts) {
- if (cuts.length != this.cutNames.length) {
- throw new IllegalArgumentException("The cuts array has the wrong length: " + cuts.length);
- }
- this.cuts = cuts;
+ public NumericalCuts getCuts() {
+ return this.cuts;
}
- public double[] getCuts() {
- return cuts;
+ /**
+ * Convenience method to get the identifier helper from the ECAL subdetector.
+ * @return The identifier helper.
+ */
+ protected IIdentifierHelper getIdentifierHelper() {
+ return ecal.getDetectorElement().getIdentifierHelper();
}
-
- public double getCut(String name) {
- int index = indexFromName(name);
- if (index == -1) {
- throw new IllegalArgumentException("There is no cut called " + name + " defined by this clusterer.");
- }
- return getCut(index);
- }
-
- public double getCut(int index) {
- if (index > cuts.length || index < 0) {
- throw new IndexOutOfBoundsException("The index " + index + " is out of bounds for cuts array.");
- }
- return cuts[index];
- }
-
- public String[] getCutNames() {
- return cutNames;
- }
-
- @Override
- public void setCut(int index, double value) {
- cuts[index] = value;
- }
-
- public boolean isDefaultCuts() {
- return cuts == defaultCuts;
- }
-
- public double[] getDefaultCuts() {
- return defaultCuts;
- }
-
- public void setCut(String name, double value) {
- int index = indexFromName(name);
- cuts[index] = value;
- }
-
- protected int indexFromName(String name) {
- for (int index = 0; index < cuts.length; index++) {
- if (getCutNames()[index] == name) {
- return index;
- }
- }
- return -1;
- }
}
Modified: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterDriver.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterDriver.java (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterDriver.java Fri Dec 19 16:22:04 2014
@@ -10,6 +10,8 @@
* This is an implementation of {@link ClusterDriver} specialized for the
* {@link ClasInnerCalClusterer}. It currently implements optional
* writing of a rejected hit list to the LCSim event.
+ *
+ * @see ClasInnerCalClusterer
*
* @author Jeremy McCormick <[log in to unmask]>
*/
@@ -30,7 +32,7 @@
public void startOfData() {
if (clusterer == null) {
// Setup the Clusterer if it wasn't already initialized by a Driver argument.
- this.setClusterer("ClasInnerCalClusterer");
+ this.setClustererName("ClasInnerCalClusterer");
} else {
// Does the Clusterer have the right type if there was a custom inialization parameter?
if (!(clusterer instanceof ClasInnerCalClusterer)) {
Modified: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterer.java (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClasInnerCalClusterer.java Fri Dec 19 16:22:04 2014
@@ -24,16 +24,20 @@
* <p>
* This class creates clusters from a CalorimeterHit input collection.
* <p>
- * This clustering logic is based on that from the CLAS-Note-2005-001.
+ * This clustering logic is based on that from the
+ * <a href="https://misportal.jlab.org/ul/Physics/Hall-B/clas/viewFile.cfm/2005-001.pdf?documentId=6">CLAS-Note-2005-001</a>.
* <p>
* The analysis and position corrections are described in
* <a href="https://misportal.jlab.org/mis/physics/hps_notes/index.cfm?note_year=2014">HPS Note 2014-001</a>.
* <p>
- * The algorithm sorts hits from highest to lowest energy and build clusters around each local maximum/seed hit.
+ * The algorithm sorts hits from highest to lowest energy and builds 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. Energy
* corrections are applied separately.
+ *
+ * @see AbstractClusterer
+ * @see Clusterer
*
* @author Holly Szumila-Vance <[log in to unmask]>
* @author Kyle McCarty <[log in to unmask]>
@@ -66,64 +70,113 @@
// Make a map for quick calculation of the x-y position of crystal face
public Map<Point, double[]> correctedPositionMap = new HashMap<Point, double[]>();
-
- // 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;
+ // List of rejected hits.
List<CalorimeterHit> currentRejectedHitList;
public ClasInnerCalClusterer() {
- super(new String[] { "hitEnergyThreshold", "seedEnergyThreshold", "clusterEnergyThreshold", "minTime", "timeWindow", "timeCut" }, new double[] { 0.0075, 0.1, 0.3, 0.0, 20.0, 0. });
+ super(new String[] { "hitEnergyThreshold", "seedEnergyThreshold", "clusterEnergyThreshold", "minTime", "timeWindow", "timeCut" },
+ new double[] { 0.0075, 0.1, 0.3, 0.0, 20.0, 0. });
}
public void initialize() {
- hitEnergyThreshold = this.getCut("hitEnergyThreshold");
- seedEnergyThreshold = this.getCut("seedEnergyThreshold");
- clusterEnergyThreshold = this.getCut("clusterEnergyThreshold");
- minTime = this.getCut("minTime");
- timeWindow = this.getCut("timeWindow");
- timeCut = (this.getCut("timeCut") == 1.);
+ hitEnergyThreshold = getCuts().getValue("hitEnergyThreshold");
+ seedEnergyThreshold = getCuts().getValue("seedEnergyThreshold");
+ if (hitEnergyThreshold > seedEnergyThreshold) {
+ throw new IllegalArgumentException("hitEnergyThreshold must be <= to seedEnergyThreshold");
+ }
+ clusterEnergyThreshold = getCuts().getValue("clusterEnergyThreshold");
+ minTime = getCuts().getValue("minTime");
+ timeWindow = getCuts().getValue("timeWindow");
+ timeCut = (getCuts().getValue("timeCut") == 1.);
}
public List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hits) {
- // Create a list to store the event hits in.
- List<CalorimeterHit> hitList = new ArrayList<CalorimeterHit>();
- List<CalorimeterHit> baseList = hits;
- for (CalorimeterHit r : baseList) {
- hitList.add(r);
- }
-
- // Create a list to store the newly created clusters in.
- ArrayList<Cluster> clusterList = new ArrayList<Cluster>();
-
- // Reset the reference to the rejected hit list.
- currentRejectedHitList = new ArrayList<CalorimeterHit>();
+ // Clear corrected position map (or there is a potential memory leak!).
+ this.correctedPositionMap.clear();
+
+ // Reset the the rejected hit list.
+ resetRejectedHitList();
+
+ // Create the master list of hits to work from for clustering.
+ List<CalorimeterHit> hitList = new ArrayList<CalorimeterHit>(hits);
// Sort the list of hits by energy.
Collections.sort(hitList, ENERGY_COMP);
+ // Filter the hit list.
+ filterHitList(hitList);
+
+ // Create an ID map of the hits.
+ Map<Long, CalorimeterHit> hitMap = ClusterUtilities.createHitMap(hitList);
+
+ // 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>();
+
+ // Perform first pass calculations to find seed hits and do initial calculations on hits.
+ firstPass(hitList, hitMap, commonHits, hitSeedMap);
+
+ // Perform second pass calculations for component hits.
+ secondPass(hitList, hitMap, hitSeedMap);
+
+ // Perform third pass for calculations on common hits.
+ thirdPass(hitMap, commonHits, hitSeedMap);
+
+ // This removes common hits from the seed hit map. (Why?)
+ removeCommonHits(commonHits, hitSeedMap);
+
+ // Create map to contain the total energy of each cluster
+ Map<CalorimeterHit, Double> seedEnergy = new HashMap<CalorimeterHit, Double>();
+
+ // ?????
+ getClusterEnergies(hitList, hitSeedMap, seedEnergy);
+
+ // Create map of seeds to energies.
+ createSeedMap(hitSeedMap, seedEnergy);
+
+ // Create a map to contain final uncorrected cluster energies including common hit distributions.
+ Map<CalorimeterHit, Double> seedEnergyTot = seedEnergy;
+
+ // Distribute energies of common hits.
+ distributeCommonHitEnergies(commonHits, seedEnergy, seedEnergyTot);
+
+ //
+ // Correct the position as per HPS Note 2014-001. (See class doc for link.)
+ //
+ // Create maps with seed as key to position/centroid value.
+ Map<CalorimeterHit, double[]> rawSeedPosition = new HashMap<CalorimeterHit, double[]>();
+ Map<CalorimeterHit, double[]> corrSeedPosition = new HashMap<CalorimeterHit, double[]>();
+ correctClusterPositions(hitSeedMap, seedEnergyTot, rawSeedPosition, corrSeedPosition);
+
+ // Create the cluster collection.
+ return createClusterCollection(hitMap, commonHits, hitSeedMap, seedEnergyTot, rawSeedPosition, corrSeedPosition);
+ }
+
+ private void createSeedMap(HashMap<CalorimeterHit, CalorimeterHit> hitSeedMap, Map<CalorimeterHit, Double> seedEnergy) {
+ // 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);
+ }
+ }
+
+ private void resetRejectedHitList() {
+ if (currentRejectedHitList != null) {
+ currentRejectedHitList.clear();
+ }
+ currentRejectedHitList = new ArrayList<CalorimeterHit>();
+ }
+
+ private void filterHitList(List<CalorimeterHit> hitList) {
// Filter the hit list of any hits that fail to pass the
// designated threshold.
- filterLoop: for (int index = hitList.size() - 1; index >= 0; index--) {
+ 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)))) {
currentRejectedHitList.add(hitList.get(index));
@@ -138,21 +191,287 @@
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.
+ }
+
+ private void getClusterEnergies(List<CalorimeterHit> hitList, Map<CalorimeterHit, CalorimeterHit> hitSeedMap, Map<CalorimeterHit, Double> seedEnergy) {
+ // Get energy of each cluster, excluding common hits
+ for (CalorimeterHit iSeed : hitList) {
+ if (hitSeedMap.get(iSeed) == iSeed) {
+ seedEnergy.put(iSeed, 0.0);
+ }
+ }
+ }
+
+ /*
+ * Outputs results to cluster collection.
+ */
+ private List<Cluster> createClusterCollection(Map<Long, CalorimeterHit> hitMap, Map<CalorimeterHit, List<CalorimeterHit>> commonHits, HashMap<CalorimeterHit, CalorimeterHit> hitSeedMap, Map<CalorimeterHit, Double> seedEnergyTot, Map<CalorimeterHit, double[]> rawSeedPosition, Map<CalorimeterHit, double[]> corrSeedPosition) {
+
+ // Create a list to store the newly created clusters.
+ ArrayList<Cluster> clusterList = new ArrayList<Cluster>();
+
+ // 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 (seedEnergyTot.get(entry2.getKey()) < clusterEnergyThreshold) {
+ // Not clustered for not passing cuts
+ currentRejectedHitList.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 (currentRejectedHitList.contains(entry2.getValue())) {
+ currentRejectedHitList.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 uncorrected cluster energies
+ if (seedEnergyTot.values().size() > 0) {
+ cluster.setEnergy(seedEnergyTot.get(entry2.getKey()));
+ cluster.setUncorrectedEnergy(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.
+ return clusterList;
+ }
+
+ private void correctClusterPositions(HashMap<CalorimeterHit, CalorimeterHit> hitSeedMap, Map<CalorimeterHit, Double> seedEnergyTot, Map<CalorimeterHit, double[]> rawSeedPosition, Map<CalorimeterHit, double[]> corrSeedPosition) {
+ // 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 = correctPosition(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
+ }
+
+ // Distribute common hit energies with clusters
+ private void distributeCommonHitEnergies(Map<CalorimeterHit, List<CalorimeterHit>> commonHits, Map<CalorimeterHit, Double> seedEnergy, Map<CalorimeterHit, Double> seedEnergyTot) {
+ 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);
+ }
+ }
+
+ // Remove any common hits from the clustered hits collection.
+ private void removeCommonHits(Map<CalorimeterHit, List<CalorimeterHit>> commonHits, HashMap<CalorimeterHit, CalorimeterHit> hitSeedMap) {
+ for (CalorimeterHit commonHit : commonHits.keySet()) {
+ hitSeedMap.remove(commonHit);
+ }
+ }
+
+ // Performs second pass calculations for common hits.
+ private void thirdPass(Map<Long, CalorimeterHit> hitMap, Map<CalorimeterHit, List<CalorimeterHit>> commonHits, Map<CalorimeterHit, CalorimeterHit> hitSeedMap) {
+ 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.
+ }
+
+ private void secondPass(List<CalorimeterHit> hitList, Map<Long, CalorimeterHit> hitMap, Map<CalorimeterHit, CalorimeterHit> hitSeedMap) {
+ // 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.
+ }
+
+ // Loop through all calorimeter hits to locate seeds and perform first pass calculations for component and common hits.
+ private void firstPass(List<CalorimeterHit> hitList, Map<Long, CalorimeterHit> hitMap, Map<CalorimeterHit, List<CalorimeterHit>> commonHits, Map<CalorimeterHit, CalorimeterHit> hitSeedMap) {
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.
@@ -240,310 +559,24 @@
}
}
} // 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);
- }
-
- // 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 (seedEnergyTot.get(entry2.getKey()) < clusterEnergyThreshold) {
- // Not clustered for not passing cuts
- currentRejectedHitList.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 (currentRejectedHitList.contains(entry2.getValue())) {
- currentRejectedHitList.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 uncorrected cluster energies
- if (seedEnergyTot.values().size() > 0) {
- cluster.setEnergy(seedEnergyTot.get(entry2.getKey()));
- cluster.setUncorrectedEnergy(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.
- return clusterList;
}
/**
* Get the list of rejected hits that was made from processing the last event.
* @return The list of rejected hit.
*/
- public List<CalorimeterHit> getRejectedHitList() {
+ List<CalorimeterHit> getRejectedHitList() {
return this.currentRejectedHitList;
}
/**
- * Handles pathological case where multiple neighboring crystals have EXACTLY the same energy.
- *
+ * Handles pathological case where multiple neighboring crystals have EXACTLY the same energy.
* @param hit
- * @param neighbor
- * Neighbor to 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");
@@ -559,54 +592,8 @@
}
return isSeed;
}
-
- /**
- * 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;
- }
-
+
+ // FIXME: Is this really necessary?
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
@@ -674,4 +661,70 @@
}
}
}
+
+ // 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;
+
+ /**
+ * 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
+ */
+ static double correctPosition(int pdg, double xPos, double rawEnergy) {
+ double xCl = xPos / 10.0;// convert to mm
+ if (pdg == 11) { // Particle is electron
+ double xCorr = correctPosition(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 = correctPosition(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 = correctPosition(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
+ */
+ static double correctPosition(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;
+ }
+
}
Modified: 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 (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterDriver.java Fri Dec 19 16:22:04 2014
@@ -3,26 +3,29 @@
import java.util.List;
import java.util.logging.Logger;
-import org.lcsim.detector.identifier.IIdentifierHelper;
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.geometry.subdetector.HPSEcal3.NeighborMap;
import org.lcsim.lcio.LCIOConstants;
import org.lcsim.util.Driver;
+import org.lcsim.util.log.BasicFormatter;
import org.lcsim.util.log.LogUtil;
-import org.lcsim.util.log.BasicFormatter;
/**
* <p>
* This is a basic Driver that creates ECAL <code>Cluster</code> collections
* through the {@link Clusterer} interface.
* <p>
- * A specific clustering engine can be created with the {@link #setClusterer(String)} method
+ * A specific clustering engine can be created with the {@link #setClustererName(String)} method
* which will use a factory to create it by name. The cuts of the {@link Clusterer}
* can be set generically with the {@link #setCuts(double[])} method.
+ *
+ * @see Clusterer
+ * @see org.lcsim.event.Cluster
*
* @author Jeremy McCormick <[log in to unmask]>
*/
@@ -30,7 +33,7 @@
protected String ecalName = "Ecal";
protected HPSEcal3 ecal;
- protected IIdentifierHelper helper;
+ protected NeighborMap neighborMap;
protected String outputClusterCollectionName = "EcalClusters";
protected String inputHitCollectionName = "EcalCalHits";
protected Clusterer clusterer;
@@ -39,8 +42,8 @@
protected boolean skipNoClusterEvents = false;
protected boolean writeClusterCollection = true;
protected boolean storeHits = true;
- protected double[] cuts;
- protected Logger logger = LogUtil.create(ClusterDriver.class, new BasicFormatter(ClusterDriver.class.getSimpleName()));
+ protected double[] cutValues;
+ protected static Logger logger = LogUtil.create(ClusterDriver.class, new BasicFormatter(ClusterDriver.class.getSimpleName()));
/**
* No arg constructor.
@@ -118,7 +121,7 @@
* class name and try to instantiate it using the Class API.
* @param The name or canonical class name of the Clusterer.
*/
- public void setClusterer(String name) {
+ public void setClustererName(String name) {
clusterer = ClustererFactory.create(name);
this.getLogger().config("Clusterer was set to " + this.clusterer.getClass().getSimpleName());
}
@@ -145,7 +148,7 @@
* @param cuts The numerical cuts.
*/
public void setCuts(double[] cuts) {
- this.cuts = cuts;
+ this.cutValues = cuts;
}
/**
@@ -161,7 +164,7 @@
throw new RuntimeException("Ther subdetector " + ecalName + " does not have the right type.");
}
ecal = (HPSEcal3) subdetector;
- helper = ecal.getDetectorElement().getIdentifierHelper();
+ neighborMap = ecal.getNeighborMap();
}
/**
@@ -172,13 +175,14 @@
if (this.clusterer == null) {
throw new RuntimeException("The clusterer was never initialized.");
}
- if (this.cuts != null) {
+ if (this.cutValues != null) {
logger.config("setting cuts on clusterer");
- this.clusterer.setCuts(cuts);
- for (int cutIndex = 0; cutIndex < clusterer.getCuts().length; cutIndex++) {
- logger.config(" " + this.clusterer.getCutNames()[cutIndex] + " = " + this.clusterer.getCut(cutIndex));
- }
+ this.clusterer.getCuts().setValues(cutValues);
}
+ logger.config("Clusterer has the following cuts ...");
+ for (int cutIndex = 0; cutIndex < clusterer.getCuts().getValues().length; cutIndex++) {
+ logger.config(" " + this.clusterer.getCuts().getNames()[cutIndex] + " = " + this.clusterer.getCuts().getValue(cutIndex));
+ }
logger.config("initializing clusterer");
this.clusterer.initialize();
}
@@ -190,13 +194,13 @@
this.getLogger().fine("processing event #" + event.getEventNumber());
if (event.hasCollection(CalorimeterHit.class, inputHitCollectionName)) {
List<CalorimeterHit> hits = event.get(CalorimeterHit.class, inputHitCollectionName);
- this.getLogger().fine("Input hit collection " + inputHitCollectionName + " has " + hits.size() + " hits.");
+ logger.fine("input hit collection " + inputHitCollectionName + " has " + hits.size() + " hits");
List<Cluster> clusters = clusterer.createClusters(event, hits);
if (clusters == null) {
- throw new RuntimeException("The clusterer returned null from its createClusters method.");
+ throw new RuntimeException("The clusterer returned a null list from its createClusters method.");
}
if (clusters.isEmpty() && this.skipNoClusterEvents) {
- logger.finer("Skipping event with no clusters.");
+ logger.finer("skipping event with no clusters");
throw new NextEventException();
}
if (event.hasCollection(Cluster.class, this.outputClusterCollectionName)) {
@@ -237,7 +241,7 @@
*/
@SuppressWarnings("unchecked")
<ClustererType extends Clusterer> ClustererType getClusterer() {
- // Return the Clusterer casting to the right type, which should always work because ClustererType must extend Clusterer.
+ // Return the Clusterer and cast it to the type provided by the caller.
return (ClustererType) clusterer;
}
}
Modified: 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 (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterUtilities.java Fri Dec 19 16:22:04 2014
@@ -1,5 +1,7 @@
package org.hps.recon.ecal.cluster;
+import java.util.ArrayList;
+import java.util.Collections;
import java.util.Comparator;
import java.util.HashSet;
import java.util.LinkedHashMap;
@@ -13,7 +15,10 @@
import org.lcsim.geometry.subdetector.HPSEcal3;
/**
- * This is a set of simple clustering utility methods.
+ * This is a set of simple utility methods for clustering algorithms.
+ *
+ * @see org.lcsim.event.Cluster
+ * @see org.lcsim.event.base.BaseCluster
*/
public final class ClusterUtilities {
@@ -64,14 +69,52 @@
cluster.setEnergy(totalEnergy);
return cluster;
}
+
+ /**
+ * Compute the raw energy of a cluster which is just the
+ * sum of all its hit energies.
+ * @param cluster The input cluster.
+ * @return The total raw energy.
+ */
+ public static double computeRawEnergy(Cluster cluster) {
+ double uncorrectedEnergy = 0;
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ uncorrectedEnergy += hit.getCorrectedEnergy();
+ }
+ return uncorrectedEnergy;
+ }
/**
- * Compare CalorimeterHit objects by their energy using default double comparison strategy.
+ * Find the hit with the highest energy value.
+ * @param cluster The input cluster.
+ * @return The hit with the highest energy value.
*/
- public static class HitEnergyComparator implements Comparator<CalorimeterHit> {
- @Override
- public int compare(CalorimeterHit o1, CalorimeterHit o2) {
- return Double.compare(o1.getCorrectedEnergy(), o2.getCorrectedEnergy());
+ public static CalorimeterHit getHighestEnergyHit(Cluster cluster) {
+ double maxEnergy = Double.MIN_VALUE;
+ CalorimeterHit highestEnergyHit = null;
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ if (hit.getCorrectedEnergy() > maxEnergy) {
+ highestEnergyHit = hit;
+ }
}
- }
+ return highestEnergyHit;
+ }
+
+ /**
+ * Sort the hits in the cluster using a <code>Comparator</code>.
+ * This method will not change the hits in place. It returns a new list.
+ * @param cluster The input cluster.
+ * @param comparator The Comparator to use for sorting.
+ * @param reverseOrder True to use reverse rather than default ordering.
+ * @return The sorted list of hits.
+ */
+ public static List<CalorimeterHit> sortedHits(Cluster cluster, Comparator<CalorimeterHit> comparator, boolean reverseOrder) {
+ List<CalorimeterHit> sortedHits = new ArrayList<CalorimeterHit>(cluster.getCalorimeterHits());
+ Comparator<CalorimeterHit>sortComparator = comparator;
+ if (reverseOrder) {
+ sortComparator = Collections.reverseOrder(comparator);
+ }
+ Collections.sort(sortedHits, sortComparator);
+ return sortedHits;
+ }
}
Modified: 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 (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/Clusterer.java Fri Dec 19 16:22:04 2014
@@ -10,6 +10,11 @@
/**
* This is an interface for creating clusters and providing cut values
* to the clustering algorithms in a generic fashion.
+ *
+ * @see org.lcsim.event.Cluster
+ * @see org.lcsim.event.CalorimeterHit
+ * @see org.lcsim.event.EventHeader
+ * @see NumericalCuts
*/
public interface Clusterer extends ConditionsListener {
@@ -22,65 +27,28 @@
List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hits);
/**
- * Perform start of job intialization on this object.
+ * <p>
+ * Perform start of job initialization on this object.
+ * <p>
+ * This method would typically be used to cache cluster cut values
+ * from the {@link NumericalCuts} into instance variables for
+ * convenience and runtime performance purposes. If the cuts
+ * have certain constraints on their reasonable values for the
+ * algorithm, then this method should throw an <code>IllegalArgumentException</code>
+ * if the parameter value is invalid.
+ * <p>
+ * The Detector object from LCSim is not available yet when this
+ * method is typically called, so the conditions system should not
+ * be used. Instead, the inherited callback method
+ * {@link #conditionsChanged(org.lcsim.conditions.ConditionsEvent)}
+ * can be used to configure the class depending on the available
+ * conditions when they are available.
*/
void initialize();
/**
- * Get the list of numerical cuts.
- * @return The list of numerical cuts.
+ * Get numerical cut settings.
+ * @return The numerical cut settings.
*/
- double[] getCuts();
-
- /**
- * Get the default cut values.
- * @return The default cut values.
- */
- double[] getDefaultCuts();
-
- /**
- * True if algorithm is using its default cuts.
- * @return True if using the default cuts.
- */
- boolean isDefaultCuts();
-
- /**
- * Set numerical cuts array.
- * @param cuts The numerical cuts.
- */
- void setCuts(double[] cuts);
-
- /**
- * Get a cut value by its index.
- * @param index The index of the cut.
- * @return The cut value at index.
- */
- double getCut(int index);
-
- /**
- * Get a cut value by name.
- * @param name The name of the cut.
- * @return The named cut.
- */
- double getCut(String name);
-
- /**
- * Set a cut value by name.
- * @param name The name of the cut.
- * @param value The value of the cut.
- */
- void setCut(String name, double value);
-
- /**
- * Set a cut value by index.
- * @param index The index of the cut.
- * @param value The value of the cut.
- */
- void setCut(int index, double value);
-
- /**
- * Get the names of the cuts.
- * @return The names of the cuts.
- */
- String[] getCutNames();
-}
+ NumericalCuts getCuts();
+}
Modified: 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 (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClustererFactory.java Fri Dec 19 16:22:04 2014
@@ -1,18 +1,24 @@
package org.hps.recon.ecal.cluster;
-import org.lcsim.conditions.ConditionsListener;
import org.lcsim.conditions.ConditionsManager;
/**
* <p>
- * This is a convenience class for creating different kinds of clustering algorithms via their name.
+ * This is a convenience class for creating specific clustering algorithms via their name
+ * in the package <code>org.hps.recon.ecal.cluster</code>. They must implement the {@link Clusterer}
+ * interface.
* <p>
- * The currently available types include:
- * <ul>
- * <li>{@link LegacyClusterer}</li>
- * <li>{@link SimpleClasInnerCalClusterer}</li>
- * <li>{@link ClasInnerCalClusterer</li>
- * </ul>
+ * If the name does not match one in the clustering package, the factory will attempt to create
+ * a new class instance, assuming that the string is a canonical class name. It then checks if
+ * this class implements the {@link Clusterer} interface and will throw an error if it does not.
+ *
+ * @see Clusterer
+ * @see ClasInnerCalClusterer
+ * @see DualThresholdCosmicClusterer
+ * @see LegacyClusterer
+ * @see NearestNeighborClusterer
+ * @see SimpleClasInnerCalClusterer
+ * @see SimpleCosmicClusterer
*
* @author Jeremy McCormick <[log in to unmask]>
*/
@@ -34,30 +40,36 @@
public static Clusterer create(String name, double[] cuts) {
Clusterer clusterer;
if (LegacyClusterer.class.getSimpleName().equals(name)) {
- // Test Run legacy clusterer.
clusterer = new LegacyClusterer();
} else if (SimpleClasInnerCalClusterer.class.getSimpleName().equals(name)) {
- // Simple IC clusterer.
clusterer = new SimpleClasInnerCalClusterer();
} else if (ClasInnerCalClusterer.class.getSimpleName().equals(name)) {
- // Full IC algorithm.
clusterer = new ClasInnerCalClusterer();
+ } else if (NearestNeighborClusterer.class.getSimpleName().equals(name)) {
+ clusterer = new NearestNeighborClusterer();
+ } else if (DualThresholdCosmicClusterer.class.getSimpleName().equals(name)) {
+ clusterer = new DualThresholdCosmicClusterer();
+ } else if (SimpleCosmicClusterer.class.getSimpleName().equals(name)) {
+ clusterer = new SimpleCosmicClusterer();
} else {
- // Try to instantiate the class from the name argument, assuming it is a canonical class name.
+ // Try to instantiate a Clusterer object from the name argument, assuming it is a canonical class name.
try {
clusterer = fromCanonicalClassName(name);
+ if (!clusterer.getClass().isAssignableFrom(Clusterer.class)) {
+ throw new IllegalArgumentException("The clsas " + name + " does not implement the Clusterer interface.");
+ }
} catch (Exception e) {
// Okay nothing worked, so we have a problem!
throw new IllegalArgumentException("Unknown Clusterer type " + name + " cannot be instantiated.", e);
}
}
- // Add the Clusterer as a conditions listener so it can set itself up via the conditions system.
- if (clusterer instanceof ConditionsListener) {
- ConditionsManager.defaultInstance().addConditionsListener((ConditionsListener) clusterer);
- }
+
+ // Register the Clusterer for notification when conditions change.
+ ConditionsManager.defaultInstance().addConditionsListener(clusterer);
+
// Set cuts if they were provided.
if (cuts != null) {
- clusterer.setCuts(cuts);
+ clusterer.getCuts().setValues(cuts);
}
return clusterer;
}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterDriver.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterDriver.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterDriver.java Fri Dec 19 16:22:04 2014
@@ -0,0 +1,75 @@
+package org.hps.recon.ecal.cluster;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.RawTrackerHit;
+
+
+/**
+ * This is a specialized Driver to use the dual threshold cosmic clustering algorithm.
+ * It really only sets up the tight hit collection name for the clustering algorithm.
+ *
+ * @see ClusterDriver
+ * @see DualThresholdCosmicClusterer
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public class DualThresholdCosmicClusterDriver extends ClusterDriver {
+
+ String inputTightHitCollectionName = "TightEcalCosmicReadoutHits";
+
+ DualThresholdCosmicClusterer cosmicClusterer;
+
+ int minClusterSize;
+ int minRows;
+
+ public DualThresholdCosmicClusterDriver() {
+ }
+
+ public void setMinClusterSize(int minClusterSize) {
+ this.minClusterSize = minClusterSize;
+ }
+
+ public void setMinRows(int minRows) {
+ this.minRows = minRows;
+ }
+
+ public void setInputTightHitCollectionName(String inputTightHitCollectionName) {
+ this.inputTightHitCollectionName = inputTightHitCollectionName;
+ }
+
+ public void initialize() {
+ clusterer.getCuts().setValue("minClusterSize", minClusterSize);
+ clusterer.getCuts().setValue("minRows", minRows);
+ }
+
+ /**
+ * Perform job initialization.
+ */
+ public void startOfData() {
+ if (clusterer == null) {
+ // Setup the Clusterer if it wasn't already initialized by a Driver argument.
+ this.setClustererName("DualThresholdCosmicClusterer");
+ } else {
+ // Does the Clusterer have the right type if there was a custom initialization parameter?
+ if (!(clusterer instanceof DualThresholdCosmicClusterer)) {
+ // The Clusterer does not appear to have the right type for this Driver!
+ throw new IllegalArgumentException("The Clusterer " + this.clusterer.getClass().getCanonicalName() + " does not have the right type.");
+ }
+ }
+
+ // Perform standard start of data initialization from super-class.
+ super.startOfData();
+
+ // Set a reference to the specific type of Clusterer.
+ cosmicClusterer = getClusterer();
+ cosmicClusterer.setInputTightHitCollectionName(inputTightHitCollectionName);
+ }
+
+ public void process(EventHeader event) {
+ if (event.hasCollection(RawTrackerHit.class, inputTightHitCollectionName)) {
+ getLogger().info("tight hit collection " + inputTightHitCollectionName + " has " +
+ event.get(RawTrackerHit.class, inputTightHitCollectionName).size() + " hits");
+ }
+ super.process(event);
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/DualThresholdCosmicClusterer.java Fri Dec 19 16:22:04 2014
@@ -0,0 +1,117 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+import java.util.HashSet;
+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.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.RawTrackerHit;
+
+/**
+ * This is a more complicated version of the {@link SimpleCosmicClusterer}
+ * which uses a list of tight hits for seeding.
+ *
+ * @see SimpleCosmicClusterer
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public class DualThresholdCosmicClusterer extends SimpleCosmicClusterer {
+
+ String inputTightHitCollectionName;
+
+ void setInputTightHitCollectionName(String inputTightHitCollectionName) {
+ this.inputTightHitCollectionName = inputTightHitCollectionName;
+ }
+
+ /**
+ * Create and return clusters.
+ */
+ public List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hitList) {
+
+ if (inputTightHitCollectionName == null) {
+ throw new RuntimeException("The inputTightHitCollectionName was never set.");
+ }
+
+ // Get tight hits collection of RawTrackerHit (raw data).
+ if (!event.hasCollection(RawTrackerHit.class, inputTightHitCollectionName)) {
+ throw new RuntimeException("The hit collection " + inputTightHitCollectionName + " does not exist.");
+ }
+ List<RawTrackerHit> currentTightHits = event.get(RawTrackerHit.class, inputTightHitCollectionName);
+ Set<Long> tightIDs = this.createCellIDSet(currentTightHits);
+
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ // Create map of IDs to hits for convenience.
+ Map<Long, CalorimeterHit> hitMap = ClusterUtilities.createHitMap(hitList);
+
+ // Create list of CalorimeterHits that are clusterable, which is initially all of them.
+ Set<CalorimeterHit> clusterable = new HashSet<CalorimeterHit>();
+ clusterable.addAll(hitList);
+
+ // Loop over all hits in the map.
+ for (CalorimeterHit currentHit : hitList) {
+
+ // Is hit not clusterable (e.g. already used) or not in tight hit list?
+ if (!clusterable.contains(currentHit) || !tightIDs.contains(currentHit.getCellID())) {
+ // Continue to the next hit.
+ continue;
+ }
+
+ // Create list for clustering this hit.
+ List<CalorimeterHit> clusterHits = new ArrayList<CalorimeterHit>();
+
+ // Set of hits whose neighbors have not been checked yet.
+ LinkedList<CalorimeterHit> uncheckedHits = new LinkedList<CalorimeterHit>();
+ uncheckedHits.add(currentHit);
+
+ // While there are still unchecked hits.
+ while (uncheckedHits.size() > 0) {
+
+ // Get the first hit.
+ CalorimeterHit clusterHit = uncheckedHits.removeFirst();
+
+ // Add hit to the cluster.
+ clusterHits.add(clusterHit);
+
+ // Remove the hit from the clusterable list as we have used it in a cluster.
+ clusterable.remove(clusterHit);
+
+ // Loop over the neighbors and to the unchecked list neighboring hits.
+ for (Long neighborHitID : ClusterUtilities.findNeighborHitIDs(ecal, clusterHit, hitMap)) {
+ CalorimeterHit neighborHit = hitMap.get(neighborHitID);
+ if (clusterable.contains(neighborHit)) {
+ uncheckedHits.add(neighborHit);
+ }
+ }
+ }
+
+ // Are there enough hits in the cluster?
+ if (clusterHits.size() >= this.minClusterSize) {
+ // Create cluster and add it to the output list.
+ clusters.add(ClusterUtilities.createBasicCluster(clusterHits));
+ }
+ }
+
+ // Return the cluster list, applying cosmic clustering cuts (from super-class).
+ return applyCuts(clusters);
+ }
+
+ /**
+ * Create a list of cell IDs from a set of RawTrackerHit (raw data) objects.
+ * @param rawHits The input hits.
+ * @return The ID set.
+ */
+ Set<Long> createCellIDSet(List<RawTrackerHit> rawHits) {
+ Set<Long> set = new HashSet<Long>();
+ for (RawTrackerHit hit : rawHits) {
+ set.add(hit.getCellID());
+ }
+ return set;
+ }
+
+}
Modified: 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 (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/LegacyClusterer.java Fri Dec 19 16:22:04 2014
@@ -17,8 +17,6 @@
* <p>
* The clustering algorithm is implemented according to the description in pages 83 and 84 of the
* <a href="https://confluence.slac.stanford.edu/download/attachments/86676777/HPSProposal-FINAL_Rev2.pdf">HPS Proposal document</a>.
- * <p>
- * This is a simple algorithm that is obsolete! The current IC or hardware algorithm clustering algorithms should generally be used instead.
*
* @author Jeremy McCormick <[log in to unmask]>
* @author Tim Nelson <[log in to unmask]>
@@ -33,8 +31,8 @@
}
public void initialize() {
- minClusterSeedEnergy = this.getCut("minClusterSeedEnergy");
- minHitEnergy = this.getCut("minHitEnergy");
+ minClusterSeedEnergy = getCuts().getValue("minClusterSeedEnergy");
+ minHitEnergy = getCuts().getValue("minHitEnergy");
}
public List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hits) {
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NearestNeighborClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NearestNeighborClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NearestNeighborClusterer.java Fri Dec 19 16:22:04 2014
@@ -0,0 +1,91 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+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.event.base.BaseCluster;
+
+/**
+ * This is a simple (example) nearest-neighbor clustering algorithm.
+ *
+ * @author Norman A. Graf
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public class NearestNeighborClusterer extends AbstractClusterer {
+
+ double minHitEnergy = 0;
+ double minSize = 2;
+
+ public NearestNeighborClusterer() {
+ super(new String[] {"minHitEnergy", "minSize" }, new double[] { 0.0, 2.0 });
+ }
+
+ public void initialize() {
+ minHitEnergy = getCuts().getValue("minHitEnergy");
+ minSize = getCuts().getValue("minSize");
+ }
+
+ @Override
+ public List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hits) {
+ Map<Long, CalorimeterHit> hitMap = ClusterUtilities.createHitMap(hits);
+ // New Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ while (!hitMap.isEmpty()) {
+ Long k = hitMap.keySet().iterator().next();
+ CalorimeterHit hit = hitMap.get(k);
+ BaseCluster cluster = this.createCluster(hitMap, hit);
+ //done with this cluster, let's compute some shape properties
+ if (cluster.getSize() >= minSize) {
+ cluster.calculateProperties();
+ clusters.add(cluster);
+ }
+ }
+ return clusters;
+ }
+
+ /**
+ *
+ * @param map A Map of CalorimeterHit keyed on CellID, Hits added to this
+ * cluster are removed from the map
+ * @param hit A CalorimeterHit representing a hit crystal
+ * @param threshold The energy threshold below which a hit crystal should
+ * NOT be added to the cluster
+ * @param neighborMap The map of HPS crystal neighbors.
+ */
+ BaseCluster createCluster(Map<Long, CalorimeterHit> map, CalorimeterHit hit) {
+ BaseCluster cluster = new BaseCluster();
+ // start by adding this hit to the cluster
+ cluster.addHit(hit);
+ // remove this hit from the map so it can't be used again
+ map.remove(hit.getCellID());
+ List<CalorimeterHit> hits = cluster.getCalorimeterHits();
+ // loop over the hits in the cluster and add all its neighbors
+ // note that hits.size() grows as we add cells, so we recursively find neighbors of neighbors
+ for (int i = 0; i < hits.size(); ++i) {
+ CalorimeterHit c = hits.get(i);
+ // Get neighbor crystal IDs.
+ Set<Long> neighbors = neighborMap.get(c.getCellID());
+ // loop over all neighboring cell Ids
+ for (Long neighborId : neighbors) {
+ // Find the neighbor hit in the event if it exists.
+ CalorimeterHit neighborHit = map.get(neighborId);
+ // Was this neighbor cell hit?
+ if (neighborHit != null) {
+ // if so, does it meet or exceed threshold?
+ if (neighborHit.getRawEnergy() >= minHitEnergy) {
+ cluster.addHit(neighborHit);
+ }
+ //remove this hit from the map so it can't be used again
+ map.remove(neighborHit.getCellID());
+ }
+ }
+ }
+ return cluster;
+ }
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCuts.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCuts.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCuts.java Fri Dec 19 16:22:04 2014
@@ -0,0 +1,75 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.Map;
+
+/**
+ * This is an interface for accessing numerical cut values
+ * in a clustering algorithm by index or name.
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public interface NumericalCuts {
+
+ /**
+ * Set all the cut values from an array.
+ * @param values The cut values.
+ */
+ void setValues(double[] values);
+
+ /**
+ * Get the cut values.
+ * @return The cut values array.
+ */
+ double[] getValues();
+
+ /**
+ * Get a cut setting by name.
+ * @param name The name of the cut.
+ * @return The cut value.
+ */
+ double getValue(String name);
+
+ /**
+ * Get a cut value by index.
+ * @param index The index of the cut.
+ * @return The cut value from the index.
+ */
+ double getValue(int index);
+
+ /**
+ * Get the names of the cuts.
+ * @return The names of the cuts.
+ */
+ String[] getNames();
+
+ /**
+ * Set the value of a cut by index.
+ * @param index The index of the cut.
+ * @param value The value of the cut.
+ */
+ void setValue(int index, double value);
+
+ /**
+ * Set a cut value by name.
+ * @param name The name of the cut.
+ * @param value The value of the cut.
+ */
+ void setValue(String name, double value);
+
+ /**
+ * True if using the default cuts.
+ * @return True if using default cuts.
+ */
+ boolean isDefaultValues();
+
+ /**
+ * Get the default cuts.
+ * @return The default cuts.
+ */
+ double[] getDefaultValues();
+
+ /**
+ * Set the cut values from a map of keys to values.
+ * @param valueMap The cut values as a map.
+ */
+ void setValues(Map<String, Double> valueMap);
+}
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCutsImpl.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCutsImpl.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/NumericalCutsImpl.java Fri Dec 19 16:22:04 2014
@@ -0,0 +1,101 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.Map;
+import java.util.Map.Entry;
+
+/**
+ * This is the basic implementation of the {@link NumericalCuts} interface.
+ *
+ * @see NumericalCuts
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ */
+public class NumericalCutsImpl implements NumericalCuts {
+
+ String[] names;
+ double[] values;
+ double[] defaultValues;
+
+ public NumericalCutsImpl(String[] names, double[] defaultValues) {
+ this.names = names;
+ this.defaultValues = defaultValues;
+ this.values = this.defaultValues;
+ }
+
+ @Override
+ public void setValues(double[] values) {
+ if (values.length != this.names.length) {
+ throw new IllegalArgumentException("The values array has the wrong length: " + values.length);
+ }
+ this.values = values;
+ }
+
+ @Override
+ public double[] getValues() {
+ return values;
+ }
+
+ @Override
+ public double getValue(String name) {
+ int index = indexFromName(name);
+ if (index == -1) {
+ throw new IllegalArgumentException("There is no cut called " + name + " defined by this clusterer.");
+ }
+ return getValue(index);
+ }
+
+ @Override
+ public double getValue(int index) {
+ if (index > values.length || index < 0) {
+ throw new IndexOutOfBoundsException("The index " + index + " is out of bounds for cuts array.");
+ }
+ return values[index];
+ }
+
+ @Override
+ public String[] getNames() {
+ return names;
+ }
+
+ @Override
+ public void setValue(int index, double value) {
+ values[index] = value;
+ }
+
+ @Override
+ public boolean isDefaultValues() {
+ return values == defaultValues;
+ }
+
+ @Override
+ public double[] getDefaultValues() {
+ return defaultValues;
+ }
+
+ @Override
+ public void setValue(String name, double value) {
+ int index = indexFromName(name);
+ values[index] = value;
+ }
+
+ @Override
+ public void setValues(Map<String, Double> valueMap) {
+ for (Entry<String, Double> entry : valueMap.entrySet()) {
+ this.setValue(entry.getKey(), entry.getValue());
+ }
+ }
+
+ /**
+ * Get the index of a cut from its name.
+ * @param name The name of the cut.
+ * @return The index of the cut from the name.
+ */
+ protected int indexFromName(String name) {
+ for (int index = 0; index < values.length; index++) {
+ if (getNames()[index] == name) {
+ return index;
+ }
+ }
+ return -1;
+ }
+}
Modified: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleClasInnerCalClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleClasInnerCalClusterer.java (original)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleClasInnerCalClusterer.java Fri Dec 19 16:22:04 2014
@@ -37,15 +37,15 @@
* Initialize the algorithm with default cuts.
*/
public SimpleClasInnerCalClusterer() {
- super(new String[] { "minEnergy", "minTime", "timeWindow", "timeCut" }, new double[] { 0., 0.001, 0.0, 20.0 });
+ super(new String[] { "minEnergy", "minTime", "timeWindow", "timeCut" }, new double[] { 0.001, 0.0, 20.0, 0. });
}
public void initialize() {
// Setup class variables from cuts.
- timeCut = (this.getCut("timeCut") == 1.0);
- minEnergy = this.getCut("minEnergy");
- minTime = this.getCut("minTime");
- timeWindow = this.getCut("timeWindow");
+ timeCut = (getCuts().getValue("timeCut") == 1.0);
+ minEnergy = getCuts().getValue("minEnergy");
+ minTime = getCuts().getValue("minTime");
+ timeWindow = getCuts().getValue("timeWindow");
}
public List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hitCollection) {
Added: java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleCosmicClusterer.java
=============================================================================
--- java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleCosmicClusterer.java (added)
+++ java/trunk/ecal-recon/src/main/java/org/hps/recon/ecal/cluster/SimpleCosmicClusterer.java Fri Dec 19 16:22:04 2014
@@ -0,0 +1,159 @@
+package org.hps.recon.ecal.cluster;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.LinkedList;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import java.util.Map.Entry;
+
+import org.lcsim.detector.identifier.Identifier;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+
+/**
+ * <p>
+ * This algorithm clusters an input list of CalorimeterHits using an iterative nearest neighbor algorithm.
+ * <p>
+ * There is a set of cuts applied on the initial cluster list that includes the following defaults:
+ * <ul>
+ * <li>must have at least 3 hits total in the cluster</li>
+ * <li>must have at least three rows of crystals with hits</li>
+ * <li>must have no more than 2 hits in each row of crystals</li>
+ * </ul>
+ * <p>
+ * This algorithm does not cluster across the beam gap. Separate clusters are made for the top
+ * and bottom sets of crystals.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Tim "THammer" Nelson <[log in to unmask]>
+ */
+public class SimpleCosmicClusterer extends AbstractClusterer {
+
+ protected int minClusterSize;
+ protected int minRows;
+ protected int maxHitsPerRow;
+
+ /**
+ * Default constructor that sets cut names and default values.
+ */
+ public SimpleCosmicClusterer() {
+ super(new String[] { "minClusterSize", "minRows", "maxHitsPerRow" }, new double[] { 3, 3, 2 });
+ }
+
+ /**
+ * Initialize some instance variables from cut settings.
+ */
+ public void initialize() {
+ this.minClusterSize = (int) getCuts().getValue("minClusterSize");
+ this.minRows = (int) getCuts().getValue("minRows");
+ this.maxHitsPerRow = (int) getCuts().getValue("maxHitsPerRow");
+ }
+
+ /**
+ * This method implements the cosmic clustering algorithm which uses a iterative
+ * nearest neighbor algorithm. It uses each hit as a seed and tries to then
+ * extend the cluster by adding neighboring hits. Clustered hits are added to a list
+ * which is checked so that they are not reused.
+ * @param hitList The input hit list of hits.
+ * @param event The current LCSim event.
+ * @return A list of calorimeter hits that can be turned into clusters.
+ */
+ public List<Cluster> createClusters(EventHeader event, List<CalorimeterHit> hitList) {
+
+ // Create empty list of clusters which are just lists of hits.
+ //List<List<CalorimeterHit>> clusterList = new ArrayList<List<CalorimeterHit>>();
+ List<Cluster> clusterList = new ArrayList<Cluster>();
+
+ // Create map of IDs to hits for convenience.
+ Map<Long, CalorimeterHit> hitMap = ClusterUtilities.createHitMap(hitList);
+
+ // Create list of hits that are clusterable, which is initially all of them.
+ Set<CalorimeterHit> clusterable = new HashSet<CalorimeterHit>();
+ clusterable.addAll(hitList);
+
+ // Loop over all hits in the map.
+ for (CalorimeterHit hit : hitList) {
+
+ // Is hit clusterable?
+ if (!clusterable.contains(hit)) {
+ // Continue to next hit.
+ continue;
+ }
+
+ // Create list for clustering this hit.
+ List<CalorimeterHit> clusterHits = new ArrayList<CalorimeterHit>();
+
+ // Set of hits whose neighbors have not been checked yet.
+ LinkedList<CalorimeterHit> uncheckedHits = new LinkedList<CalorimeterHit>();
+ uncheckedHits.add(hit);
+
+ // While there are still unchecked hits in the cluster.
+ while (uncheckedHits.size() > 0) {
+
+ // Get the first hit and add it to the cluster.
+ CalorimeterHit clusterHit = uncheckedHits.removeFirst();
+
+ // Add hit to the cluster.
+ clusterHits.add(clusterHit);
+
+ // Remove the hit from the clusterable list.
+ clusterable.remove(clusterHit);
+
+ // Loop over the neighbors and add IDs with hits to the unchecked list.
+ for (Long neighborHitID : ClusterUtilities.findNeighborHitIDs(ecal, clusterHit, hitMap)) {
+ CalorimeterHit neighborHit = hitMap.get(neighborHitID);
+ if (clusterable.contains(neighborHit)) {
+ uncheckedHits.add(neighborHit);
+ }
+ }
+ }
+
+ if (clusterHits.size() >= this.minClusterSize) {
+ Cluster cluster = ClusterUtilities.createBasicCluster(clusterHits);
+ clusterList.add(cluster);
+ }
+ }
+
+ // Apply cuts on the cluster list that was generated.
+ List<Cluster> selectedClusters = applyCuts(clusterList);
+
+ return selectedClusters;
+ }
+
+ /**
+ * This method takes a list of potential cluster hits and applies selection cuts,
+ * returning a new list that has the hit lists which did not pass the cuts removed.
+ * @param clusteredHitLists The input hit lists.
+ * @return The hit lists that passed the cuts.
+ */
+ protected List<Cluster> applyCuts(List<Cluster> clusterList) {
+ List<Cluster> selectedClusters = new ArrayList<Cluster>();
+ for (Cluster cluster : clusterList) {
+ Map<Integer, Set<CalorimeterHit>> rowMap = new HashMap<Integer, Set<CalorimeterHit>>();
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ int row = ecal.getDetectorElement().getIdentifierHelper().getValue(new Identifier(hit.getCellID()), "iy");
+ if (rowMap.get(row) == null) {
+ rowMap.put(row, new HashSet<CalorimeterHit>());
+ }
+ rowMap.get(row).add(hit);
+ }
+ if (rowMap.size() >= this.minRows) {
+ boolean okay = true;
+ rowMapLoop: for (Entry<Integer, Set<CalorimeterHit>> entries : rowMap.entrySet()) {
+ if (entries.getValue().size() > this.maxHitsPerRow) {
+ okay = false;
+ break rowMapLoop;
+ }
+ }
+ if (okay) {
+ selectedClusters.add(cluster);
+ }
+ }
+ }
+ return selectedClusters;
+ }
+}
|
