lcsim/src/org/lcsim/recon/cluster/mipfinder
diff -N FlexibleMIPFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FlexibleMIPFinder.java 18 Oct 2007 03:08:41 -0000 1.1
@@ -0,0 +1,343 @@
+package org.lcsim.recon.cluster.mipfinder;
+
+import hep.physics.vec.*;
+
+import java.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.*;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+/**
+ * Find MIP-like segments
+ *
+ * Allow for curved MIPs and movement at an angle of > 45 degrees
+ * by not requiring one hit per layer.
+ *
+ * Instead, require that most hits be locally isolated,
+ * i.e. only 1 or 2 other hits in a 3x3x3 grid centered on the hit.
+ * Allow instances of 3 other hits to take into account delta rays etc.
+ * Also allow for hits that are have nothing else within a 3x3x3 grid
+ * and have 1 or 2 other hits within a 5x5x3 grid
+ *
+ * @author Mat Charles
+ * @version $Id: FlexibleMIPFinder.java,v 1.1 2007/10/18 03:08:41 mcharles Exp $
+ */
+
+public class FlexibleMIPFinder extends Driver
+{
+ public FlexibleMIPFinder(String inputHitMapName, String outputClusterListName, String outputHitMapName) {
+ m_inputHitMapName = inputHitMapName;
+ m_outputClusterListName = outputClusterListName;
+ m_outputHitMapName = outputHitMapName;
+ }
+
+ String m_inputHitMapName;
+ String m_outputClusterListName;
+ String m_outputHitMapName;
+
+ /**
+ * Find track segments in one event.
+ */
+ public void process(EventHeader event)
+ {
+ // Fetch the hit map from the event:
+ HitMap inputHitMap = (HitMap) (event.get(m_inputHitMapName));
+
+ // Use this to record calorimeter hits that have only 1 or 2 other hits in 3x3 grid:
+ Set<CalorimeterHit> allLocallyIsolatedHits_12 = new HashSet<CalorimeterHit>();
+
+ // Use this to record calorimeter hits that have only 3 other hits in 3x3 grid:
+ Set<CalorimeterHit> allLocallyIsolatedHits_3 = new HashSet<CalorimeterHit>();
+
+ // Cache neighbour information to save lookup time
+ Map<CalorimeterHit,List<CalorimeterHit>> cacheNeighbours = new HashMap<CalorimeterHit,List<CalorimeterHit>>();
+
+ // Debug stuff
+ List<CalorimeterHit> listOfIsolatedHitsForEvent = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> listOfBusyHitsForEvent = new Vector<CalorimeterHit>();
+
+ // There is a bug -- neighbour relationships should be transient but aren't always.
+ // We'll keep a note of these.
+ Map<Long,List<Long>> buggedNeighbourRelationships = new HashMap<Long,List<Long>>();
+
+ Map<CalorimeterHit, List<CalorimeterHit>> cacheNeighboursOfHitsThatNeedWideWindow = new HashMap<CalorimeterHit,List<CalorimeterHit>>();
+ for (Long hitID : inputHitMap.keySet()) {
+ // Find all neighbouring cells
+ int dLayer = 1;
+ int dU = 1;
+ int dV = 1;
+ CalorimeterHit hit = inputHitMap.get(hitID);
+ if (hit.getCellID() != hitID) { throw new AssertionError("Internal consistency failure"); }
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hitID);
+ long[] neighbours = id.getNeighbourIDs(dLayer, dU, dV);
+
+ // Check which have a hit
+ List<CalorimeterHit> listOfNeighbours = new ArrayList<CalorimeterHit>(26);
+ for (long neighbourID : neighbours) {
+ if (inputHitMap.keySet().contains(neighbourID)) {
+ CalorimeterHit neighbourHit = inputHitMap.get(neighbourID);
+ if (neighbourHit == null) { throw new AssertionError("Internal consistency failure"); }
+ listOfNeighbours.add(neighbourHit);
+
+ // Check whether I am a neighbour of my neighbours:
+ id.setID(neighbourID);
+ long[] neighboursOfNeighbour = id.getNeighbourIDs(dLayer, dU, dV);
+ boolean foundOK = false;
+ for (long idOfNeighbourOfNeighbour : neighboursOfNeighbour) {
+ if (idOfNeighbourOfNeighbour == hitID) { foundOK = true; break; }
+ }
+ if (!foundOK) {
+ // Here's an example of the bug
+ List<Long> buggedLinks = buggedNeighbourRelationships.get(hitID);
+ if (buggedLinks == null) {
+ buggedLinks = new Vector<Long>();
+ buggedNeighbourRelationships.put(hitID, buggedLinks);
+ }
+ buggedLinks.add(new Long(neighbourID));
+ }
+ }
+ }
+
+ cacheNeighbours.put(hit, listOfNeighbours);
+
+ if (listOfNeighbours.size()==0) {
+ // This hit looks isolated. Check again in
+ // 2x2x1 window in case it's a low-angle track
+ id.setID(hitID);
+ long[] neighboursInWideWindow = id.getNeighbourIDs(dLayer, dU+1, dV+1);
+ List<CalorimeterHit> listOfNeighboursInWideWindow = new ArrayList<CalorimeterHit>(26);
+ for (long neighbourID : neighboursInWideWindow) {
+ if (inputHitMap.keySet().contains(neighbourID)) {
+ CalorimeterHit neighbourHit = inputHitMap.get(neighbourID);
+ if (neighbourHit == null) { throw new AssertionError("Internal consistency failure"); }
+ listOfNeighboursInWideWindow.add(neighbourHit);
+ }
+ }
+ if (listOfNeighboursInWideWindow.size()>0 && listOfNeighboursInWideWindow.size()<3) {
+ // Accept it for now
+ cacheNeighboursOfHitsThatNeedWideWindow.put(hit, listOfNeighboursInWideWindow);
+ }
+ }
+ }
+
+ // Now, what about those hits that needed a wide window?
+ // We may need to go back and reverse-add links to them.
+ for (CalorimeterHit hit : cacheNeighboursOfHitsThatNeedWideWindow.keySet()) {
+ List<CalorimeterHit> mainNeighbourList = cacheNeighbours.get(hit);
+ for (CalorimeterHit neighbour : cacheNeighboursOfHitsThatNeedWideWindow.get(hit)) {
+ if (cacheNeighboursOfHitsThatNeedWideWindow.keySet().contains(neighbour)) {
+ // This is another wide-window hit. Make sure they link reciprocally.
+ // Put it in the primary neighbour list:
+ mainNeighbourList.add(neighbour);
+ // Crosscheck of the secondary neighbour list:
+ List<CalorimeterHit> wideNeighboursOfNeighbour = cacheNeighboursOfHitsThatNeedWideWindow.get(neighbour);
+ if (!(wideNeighboursOfNeighbour.contains(hit))) { throw new AssertionError("Book-keeping error"); }
+ } else {
+ // This is not a wide-window hit.
+ List<CalorimeterHit> neighboursOfNeighbour = cacheNeighbours.get(neighbour);
+ if (neighboursOfNeighbour.contains(hit)) { throw new AssertionError("Book-keeping error"); }
+ if (neighboursOfNeighbour.size() > 2) {
+ // This is a DoubleHit so we can't add any more to it
+ } else {
+ neighboursOfNeighbour.add(hit);
+ mainNeighbourList.add(neighbour);
+ }
+ }
+ }
+ }
+
+ // Verify book-keeping
+ for (CalorimeterHit hit : inputHitMap.values()) {
+ List<CalorimeterHit> listOfNeighbours = cacheNeighbours.get(hit);
+ for (CalorimeterHit neighbour : listOfNeighbours) {
+ List<CalorimeterHit> listofNeighboursOfNeighbour = cacheNeighbours.get(neighbour);
+ if (!(listofNeighboursOfNeighbour.contains(hit))) {
+ // This could be the known bug
+ List<Long> checkBug = buggedNeighbourRelationships.get(new Long(hit.getCellID()));
+ if (checkBug != null && checkBug.contains(neighbour.getCellID())) {
+ // OK, no problem
+ } else {
+ System.out.println("Looking at hit ["+hit.getCellID()+"] which has "+listOfNeighbours.size()+" neighbours");
+ if (cacheNeighboursOfHitsThatNeedWideWindow.keySet().contains(hit)) {
+ System.out.println("Hit needed a wide window, and had "+cacheNeighboursOfHitsThatNeedWideWindow.get(hit).size()+" wide-window neighbours.");
+ }
+ System.out.println("Looking at its neighbour ["+neighbour.getCellID()+"] which has "+listofNeighboursOfNeighbour.size()+" neighbours");
+ if (cacheNeighboursOfHitsThatNeedWideWindow.keySet().contains(neighbour)) {
+ System.out.println("Neighbour needed a wide window, and had "+cacheNeighboursOfHitsThatNeedWideWindow.get(neighbour).size()+" wide-window neighbours.");
+ }
+ throw new AssertionError("Book-keeping error");
+ }
+ }
+ }
+ }
+
+ // Categorise hits:
+ for (CalorimeterHit hit : inputHitMap.values()) {
+ List listOfNeighbours = cacheNeighbours.get(hit);
+ if (listOfNeighbours.size()==1 || listOfNeighbours.size()==2) {
+ // Hit looks like a candidate for a MIP
+ allLocallyIsolatedHits_12.add(hit);
+ } else if (listOfNeighbours.size()==3) {
+ // 3 neighbours... not really isolated, but could be a one-off
+ allLocallyIsolatedHits_3.add(hit);
+ } else if (listOfNeighbours.size()==0) {
+ listOfIsolatedHitsForEvent.add(hit);
+ } else {
+ listOfBusyHitsForEvent.add(hit);
+ }
+ }
+
+ List<CalorimeterHit> listOfSingleHitsForEvent = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> listOfDoubleHitsForEvent = new Vector<CalorimeterHit>();
+ listOfSingleHitsForEvent.addAll(allLocallyIsolatedHits_12);
+ listOfDoubleHitsForEvent.addAll(allLocallyIsolatedHits_3);
+ String singleName = m_inputHitMapName + "SingleHits";
+ String doubleName = m_inputHitMapName + "DoubleHits";
+ String isolatedName = m_inputHitMapName + "IsolatedHits";
+ String busyName = m_inputHitMapName + "BusyHits";
+ event.put(singleName, listOfSingleHitsForEvent);
+ event.put(doubleName, listOfDoubleHitsForEvent);
+ event.put(isolatedName, listOfIsolatedHitsForEvent);
+ event.put(busyName, listOfBusyHitsForEvent);
+
+ List<CalorimeterHit> listOfHitsWithWideWindow = new Vector<CalorimeterHit>();
+ listOfHitsWithWideWindow.addAll(cacheNeighboursOfHitsThatNeedWideWindow.keySet());
+ String wideName = m_inputHitMapName + "WideHits";
+ event.put(wideName, listOfHitsWithWideWindow);
+
+ // OK. Now we try to build the MIPs.
+ // We need to avoid forking. Do this with 2+1 rules:
+ //
+ // Rule 1: A double-hit cannot connect to >2 single-hits
+ // Rule 2: A contiguous block of double-hits cannot connect to >2 single-hits
+ // Rule 3: A mip cluster must have at least 2 SingleHits in it.
+ //
+ // where a double-hit is a member of locallyIsolatedHits_3
+ // and a single-hit is a member of locallyIsolatedHits_12.
+ //
+ // (Really, rule 2 is an instance of rule 1, but never mind.)
+
+ Set<CalorimeterHit> unusedLocallyIsolatedHits_12 = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> unusedLocallyIsolatedHits_3 = new HashSet<CalorimeterHit>();
+ unusedLocallyIsolatedHits_12.addAll(allLocallyIsolatedHits_12);
+ unusedLocallyIsolatedHits_3.addAll(allLocallyIsolatedHits_3);
+
+ List<Cluster> outputList = new Vector<Cluster>();
+ while (unusedLocallyIsolatedHits_12.size() > 0) {
+ List<CalorimeterHit> protoCluster = new Vector<CalorimeterHit>();
+ CalorimeterHit seed = unusedLocallyIsolatedHits_12.iterator().next();
+ recursivelyAddSingleHit(seed, protoCluster, unusedLocallyIsolatedHits_12, unusedLocallyIsolatedHits_3, cacheNeighbours, allLocallyIsolatedHits_12, allLocallyIsolatedHits_3);
+ if (protoCluster.size() > 3) {
+ // Enough to make a real cluster. Check for single hits...
+ int countSingleHits = 0;
+ for (CalorimeterHit hit : protoCluster) {
+ if (allLocallyIsolatedHits_12.contains(hit)) {
+ countSingleHits++;
+ }
+ }
+ if (countSingleHits >= 2) {
+ // Great -- passes all cuts
+ BasicCluster newCluster = new BasicCluster();
+ for (CalorimeterHit hit : protoCluster) {
+ newCluster.addHit(hit);
+ }
+ outputList.add(newCluster);
+ } else {
+ }
+ }
+ }
+
+ // OK, done.
+ event.put(m_outputClusterListName, outputList);
+ }
+
+ void recursivelyAddSingleHit(CalorimeterHit hit, List<CalorimeterHit> protoCluster, Set<CalorimeterHit> unusedLocallyIsolatedHits_12, Set<CalorimeterHit> unusedLocallyIsolatedHits_3, Map<CalorimeterHit,List<CalorimeterHit>> cacheNeighbours, Set<CalorimeterHit> allLocallyIsolatedHits_12, Set<CalorimeterHit> allLocallyIsolatedHits_3)
+ {
+ // Remove the hit, since we won't use it in another cluster.
+ boolean removedOK = unusedLocallyIsolatedHits_12.remove(hit);
+ if (!removedOK) { throw new AssertionError("Internal consistency failure"); }
+ // Single hits are always safe to add => do so
+ protoCluster.add(hit);
+ // Find the neighbours
+ List<CalorimeterHit> listOfNeighbours = cacheNeighbours.get(hit);
+ for (CalorimeterHit neighbour : listOfNeighbours) {
+ if (unusedLocallyIsolatedHits_12.contains(neighbour)) {
+ recursivelyAddSingleHit(neighbour, protoCluster, unusedLocallyIsolatedHits_12, unusedLocallyIsolatedHits_3, cacheNeighbours, allLocallyIsolatedHits_12, allLocallyIsolatedHits_3);
+ } else if (unusedLocallyIsolatedHits_3.contains(neighbour)) {
+ recursivelyAddDoubleHit(neighbour, protoCluster, unusedLocallyIsolatedHits_12, unusedLocallyIsolatedHits_3, cacheNeighbours, allLocallyIsolatedHits_12, allLocallyIsolatedHits_3);
+ } else {
+ boolean neighbourIsSingleHit = (allLocallyIsolatedHits_12.contains(neighbour));
+ boolean neighbourIsUnusedSingleHit = (unusedLocallyIsolatedHits_12.contains(neighbour));
+ boolean neighbourIsDoubleHit = (allLocallyIsolatedHits_3.contains(neighbour));
+ boolean neighbourIsUnusedDoubleHit = (unusedLocallyIsolatedHits_3.contains(neighbour));
+ }
+ }
+ return;
+ }
+
+ // This ONLY gets called from within recursivelyAddSingleHit()
+ void recursivelyAddDoubleHit(CalorimeterHit hit, List<CalorimeterHit> protoCluster, Set<CalorimeterHit> unusedLocallyIsolatedHits_12, Set<CalorimeterHit> unusedLocallyIsolatedHits_3, Map<CalorimeterHit,List<CalorimeterHit>> cacheNeighbours, Set<CalorimeterHit> allLocallyIsolatedHits_12, Set<CalorimeterHit> allLocallyIsolatedHits_3)
+ {
+ // Is it safe to add this hit + contiguous DoubleHits?
+ // Start by finding the block of contiguous DoubleHits.
+ List<CalorimeterHit> contiguousDoubleHits = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> scannedDoubleHits = new Vector<CalorimeterHit>();
+ List<CalorimeterHit> unscannedDoubleHits = new Vector<CalorimeterHit>();
+ contiguousDoubleHits.add(hit);
+ unscannedDoubleHits.add(hit);
+ while (scannedDoubleHits.size() < contiguousDoubleHits.size()) {
+ if (scannedDoubleHits.size() + unscannedDoubleHits.size() != contiguousDoubleHits.size()) { throw new AssertionError("book-keeping failure"); }
+ CalorimeterHit currentHit = unscannedDoubleHits.get(0);
+ unscannedDoubleHits.remove(currentHit);
+ scannedDoubleHits.add(currentHit);
+ List<CalorimeterHit> listOfNeighbours = cacheNeighbours.get(currentHit);
+ for (CalorimeterHit neighbourOfCurrentHit : listOfNeighbours) {
+ boolean isDoubleHit = unusedLocallyIsolatedHits_3.contains(neighbourOfCurrentHit);
+ boolean unscanned = !(scannedDoubleHits.contains(neighbourOfCurrentHit));
+ boolean notPending = !(unscannedDoubleHits.contains(neighbourOfCurrentHit));
+ if (isDoubleHit && unscanned && notPending) {
+ unscannedDoubleHits.add(neighbourOfCurrentHit);
+ contiguousDoubleHits.add(neighbourOfCurrentHit);
+ }
+ }
+ }
+ if (scannedDoubleHits.size() + unscannedDoubleHits.size() != contiguousDoubleHits.size()) { throw new AssertionError("book-keeping failure"); }
+ if (unscannedDoubleHits.size() != 0) { throw new AssertionError("book-keeping failure"); }
+
+ // OK. We have our block of contiguous DoubleHits.
+ // Check that between them, they link to exactly two SingleHits.
+ Set<CalorimeterHit> linkedSingleHits = new HashSet<CalorimeterHit>();
+ for (CalorimeterHit doubleHit : contiguousDoubleHits) {
+ List<CalorimeterHit> listOfNeighbours = cacheNeighbours.get(doubleHit);
+ for (CalorimeterHit neighbourOfDoubleHit : listOfNeighbours) {
+ if (allLocallyIsolatedHits_12.contains(neighbourOfDoubleHit)) {
+ linkedSingleHits.add(neighbourOfDoubleHit);
+ }
+ }
+ }
+ if (linkedSingleHits.size() > 2) {
+ // Forbidden -- this implies a fork
+ // Return without adding any hits to cluster and remove
+ // the DoubleHits so that they don't get used again.
+ for (CalorimeterHit doubleHit : contiguousDoubleHits) {
+ unusedLocallyIsolatedHits_3.remove(doubleHit);
+ }
+ return;
+ } else {
+ // Allowed.
+ // Add the whole contiguous block, then add the SingleHits.
+ for (CalorimeterHit doubleHit : contiguousDoubleHits) {
+ unusedLocallyIsolatedHits_3.remove(doubleHit);
+ protoCluster.add(doubleHit);
+ }
+ for (CalorimeterHit singleHit : linkedSingleHits) {
+ if (unusedLocallyIsolatedHits_12.contains(singleHit)) {
+ recursivelyAddSingleHit(singleHit, protoCluster, unusedLocallyIsolatedHits_12, unusedLocallyIsolatedHits_3, cacheNeighbours, allLocallyIsolatedHits_12, allLocallyIsolatedHits_3);
+ }
+ }
+ }
+ }
+}