diff -u -r1.35 -r1.36
--- ReclusterDTreeDriver.java	23 Jul 2008 17:41:06 -0000	1.35
+++ ReclusterDTreeDriver.java	4 Aug 2008 17:38:17 -0000	1.36
@@ -34,7 +34,7 @@

   * in this package, which uses the implementation in
   * org.lcsim.recon.cluster.directedtree developed by NIU).
   *

-  * @version $Id: ReclusterDTreeDriver.java,v 1.35 2008/07/23 17:41:06 mcharles Exp $

+  * @version $Id: ReclusterDTreeDriver.java,v 1.36 2008/08/04 17:38:17 mcharles Exp $

   * @author Mat Charles <[log in to unmask]>
   */
 

@@ -81,6 +81,8 @@

     protected boolean m_fixSingleTracksWithCone = true;
     protected boolean m_fixJetsWithCone = true;
     protected boolean m_useSimpleConeForReassignment = true;

+    
+    protected boolean m_debugSeedSplitting = false;

 
     public void writeExtraEventOutput(boolean writeExtra) {
 	m_writeExtraEventOutput = writeExtra;

@@ -204,6 +206,10 @@

 	Map<Cluster, Cluster> treeOfClump = new HashMap<Cluster,Cluster>();
 	Map<Cluster, Cluster> treeOfLeftoverHits = new HashMap<Cluster,Cluster>();
 

+	// Identify the start point of showers
+	org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator findCluster = new org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator();
+	findCluster.process(m_event); // picks up geometry
+

 	// Match tracks
 	// ------------
 	Map<Track,Cluster> tracksMatchedToClusters = new HashMap<Track,Cluster>();

@@ -235,6 +241,89 @@

 	    }
 	    Cluster matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, allMatchableClusters);
 	    if (matchedCluster != null) {

+ 		if (leftoverHitClusters.contains(matchedCluster)) {
+		    if (m_debugSeedSplitting) {
+			// Debug printout
+			Hep3Vector interceptPoint = findCluster.performExtrapolation(tr);
+			if (interceptPoint != null) {
+			    double primaryDist = proximity(matchedCluster, interceptPoint);
+			    int innermostLayerOfMatchedCluster = 99;
+			    for (CalorimeterHit hit : matchedCluster.getCalorimeterHits()) {
+				int layer = getLayer(hit);
+				if (layer < innermostLayerOfMatchedCluster) { innermostLayerOfMatchedCluster = layer; }
+			    }
+			    System.out.println("DEBUG: Matched track with p="+momentum(tr).magnitude()+" to leftoverHitCluster with "+matchedCluster.getCalorimeterHits().size()+" hits at proximity="+primaryDist+" with innermost layer "+innermostLayerOfMatchedCluster);
+			}
+			Map<MCParticle, List<SimCalorimeterHit>> truthInfo = truthFromCluster(matchedCluster);
+			for (MCParticle part : truthInfo.keySet()) {
+			    List<SimCalorimeterHit> truthHits = truthInfo.get(part);
+			    System.out.println("DEBUG:     -> contribution of "+truthHits.size()+" hits from particle with p="+part.getMomentum().magnitude());
+			}
+		    }
+		    Cluster tree = treeOfSharedCluster.get(matchedCluster);
+		    if (tree != null) {
+			List<Cluster> targets = targetsInTree.get(tree);
+			if (m_debugSeedSplitting) {
+			    // Debug printout
+			    for (Cluster target : targets) {
+				String idname = new String("UNKNOWN");
+				if (mipsOld.contains(target)) {
+				    idname = new String("MipOld");
+				} else if (mipsNew.contains(target)) {
+				    idname = new String("MipNew");
+				} else if (clumps.contains(target)) {
+				    idname = new String("Clump");
+				}
+				MCParticle targetTruth = quoteDominantParticle(target);
+				double targetTruthMom = targetTruth.getMomentum().magnitude();
+				System.out.println("DEBUG:     -> structure found: "+idname+" with "+target.getCalorimeterHits().size()+" hits from particle with p="+targetTruthMom);
+			    }
+			}
+			if (targets.size()>0) {
+			    Cluster rematchedSubCluster = m_trackClusterMatcher.matchTrackToCluster(tr, targets);
+			    if (rematchedSubCluster != null) {
+				matchedCluster = rematchedSubCluster; // CHANGES BEHAVIOUR!
+				if (m_debugSeedSplitting) {
+				    // Debug printout
+				    String idname = new String("UNKNOWN");
+				    if (mipsOld.contains(rematchedSubCluster)) {
+					idname = new String("MipOld");
+				    } else if (mipsNew.contains(rematchedSubCluster)) {
+					idname = new String("MipNew");
+				    } else if (clumps.contains(rematchedSubCluster)) {
+					idname = new String("Clump");
+				    }
+				    MCParticle targetTruth = quoteDominantParticle(rematchedSubCluster);
+				    double targetTruthMom = targetTruth.getMomentum().magnitude();
+				    if (rematchedSubCluster != null) {
+					System.out.println("DEBUG:     -> Rematch to: "+idname+" with "+rematchedSubCluster.getCalorimeterHits().size()+" hits from particle with p="+targetTruthMom);
+				    }
+				}
+			    } else {
+				if (m_debugSeedSplitting) {
+				    // Debug printout
+				    System.out.println("DEBUG:     -> Rematch failed (no targets passed track-matching)");
+				    Hep3Vector interceptPoint = findCluster.performExtrapolation(tr);
+				    if (interceptPoint != null) {
+					for (Cluster target : targets) {
+					    double dist = proximity(target, interceptPoint);
+					    int innermostLayerOfTarget = 99;
+					    for (CalorimeterHit hit : target.getCalorimeterHits()) {
+						int layer = getLayer(hit);
+						if (layer < innermostLayerOfTarget) { innermostLayerOfTarget = layer; }
+					    }
+					    System.out.println("DEBUG:           * Target with "+target.getCalorimeterHits().size()+" hits at dist="+dist+" with innermost layer="+innermostLayerOfTarget);
+					}
+				    }
+				}
+			    }
+			}
+		    }
+ 		} else if (treesWithNoStructure.contains(matchedCluster)) {
+		    if (m_debugSeedSplitting) {
+			System.out.println("DEBUG: Matched track with p="+momentum(tr).magnitude()+" to block with "+matchedCluster.getCalorimeterHits().size()+" hits.");
+		    }
+ 		}

 		tracksMatchedToClusters.put(tr, matchedCluster);
 		List<Track> clusTrList = clustersMatchedToTracks.get(matchedCluster);
 		if (clusTrList == null) { 

@@ -281,46 +370,184 @@

 	    }
 	}
 

+	ShowerPointFinder showerFinder = new ShowerPointFinder(findCluster, allHits, tracksMatchedToClusters);
+	Map<Track,BasicCluster> MapTrkToMIP = showerFinder.findMips(); 
+	event.put("ShowerFinderMapTrackToMip", MapTrkToMIP);
+	List<Cluster> preShowerMips = new Vector<Cluster>();
+	preShowerMips.addAll(MapTrkToMIP.values());
+	event.put("ShowerFinderMips", preShowerMips);
+

 	// Figure out whether tracks were uniquely matched or not:
 	Set<Track> uniquelyMatchedTracks = new HashSet<Track>();
 	Set<Track> ambiguouslyMatchedTracks = new HashSet<Track>();

-	List<Set<Track>> ambiguouslyMatchedTrackGroups = new Vector<Set<Track>>();
-	for (Cluster clus : clustersMatchedToTracks.keySet()) {
-	    List<Track> tracksMatchedToThisCluster = clustersMatchedToTracks.get(clus);
-	    if (tracksMatchedToThisCluster.size() > 1) {
-		ambiguouslyMatchedTracks.addAll(tracksMatchedToThisCluster);
-		Set<Track> newSet = new HashSet<Track>();
-		newSet.addAll(tracksMatchedToThisCluster);
-		ambiguouslyMatchedTrackGroups.add(newSet);
-	    } else {
-		uniquelyMatchedTracks.addAll(tracksMatchedToThisCluster);
-	    }
-	}
-	if (uniquelyMatchedTracks.size() + ambiguouslyMatchedTracks.size() != tracksMatchedToClusters.keySet().size()) { throw new AssertionError("Book-keeping error: "+uniquelyMatchedTracks.size()+" + "+ambiguouslyMatchedTracks.size()+" != "+tracksMatchedToClusters.keySet().size()); }
-	// Redo ambiguous tracks hackishly:
-	// Similarly, redo track -> cluster map to use tweaked tracks

 	Map<Track, Cluster> tweakedTracksMatchedToClusters = new HashMap<Track,Cluster>();
 	Map<Cluster, Track> clustersMatchedToTweakedTracks = new HashMap<Cluster, Track>();
 	List<Track> tweakedTracks = new Vector<Track>();
 	Map<Track, Track> mapOrigTrackToTweakedTrack = new HashMap<Track,Track>();

-	tweakedTracks.addAll(uniquelyMatchedTracks);
-	for (Track tr : uniquelyMatchedTracks) {
-	    mapOrigTrackToTweakedTrack.put(tr,tr);
-	}
-	for (Set<Track> trackSet : ambiguouslyMatchedTrackGroups) {
-	    MultipleTrackTrack tweakedTrack = new MultipleTrackTrack(trackSet);
-	    tweakedTracks.add(tweakedTrack);
-	    for (Track origTrack : trackSet) {
-		Cluster clus = tracksMatchedToClusters.get(origTrack);
-		tweakedTracksMatchedToClusters.put(tweakedTrack, clus); // over-writes sometimes, but that's OK.
-		clustersMatchedToTweakedTracks.put(clus, tweakedTrack);
-		mapOrigTrackToTweakedTrack.put(origTrack, tweakedTrack);

+	handleTrackMatchingAmbiguities(clustersMatchedToTracks, tracksMatchedToClusters, uniquelyMatchedTracks, ambiguouslyMatchedTracks, tweakedTracksMatchedToClusters, clustersMatchedToTweakedTracks, tweakedTracks, mapOrigTrackToTweakedTrack);
+
+	// Check for cases where seed cluster is just plain too big
+	List<Track> flaggedBadTracks = new Vector<Track>();
+	for (Track tr : tweakedTracks) {
+	    if (electronTracks.contains(tr)) {
+		continue; // Electron clusters are not too big by construction (E/p ~ 1)
+	    }
+	    Cluster seed = tweakedTracksMatchedToClusters.get(tr);
+	    double seedEnergy = energy(seed);
+	    double tolerance = 3.0;
+	    boolean testValidEoverP = testEoverP_oneSided(seedEnergy, tr, tolerance);
+	    if (!testValidEoverP) {
+		// More than 3sigma off -- broken.
+		flaggedBadTracks.add(tr);

 	    }
 	}

-	for (Track tr : uniquelyMatchedTracks) {
-	    Cluster clus = tracksMatchedToClusters.get(tr);
-	    tweakedTracksMatchedToClusters.put(tr, clus);
-	    clustersMatchedToTweakedTracks.put(clus, tr);

+	for (Track tr : flaggedBadTracks) {
+	    // Track's seed has E>>p -- split it up.
+	    // First find the seed using the parent track...
+	    Cluster seed = tweakedTracksMatchedToClusters.get(tr);
+
+	    // We'll be redoing the matching, so if we have several tracks
+	    // bundled into one, we should split them up:
+	    List<Track> individualTracks = new Vector<Track>();
+	    if (tr instanceof MultipleTrackTrack) {
+		individualTracks.addAll(tr.getTracks());
+	    } else {
+		individualTracks.add(tr);
+	    }
+
+	    // Split the seed up
+	    Map<Track,Cluster> splitOutputMap = new HashMap<Track,Cluster>();
+	    List<Cluster> splitClusters = splitPhoton(seed, individualTracks, splitOutputMap, 10.0, false);
+	    if (m_debugSeedSplitting) {
+		System.out.println("DEBUG: Split big cluster with "+seed.getCalorimeterHits().size()+" hits into pieces...");
+		for (Track tmpTr : splitOutputMap.keySet()) {
+		    double p = (new BasicHep3Vector(tmpTr.getMomentum())).magnitude();
+		    Cluster clus = splitOutputMap.get(tmpTr);
+		    System.out.println(" * Track with p="+p+" mapped to piece with "+clus.getCalorimeterHits().size()+" hits.");
+		}
+	    }
+
+	    // Remove old listings / book-keeping info
+	    mipsOld.remove(seed);
+	    mipsNew.remove(seed);
+	    mips.remove(seed);
+	    clumps.remove(seed);
+	    leftoverHitClusters.remove(seed);
+	    treesWithNoStructure.remove(seed);
+	    clustersMatchedToTracks.remove(seed);
+	    tweakedTracksMatchedToClusters.remove(tr);
+	    clustersMatchedToTweakedTracks.remove(seed);
+	    tweakedTracks.remove(tr);
+	    for (Track dauTr : individualTracks) { 
+		tracksMatchedToClusters.remove(dauTr); 
+		uniquelyMatchedTracks.remove(dauTr);
+		ambiguouslyMatchedTracks.remove(dauTr);
+		mapOrigTrackToTweakedTrack.remove(dauTr);
+	    }
+	    List<Cluster> treesContainingOldSeed = new Vector<Cluster>();
+	    for (Cluster tree : targetsInTree.keySet()) {
+		List<Cluster> targets = targetsInTree.get(tree);
+		if (targets.contains(seed)) {
+		    treesContainingOldSeed.add(tree);
+		    targets.remove(seed);
+		}
+	    }
+
+	    // Some tracks in the bundle may not have a seed cluster any more. They will be jettisoned.
+	    // This is rather brutal -- really we should refine the matching somehow (iteratively?)
+	    // For the rest, see if they're now ambiguous or unique in the matching to newly broken-up seeds
+	    Map<Cluster, List<Track>> local_clustersMatchedToTracks = new HashMap<Cluster, List<Track>>();
+	    Map<Track,Cluster> local_tracksMatchedToClusters = new HashMap<Track,Cluster>();
+	    Set<Track> local_uniquelyMatchedTracks = new HashSet<Track>();
+	    Set<Track> local_ambiguouslyMatchedTracks = new HashSet<Track>();
+	    Map<Track, Cluster> local_tweakedTracksMatchedToClusters = new HashMap<Track, Cluster>();
+	    Map<Cluster, Track> local_clustersMatchedToTweakedTracks = new HashMap<Cluster, Track>();
+	    List<Track> local_tweakedTracks = new Vector<Track>();
+	    Map<Track, Track> local_mapOrigTrackToTweakedTrack = new HashMap<Track, Track>();
+	    for (Track dauTr : individualTracks) {
+		Cluster newSeed = splitOutputMap.get(dauTr);
+		if (newSeed != null) {
+		    List<Track> tracksMatchedToSeed = local_clustersMatchedToTracks.get(newSeed);
+		    if (tracksMatchedToSeed == null) { tracksMatchedToSeed = new Vector<Track>(); local_clustersMatchedToTracks.put(newSeed, tracksMatchedToSeed); }
+		    tracksMatchedToSeed.add(dauTr);
+		    local_tracksMatchedToClusters.put(dauTr, newSeed);
+		}
+	    }
+	    handleTrackMatchingAmbiguities(local_clustersMatchedToTracks, local_tracksMatchedToClusters, local_uniquelyMatchedTracks, local_ambiguouslyMatchedTracks, local_tweakedTracksMatchedToClusters, local_clustersMatchedToTweakedTracks, local_tweakedTracks, local_mapOrigTrackToTweakedTrack);
+
+	    // Book-keeping
+	    Set<Cluster> newSeeds = local_clustersMatchedToTracks.keySet();
+	    Set<Cluster> splitOutputMapSeeds = new HashSet<Cluster>();
+	    splitOutputMapSeeds.addAll(splitOutputMap.values());
+
+	    // Any piece not matched to a track becomes a photon
+	    for (Cluster piece : splitClusters) {
+		if (!newSeeds.contains(piece)) {
+		    chargedHadronLikePhotons.add(piece);
+		    modifiedPhotonClusters.add(piece);
+		}
+	    }
+
+	    // Any track not matched to anything is junked
+	    for (Track subTr : individualTracks) {
+		if (local_tracksMatchedToClusters.get(subTr) == null) {
+		    if (m_debugSeedSplitting) { System.out.println("DEBUG: After splitting cluster, couldn't find a match for track -- junking"); }
+		    unmatchedTracks.add(subTr);
+		}
+	    }
+
+	    // Now, for each track matched to a seed, re-test E/p. If it passes, keep it.
+	    for (Track newTweakedTrack : local_tweakedTracks) {
+		Cluster newSeed = local_tweakedTracksMatchedToClusters.get(newTweakedTrack);
+		if (newSeed == null) { throw new AssertionError("Seed is null, but it shouldn't be -- book-keeping error"); }
+		if (!newSeeds.contains(newSeed)) { throw new AssertionError("New seed is not in newSeeds!"); }
+		if (unmatchedTracks.contains(newTweakedTrack)) { throw new AssertionError("Book-keeping error!"); }
+		double newSeedEnergy = energy(newSeed);
+		double tolerance = 3.0;
+		boolean testValidEoverP_new = testEoverP_oneSided(newSeedEnergy, newTweakedTrack, tolerance);
+
+		if (testValidEoverP_new) {
+		    // REPLACE -- we managed to break off a piece of acceptable size.
+		    if (m_debugSeedSplitting) { System.out.println("DEBUG: Will replace old seed of "+seed.getCalorimeterHits().size()+" hits with new seed of "+newSeed.getCalorimeterHits().size()+" hits plus photon pieces"); }
+		    // Check if this is one track or several bundled together:
+		    List<Track> individualTracksOfNewTweakedTrack = new Vector<Track>();
+		    if (newTweakedTrack instanceof MultipleTrackTrack) {
+			individualTracksOfNewTweakedTrack.addAll(newTweakedTrack.getTracks());
+			ambiguouslyMatchedTracks.addAll(newTweakedTrack.getTracks());
+			for (Track subTr : newTweakedTrack.getTracks()) {
+			    mapOrigTrackToTweakedTrack.put(subTr, newTweakedTrack);
+			}
+		    } else {
+			uniquelyMatchedTracks.add(newTweakedTrack);
+			individualTracksOfNewTweakedTrack.add(newTweakedTrack);
+			mapOrigTrackToTweakedTrack.put(newTweakedTrack, newTweakedTrack);
+		    }
+		    // Add new listings
+		    clumps.add(newSeed);
+		    clustersMatchedToTracks.put(newSeed, individualTracksOfNewTweakedTrack);
+		    for (Track subTr : individualTracksOfNewTweakedTrack) {
+			tracksMatchedToClusters.put(subTr, newSeed);
+		    }
+		    tweakedTracks.add(newTweakedTrack);
+		    tweakedTracksMatchedToClusters.put(newTweakedTrack, newSeed);
+		    clustersMatchedToTweakedTracks.put(newSeed, newTweakedTrack);
+		    for (Cluster tree : treesContainingOldSeed) {
+			List<Cluster> targets = targetsInTree.get(tree);
+			targets.add(newSeed);
+		    }
+		} else {
+		    // FAIL ENTIRELY
+		    if (m_debugSeedSplitting) {
+			System.out.println("DEBUG: Unable to replace old seed of "+seed.getCalorimeterHits().size()+" hits with smaller seed. Will turn into photon and remove track from consideration.");
+			System.out.println("DEBUG: Was working with "+local_tweakedTracks.size()+" local tweaked tracks from "+individualTracks.size()+" individual tracks");
+			if (clustersMatchedToTracks.get(seed) != null) { System.out.println("OLD SEED MATCHED TO NON-NULL TRACK"); }
+		    }
+		    // Add new listings as a photon
+		    chargedHadronLikePhotons.add(newSeed);
+		    modifiedPhotonClusters.add(newSeed);
+		    unmatchedTracks.add(newTweakedTrack);
+		}
+	    }

 	}
 
 	// Track seeds

@@ -585,7 +812,13 @@

 	    sharedLeftoverHitClustersHCAL.rebuildHints();
 	    allSharedClusters.add(sharedLeftoverHitClustersHCAL);
 	} else if (steveMipShareForHaloHCAL) {

-	    SteveMIPReassignmentAlgorithm tmpMipAlg = new SteveMIPReassignmentAlgorithm(event, 1.0);

+	    String mapName;
+	    if (m_useSteveMipsForConeScoring) {
+		mapName = "TrackMipMap";
+	    } else {
+		mapName = "ShowerFinderMapTrackToMip";
+	    }
+	    SteveMIPReassignmentAlgorithm tmpMipAlg = new SteveMIPReassignmentAlgorithm(event, 1.0, mapName);

 	    DownstreamTrackClusterSharingAlgorithm coneSharingAlgHCAL  = new DownstreamTrackClusterSharingAlgorithm(clustersMatchedToTweakedTracks, tmpMipAlg); // TEST
 	    SharedClusterGroup sharedLeftoverHitClustersHCAL = new SharedClusterGroup(leftoverHitClustersToShareHCAL, coneSharingAlgHCAL); // TEST
 	    List<Cluster> tmpSeedList = new Vector<Cluster>(); // TEST

@@ -594,8 +827,6 @@

 	    sharedLeftoverHitClustersHCAL.rebuildHints();
 	    allSharedClusters.add(sharedLeftoverHitClustersHCAL);
 	} else if (taejeongMipShareForHaloHCAL) {

-	    LocalHelixExtrapolator tmpFindCluster = new LocalHelixExtrapolator(); // TEST    
-	    tmpFindCluster.process(m_event); // TEST

 	    Map<Track, Set<Cluster>> tmpMap = new HashMap<Track, Set<Cluster>>(); // TEST
 	    for (Cluster clus : clustersMatchedToTweakedTracks.keySet()) {
 		Track tr = clustersMatchedToTweakedTracks.get(clus);

@@ -606,7 +837,6 @@

 	    }
 	    System.out.println("DEBUG: [seeds] contains "+seeds.size()+" clusters.");
 	    System.out.println("DEBUG: [clustersMatchedToTweakedTracks] contains "+clustersMatchedToTweakedTracks.keySet().size()+" clusters.");

-	    //CachedMIPReassignmentAlgorithm tmpMipAlg = new CachedMIPReassignmentAlgorithm(1.0, tmpFindCluster, tmpMap, mips, allHits); // TEST

 	    MIPReassignmentAlgorithm tmpMipAlg = null; // FIXME: Use actual constructor
 	    DownstreamTrackClusterSharingAlgorithm coneSharingAlgHCAL  = new DownstreamTrackClusterSharingAlgorithm(clustersMatchedToTweakedTracks, tmpMipAlg); // TEST
 	    SharedClusterGroup sharedLeftoverHitClustersHCAL = new SharedClusterGroup(leftoverHitClustersToShareHCAL, coneSharingAlgHCAL); // TEST

@@ -763,20 +993,11 @@

 	    }
 	}
 

-	LocalHelixExtrapolator findCluster = new LocalHelixExtrapolator();
-	findCluster.process(m_event); // picks up geometry
-	ShowerPointFinder showerFinder = new ShowerPointFinder(findCluster, allHits, tracksMatchedToClusters);
-	Map<Track,BasicCluster> MapTrkToMIP = showerFinder.findMips(); 
-	event.put("ShowerFinderMapTrackToMip", MapTrkToMIP);
-	List<Cluster> preShowerMips = new Vector<Cluster>();
-	preShowerMips.addAll(MapTrkToMIP.values());
-	event.put("ShowerFinderMips", preShowerMips);
-

 	ReassignClustersAlgorithm algorithm = null;
 	if (m_useSimpleConeForReassignment) {
 	    algorithm = new ConeReassignmentAlgorithm(1.0, findCluster);
 	} else {

-	    algorithm = new MIPReassignmentAlgorithm(0.79, findCluster, MapTrkToMIP);

+	    algorithm = new MIPReassignmentAlgorithm(1.05, findCluster, MapTrkToMIP, 7);

 	}
 	if (m_fixSingleTracksWithCone) {
 	    // First, try to fix the simplest case: single tracks with E/p < 1

@@ -817,6 +1038,9 @@

 	m_event = null;
     }
 

+    Vector<Double> m_tmp_calib = new Vector<Double>();
+    Vector<Double> m_tmp_calibErr = new Vector<Double>();
+

     boolean checkIfReassignmentNeeded(Set<Track> jet, Set<Cluster> showerComponents, List<SharedClusterGroup> allSharedClusters, double tolerance)
     {
 	double jetMomentum = jetScalarMomentum(jet);

@@ -1327,7 +1551,82 @@

 	}
     }
 

+    private void findTrackSeedInFirstLayersOnly(Track tr, Map<Integer, Set<CalorimeterHit>> unusedHitsByLayer, List<CalorimeterHit> trackSeedHits, Map<CalorimeterHit, Hep3Vector> mapTrackSeedToTangent, Map<Track, CalorimeterHit> mapTrackToTrackSeed, double cutTrackSeedDist) {
+	LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
+	debugTrackMatch.process(m_event);
+	List<Cluster> tmpListLayer0 = new Vector<Cluster>();
+	List<Cluster> tmpListLayer1 = new Vector<Cluster>();
+	if (unusedHitsByLayer.get(0) != null) {
+	    for (CalorimeterHit hit : unusedHitsByLayer.get(0)) {
+		if (!trackSeedHits.contains(hit)) {
+		    BasicCluster tmpClus = new BasicCluster();
+		    tmpClus.addHit(hit);
+		    tmpListLayer0.add(tmpClus);
+		}
+	    }
+	}
+	if (unusedHitsByLayer.get(1) != null) {
+	    for (CalorimeterHit hit : unusedHitsByLayer.get(1)) {
+		if (!trackSeedHits.contains(hit)) {
+		    BasicCluster tmpClus = new BasicCluster();
+		    tmpClus.addHit(hit);
+		    tmpListLayer1.add(tmpClus);
+		}
+	    }
+	}
+	Cluster seedClusToUse = null;
+	Cluster tmpMatchedClusterLayer0 = debugTrackMatch.matchTrackToCluster(tr, tmpListLayer0);
+	Hep3Vector interceptPointLayer0 = debugTrackMatch.getExtrapolator().extendToECALLayer(0);
+	Cluster tmpMatchedClusterLayer1 = debugTrackMatch.matchTrackToCluster(tr, tmpListLayer1);
+	Hep3Vector interceptPointLayer1 = debugTrackMatch.getExtrapolator().extendToECALLayer(1);
+	if (tmpMatchedClusterLayer0 == null && tmpMatchedClusterLayer1 == null) {
+	    // No seed found for this track -- maybe extrapolation is too imprecise.
+	    // Handle it the usual way.
+	} else if (tmpMatchedClusterLayer0 == null && tmpMatchedClusterLayer1 != null) {
+	    seedClusToUse = tmpMatchedClusterLayer1;
+	} else if (tmpMatchedClusterLayer0 != null && tmpMatchedClusterLayer1 == null) {
+	    seedClusToUse = tmpMatchedClusterLayer0;
+	} else {
+	    // Hits found in both layers -- need to pick one.
+	    // Look at distance from track intercept point IN LAYER 0 to hit.
+	    // This favours the hits in the innermost layer... which is what we
+	    // want in general.
+	    Hep3Vector positionOfMatchedClusterLayer0 = new BasicHep3Vector(tmpMatchedClusterLayer0.getCalorimeterHits().iterator().next().getPosition());
+	    Hep3Vector positionOfMatchedClusterLayer1 = new BasicHep3Vector(tmpMatchedClusterLayer1.getCalorimeterHits().iterator().next().getPosition());
+	    double distFromFaceToHit0 = VecOp.sub(positionOfMatchedClusterLayer0, interceptPointLayer0).magnitude();
+	    double distFromFaceToHit1 = VecOp.sub(positionOfMatchedClusterLayer1, interceptPointLayer0).magnitude();
+	    if (distFromFaceToHit0 > distFromFaceToHit1) {
+		seedClusToUse = tmpMatchedClusterLayer1;
+	    } else {
+		seedClusToUse = tmpMatchedClusterLayer0;
+	    }
+	}
+	CalorimeterHit seedHitToUse = null;
+	if (seedClusToUse != null) {
+	    seedHitToUse = seedClusToUse.getCalorimeterHits().iterator().next();
+	    // Require within a certain distance (sanity check)
+	    Hep3Vector positionOfSeedHit = new BasicHep3Vector(seedHitToUse.getPosition());
+	    Hep3Vector trackExtrapPointInLayer = debugTrackMatch.getExtrapolator().extendToECALLayer(getLayer(seedHitToUse));
+	    // It's possible (but rare) for trackExtrapPointInLayer to be null -- e.g. if track just clipped the calorimeter
+	    // and never entered layer 2. Watch for that case.
+	    if (trackExtrapPointInLayer != null) {
+		double distForCut = VecOp.sub(positionOfSeedHit, trackExtrapPointInLayer).magnitude();
+		if (distForCut < cutTrackSeedDist) {
+		    // Within 1cm
+		    trackSeedHits.add(seedHitToUse);
+		    Hep3Vector tangent = VecOp.unit(VecOp.sub(interceptPointLayer1, interceptPointLayer0));
+		    mapTrackSeedToTangent.put(seedHitToUse, tangent);
+		    mapTrackToTrackSeed.put(tr, seedHitToUse);
+		}
+	    }
+	}
+    }
+

     protected List<Cluster> splitPhoton(Cluster photon, List<Track> tracksOfThisPhoton, Map<Track,Cluster> mapTrackToCluster) {

+	return splitPhoton(photon, tracksOfThisPhoton, mapTrackToCluster, 10.0, true);
+    }
+
+    protected List<Cluster> splitPhoton(Cluster photon, List<Track> tracksOfThisPhoton, Map<Track,Cluster> mapTrackToCluster, double cutTrackSeedDist, boolean trackSeedsInFirstLayersOnly) {

 	// Output
 	List<Cluster> output = new Vector<Cluster>();
 	if (mapTrackToCluster == null || mapTrackToCluster.size() != 0) { throw new AssertionError("Must supply initialized, blank map"); }

@@ -1353,77 +1652,42 @@

 	Map<CalorimeterHit, Hep3Vector> mapTrackSeedToTangent = new HashMap<CalorimeterHit, Hep3Vector>();
 	Map<Track, CalorimeterHit> mapTrackToTrackSeed = new HashMap<Track, CalorimeterHit>();
 	for (Track tr : tracksOfThisPhoton) {

-	    LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
-	    debugTrackMatch.process(m_event);
-	    List<Cluster> tmpListLayer0 = new Vector<Cluster>();
-	    List<Cluster> tmpListLayer1 = new Vector<Cluster>();
-	    if (unusedHitsByLayer.get(0) != null) {
-		for (CalorimeterHit hit : unusedHitsByLayer.get(0)) {
-		    if (!trackSeedHits.contains(hit)) {
-			BasicCluster tmpClus = new BasicCluster();
-			tmpClus.addHit(hit);
-			tmpListLayer0.add(tmpClus);
-		    }
-		}
-	    }
-	    if (unusedHitsByLayer.get(1) != null) {
-		for (CalorimeterHit hit : unusedHitsByLayer.get(1)) {
-		    if (!trackSeedHits.contains(hit)) {
-			BasicCluster tmpClus = new BasicCluster();
-			tmpClus.addHit(hit);
-			tmpListLayer1.add(tmpClus);
-		    }
-		}
-	    }
-	    Cluster seedClusToUse = null;
-	    Cluster tmpMatchedClusterLayer0 = debugTrackMatch.matchTrackToCluster(tr, tmpListLayer0);
-	    Hep3Vector interceptPointLayer0 = debugTrackMatch.getExtrapolator().extendToECALLayer(0);
-	    Cluster tmpMatchedClusterLayer1 = debugTrackMatch.matchTrackToCluster(tr, tmpListLayer1);
-	    Hep3Vector interceptPointLayer1 = debugTrackMatch.getExtrapolator().extendToECALLayer(1);
-	    if (tmpMatchedClusterLayer0 == null && tmpMatchedClusterLayer1 == null) {
-		// No seed found for this track -- maybe extrapolation is too imprecise.
-		// Handle it the usual way.
-	    } else if (tmpMatchedClusterLayer0 == null && tmpMatchedClusterLayer1 != null) {
-		seedClusToUse = tmpMatchedClusterLayer1;
-	    } else if (tmpMatchedClusterLayer0 != null && tmpMatchedClusterLayer1 == null) {
-		seedClusToUse = tmpMatchedClusterLayer0;

+	    if (trackSeedsInFirstLayersOnly) {
+		findTrackSeedInFirstLayersOnly(tr, unusedHitsByLayer, trackSeedHits, mapTrackSeedToTangent, mapTrackToTrackSeed, cutTrackSeedDist);

 	    } else {

-		// Hits found in both layers -- need to pick one.
-		// Look at distance from track intercept point IN LAYER 0 to hit.
-		// This favours the hits in the innermost layer... which is what we
-		// want in general.
-		Hep3Vector positionOfMatchedClusterLayer0 = new BasicHep3Vector(tmpMatchedClusterLayer0.getCalorimeterHits().iterator().next().getPosition());
-		Hep3Vector positionOfMatchedClusterLayer1 = new BasicHep3Vector(tmpMatchedClusterLayer1.getCalorimeterHits().iterator().next().getPosition());
-		double distFromFaceToHit0 = VecOp.sub(positionOfMatchedClusterLayer0, interceptPointLayer0).magnitude();
-		double distFromFaceToHit1 = VecOp.sub(positionOfMatchedClusterLayer1, interceptPointLayer0).magnitude();
-		if (distFromFaceToHit0 > distFromFaceToHit1) {
-		    seedClusToUse = tmpMatchedClusterLayer1;
-		} else {
-		    seedClusToUse = tmpMatchedClusterLayer0;
-		}
-	    }
-	    CalorimeterHit seedHitToUse = null;
-	    if (seedClusToUse != null) {
-		seedHitToUse = seedClusToUse.getCalorimeterHits().iterator().next();
-		// Require within a certain distance (sanity check)
-		Hep3Vector positionOfSeedHit = new BasicHep3Vector(seedHitToUse.getPosition());
-		Hep3Vector trackExtrapPointInLayer = debugTrackMatch.getExtrapolator().extendToECALLayer(getLayer(seedHitToUse));
-		// It's possible (but rare) for trackExtrapPointInLayer to be null -- e.g. if track just clipped the calorimeter
-		// and never entered layer 2. Watch for that case.
-		if (trackExtrapPointInLayer != null) {
-		    double distForCut = VecOp.sub(positionOfSeedHit, trackExtrapPointInLayer).magnitude();
-		    if (distForCut < 10.0) {
-		    	// Within 1cm
-		    	trackSeedHits.add(seedHitToUse);
-		    	Hep3Vector tangent = VecOp.unit(VecOp.sub(interceptPointLayer1, interceptPointLayer0));
-		    	mapTrackSeedToTangent.put(seedHitToUse, tangent);
-		    	mapTrackToTrackSeed.put(tr, seedHitToUse);

+		LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
+		debugTrackMatch.process(m_event);
+		for (int iLayer=minLayer; iLayer<maxLayer; iLayer++) {
+		    Set<CalorimeterHit> unusedHitsInLayer = unusedHitsByLayer.get(iLayer);
+		    List<Cluster> tmpClusterList = new Vector<Cluster>();
+		    for (CalorimeterHit hit : unusedHitsInLayer) {
+			if (!trackSeedHits.contains(hit)) {
+			    BasicCluster tmpClus = new BasicCluster();
+			    tmpClus.addHit(hit);
+			    tmpClusterList.add(tmpClus);
+			}
+		    }
+		    Cluster bestClusterMatchInLayer = debugTrackMatch.matchTrackToCluster(tr, tmpClusterList);
+		    Hep3Vector interceptPointInLayer = debugTrackMatch.getExtrapolator().extendToECALLayer(iLayer);
+		    if (bestClusterMatchInLayer != null) {
+			CalorimeterHit seedHit = bestClusterMatchInLayer.getCalorimeterHits().get(0);
+			Hep3Vector positionOfMatchedHit = new BasicHep3Vector(seedHit.getPosition());
+			double transverseDistance = VecOp.sub(positionOfMatchedHit, interceptPointInLayer).magnitude();
+			if (transverseDistance < cutTrackSeedDist) {
+			    // Within 1cm => OK
+			    trackSeedHits.add(seedHit);
+			    Hep3Vector tangent = debugTrackMatch.getExtrapolator().getTangent(interceptPointInLayer, tr);
+			    mapTrackSeedToTangent.put(seedHit, tangent);
+			    mapTrackToTrackSeed.put(tr, seedHit);
+			    break; // stop looping over layers
+			}

 		    }
 		}
 	    }
 	}
 
 	// Start building from track seeds:

+	double dotProductCutForTracks = 0.7; // FIXME: Use 0.85?

 	for (Track tr : tracksOfThisPhoton) {
 	    if (m_photonSplitDebug) { System.out.println("DEBUG: Attempting to build cluster for track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()); }
 	    CalorimeterHit seed = mapTrackToTrackSeed.get(tr);

@@ -1431,15 +1695,20 @@

 		if (unusedHitsByLayer.get(getLayer(seed)).contains(seed)) {
 		    // Not yet used => build a cluster from it
 		    Hep3Vector tangent = mapTrackSeedToTangent.get(seed);

-		    Cluster protoCluster = buildProtoCluster(seed, unusedHitsByLayer, tangent);

+		    Cluster protoCluster = buildProtoCluster(seed, unusedHitsByLayer, tangent, dotProductCutForTracks);

 		    output.add(protoCluster);
 		    mapTrackToCluster.put(tr, protoCluster);
 		    if (m_photonSplitDebug) { debugPrintProtoCluster(protoCluster, photon); }

+		} else {
+		    if (m_photonSplitDebug) { System.out.println("DEBUG: Unable to build cluster for track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" because seed was not unused"); }

 		}

+	    } else {
+		if (m_photonSplitDebug) { System.out.println("DEBUG: Unable to build cluster for track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" because seed was null"); }

 	    }
 	}
 
 	// Start building from non-track seeds:

+	double dotProductCutForNonTracks = 0.7;

 	if (m_photonSplitDebug) { System.out.println("DEBUG: Attempting to build clusters from non-track seeds:"); }
 	for (int iLayer=minLayer; iLayer<=maxLayer; iLayer++) {
 	    Set<CalorimeterHit> unusedHitsInCurrentLayer = unusedHitsByLayer.get(iLayer);

@@ -1447,7 +1716,7 @@

 		CalorimeterHit seed = unusedHitsInCurrentLayer.iterator().next();
 		Hep3Vector seedPosition = new BasicHep3Vector(seed.getPosition());
 		Hep3Vector tangentAtSeed = VecOp.unit(seedPosition);

-		Cluster protoCluster = buildProtoCluster(seed, unusedHitsByLayer, tangentAtSeed);

+		Cluster protoCluster = buildProtoCluster(seed, unusedHitsByLayer, tangentAtSeed, dotProductCutForNonTracks);

 		output.add(protoCluster);
 		if (m_photonSplitDebug) { debugPrintProtoCluster(protoCluster, photon); }
 	    }

@@ -1474,7 +1743,7 @@

 	}
     }
 

-    Cluster buildProtoCluster(CalorimeterHit seed, Map<Integer, Set<CalorimeterHit>> unusedHitsByLayer, Hep3Vector tangentAtSeed) {

+    Cluster buildProtoCluster(CalorimeterHit seed, Map<Integer, Set<CalorimeterHit>> unusedHitsByLayer, Hep3Vector tangentAtSeed, double dotProductCut) {

 	int iLayer = getLayer(seed);
 	Set<CalorimeterHit> protoCluster = new HashSet<CalorimeterHit>();
 	protoCluster.add(seed);

@@ -1493,7 +1762,7 @@

 	unusedHitsByLayer.get(iLayer).removeAll(protoCluster);
 	int currentProtoClusterExpansionLayer = iLayer;
 	while (true) {

-	    Set<CalorimeterHit> extraProtoClusterHits = expandProtoCluster(protoCluster, currentProtoClusterExpansionLayer, unusedHitsByLayer, tangentAtSeed);

+	    Set<CalorimeterHit> extraProtoClusterHits = expandProtoCluster(protoCluster, currentProtoClusterExpansionLayer, unusedHitsByLayer, tangentAtSeed, dotProductCut);

 	    for (CalorimeterHit hit : extraProtoClusterHits) {
 		int hitLayer = getLayer(hit);
 		unusedHitsByLayer.get(hitLayer).remove(hit);

@@ -1512,7 +1781,7 @@

 	return newClus;
     }
 

-    Set<CalorimeterHit> expandProtoCluster(Set<CalorimeterHit> protoCluster, int currentProtoClusterExpansionLayer, Map<Integer, Set<CalorimeterHit>> unusedHitsByLayer, Hep3Vector tangentAtSeed) {

+    Set<CalorimeterHit> expandProtoCluster(Set<CalorimeterHit> protoCluster, int currentProtoClusterExpansionLayer, Map<Integer, Set<CalorimeterHit>> unusedHitsByLayer, Hep3Vector tangentAtSeed, double dotProductCut) {

 	Set<CalorimeterHit> output = new HashSet<CalorimeterHit>();
 	List<CalorimeterHit> protoClusterHitsInLayer = new Vector<CalorimeterHit>();
 	for (CalorimeterHit hit : protoCluster) {

@@ -1525,7 +1794,6 @@

 	    //   a) Within n layers (moving outward)
 	    //   b) Inside a cone of angle theta
 	    int n = 2;

-	    double dotProductCut = 0.7; // Corresponds to theta = 20 degrees

 	    Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
 	    for (int targetLayer = currentProtoClusterExpansionLayer+1; targetLayer <= currentProtoClusterExpansionLayer+n; targetLayer++) {
 		Set<CalorimeterHit> unusedHitsInTargetLayer = unusedHitsByLayer.get(targetLayer);

@@ -2934,10 +3202,13 @@

 	Set<Cluster> allUsedCores = new HashSet<Cluster>();
 	allUsedCores    .addAll(neutralClusterCoresUsed);
 	allUsedCoresList.addAll(neutralClusterCoresUsed);

+	testForHitUniqueness(allUsedCores);

 	if (allUsedCores.size() != allUsedCoresList.size()) { throw new AssertionError("Mis-counting of clusters: "+allUsedCores.size()+" vs "+allUsedCoresList.size()); }
 	allUsedCores    .addAll(photonsUsed);
 	allUsedCoresList.addAll(photonsUsed);
 	if (allUsedCores.size() != allUsedCoresList.size()) { throw new AssertionError("Mis-counting of clusters: "+allUsedCores.size()+" vs "+allUsedCoresList.size()); }

+	testForHitUniqueness(photonsUsed);
+	testForHitUniqueness(allUsedCores);

 	for (Track tr : tracksSortedByMomentum) {
 	    Set<Track> jet = mapTrackToJet.get(tr);
 	    boolean useJetInstead = m_clusterAsJets && jet != null;

@@ -2953,6 +3224,7 @@

 	    }
 	}
 	if (allUsedCores.size() != allUsedCoresList.size()) { throw new AssertionError("Mis-counting of clusters: "+allUsedCores.size()+" vs "+allUsedCoresList.size()); }

+	testForHitUniqueness(allUsedCores);

 	if (m_clusterAsJets) {
 	    for (Set<Track> jet : jets) {
 		for (Cluster clus : newMapJetToShowerComponents.get(jet)) {

@@ -3052,6 +3324,16 @@

 	}
     }
 

+    void testForHitUniqueness(Collection<Cluster> clusters) {
+	Set<CalorimeterHit> hits = new HashSet<CalorimeterHit>();
+	List<CalorimeterHit> hitsList = new Vector<CalorimeterHit>();
+	for (Cluster clus : clusters) {
+	    hits    .addAll(clus.getCalorimeterHits());
+	    hitsList.addAll(clus.getCalorimeterHits());
+	}
+	if (hits.size() != hitsList.size()) { throw new AssertionError("Mis-counting of hits: "+hits.size()+" vs "+hitsList.size()); }
+    }
+

     void clusteringIteration(List<Track> tracksSortedByMomentum, Map<Track, Cluster> tweakedTracksMatchedToClusters, List<Track> electronTracks, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, Map<Cluster, Track> newMapShowerComponentToTrack, Map<Track, Set<Cluster>> newMapTrackToVetoedAdditions, Map<Track, Set<Cluster>> mapTrackToVetoedAdditionsDueToTrackSeed, List<SharedClusterGroup> allSharedClusters, Map<Track,Double> newMapTrackToThreshold, Map<Track,Double> newMapTrackToTolerance) {
 	// Pop the seeds in first:
 	for (Track tr : tracksSortedByMomentum) {

@@ -3589,6 +3871,56 @@

 	m_event.put("DebugMipsRejected", debugList, Cluster.class, flag);
     }
 

-

+    void handleTrackMatchingAmbiguities(Map<Cluster, List<Track>> clustersMatchedToTracks,
+					Map<Track,Cluster> tracksMatchedToClusters,
+					Set<Track> uniquelyMatchedTracks,
+					Set<Track> ambiguouslyMatchedTracks,
+					Map<Track, Cluster> tweakedTracksMatchedToClusters,
+					Map<Cluster, Track> clustersMatchedToTweakedTracks,
+					List<Track> tweakedTracks,
+					Map<Track, Track> mapOrigTrackToTweakedTrack
+					)
+    {
+	if (uniquelyMatchedTracks == null || uniquelyMatchedTracks.size()!=0) { throw new AssertionError("Output collection not in correct format"); }
+	if (ambiguouslyMatchedTracks == null || ambiguouslyMatchedTracks.size()!=0) { throw new AssertionError("Output collection not in correct format"); }
+	if (tweakedTracksMatchedToClusters == null || tweakedTracksMatchedToClusters.size()!=0) { throw new AssertionError("Output collection not in correct format"); }
+	if (clustersMatchedToTweakedTracks == null || clustersMatchedToTweakedTracks.size()!=0) { throw new AssertionError("Output collection not in correct format"); }
+	if (tweakedTracks == null || tweakedTracks.size()!=0) { throw new AssertionError("Output collection not in correct format"); }
+	if (mapOrigTrackToTweakedTrack == null || mapOrigTrackToTweakedTrack.size()!=0)  { throw new AssertionError("Output collection not in correct format"); }
+	List<Set<Track>> ambiguouslyMatchedTrackGroups = new Vector<Set<Track>>();
+	for (Cluster clus : clustersMatchedToTracks.keySet()) {
+	    List<Track> tracksMatchedToThisCluster = clustersMatchedToTracks.get(clus);
+	    if (tracksMatchedToThisCluster.size() > 1) {
+		ambiguouslyMatchedTracks.addAll(tracksMatchedToThisCluster);
+		Set<Track> newSet = new HashSet<Track>();
+		newSet.addAll(tracksMatchedToThisCluster);
+		ambiguouslyMatchedTrackGroups.add(newSet);
+	    } else {
+		uniquelyMatchedTracks.addAll(tracksMatchedToThisCluster);
+	    }
+	}
+	if (uniquelyMatchedTracks.size() + ambiguouslyMatchedTracks.size() != tracksMatchedToClusters.keySet().size()) { throw new AssertionError("Book-keeping error: "+uniquelyMatchedTracks.size()+" + "+ambiguouslyMatchedTracks.size()+" != "+tracksMatchedToClusters.keySet().size()); }	
+	// Redo ambiguous tracks hackishly:
+	// Similarly, redo track -> cluster map to use tweaked tracks
+	tweakedTracks.addAll(uniquelyMatchedTracks);
+	for (Track tr : uniquelyMatchedTracks) {
+	    mapOrigTrackToTweakedTrack.put(tr,tr);
+	}
+		for (Set<Track> trackSet : ambiguouslyMatchedTrackGroups) {
+	    MultipleTrackTrack tweakedTrack = new MultipleTrackTrack(trackSet);
+	    tweakedTracks.add(tweakedTrack);
+	    for (Track origTrack : trackSet) {
+		Cluster clus = tracksMatchedToClusters.get(origTrack);
+		tweakedTracksMatchedToClusters.put(tweakedTrack, clus); // over-writes sometimes, but that's OK.
+		clustersMatchedToTweakedTracks.put(clus, tweakedTrack);
+		mapOrigTrackToTweakedTrack.put(origTrack, tweakedTrack);
+	    }
+	}
+	for (Track tr : uniquelyMatchedTracks) {
+	    Cluster clus = tracksMatchedToClusters.get(tr);
+	    tweakedTracksMatchedToClusters.put(tr, clus);
+	    clustersMatchedToTweakedTracks.put(clus, tr);
+	}
+    }

 }