lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.6 -r1.7
--- ReclusterDriver.java 21 Nov 2007 19:48:37 -0000 1.6
+++ ReclusterDriver.java 21 Nov 2007 22:48:14 -0000 1.7
@@ -30,7 +30,7 @@
*
* This version is PRELIMINARY.
*
- * @version $Id: ReclusterDriver.java,v 1.6 2007/11/21 19:48:37 mcharles Exp $
+ * @version $Id: ReclusterDriver.java,v 1.7 2007/11/21 22:48:14 mcharles Exp $
* @author Mat Charles
*/
@@ -237,11 +237,12 @@
}
// Match tracks to clusters without regard for E/p.
- // We allow them to match the following categories of cluster:
+ // We allow them to match all categories of cluster:
// mips & clumps (from skeletonClusters)
// largeClustersWithoutSkeletons
// photonClusters (could be electrons)
// smallClusters
+ // miscellaneous hits from largeClusters that contain skeletons
List<Cluster> mixedClusters = new Vector<Cluster>();
mixedClusters.addAll(mips);
mixedClusters.addAll(clumps);
@@ -249,24 +250,41 @@
mixedClusters.addAll(photonClusters); // could match [electrons]
List<Cluster> mixedClustersAndHits = new Vector<Cluster>(mixedClusters);
mixedClustersAndHits.addAll(smallClustersWithoutSkeletons);
+ List<Cluster> allMatchableClusters = new Vector<Cluster>();
+ allMatchableClusters.addAll(mixedClusters);
+ allMatchableClusters.addAll(smallClustersWithoutSkeletons);
+ allMatchableClusters.addAll(wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons);
Map<Track,Cluster> tracksMatchedToClusters = new HashMap<Track,Cluster>();
Map<Cluster, List<Track>> clustersMatchedToTracks = new HashMap<Cluster, List<Track>>();
+ boolean tryAllClustersTogether = true;
for (Track tr : trackList) {
- // First try excluding the small clusters, which are often just single hits:
- Cluster matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, mixedClusters);
- if (matchedCluster == null) {
- // Now retry including the small clusters:
- matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, mixedClustersAndHits);
- if (matchedCluster == null) {
- // Now retry with single hits
- matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons);
- if (matchedCluster != null) {
- // Better turn that into a small cluster!
- wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.remove(matchedCluster);
- smallClustersWithoutSkeletons.add(matchedCluster);
+ Cluster matchedCluster = null;
+ if (!tryAllClustersTogether) {
+ // First try excluding the small clusters, which are often just single hits:
+ matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, mixedClusters);
+ if (matchedCluster == null) {
+ // Now retry including the small clusters:
+ matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, mixedClustersAndHits);
+ if (matchedCluster == null) {
+ // Now retry with single hits
+ matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons);
+ if (matchedCluster != null) {
+ // Better turn that into a small cluster!
+ wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.remove(matchedCluster);
+ smallClustersWithoutSkeletons.add(matchedCluster);
+ }
}
}
+ } else {
+ // Try all in one go:
+ matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, allMatchableClusters);
+ if (matchedCluster != null && wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.contains(matchedCluster)) {
+ // Better turn that into a small cluster!
+ wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.remove(matchedCluster);
+ smallClustersWithoutSkeletons.add(matchedCluster);
+ }
}
+ // If we found a match, record it.
if (matchedCluster != null) {
tracksMatchedToClusters.put(tr, matchedCluster);
List<Track> clusTrList = clustersMatchedToTracks.get(matchedCluster);
@@ -407,6 +425,35 @@
printme += "r="+r+", z="+z;
}
+ // Check in case incorrect track match
+ if (maxParticle != mcTruth) {
+ printme += " -- MISTAKE in track-cluster matching! Nearest correct hit is: ";
+ double minDistSq = -1.0;
+ HitMap hitsEcal = ((HitMap)(m_event.get("inputHitMapEcal")));
+ for (CalorimeterHit hit : hitsEcal.values()) {
+ boolean truthMatch = false;
+ SimCalorimeterHit simhit = ((SimCalorimeterHit)(hit));
+ for (int i=0; i<simhit.getMCParticleCount(); i++) {
+ if (simhit.getMCParticle(i) == mcTruth) {
+ truthMatch = true;
+ break;
+ }
+ }
+ if (truthMatch) {
+ double[] hitPos = hit.getPosition();
+ double distSq = (hitPos[0]-point.x())*(hitPos[0]-point.x()) + (hitPos[1]-point.y())*(hitPos[1]-point.y()) + (hitPos[2]-point.z())*(hitPos[2]-point.z());
+ if (distSq < minDistSq || minDistSq < 0.0) { minDistSq = distSq; }
+ }
+ }
+ printme += Math.sqrt(minDistSq)+" away from point (vs "+proximity(clus,point)+" for found cluster). By category, nearby correct clusters are:";
+ double window = 50.0; // somewhat arbitrary
+ for (Cluster tmpclus : mips) { if (proximity(tmpclus,point) < window) { printme += "\n MIP["+tmpclus.getCalorimeterHits().size()+"] at "+proximity(tmpclus,point)+" ;"; if (quoteDominantParticle(tmpclus)==mcTruth) { printme += " (TM)"; } else { printme += " (not TM)"; } } }
+ for (Cluster tmpclus : clumps) { if (proximity(tmpclus,point) < window) { printme += "\n Clump["+tmpclus.getCalorimeterHits().size()+"] at "+proximity(tmpclus,point)+" ;"; if (quoteDominantParticle(tmpclus)==mcTruth) { printme += " (TM)"; } else { printme += " (not TM)"; } } }
+ for (Cluster tmpclus : largeClustersWithoutSkeletons) { if (proximity(tmpclus,point) < window) { printme += "\n Large["+tmpclus.getCalorimeterHits().size()+"] at "+proximity(tmpclus,point)+" ;"; if (quoteDominantParticle(tmpclus)==mcTruth) { printme += " (TM)"; } else { printme += " (not TM)"; } } }
+ for (Cluster tmpclus : smallClustersWithoutSkeletons) { if (proximity(tmpclus,point) < window) { printme += "\n Small["+tmpclus.getCalorimeterHits().size()+"] at "+proximity(tmpclus,point)+" ;"; if (quoteDominantParticle(tmpclus)==mcTruth) { printme += " (TM)"; } else { printme += " (not TM)"; } } }
+ for (Cluster tmpclus : wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons) { if (proximity(tmpclus,point) < window) { printme += "\n HitFromLarge["+tmpclus.getCalorimeterHits().size()+"] at "+proximity(tmpclus,point)+" ;"; if (quoteDominantParticle(tmpclus)==mcTruth) { printme += " (TM)"; } else { printme += " (not TM)"; } } }
+ }
+
System.out.println(printme);
}
@@ -433,25 +480,6 @@
}
Collections.sort(tracksSortedByMomentum, new MomentumSort());
- // Sanity check: Total cluster energy in event
- {
- Set<Cluster> allClustersInEvent = new HashSet<Cluster>();
- Set<CalorimeterHit> allHitsInEvent = new HashSet<CalorimeterHit>();
- for (ReconstructedParticle p : muonParticles) { allClustersInEvent.addAll(p.getClusters()); }
- allClustersInEvent.addAll(photonClusters);
- allClustersInEvent.addAll(skeletonClusters);
- allClustersInEvent.addAll(smallClustersWithoutSkeletons);
- allClustersInEvent.addAll(largeClusters);
- allClustersInEvent.addAll(mips);
- allClustersInEvent.addAll(clumps);
- allClustersInEvent.addAll(splitSkeletonClusters);
- for (Cluster c : allClustersInEvent) { allHitsInEvent.addAll(c.getCalorimeterHits()); }
- allHitsInEvent.addAll(unusedHits.values());
- BasicCluster tmpClus = new BasicCluster();
- for (CalorimeterHit h : allHitsInEvent) { tmpClus.addHit(h); }
- System.out.println("DEBUG: Total calorimeter energy in event = "+energy(tmpClus));
- }
-
////////////////////////////////////////////////////////////////////////
List<SharedClusterGroup> allSharedClusters = new Vector<SharedClusterGroup>();
// We can't share seeds... on occasion, a small cluster will be a seed.
@@ -1120,6 +1148,7 @@
)
{
System.out.println(desc);
+ double totalChiSquared = 0.0;
for (Track tr : tracksSortedByMomentum) {
Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
double clusterEnergy = energy(showerComponents, allSharedClusters);
@@ -1137,7 +1166,10 @@
+", E/p is "+clusterEnergy+" / "+trackMomentum+" => NormResid = "+normalizedResidual
+". Core: effic="+coreEffic+", purity="+corePurity
+". Real: effic="+realEffic+", purity="+realPurity);
+ double trackChiSquared = normalizedResidual * normalizedResidual;
+ totalChiSquared += trackChiSquared;
}
+ System.out.println("Total CHI^2 = "+totalChiSquared+" with NDF ~ "+tracksSortedByMomentum.size());
for (Cluster clus : photonClusters) {
debugPrint("Photon", clus,tracksSortedByMomentum,newMapTrackToShowerComponents,newMapShowerComponentToTrack,allSharedClusters,newMapTrackToVetoedAdditions);
debugPrint(clus, newMapShowerComponentToTrack, photonClusters, mips, clumps,largeClustersWithoutSkeletons);
@@ -1200,6 +1232,20 @@
return output;
}
+ protected double proximity(Cluster clus, Hep3Vector point) {
+ if (clus.getCalorimeterHits().size()<1) { throw new AssertionError("Empty cluster"); }
+ double minDist = 0;
+ boolean found = false;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ if (distance<minDist || found==false) {
+ found = true;
+ minDist = distance;
+ }
+ }
+ return minDist;
+ }
protected double proximity(Cluster clus1, Cluster clus2) {
if (clus1.getCalorimeterHits().size()<1) { throw new AssertionError("Empty cluster"); }
if (clus2.getCalorimeterHits().size()<1) { throw new AssertionError("Empty cluster"); }
@@ -1207,13 +1253,10 @@
boolean found = false;
for (CalorimeterHit hit1 : clus1.getCalorimeterHits()) {
Hep3Vector hitPosition1 = new BasicHep3Vector(hit1.getPosition());
- for (CalorimeterHit hit2 : clus2.getCalorimeterHits()) {
- Hep3Vector hitPosition2 = new BasicHep3Vector(hit2.getPosition());
- double distance = VecOp.sub(hitPosition1,hitPosition2).magnitude();
- if (distance<minDist || found==false) {
- found = true;
- minDist = distance;
- }
+ double currentHitProximity = proximity(clus2, hitPosition1);
+ if (currentHitProximity<minDist || found==false) {
+ found = true;
+ minDist = currentHitProximity;
}
}
return minDist;
