diff -N SteveMIPReassignmentAlgorithm.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ SteveMIPReassignmentAlgorithm.java	8 Jul 2008 16:36:23 -0000	1.1
@@ -0,0 +1,306 @@

+package org.lcsim.contrib.uiowa;
+
+import java.util.*; 
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import hep.physics.vec.*;
+import org.lcsim.geometry.Subdetector;
+import org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+public class SteveMIPReassignmentAlgorithm implements ReassignClustersAlgorithm {
+
+    protected LocalHelixExtrapolator m_extrap;
+    protected EventHeader m_event;
+    protected double m_limit;
+    protected Map<Track, Hep3Vector> m_cachePoint;
+    protected Map<Track, Hep3Vector> m_cacheTangent;
+
+    public SteveMIPReassignmentAlgorithm(EventHeader event, double limit) {
+	m_event = event;
+	m_extrap = new LocalHelixExtrapolator();
+	m_extrap.process(event); // pick up geometry info
+	m_limit = limit;
+	m_cachePoint = new HashMap<Track, Hep3Vector>();
+	m_cacheTangent = new HashMap<Track, Hep3Vector>();
+    }
+
+    public Double computeFigureOfMerit(Track tr, Cluster clus) {
+	// Handle possibility of multiple-track-tracks
+	if (tr.getTracks().size() != 0) {
+	    // There's more than one daughter -- is this a known case?
+	    if (tr instanceof ReclusterDriver.MultipleTrackTrack) {
+		// OK, this is understood & happens when multiple tracks point
+		// to a single seed. Check extrapolation for all the tracks and
+		// take the best score among them.
+		Double bestAmongDaughters = null;
+		for (Track daughterTrack : tr.getTracks()) {
+		    Double daughterFigureOfMerit = this.computeFigureOfMerit(daughterTrack, clus);
+		    if (daughterFigureOfMerit == null) {
+			// No valid score for this daughter
+			continue;
+		    } else if (bestAmongDaughters == null) {
+			// No previous best score => this is the new best
+			bestAmongDaughters = daughterFigureOfMerit;
+		    } else if (daughterFigureOfMerit < bestAmongDaughters) {
+			// This is the new best (lowest) score
+			bestAmongDaughters = daughterFigureOfMerit;
+		    }
+		}
+		// All done
+		return bestAmongDaughters;
+	    } else {
+		// Don't know about this case -- complain.
+		throw new AssertionError("WATCH OUT! Unhandled case where track of class "+tr.getClass().getName()+" has "+tr.getTracks().size()+" daughter tracks!");
+	    }
+	} else {
+	    // Only one track -- find track helix parameters
+	    m_extrap.performExtrapolation(tr);
+	}
+
+	BasicHep3Vector endOfMipPoint = new BasicHep3Vector();
+	BasicHep3Vector endOfMipTangentUnit = new BasicHep3Vector();
+	try {
+	    findPointAndTangent(tr, endOfMipPoint, endOfMipTangentUnit);
+	} catch (ExtrapolationFailureException x) {
+	    // Exception due to failure to extrapolate track to ECAL
+	    return null;
+	}
+
+	Hep3Vector clusterPosition = new BasicHep3Vector(clus.getPosition());
+	Hep3Vector displacement = VecOp.sub(clusterPosition, endOfMipPoint);
+        Hep3Vector displacementUnit = VecOp.unit(displacement);
+	double cosTheta = VecOp.dot(endOfMipTangentUnit, displacementUnit);
+	double angle=Math.acos(cosTheta);
+
+	if (angle < m_limit) {
+	    return new Double(angle);
+	} else {
+	    return null;
+	}
+    }
+
+    protected void findPointAndTangent(Track tr, BasicHep3Vector outputPoint, BasicHep3Vector outputTangentUnit) throws ExtrapolationFailureException {
+	Hep3Vector point = m_cachePoint.get(tr);
+	Hep3Vector tangent = m_cacheTangent.get(tr);
+	if (point == null && tangent == null) {
+	    BasicHep3Vector tmpPoint = new BasicHep3Vector();
+	    BasicHep3Vector tmpTangent = new BasicHep3Vector();
+	    try {
+		findPointAndTangentNoCache(tr, tmpPoint, tmpTangent);
+		m_cachePoint.put(tr, tmpPoint);
+		m_cacheTangent.put(tr, tmpTangent);
+		point = tmpPoint;
+		tangent = tmpTangent;
+	    } catch (ExtrapolationFailureException x) {
+		// Exception due to failure to extrapolate track to ECAL
+		throw x;
+	    }
+	} else if (point!=null && tangent!=null) {
+	    // retrieved from cache OK
+	} else {
+	    throw new AssertionError("Cache error");
+	}
+
+	outputPoint.setV(point.x(), point.y(), point.z());
+	outputTangentUnit.setV(tangent.x(), tangent.y(), tangent.z());
+	return;
+    }
+
+    protected void findPointAndTangentNoCache(Track tr, BasicHep3Vector outputPoint, BasicHep3Vector outputTangentUnit) throws ExtrapolationFailureException {
+	// Find the MIP trace for the track
+	Object objTrackMipMap = m_event.get("TrackMipMap");
+	Object objMipClusILMap = m_event.get("MipClusILMap");
+	Map<Track, BasicCluster> trkmipmap = ( Map<Track, BasicCluster> ) (objTrackMipMap);
+	Map<BasicCluster, Integer> clusILmap = ( Map<BasicCluster, Integer> ) (objMipClusILMap);
+
+	BasicCluster mip = trkmipmap.get(tr);
+	if (mip == null) {
+	    throw new AssertionError("Track of class "+tr.getClass().getName()+" with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" has null mip!");
+	} else if (mip.getCalorimeterHits().size()==0) {
+	    Hep3Vector interceptPoint = m_extrap.getLastInterceptPoint();
+	    if (interceptPoint == null) {
+		throw new ExtrapolationFailureException("Failure to extrapolate");
+	    } else {
+		outputPoint.setV(interceptPoint.x(), interceptPoint.y(), interceptPoint.z());
+	    }
+	    Hep3Vector tangent = m_extrap.getTangent();
+	    outputTangentUnit.setV(tangent.x(), tangent.y(), tangent.z());
+	    return;
+	}
+
+	// Find the outermost hit(s)
+
+	// 1) Split hits up by subdetector
+	Map<Subdetector, List<CalorimeterHit>> hitsInEachSubdetector = splitHitsBySubdetector(mip);
+	// 2) In each subdetector, find outermost layer hit & hits in it
+	Map<Subdetector, List<CalorimeterHit>> outermostHitsInEachSubdetector = findOutermostHitsInEachSubdetector(hitsInEachSubdetector);
+	// 3) Find which subdetector has the furthest-out hits (checked with |z|)
+	Subdetector outermostSubdet = findSubdetectorWithOutermostHits(outermostHitsInEachSubdetector);
+	// 4) For the hits in the furthest-out layer in that subdetector, find the mean position:
+	Hep3Vector outermostPosition = findMeanPositionOfHits(outermostHitsInEachSubdetector.get(outermostSubdet));
+	outputPoint.setV(outermostPosition.x(), outermostPosition.y(), outermostPosition.z());
+
+	// Check layer index:
+	int lastLayer = getLayer(outermostHitsInEachSubdetector.get(outermostSubdet).get(0));
+	for (CalorimeterHit hit : outermostHitsInEachSubdetector.get(outermostSubdet)) {
+	    if (getLayer(hit) != lastLayer) { throw new AssertionError("Book-keeping error: layer mis-match"); }
+	    if (hit.getSubdetector() != outermostSubdet) { throw new AssertionError("Book-keeping error"); }
+	}
+
+	// OK, now we have the position & layer. Next, we need the tangent at/near that point.
+	try {
+	    Hep3Vector tangentUnit = getTangentUnit(lastLayer, outermostSubdet);
+	    outputTangentUnit.setV(tangentUnit.x(), tangentUnit.y(), tangentUnit.z());
+	} catch (ExtrapolationFailureException x) {
+	    // Failure...
+	    Hep3Vector tangentUnit = m_extrap.getTangent();
+	    if (tangentUnit != null) {
+		// Use tangent at ECAL front surface instead (iffy...)
+		outputTangentUnit.setV(tangentUnit.x(), tangentUnit.y(), tangentUnit.z());
+	    } else {
+		// Complete failure
+		throw new ExtrapolationFailureException("Failure to compute tangent");
+	    }
+	}
+
+	// All done
+	return;
+    }
+
+    protected Hep3Vector getTangentUnit(int lastLayer, Subdetector outermostSubdet) throws ExtrapolationFailureException {
+	// We need the tangent in the layer of interest.
+	// To do this, we extrapolate the track to the corresponding layer, and the previous one.
+	// (Or, if this is the first layer, use the next one).
+	Hep3Vector trackPointInLayer_N = null;
+	Hep3Vector trackPointInLayer_NminusOne = null;
+	int layerN = lastLayer;
+	if (lastLayer == 0) {
+	    layerN = lastLayer + 1;
+	}
+
+	// To do the extrapolation, we need to know if we're looking at a barrel or an endcap subdetector.
+	String calName = outermostSubdet.getName();
+	if (calName.compareTo("HADBarrel")==0) {
+	    trackPointInLayer_N = m_extrap.extendToHCALBarrelLayer(layerN);
+	    trackPointInLayer_NminusOne = m_extrap.extendToHCALBarrelLayer(layerN-1);
+	} else if (calName.compareTo("HADEndcap")==0) {
+	    trackPointInLayer_N = m_extrap.extendToHCALEndcapLayer(layerN);
+	    trackPointInLayer_NminusOne = m_extrap.extendToHCALEndcapLayer(layerN-1);
+	} else if (calName.compareTo("EMBarrel")==0) {
+	    trackPointInLayer_N = m_extrap.extendToECALBarrelLayer(layerN);
+	    trackPointInLayer_NminusOne = m_extrap.extendToECALBarrelLayer(layerN-1);
+	} else if (calName.compareTo("EMEndcap")==0) {
+	    trackPointInLayer_N = m_extrap.extendToECALEndcapLayer(layerN);
+	    trackPointInLayer_NminusOne = m_extrap.extendToECALEndcapLayer(layerN-1);
+	} else {
+	    throw new AssertionError("Calorimeter component "+calName+" not recognized!");
+	}
+
+	if (trackPointInLayer_N == null || trackPointInLayer_NminusOne==null) {
+	    System.out.println("ERROR: Extrapolation failure when computing tangent for lastLayer="+lastLayer+" and subdet="+outermostSubdet.getClass().getName());
+	    System.out.println("Tried to extrapolate to layer "+layerN+" and layer "+(layerN-1)+" of "+calName);
+	    if (trackPointInLayer_N == null) {
+		System.out.println("  -- point in layer "+layerN+" was NULL!");
+	    } else {
+		System.out.println("  -- point in layer "+layerN+" found OK: ("+trackPointInLayer_N.x()+", "+trackPointInLayer_N.y()+", "+trackPointInLayer_N.z()+")");
+	    }
+	    if (trackPointInLayer_NminusOne == null) {
+		System.out.println("  -- point in layer "+(layerN-1)+" was NULL!");
+	    } else {
+		System.out.println("  -- point in layer "+(layerN-1)+" found OK: ("+trackPointInLayer_NminusOne.x()+", "+trackPointInLayer_NminusOne.y()+", "+trackPointInLayer_NminusOne.z()+")");
+	    }
+	    throw new ExtrapolationFailureException("Unable to compute tangent");
+	}
+	    
+
+	Hep3Vector tangentUnit = VecOp.unit(VecOp.sub(trackPointInLayer_N, trackPointInLayer_NminusOne));
+	return tangentUnit;
+    }
+
+    // Split hits up by subdetector
+    protected Map<Subdetector, List<CalorimeterHit>> splitHitsBySubdetector(Cluster mip) {
+	if (mip.getCalorimeterHits().size()==0) { throw new AssertionError("MIP has no hits!"); }
+	Map<Subdetector, List<CalorimeterHit>> hitsInEachSubdetector = new HashMap<Subdetector, List<CalorimeterHit>>();
+	for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+	    Subdetector det = hit.getSubdetector();
+	    List<CalorimeterHit> hitList = hitsInEachSubdetector.get(det);
+	    if (hitList == null) {
+		hitList = new Vector<CalorimeterHit>();
+		hitsInEachSubdetector.put(det, hitList);
+	    }
+	    hitList.add(hit);
+	}
+	return hitsInEachSubdetector;
+    }
+
+    // In each subdetector, find outermost layer hit & hits in it
+    Map<Subdetector, List<CalorimeterHit>> findOutermostHitsInEachSubdetector(Map<Subdetector, List<CalorimeterHit>> hitsInEachSubdetector) {
+	Map<Subdetector, List<CalorimeterHit>> outermostHitsInEachSubdetector = new HashMap<Subdetector, List<CalorimeterHit>>();
+	if (hitsInEachSubdetector.keySet().size()==0) { throw new AssertionError("No subdetectors supplied!"); }
+	for (Subdetector det : hitsInEachSubdetector.keySet()) {
+	    List<CalorimeterHit> hits = hitsInEachSubdetector.get(det);
+	    if (hits.size()==0) { throw new AssertionError("Book-keeping error"); }
+	    int layerMax = -1;
+	    List<CalorimeterHit> hitsInMaxLayer = new Vector<CalorimeterHit>();
+	    for (CalorimeterHit hit : hits) {
+		int layer = getLayer(hit);
+		if (layer > layerMax) {
+		    layerMax = layer;
+		    hitsInMaxLayer.clear();
+		}
+		if (layer == layerMax) {
+		    hitsInMaxLayer.add(hit);
+		}
+	    }
+	    if (layerMax<0) { throw new AssertionError("Book-keeping error"); }
+	    if (hitsInMaxLayer.size()==0) { throw new AssertionError("Book-keeping error"); }
+	    outermostHitsInEachSubdetector.put(det, hitsInMaxLayer);
+	}
+	return outermostHitsInEachSubdetector;
+    }
+
+    // Find which subdetector has the furthest-out hits (checked with |z|)
+    Subdetector findSubdetectorWithOutermostHits(Map<Subdetector, List<CalorimeterHit>> outermostHitsInEachSubdetector) {
+	if (outermostHitsInEachSubdetector.keySet().size()==0) { throw new AssertionError("No subdetectors supplied!"); }
+	Subdetector subdetectorWithMaxZ = null;
+	double maxZ = 0;
+	for (Subdetector det : outermostHitsInEachSubdetector.keySet()) {
+	    List<CalorimeterHit> hits = outermostHitsInEachSubdetector.get(det);
+	    if (hits.size()==0) { throw new AssertionError("No hits supplied for this subdetector!"); }
+	    for (CalorimeterHit hit : hits) {
+		double modZ = Math.abs(hit.getPosition()[2]);
+		if (modZ > maxZ) {
+		    maxZ = modZ;
+		    subdetectorWithMaxZ = det;
+		}
+	    }
+	}
+	if (subdetectorWithMaxZ == null) { throw new AssertionError("Book-keeping error when scanning over "+outermostHitsInEachSubdetector.keySet().size()+" subdetectors (maxZ="+maxZ+")"); }
+	return subdetectorWithMaxZ;
+    }
+
+    // For the hits in the furthest-out layer in that subdetector, find the mean position:
+    Hep3Vector findMeanPositionOfHits(Collection<CalorimeterHit> hitsInOutermostLayer) {
+	Hep3Vector meanPosition = new BasicHep3Vector(0.0, 0.0, 0.0);
+	for (CalorimeterHit hit : hitsInOutermostLayer) {
+	    Hep3Vector position = new BasicHep3Vector(hit.getPosition());
+	    double scaleFactor = 1.0 / ((double)(hitsInOutermostLayer.size()));
+	    Hep3Vector contribution = VecOp.mult(scaleFactor, position);
+	    meanPosition = VecOp.add(meanPosition, contribution);
+	}
+	return meanPosition;
+    }
+
+    protected int getLayer(CalorimeterHit hit) {
+        org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+        id.setID(hit.getCellID());
+        int layer = id.getLayer();
+        return layer;
+    }
+
+    private class ExtrapolationFailureException extends java.lang.Exception {
+	public ExtrapolationFailureException(String m) { super(m); }
+    }
+}

diff -N SteveMipWrapper.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ SteveMipWrapper.java	8 Jul 2008 16:36:23 -0000	1.1
@@ -0,0 +1,55 @@

+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.event.*;
+
+public class SteveMipWrapper extends Driver
+{
+    public SteveMipWrapper() 
+    {
+	// Set up Hitmaps
+	HitListToHitMapDriver hitmapEcalBarrel = new HitListToHitMapDriver();
+	HitListToHitMapDriver hitmapHcalBarrel = new HitListToHitMapDriver();
+	HitListToHitMapDriver hitmapEcalEndcap = new HitListToHitMapDriver();
+	HitListToHitMapDriver hitmapHcalEndcap = new HitListToHitMapDriver();
+	hitmapEcalBarrel.addInputList("EcalBarrDigiHits");
+	hitmapEcalEndcap.addInputList("EcalEndcapDigiHits");
+	hitmapHcalBarrel.addInputList("HcalBarrDigiHits");
+	hitmapHcalEndcap.addInputList("HcalEndcapDigiHits");
+	hitmapEcalBarrel.setOutput("EMBarrhitmap");
+	hitmapEcalEndcap.setOutput("EMEndcaphitmap");
+	hitmapHcalBarrel.setOutput("HADBarrhitmap");
+	hitmapHcalEndcap.setOutput("HADEndcaphitmap");
+	add(hitmapEcalBarrel);
+	add(hitmapEcalEndcap);
+	add(hitmapHcalBarrel);
+	add(hitmapHcalEndcap);
+
+	// Steve's MIP code needs an input track list called "PerfectTracks".
+	// Set that up by cloning FSReconTracks:
+	DecisionMakerSingle<Track> dummy = new DummyDecisionMakerSingle<Track>();
+	add(new ListFilterDriver(dummy, "FSReconTracks", "PerfectTracks", Track.class));
+
+	// Steve's MIP code
+	// ----------------
+
+	// Four parameter values below taken from his driver routine:
+	double dminE = 0.005;  // hard-coded now - change to detector independent value
+	double dminH = 3*dminE;  //  same here
+	double dminTrC = 0.01;  // default range is 0.0175 - 0.0225 tuned with muons
+	double mincd = 3.5;  // min density to make mip - default is 3.5 tuned with muons
+	// Set up & turn on:
+	org.lcsim.contrib.SteveMagill.TrackMipClusterDriver TMdriver = new org.lcsim.contrib.SteveMagill.TrackMipClusterDriver(dminE, dminH, dminTrC, mincd);
+	String[] tmcollnames = {"EcalBarrDigiHits","HcalBarrDigiHits","EcalEndcapDigiHits","HcalEndcapDigiHits"};
+	TMdriver.setCollectionNames(tmcollnames);
+	TMdriver.setClusterNameExtension("TMClusters");
+	add(TMdriver);
+    }
+
+    public void process(EventHeader event) {
+	super.process(event);
+    }
+}