diff -N MIPChargedParticleMaker.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ MIPChargedParticleMaker.java	10 Aug 2006 22:57:56 -0000	1.1
@@ -0,0 +1,218 @@

+package org.lcsim.recon.pfa.identifier;
+
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.particle.Particle;
+import org.lcsim.event.Track;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Cluster;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.base.BaseReconstructedParticle;
+
+/**
+ * Given a list of MIP clusters and a list of tracks,
+ * try to connect tracks to MIP segments and
+ * make a list of charged ReconstructedParticles.
+ *
+ * Each track is matched to exactly zero or one MIP,
+ * and appears in zero or one ReconstructedParticles.
+ * But a MIP/particle may be associated with more than
+ * one track.
+ *
+ * Optionally, the user may supply further list of clusters.
+ * If the MIP is part of one of these clusters (via Cluster.getClusters),
+ * then the entire cluster is added to the ReconstructedParticle instead.
+ * The parent must be unique.
+ *
+ * @version $Id: MIPChargedParticleMaker.java,v 1.1 2006/08/10 22:57:56 mcharles Exp $
+ */
+
+public class MIPChargedParticleMaker
+ extends Driver
+{
+    /** Simple constructor. */
+    public MIPChargedParticleMaker() {
+	m_clusterLists = new HashMap<String,String>();
+    }
+
+    public void setInputTrackList(String name) { m_inputTrackListName = name; }
+    public void setOutputTrackList(String name){ m_outputTrackListName = name; }
+    public void setInputMIPList(String name){ m_inputMIPListName = name; }
+    public void setOutputMIPList(String name){ m_outputMIPListName = name; }
+    public void setOutputParticleList(String name){ m_outputParticleListName = name; }
+    public void setTrackMatcher(TrackClusterMatcher matcher) { m_matcher = matcher; }
+   
+    public void addClusterList(String inputName, String outputName) {
+	m_clusterLists.put(inputName, outputName);
+    }
+
+    public void process(EventHeader event)
+    {
+	super.process(event);
+	
+	// Inputs:
+	List<Track> inputTrackList = event.get(Track.class, m_inputTrackListName);
+	List<Cluster> inputMIPList = event.get(Cluster.class, m_inputMIPListName);
+
+	// Outputs
+	// Initially, all tracks and MIPs are unmatched.
+	Map<Track,Cluster> matchedTracks = new HashMap<Track,Cluster>();
+	List<Track> unmatchedTracks = new Vector<Track>(inputTrackList);
+	List<Cluster> unmatchedMIPs = new Vector<Cluster>(inputMIPList);
+	List<ReconstructedParticle> outputParticleList = new Vector<ReconstructedParticle>();
+
+	// Optional inputs and outputs
+	// Output lists are initially identical to the input lists; we will
+	// remove clusters as they are matched.
+	Map<String,List<Cluster>> inputClusterLists = new HashMap<String,List<Cluster>>();
+	Map<String,List<Cluster>> outputClusterLists = new HashMap<String,List<Cluster>>();
+	for (String inputName : m_clusterLists.keySet()) {
+	    String outputName = m_clusterLists.get(inputName);
+	    List<Cluster> inputList = event.get(Cluster.class, inputName);
+	    List<Cluster> outputList = new Vector<Cluster>(inputList);
+	    inputClusterLists.put(inputName, inputList);
+	    outputClusterLists.put(outputName, outputList);
+	}
+
+	// Try to match each track to a MIP segment
+	Map<Cluster,List<Track>> matchedMIPs = new HashMap<Cluster,List<Track>>();
+	for (Track tr : inputTrackList) {
+	    Cluster matchedMIP = m_matcher.matchTrackToCluster(tr, inputMIPList);
+	    if (matchedMIP != null) {
+		// Verify that the returned MIP is a member of inputMIPList
+		if ( !(inputMIPList.contains(matchedMIP)) ) {
+		    throw new AssertionError("Book-keeping error: MIP Matcher must return a member of the input list or null");
+		}
+		if ( !(matchedMIPs.keySet().contains(matchedMIP)) ) {
+		    // First time we've seen this MIP => make its track list
+		    matchedMIPs.put(matchedMIP, new Vector<Track>());
+		}
+		matchedMIPs.get(matchedMIP).add(tr);
+		unmatchedMIPs.remove(matchedMIP);
+		unmatchedTracks.remove(tr);
+		matchedTracks.put(tr, matchedMIP);
+	    }
+	}
+
+	// Optional: Look for each MIP's parent cluster
+	for (Cluster matchedMIP : matchedMIPs.keySet()) {
+	    Cluster uniqueParent = null;
+	    for (String inputListName : m_clusterLists.keySet()) {
+		String outputListName = m_clusterLists.get(inputListName);
+		List<Cluster> inputList = inputClusterLists.get(inputListName);
+		List<Cluster> outputList = outputClusterLists.get(outputListName);
+		
+		for (Cluster clus : inputList) {
+		    List<Cluster> daughters = recursivelyFindSubClusters(clus);
+		    if (daughters.contains(matchedMIP)) {
+			// Found a parent containing this MIP.
+			if (uniqueParent != null) {
+			    throw new AssertionError("Book-keeping error: Non-unique parent of MIP");
+			}
+			uniqueParent = clus;
+		    }
+		}
+		
+		// If we found a unique parent in this list, do book-keeping:
+		//    1) Remove parent from output list of unmatched clusters
+		//    2) Make associated track(s) point to parent cluster
+		if (uniqueParent != null) {
+		    outputList.remove(uniqueParent);
+		    for (Track tr : matchedMIPs.get(matchedMIP)) {
+			matchedTracks.put(tr, uniqueParent); // over-write previous mapping
+		    }
+		}
+	    }
+	}
+
+	// Now that we have the track:cluster association, make output
+	// particles. We have to watch for the special case where >1 track
+	// is matched to a cluster.
+	Map<Cluster,LocalReconstructedParticle> matchedClusters = new HashMap<Cluster,LocalReconstructedParticle> ();
+	for (Track tr : matchedTracks.keySet()) {
+	    Cluster clus = matchedTracks.get(tr);
+	    if (matchedClusters.keySet().contains(clus)) {
+		// Already used this cluster in a particle => just add the track to that particle
+		LocalReconstructedParticle part = matchedClusters.get(clus);
+		part.addTrack(tr);
+		part.recomputeKinematics();
+	    } else {
+		// This cluster hasn't been used yet -- initialize its particle
+		LocalReconstructedParticle part = new LocalReconstructedParticle();
+		part.addTrack(tr);
+		part.addCluster(clus);
+		part.recomputeKinematics();
+		matchedClusters.put(clus, part);
+	    }
+	}
+	outputParticleList.addAll(matchedClusters.values());
+
+	// Write out
+	event.put(m_outputTrackListName, unmatchedTracks);
+	event.put(m_outputMIPListName, unmatchedMIPs);
+	event.put(m_outputParticleListName, outputParticleList);
+	for (String inputName : m_clusterLists.keySet()) {
+	    String outputName = m_clusterLists.get(inputName);
+	    List<Cluster> outputList = outputClusterLists.get(outputName);
+	    event.put(outputName, outputList);
+	}
+    }
+
+    /**
+     * Internal utility routine
+     */
+    protected List<Cluster> recursivelyFindSubClusters(Cluster clus) 
+    {
+	List<Cluster> output = new Vector<Cluster>();
+	for (Cluster dau : clus.getClusters()) {
+	    output.addAll(recursivelyFindSubClusters(dau));
+	}
+	output.add(clus);
+	return output;
+    }
+
+    protected TrackClusterMatcher m_matcher;
+    protected String m_inputTrackListName;
+    protected String m_outputTrackListName;
+    protected String m_inputMIPListName;
+    protected String m_outputMIPListName;
+    protected String m_outputParticleListName;
+    protected Map<String,String> m_clusterLists;
+
+    // Really, this belongs in a central implementation.
+    // It may change at any time (hence private).
+    private class LocalReconstructedParticle extends BaseReconstructedParticle
+    {
+	public void setEnergy(double e) {
+	    _fourVec.setT(e);
+	}
+	public void setMomentum(Hep3Vector p3) {
+	    _fourVec.setV3(_fourVec.t(), p3);
+	}
+	public void setReferencePoint(Hep3Vector p3) {
+	    _referencePoint = new BasicHep3Vector(p3.x(), p3.y(), p3.z());
+	}
+	public void setCharge(double q) {
+	    _charge = q;
+	}
+	public void recomputeKinematics() {
+	    double energy = 0.0;
+	    for (Track tr : this.getTracks()) {
+		double[] trackMomentum = tr.getMomentum();
+		double trackMomentumMagSq = (trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
+		double mass = 0.140;
+		if (tr instanceof ReconTrack) {
+		    Particle truthParticle = ((ReconTrack)(tr)).getMCParticle();
+		    if (truthParticle != null) {
+			mass = truthParticle.getMass();
+		    }
+		}
+		double trackEnergy = Math.sqrt(trackMomentumMagSq + mass*mass);
+		energy += trackEnergy;
+	    }
+	    this.setEnergy(energy);
+	}
+    }
+}