lcsim/src/org/lcsim/recon/cluster/mipfinder
diff -N MipFinderCrossingBarrelEndcapBorder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MipFinderCrossingBarrelEndcapBorder.java 5 Jun 2008 17:01:36 -0000 1.1
@@ -0,0 +1,408 @@
+package org.lcsim.recon.cluster.mipfinder;
+
+import java.util.*;
+import java.io.IOException;
+
+import hep.physics.vec.*;
+
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.util.decision.*;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.recon.cluster.util.ClusterSizeDecision;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.recon.cluster.util.Clusterer;
+import org.lcsim.recon.cluster.mipfinder.*;
+
+/**
+ * A version of the MIP finding algorithm designed to be
+ * able to move from the barrel to the endcap (or vice versa).
+ *
+ * @author [log in to unmask]
+ * @version $Id: MipFinderCrossingBarrelEndcapBorder.java,v 1.1 2008/06/05 17:01:36 mcharles Exp $
+ */
+
+public class MipFinderCrossingBarrelEndcapBorder extends TrackClusterDriver
+{
+ /** Constructor, suitable for stand-alone use. */
+ public MipFinderCrossingBarrelEndcapBorder(Collection<CalorimeterHit> barrelHits, Collection<CalorimeterHit> barrelHitsNearBoundary,
+ Collection<CalorimeterHit> endcapHits, Collection<CalorimeterHit> endcapHitsNearBoundary)
+ {
+ super();
+ setInputBarrelHits(barrelHits, barrelHitsNearBoundary);
+ setInputEndcapHits(endcapHits, endcapHitsNearBoundary);
+ }
+
+ /** Constructor, suitable for Driver use. */
+ public MipFinderCrossingBarrelEndcapBorder(String inputHitMapName, String outputClusterListName, String outputHitMapName,
+ String inputBarrelHitList, String inputBarrelHitListNearBoundary,
+ String inputEndcapHitList, String inputEndcapHitListNearBoundary)
+ {
+ super(inputHitMapName, outputClusterListName, outputHitMapName);
+ m_inputBarrelHitListName = inputBarrelHitList;
+ m_inputEndcapHitListName = inputEndcapHitList;
+ m_inputBarrelHitListNearBoundaryName = inputBarrelHitListNearBoundary;
+ m_inputEndcapHitListNearBoundaryName = inputEndcapHitListNearBoundary;
+ }
+ /** Process routine for Driver use.*/
+ public void process(EventHeader event) {
+ List<CalorimeterHit> barrelHits = event.get(CalorimeterHit.class, m_inputBarrelHitListName);
+ List<CalorimeterHit> endcapHits = event.get(CalorimeterHit.class, m_inputEndcapHitListName);
+ List<CalorimeterHit> barrelHitsNearBoundary = event.get(CalorimeterHit.class, m_inputBarrelHitListNearBoundaryName);
+ List<CalorimeterHit> endcapHitsNearBoundary = event.get(CalorimeterHit.class, m_inputEndcapHitListNearBoundaryName);
+ setInputBarrelHits(barrelHits, barrelHitsNearBoundary);
+ setInputEndcapHits(endcapHits, endcapHitsNearBoundary);
+ super.process(event);
+ }
+
+ String m_inputBarrelHitListName;
+ String m_inputBarrelHitListNearBoundaryName;
+ String m_inputEndcapHitListName;
+ String m_inputEndcapHitListNearBoundaryName;
+ HitMap m_barrelHits;
+ HitMap m_endcapHits;
+ HitMap m_barrelHitsNearBoundary;
+ HitMap m_endcapHitsNearBoundary;
+
+ public void setInputBarrelHits(Collection<CalorimeterHit> barrelHits, Collection<CalorimeterHit> barrelHitsNearBoundary) {
+ m_barrelHits = new HitMap(barrelHits);
+ m_barrelHitsNearBoundary = new HitMap(barrelHitsNearBoundary);
+ }
+
+ public void setInputEndcapHits(Collection<CalorimeterHit> endcapHits, Collection<CalorimeterHit> endcapHitsNearBoundary) {
+ m_endcapHits = new HitMap(endcapHits);
+ m_endcapHitsNearBoundary = new HitMap(endcapHitsNearBoundary);
+ }
+
+ /** Take a collection which should contain both barrel and endcap hits, and attempt to find MIPs. */
+ public List<Cluster> createClusters(List<CalorimeterHit> hits) {
+ // Split into barrel & endcap collections
+ List<CalorimeterHit> hitsBarrel = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> hitsEndcap = new Vector<CalorimeterHit>();
+ Set<Long> keySetBarrel = m_barrelHits.keySet();
+ Set<Long> keySetEndcap = m_endcapHits.keySet();
+ Set<Long> keySetBarrelHitsNearBoundary = m_barrelHitsNearBoundary.keySet();
+ Set<Long> keySetEndcapHitsNearBoundary = m_endcapHitsNearBoundary.keySet();
+ for (CalorimeterHit hit : hits) {
+ long id = hit.getCellID();
+ boolean inBarrel = keySetBarrel.contains(id);
+ boolean inEndcap = keySetEndcap.contains(id);
+ if (inBarrel && inEndcap) { throw new AssertionError("Book-keeping error"); }
+ if (!inBarrel && !inEndcap) { throw new AssertionError("Book-keeping error"); }
+ if (inBarrel) {
+ hitsBarrel.add(hit);
+ } else {
+ hitsEndcap.add(hit);
+ }
+ }
+
+ // Setup
+ List<AbstractHitType> hitTypes = defineHitTypes();
+ int nLayersBarrel = countLayersInCalorimeter(hitsBarrel);
+ int nLayersEndcap = countLayersInCalorimeter(hitsEndcap);
+ Map<Integer, List<CalorimeterHit>> hitsPerLayerBarrel = sortHitsByLayer(hitsBarrel);
+ Map<Integer, List<CalorimeterHit>> hitsPerLayerEndcap = sortHitsByLayer(hitsEndcap);
+ List<Cluster> outputClusterList = new Vector<Cluster>();
+
+ MIPClusterBuilder clusterBuilderEndcap = new MIPClusterBuilder(hitsEndcap, hitTypes);
+ List<Cluster> trackSegmentsEndcap = subBuildClusters(nLayersEndcap, clusterBuilderEndcap, hitsPerLayerEndcap);
+ MIPClusterBuilder clusterBuilderBarrel = new MIPClusterBuilder(hitsBarrel, hitTypes);
+ List<Cluster> trackSegmentsBarrel = subBuildClusters(nLayersBarrel, clusterBuilderBarrel, hitsPerLayerBarrel);
+
+ // Now, look for merges we should make. Specifically, look for this case:
+ // * MIP segment in one component stops in border region
+ // * MIP segment in other component stops in border region
+ // * The two segments stop within a short distance of one another
+ // * The two segments have a compatible local direction / point to one another
+
+ // First, find the ones that stop in the border region:
+ List<Cluster> trackSegmentsEndcapStoppingNearBoundary = findStubsThatStopInBoundary(trackSegmentsEndcap, keySetEndcapHitsNearBoundary);
+ List<Cluster> trackSegmentsBarrelStoppingNearBoundary = findStubsThatStopInBoundary(trackSegmentsBarrel, keySetBarrelHitsNearBoundary);
+
+ // Next, look at combinations and see if they can be joined:
+ class Merge extends ArrayList<Cluster> { public Merge() { super(); } }
+ Map<Cluster, List<Merge>> mapStubToMergeCands = new HashMap<Cluster, List<Merge>>();
+ Map<Merge, Double> mapMergeCandToScore = new HashMap<Merge, Double>();
+ for (Cluster stubEndcap : trackSegmentsEndcapStoppingNearBoundary) {
+ // Find the innermost & outermost layers
+ List<CalorimeterHit> hitsInInnermostLayerEndcap = findHitsByLayerFromInnermost(stubEndcap, 0);
+ List<CalorimeterHit> hitsInOutermostLayerEndcap = findHitsByLayerFromOutermost(stubEndcap, 0);
+ for (Cluster stubBarrel : trackSegmentsBarrelStoppingNearBoundary) {
+ // Find the innermost & outermost layers
+ List<CalorimeterHit> hitsInInnermostLayerBarrel = findHitsByLayerFromInnermost(stubBarrel, 0);
+ List<CalorimeterHit> hitsInOutermostLayerBarrel = findHitsByLayerFromOutermost(stubBarrel, 0);
+ // Now see if they're compatible...
+ Hep3Vector minDisplacement1 = minSeparation(hitsInInnermostLayerEndcap, hitsInInnermostLayerBarrel);
+ Hep3Vector minDisplacement2 = minSeparation(hitsInInnermostLayerEndcap, hitsInOutermostLayerBarrel);
+ Hep3Vector minDisplacement3 = minSeparation(hitsInOutermostLayerEndcap, hitsInInnermostLayerBarrel);
+ Hep3Vector minDisplacement4 = minSeparation(hitsInOutermostLayerEndcap, hitsInOutermostLayerBarrel);
+ double minDist1 = minDisplacement1.magnitude();
+ double minDist2 = minDisplacement2.magnitude();
+ double minDist3 = minDisplacement3.magnitude();
+ double minDist4 = minDisplacement4.magnitude();
+ // Figure out which is the best combination (e.g. barrel outer and endcap inner, ...)
+ boolean bestEndcapInner, bestBarrelInner;
+ double minDist = -1.0;
+ double distCutOff = 30.0; // 3cm -- FIXME: Don't hard-code
+ if (minDist1 <= minDist2 && minDist1 <= minDist3 && minDist1 <= minDist4) {
+ bestEndcapInner = true; bestBarrelInner = true;
+ minDist = minDist1;
+ } else if (minDist2 <= minDist1 && minDist2 <= minDist3 && minDist3 <= minDist4) {
+ bestEndcapInner = true; bestBarrelInner = false;
+ minDist = minDist2;
+ } else if (minDist3 <= minDist1 && minDist3 <= minDist2 && minDist3 <= minDist4) {
+ bestEndcapInner = false; bestBarrelInner = true;
+ minDist = minDist3;
+ } else if (minDist4 <= minDist1 && minDist4 <= minDist2 && minDist4 <= minDist3) {
+ bestEndcapInner = false; bestBarrelInner = false;
+ minDist = minDist4;
+ } else {
+ throw new AssertionError("Unhandled case");
+ }
+ boolean passesDistanceCut = (minDist < distCutOff);
+ System.out.println("DEBUG: Found candidate for connecting stubs! Distances are "+minDist1+", "+minDist2+", "+minDist3+", "+minDist4);
+ if (passesDistanceCut) {
+ double debugBestDist = 0.0;
+ if (bestEndcapInner) {
+ if (bestBarrelInner) { debugBestDist = minDist1; } else { debugBestDist = minDist2; }
+ } else {
+ if (bestBarrelInner) { debugBestDist = minDist3; } else { debugBestDist = minDist4; }
+ }
+ System.out.println("DEBUG: Passes distance cut! bestEndcapInner="+bestEndcapInner+", bestBarrelInner="+bestBarrelInner+" --> bestDist="+debugBestDist);
+ // Check direction.
+ // For this, we want the hits in the previous layer too...
+ List<CalorimeterHit> hitsInBarrelLayerMinusZero;
+ List<CalorimeterHit> hitsInBarrelLayerMinusOne;
+ List<CalorimeterHit> hitsInEndcapLayerMinusZero;
+ List<CalorimeterHit> hitsInEndcapLayerMinusOne;
+ List<CalorimeterHit> hitsInEndcapLayerMinusTwo;
+ if (bestEndcapInner) {
+ hitsInEndcapLayerMinusZero = hitsInInnermostLayerEndcap;
+ hitsInEndcapLayerMinusOne = findHitsByLayerFromInnermost(stubEndcap, 1);
+ hitsInEndcapLayerMinusTwo = findHitsByLayerFromInnermost(stubEndcap, 2);
+ } else {
+ hitsInEndcapLayerMinusZero = hitsInOutermostLayerEndcap;
+ hitsInEndcapLayerMinusOne = findHitsByLayerFromOutermost(stubEndcap, 1);
+ hitsInEndcapLayerMinusTwo = findHitsByLayerFromOutermost(stubEndcap, 2);
+ }
+ if (bestBarrelInner) {
+ hitsInBarrelLayerMinusZero = hitsInInnermostLayerBarrel;
+ hitsInBarrelLayerMinusOne = findHitsByLayerFromInnermost(stubBarrel, 1);
+ } else {
+ hitsInBarrelLayerMinusZero = hitsInOutermostLayerBarrel;
+ hitsInBarrelLayerMinusOne = findHitsByLayerFromOutermost(stubBarrel, 1);
+ }
+ Hep3Vector directionEndcap = findTangentDirection(hitsInEndcapLayerMinusZero, hitsInEndcapLayerMinusOne);
+ Hep3Vector directionBarrel = findTangentDirection(hitsInBarrelLayerMinusZero, hitsInBarrelLayerMinusOne);
+ double dotProduct = VecOp.dot(directionEndcap, directionBarrel);
+ System.out.println("DEBUG: Dot product of tangents = "+dotProduct);
+ double dotProductCut = 0.7; // FIXME: Don't hard-code
+ if (Math.abs(dotProduct) > 0.7) {
+ Merge mergeCand = new Merge();
+ double score = Math.abs(dotProduct);
+ mergeCand.add(stubEndcap);
+ mergeCand.add(stubBarrel);
+ List<Merge> mergeCandsEndcap = mapStubToMergeCands.get(stubEndcap);
+ List<Merge> mergeCandsBarrel = mapStubToMergeCands.get(stubBarrel);
+ if (mergeCandsEndcap == null) { mergeCandsEndcap = new Vector<Merge>(); mapStubToMergeCands.put(stubEndcap, mergeCandsEndcap); }
+ if (mergeCandsBarrel == null) { mergeCandsBarrel = new Vector<Merge>(); mapStubToMergeCands.put(stubBarrel, mergeCandsBarrel); }
+ mergeCandsEndcap.add(mergeCand);
+ mergeCandsBarrel.add(mergeCand);
+ mapMergeCandToScore.put(mergeCand, minDist);
+ }
+ }
+ }
+ }
+
+ // Check up on merges, looking for ambiguity.
+ // In case of ambiguity, use only best merge (or none at all).
+ Map<Cluster,Cluster> bestMergeMapEndcapToBarrel = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> bestMergeMapBarrelToEndcap = new HashMap<Cluster,Cluster>();
+ for (Cluster clus : trackSegmentsEndcapStoppingNearBoundary) {
+ List<Merge> mergeCands = mapStubToMergeCands.get(clus);
+ Merge bestMerge = null;
+ if (mergeCands != null) {
+ for (Merge cand : mergeCands) {
+ double score = mapMergeCandToScore.get(cand);
+ if (bestMerge == null || score > mapMergeCandToScore.get(bestMerge)) {
+ bestMerge = cand;
+ }
+ }
+ }
+ if (bestMerge != null) {
+ Cluster endcapClus = bestMerge.get(0);
+ Cluster barrelClus = bestMerge.get(1);
+ if (endcapClus != clus) { throw new AssertionError("Book-keeping error"); }
+ if (barrelClus == clus) { throw new AssertionError("Book-keeping error"); }
+ bestMergeMapEndcapToBarrel.put(endcapClus, barrelClus);
+ }
+ }
+ for (Cluster clus : trackSegmentsBarrelStoppingNearBoundary) {
+ List<Merge> mergeCands = mapStubToMergeCands.get(clus);
+ Merge bestMerge = null;
+ if (mergeCands != null) {
+ for (Merge cand : mergeCands) {
+ double score = mapMergeCandToScore.get(cand);
+ if (bestMerge == null || score > mapMergeCandToScore.get(bestMerge)) {
+ bestMerge = cand;
+ }
+ }
+ }
+ if (bestMerge != null) {
+ Cluster endcapClus = bestMerge.get(0);
+ Cluster barrelClus = bestMerge.get(1);
+ if (endcapClus == clus) { throw new AssertionError("Book-keeping error"); }
+ if (barrelClus != clus) { throw new AssertionError("Book-keeping error"); }
+ bestMergeMapBarrelToEndcap.put(barrelClus, endcapClus);
+ }
+ }
+
+ // Make output, including merges
+ List<Cluster> outputList = new Vector<Cluster>();
+ outputList.addAll(trackSegmentsBarrel);
+ outputList.addAll(trackSegmentsEndcap);
+ for (Cluster endcapClus : bestMergeMapEndcapToBarrel.keySet()) {
+ Cluster barrelClus = bestMergeMapEndcapToBarrel.get(endcapClus);
+ Cluster targetOfBarrelClus = bestMergeMapBarrelToEndcap.get(barrelClus);
+ if (targetOfBarrelClus == endcapClus) {
+ // All self-consistent
+ System.out.println("Will do a merge of endcap MIP["+endcapClus.getCalorimeterHits().size()+"] -> barrel MIP["+barrelClus.getCalorimeterHits().size()+"]");
+ outputList.remove(endcapClus);
+ outputList.remove(barrelClus);
+ MIPCluster newClus = new MIPCluster(1,1);
+ for (CalorimeterHit hit : endcapClus.getCalorimeterHits()) { newClus.addHit(hit); }
+ for (CalorimeterHit hit : barrelClus.getCalorimeterHits()) { newClus.addHit(hit); }
+ outputList.add(newClus);
+ }
+ }
+
+ return outputList;
+ }
+
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+ protected List<Cluster> subBuildClusters(int nLayers, MIPClusterBuilder clusterBuilder, Map<Integer, List<CalorimeterHit>> hitsPerLayer) {
+ clusterBuilder.setNumberOfSeedLayers(1);
+ clusterBuilder.initialize();
+ List<Cluster> outputList = new Vector<Cluster>();
+ for (int iLayer=0; iLayer<nLayers; iLayer++) {
+ for (int iDirection=0; iDirection<2; iDirection++) {
+ if (iDirection==0) {
+ // Forward
+ clusterBuilder.setDirectionAndFirstLayer(+1, iLayer);
+ } else {
+ // Backward
+ clusterBuilder.setDirectionAndFirstLayer(-1, nLayers-iLayer-1);
+ }
+ List<CalorimeterHit> nuclei = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> hitsInLayer = hitsPerLayer.get(new Integer(iLayer));
+ if (hitsInLayer != null) {
+ nuclei.addAll(hitsInLayer);
+ }
+ clusterBuilder.provideNucleii(nuclei);
+ // Build track segments within the endcap.
+ List<Cluster> trackSegments = clusterBuilder.getMIPClusterList();
+ outputList.addAll(trackSegments);
+ }
+ }
+ return outputList;
+ }
+
+ protected Map<Integer, List<CalorimeterHit>> sortHitsByLayer(List<CalorimeterHit> hits) {
+ Map<Integer, List<CalorimeterHit>> output = new HashMap<Integer, List<CalorimeterHit>>();
+ for (CalorimeterHit hit : hits) {
+ int layer = getLayer(hit);
+ List<CalorimeterHit> hitsInLayer = output.get(layer);
+ if (hitsInLayer == null) { hitsInLayer = new Vector<CalorimeterHit>(); output.put(layer, hitsInLayer); }
+ hitsInLayer.add(hit);
+ }
+ return output;
+ }
+
+ protected List<CalorimeterHit> findHitsByLayerFromInnermost(Cluster clus, int offset) {
+ Map<Integer, List<CalorimeterHit>> sortedHits = sortHitsByLayer(clus.getCalorimeterHits());
+ Integer minLayer = Collections.min(sortedHits.keySet());
+ int targetLayer = minLayer + offset;
+ return sortedHits.get(targetLayer);
+ }
+ protected List<CalorimeterHit> findHitsByLayerFromOutermost(Cluster clus, int offset) {
+ Map<Integer, List<CalorimeterHit>> sortedHits = sortHitsByLayer(clus.getCalorimeterHits());
+ Integer maxLayer = Collections.max(sortedHits.keySet());
+ int targetLayer = maxLayer - offset;
+ return sortedHits.get(targetLayer);
+ }
+
+ protected List<Cluster> findStubsThatStopInBoundary(List<Cluster> trackSegments, Set<Long> keySetHitsNearBoundary) {
+ List<Cluster> trackSegmentsStoppingNearBoundary = new Vector<Cluster>();
+ for (Cluster clus : trackSegments) {
+ // Find the innermost & outermost layers
+ List<CalorimeterHit> hitsInInnermostLayer = findHitsByLayerFromInnermost(clus, 0);
+ List<CalorimeterHit> hitsInOutermostLayer = findHitsByLayerFromOutermost(clus, 0);
+ if (hitsInInnermostLayer.size()<1) { throw new AssertionError("Book-keeping error"); }
+ if (hitsInOutermostLayer.size()<1) { throw new AssertionError("Book-keeping error"); }
+ // Check whether any of the hits in the inner/outermost layers are in the boundary region
+ List<CalorimeterHit> boundaryHitsInInnermostLayer = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> boundaryHitsInOutermostLayer = new Vector<CalorimeterHit>();
+ for (CalorimeterHit hit : hitsInInnermostLayer) {
+ if (keySetHitsNearBoundary.contains(hit.getCellID())) {
+ boundaryHitsInInnermostLayer.add(hit);
+ }
+ }
+ for (CalorimeterHit hit : hitsInOutermostLayer) {
+ if (keySetHitsNearBoundary.contains(hit.getCellID())) {
+ boundaryHitsInOutermostLayer.add(hit);
+ }
+ }
+ if (boundaryHitsInInnermostLayer.size()>0 || boundaryHitsInOutermostLayer.size()>0) {
+ trackSegmentsStoppingNearBoundary.add(clus);
+ }
+ }
+ return trackSegmentsStoppingNearBoundary;
+ }
+
+ Hep3Vector minSeparation(List<CalorimeterHit> list1, List<CalorimeterHit> list2) {
+ Hep3Vector output = null;
+ double minDist = 0.0;
+ boolean firstPass = true;
+ for (CalorimeterHit hit1 : list1) {
+ Hep3Vector pos1 = new BasicHep3Vector(hit1.getPosition());
+ for (CalorimeterHit hit2 : list2) {
+ Hep3Vector pos2 = new BasicHep3Vector(hit2.getPosition());
+ Hep3Vector displacement = VecOp.sub(pos1,pos2);
+ double dist = displacement.magnitude();
+ if (firstPass || dist < minDist) {
+ minDist = dist;
+ output = displacement;
+ firstPass = false;
+ }
+ }
+ }
+ if (firstPass) { throw new AssertionError("Empty list or book-keeping error"); }
+ return output;
+ }
+
+ Hep3Vector getMeanPoint(List<CalorimeterHit> hits) {
+ Hep3Vector sum = new BasicHep3Vector(0.0, 0.0, 0.0);
+ for (CalorimeterHit hit : hits) {
+ sum = VecOp.add(sum, new BasicHep3Vector(hit.getPosition()));
+ }
+ double scaleFactor = 1.0 / ((double)(hits.size()));
+ Hep3Vector mean = VecOp.mult(scaleFactor, sum);
+ return mean;
+ }
+
+ Hep3Vector findTangentDirection(List<CalorimeterHit> list1, List<CalorimeterHit> list2) {
+ Hep3Vector mean1 = getMeanPoint(list1);
+ Hep3Vector mean2 = getMeanPoint(list2);
+ Hep3Vector displacement = VecOp.sub(mean2, mean1);
+ return VecOp.unit(displacement);
+ }
+}