lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.7 -r1.8
--- ReclusterDriver.java 21 Nov 2007 22:48:14 -0000 1.7
+++ ReclusterDriver.java 26 Nov 2007 00:26:46 -0000 1.8
@@ -30,7 +30,7 @@
*
* This version is PRELIMINARY.
*
- * @version $Id: ReclusterDriver.java,v 1.7 2007/11/21 22:48:14 mcharles Exp $
+ * @version $Id: ReclusterDriver.java,v 1.8 2007/11/26 00:26:46 mcharles Exp $
* @author Mat Charles
*/
@@ -39,9 +39,9 @@
List<Cluster> m_mips;
boolean m_megaDebug = false;
- boolean m_debug = true;
+ boolean m_debug = false;
boolean m_debugTiming = false;
- boolean m_debugLinkScores = true;
+ boolean m_debugLinkScores = false;
String m_outputParticleListName = "ReclusteredParticles";
String m_mcList;
String m_inputSmallClusters;
@@ -56,6 +56,7 @@
String m_inputSplitSkeletonClusters;
EventHeader m_event;
protected LikelihoodEvaluator m_eval = null;
+ protected boolean m_cheatOnScoring = false;
public ReclusterDriver(String mcList, String trackList, String muonParticles, String photonClusters, String skeletonClusters, String smallClusters, String unusedHits, String largeClusters, String mips, String clumps, String splitSkeletonClusters, LikelihoodEvaluator eval) {
m_eval = eval;
@@ -102,6 +103,11 @@
// plots
initPlots();
+
+ // Print a warning message if cheating
+ if (m_cheatOnScoring) {
+ System.out.println("WARNING: ReclusterDriver cheating on scoring");
+ }
}
public void process(EventHeader event) {
@@ -122,15 +128,17 @@
m_mips = mips;
- System.out.println("DEBUG: ReclusterDriver read in "+trackList.size()+" tracks");
- System.out.println("DEBUG: ReclusterDriver read in "+muonParticles.size()+" muons");
- System.out.println("DEBUG: ReclusterDriver read in "+photonClusters.size()+" photons");
- System.out.println("DEBUG: ReclusterDriver read in "+skeletonClusters.size()+" skeletons");
- System.out.println("DEBUG: ReclusterDriver read in "+largeClusters.size()+" large clusters");
- System.out.println("DEBUG: ReclusterDriver read in "+mips.size()+" MIP clusters");
- System.out.println("DEBUG: ReclusterDriver read in "+clumps.size()+" clump clusters");
- System.out.println("DEBUG: ReclusterDriver read in "+unusedHits.values().size()+" unused hits");
- System.out.println("DEBUG: ReclusterDriver read in "+splitSkeletonClusters.size()+" pre-split skeletons");
+ if (m_debug) {
+ System.out.println("DEBUG: ReclusterDriver read in "+trackList.size()+" tracks");
+ System.out.println("DEBUG: ReclusterDriver read in "+muonParticles.size()+" muons");
+ System.out.println("DEBUG: ReclusterDriver read in "+photonClusters.size()+" photons");
+ System.out.println("DEBUG: ReclusterDriver read in "+skeletonClusters.size()+" skeletons");
+ System.out.println("DEBUG: ReclusterDriver read in "+largeClusters.size()+" large clusters");
+ System.out.println("DEBUG: ReclusterDriver read in "+mips.size()+" MIP clusters");
+ System.out.println("DEBUG: ReclusterDriver read in "+clumps.size()+" clump clusters");
+ System.out.println("DEBUG: ReclusterDriver read in "+unusedHits.values().size()+" unused hits");
+ System.out.println("DEBUG: ReclusterDriver read in "+splitSkeletonClusters.size()+" pre-split skeletons");
+ }
Map<Cluster,Cluster> mapMIPToSkeleton = new HashMap<Cluster,Cluster>();
Map<Cluster,Cluster> mapClumpToSkeleton = new HashMap<Cluster,Cluster>();
@@ -160,7 +168,7 @@
}
- {
+ if (m_debug) {
int countSmallClusterHits = 0;
for (Cluster clus : smallClusters) {
countSmallClusterHits += clus.getCalorimeterHits().size();
@@ -213,7 +221,7 @@
}
}
- {
+ if (m_debug) {
int countLargeHitsUnused = 0;
for (Cluster clus : largeClustersWithoutSkeletons) {
countLargeHitsUnused += clus.getCalorimeterHits().size();
@@ -303,158 +311,162 @@
unmatchedTracks.removeAll(tracksMatchedToClusters.keySet());
// Some debug printout
- System.out.println("DEBUG: "+unmatchedTracks.size()+" unmatched tracks remaining:");
- for (Track tr : unmatchedTracks) {
- Hep3Vector p = new BasicHep3Vector(tr.getMomentum());
- double px = p.x();
- double py = p.y();
- double pz = p.z();
- double pmag = p.magnitude();
- double pt = Math.sqrt(px*px + py*py);
- double cosTheta = pz/pmag;
- String printme = new String("Track with p="+p.magnitude()+" (cosTheta="+cosTheta+")");
- LocalHelixExtrapolationTrackClusterMatcher extrap = new LocalHelixExtrapolationTrackClusterMatcher();
- extrap.process(m_event);
- Hep3Vector point = extrap.performExtrapolation(tr);
- if (point == null) {
- printme += " -- extrapolation point is null";
- } else {
- double x = point.x();
- double y = point.y();
- double z = point.z();
- double r = Math.sqrt(x*x + y*y);
- printme += " -- extrapolation point at x="+point.x()+", y="+point.y()+", z="+point.z()+", r="+r;
- }
- System.out.println(printme);
- if (point != null) {
- System.out.println("Dumping ECAL hits from this particle:");
- MCParticle trackTruth = ((BaseTrackMC)(tr)).getMCParticle();
- List<CalorimeterHit> EcalBarrDigiHits = event.get(CalorimeterHit.class, "EcalBarrDigiHits");
- List<CalorimeterHit> EcalEndcapDigiHits = event.get(CalorimeterHit.class, "EcalEndcapDigiHits");
- List<SimCalorimeterHit> EcalHits = new Vector<SimCalorimeterHit>();
- for (CalorimeterHit hit : EcalBarrDigiHits) { EcalHits.add((SimCalorimeterHit)(hit)); }
- for (CalorimeterHit hit : EcalEndcapDigiHits) { EcalHits.add((SimCalorimeterHit)(hit)); }
- for (SimCalorimeterHit hit : EcalHits) {
- int mcCount = hit.getMCParticleCount();
- boolean matchedTruth = false;
- for (int i=0; i<mcCount; i++) {
- MCParticle mc = hit.getMCParticle(i);
- if (mc == trackTruth) {
- matchedTruth = true;
- break;
+ if (m_debug) {
+ System.out.println("DEBUG: "+unmatchedTracks.size()+" unmatched tracks remaining:");
+ for (Track tr : unmatchedTracks) {
+ Hep3Vector p = new BasicHep3Vector(tr.getMomentum());
+ double px = p.x();
+ double py = p.y();
+ double pz = p.z();
+ double pmag = p.magnitude();
+ double pt = Math.sqrt(px*px + py*py);
+ double cosTheta = pz/pmag;
+ String printme = new String("Track with p="+p.magnitude()+" (cosTheta="+cosTheta+")");
+ LocalHelixExtrapolationTrackClusterMatcher extrap = new LocalHelixExtrapolationTrackClusterMatcher();
+ extrap.process(m_event);
+ Hep3Vector point = extrap.performExtrapolation(tr);
+ if (point == null) {
+ printme += " -- extrapolation point is null";
+ } else {
+ double x = point.x();
+ double y = point.y();
+ double z = point.z();
+ double r = Math.sqrt(x*x + y*y);
+ printme += " -- extrapolation point at x="+point.x()+", y="+point.y()+", z="+point.z()+", r="+r;
+ }
+ System.out.println(printme);
+ if (point != null) {
+ System.out.println("Dumping ECAL hits from this particle:");
+ MCParticle trackTruth = ((BaseTrackMC)(tr)).getMCParticle();
+ List<CalorimeterHit> EcalBarrDigiHits = event.get(CalorimeterHit.class, "EcalBarrDigiHits");
+ List<CalorimeterHit> EcalEndcapDigiHits = event.get(CalorimeterHit.class, "EcalEndcapDigiHits");
+ List<SimCalorimeterHit> EcalHits = new Vector<SimCalorimeterHit>();
+ for (CalorimeterHit hit : EcalBarrDigiHits) { EcalHits.add((SimCalorimeterHit)(hit)); }
+ for (CalorimeterHit hit : EcalEndcapDigiHits) { EcalHits.add((SimCalorimeterHit)(hit)); }
+ for (SimCalorimeterHit hit : EcalHits) {
+ int mcCount = hit.getMCParticleCount();
+ boolean matchedTruth = false;
+ for (int i=0; i<mcCount; i++) {
+ MCParticle mc = hit.getMCParticle(i);
+ if (mc == trackTruth) {
+ matchedTruth = true;
+ break;
+ }
+ }
+ Hep3Vector hitPos = new BasicHep3Vector(hit.getPosition());
+ Hep3Vector displacement = VecOp.sub(hitPos, point);
+ double r = Math.sqrt(hitPos.x()*hitPos.x() + hitPos.y()*hitPos.y());
+ double separation = displacement.magnitude();
+ if (separation < 50.0) {
+ String collection = "Collection=";
+ for (Cluster clus : mips) { if (clus.getCalorimeterHits().contains(hit)) { collection += "MIPs,"; break; } }
+ for (Cluster clus : clumps) { if (clus.getCalorimeterHits().contains(hit)) { collection += "Clumps,"; break; } }
+ for (Cluster clus : photonClusters) { if (clus.getCalorimeterHits().contains(hit)) { collection += "Photons,"; break; } }
+ for (Cluster clus : largeClustersWithoutSkeletons) { if (clus.getCalorimeterHits().contains(hit)) { collection += "LargeClusWithoutSkel,"; break; } }
+ for (Cluster clus : smallClustersWithoutSkeletons) { if (clus.getCalorimeterHits().contains(hit)) { collection += "SmallClusWithoutSkel,"; break; } }
+ if (unusedHits.values().contains(hit)) { collection += "UnusedHits,"; }
+ for (Cluster clus : largeClusters) { if (clus.getCalorimeterHits().contains(hit)) { collection += "LargeClus,"; break; } }
+ for (Cluster clus : skeletonClusters) { if (clus.getCalorimeterHits().contains(hit)) { collection += "Skel,"; break; } }
+ System.out.println(" hit at r="+r+", z="+hitPos.z()+", x="+hitPos.x()+", y="+hitPos.y()+" with separation="+displacement.magnitude()+" -- in list "+collection);
}
- }
- Hep3Vector hitPos = new BasicHep3Vector(hit.getPosition());
- Hep3Vector displacement = VecOp.sub(hitPos, point);
- double r = Math.sqrt(hitPos.x()*hitPos.x() + hitPos.y()*hitPos.y());
- double separation = displacement.magnitude();
- if (separation < 50.0) {
- String collection = "Collection=";
- for (Cluster clus : mips) { if (clus.getCalorimeterHits().contains(hit)) { collection += "MIPs,"; break; } }
- for (Cluster clus : clumps) { if (clus.getCalorimeterHits().contains(hit)) { collection += "Clumps,"; break; } }
- for (Cluster clus : photonClusters) { if (clus.getCalorimeterHits().contains(hit)) { collection += "Photons,"; break; } }
- for (Cluster clus : largeClustersWithoutSkeletons) { if (clus.getCalorimeterHits().contains(hit)) { collection += "LargeClusWithoutSkel,"; break; } }
- for (Cluster clus : smallClustersWithoutSkeletons) { if (clus.getCalorimeterHits().contains(hit)) { collection += "SmallClusWithoutSkel,"; break; } }
- if (unusedHits.values().contains(hit)) { collection += "UnusedHits,"; }
- for (Cluster clus : largeClusters) { if (clus.getCalorimeterHits().contains(hit)) { collection += "LargeClus,"; break; } }
- for (Cluster clus : skeletonClusters) { if (clus.getCalorimeterHits().contains(hit)) { collection += "Skel,"; break; } }
- System.out.println(" hit at r="+r+", z="+hitPos.z()+", x="+hitPos.x()+", y="+hitPos.y()+" with separation="+displacement.magnitude()+" -- in list "+collection);
}
}
}
- }
- System.out.println("DEBUG: Here are the "+tracksMatchedToClusters.keySet().size()+" matched tracks");
- for (Track tr : tracksMatchedToClusters.keySet()) {
- Cluster clus = tracksMatchedToClusters.get(tr);
- double[] p = tr.getMomentum();
- Particle mcTruth = null;
- if (tr instanceof ReconTrack) {
- mcTruth = ((ReconTrack)(tr)).getMCParticle();
- }
- if (tr instanceof BaseTrackMC) {
- mcTruth = ((BaseTrackMC)(tr)).getMCParticle();
- }
- double mom = Math.sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
- // Which list is that from?
- String printme = new String("Matched track with p=");
- printme += mom;
- if (mcTruth != null) {
- printme += " (truth type "+mcTruth.getType()+")";
- }
- printme += " to a cluster with ";
- printme += clus.getCalorimeterHits().size();
- printme += " hits and energy of ";
- printme += clus.getEnergy();
- if (mips.contains(clus)) { printme += " [from MIPS]"; }
- if (clumps.contains(clus)) { printme += " [from CLUMPS]"; }
- if (largeClustersWithoutSkeletons.contains(clus)) { printme += " [from LARGE CLUSTERS W/O SKELETONS]"; }
- if (photonClusters.contains(clus)) { printme += " [from PHOTONS]"; }
- if (smallClustersWithoutSkeletons.contains(clus)) { printme += " [from SMALL CLUSTERS]"; }
- Map<MCParticle, List<CalorimeterHit>> clusterTruth = new HashMap<MCParticle, List<CalorimeterHit>>();
- for (CalorimeterHit hit : clus.getCalorimeterHits()) {
- SimCalorimeterHit simHit = (SimCalorimeterHit)(hit);
- MCParticle mcTruthPart = simHit.getMCParticle(0);
- List<CalorimeterHit> hitList = clusterTruth.get(mcTruthPart);
- if (hitList == null) {
- hitList = new Vector<CalorimeterHit>();
- clusterTruth.put(mcTruthPart, hitList);
- }
- hitList.add(hit);
- }
- MCParticle maxParticle = null;
- for (MCParticle part : clusterTruth.keySet()) {
- if (maxParticle == null || clusterTruth.get(maxParticle).size() < clusterTruth.get(part).size()) {
- maxParticle = part;
- }
- }
- printme += " (dominant particle: "+maxParticle.getType().getPDGID()+" with "+clusterTruth.get(maxParticle).size()+" hits)";
-
- // DEBUG
- LocalHelixExtrapolationTrackClusterMatcher extrap = new LocalHelixExtrapolationTrackClusterMatcher();
- extrap.process(m_event);
- Hep3Vector point = extrap.performExtrapolation(tr);
- printme += " -- extrap point is ";
- if (point == null) {
- printme += "null";
- } else {
- double x = point.x();
- double y = point.y();
- double z = point.z();
- double r = Math.sqrt(x*x + y*y);
- 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;
+ System.out.println("DEBUG: Here are the "+tracksMatchedToClusters.keySet().size()+" matched tracks");
+ for (Track tr : tracksMatchedToClusters.keySet()) {
+ Cluster clus = tracksMatchedToClusters.get(tr);
+ double[] p = tr.getMomentum();
+ Particle mcTruth = null;
+ if (tr instanceof ReconTrack) {
+ mcTruth = ((ReconTrack)(tr)).getMCParticle();
+ }
+ if (tr instanceof BaseTrackMC) {
+ mcTruth = ((BaseTrackMC)(tr)).getMCParticle();
+ }
+ double mom = Math.sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
+ // Which list is that from?
+ String printme = new String("Matched track with p=");
+ printme += mom;
+ if (mcTruth != null) {
+ printme += " (truth type "+mcTruth.getType()+")";
+ }
+ printme += " to a cluster with ";
+ printme += clus.getCalorimeterHits().size();
+ printme += " hits and energy of ";
+ printme += clus.getEnergy();
+ if (mips.contains(clus)) { printme += " [from MIPS]"; }
+ if (clumps.contains(clus)) { printme += " [from CLUMPS]"; }
+ if (largeClustersWithoutSkeletons.contains(clus)) { printme += " [from LARGE CLUSTERS W/O SKELETONS]"; }
+ if (photonClusters.contains(clus)) { printme += " [from PHOTONS]"; }
+ if (smallClustersWithoutSkeletons.contains(clus)) { printme += " [from SMALL CLUSTERS]"; }
+ Map<MCParticle, List<CalorimeterHit>> clusterTruth = new HashMap<MCParticle, List<CalorimeterHit>>();
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ SimCalorimeterHit simHit = (SimCalorimeterHit)(hit);
+ MCParticle mcTruthPart = simHit.getMCParticle(0);
+ List<CalorimeterHit> hitList = clusterTruth.get(mcTruthPart);
+ if (hitList == null) {
+ hitList = new Vector<CalorimeterHit>();
+ clusterTruth.put(mcTruthPart, hitList);
+ }
+ hitList.add(hit);
+ }
+ MCParticle maxParticle = null;
+ for (MCParticle part : clusterTruth.keySet()) {
+ if (maxParticle == null || clusterTruth.get(maxParticle).size() < clusterTruth.get(part).size()) {
+ maxParticle = part;
+ }
+ }
+ printme += " (dominant particle: "+maxParticle.getType().getPDGID()+" with "+clusterTruth.get(maxParticle).size()+" hits)";
+
+ LocalHelixExtrapolationTrackClusterMatcher extrap = new LocalHelixExtrapolationTrackClusterMatcher();
+ extrap.process(m_event);
+ Hep3Vector point = extrap.performExtrapolation(tr);
+ printme += " -- extrap point is ";
+ if (point == null) {
+ printme += "null";
+ } else {
+ double x = point.x();
+ double y = point.y();
+ double z = point.z();
+ double r = Math.sqrt(x*x + y*y);
+ printme += "r="+r+", z="+z;
+ }
+
+ // Check in case incorrect track match
+ if (maxParticle != mcTruth) {
+ printme += " -- MISTAKE in track-cluster matching!";
+ if (point != null) {
+ printme += " 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; }
+ }
}
- }
- 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)"; } } }
}
}
- 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);
+ System.out.println(printme);
+ }
}
// OK. Now pick out those clusters that have something connected to them...
@@ -498,6 +510,7 @@
linkableClusters.addAll(clumps);
linkableClusters.addAll(photonClusters);
linkableClusters.addAll(largeClustersWithoutSkeletons);
+
// Share halo
ProximityClusterSharingAlgorithm proximityAlgForHalo = new ProximityClusterSharingAlgorithm(20.0, 100.0);
SharedClusterGroup sharedHaloHits = new SharedClusterGroup(wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons, proximityAlgForHalo);
@@ -508,12 +521,15 @@
sharedHaloHits.createShares(largeClustersWithoutSkeletons);
sharedHaloHits.rebuildHints();
allSharedClusters.add(sharedHaloHits);
- System.out.println("Checking for halo sharing... "+wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.size()+" hits -> "+sharedHaloHits.listAllSharedClusters().size()+" SharedCluster objects");
- //debugPrint(smallClusterSeeds, sharedHaloHits, "Small", "halo");
- //debugPrint(mips, sharedHaloHits, "MIP", "halo");
- //debugPrint(clumps, sharedHaloHits, "Clump", "halo");
- //debugPrint(photonClusters, sharedHaloHits, "Photon", "halo");
- //debugPrint(largeClustersWithoutSkeletons, sharedHaloHits, "Large", "halo");
+ if (m_debug) {
+ System.out.println("Checking for halo sharing... "+wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.size()+" hits -> "+sharedHaloHits.listAllSharedClusters().size()+" SharedCluster objects");
+ //debugPrint(smallClusterSeeds, sharedHaloHits, "Small", "halo");
+ //debugPrint(mips, sharedHaloHits, "MIP", "halo");
+ //debugPrint(clumps, sharedHaloHits, "Clump", "halo");
+ //debugPrint(photonClusters, sharedHaloHits, "Photon", "halo");
+ //debugPrint(largeClustersWithoutSkeletons, sharedHaloHits, "Large", "halo");
+ }
+
// Share small clusters
ProximityClusterSharingAlgorithm proximityAlgForSmallClusters = new ProximityClusterSharingAlgorithm(40.0, 250.0);
SharedClusterGroup sharedSmallClusters = new SharedClusterGroup(smallClustersExcludingSeedsEtc, proximityAlgForSmallClusters);
@@ -524,12 +540,14 @@
sharedSmallClusters.createShares(largeClustersWithoutSkeletons);
sharedSmallClusters.rebuildHints();
allSharedClusters.add(sharedSmallClusters);
- System.out.println("Checking for small cluster sharing... "+smallClustersExcludingSeedsEtc.size()+" clusters -> "+sharedSmallClusters.listAllSharedClusters().size()+" SharedCluster objects");
- //debugPrint(smallClusterSeeds, sharedSmallClusters, "Small", "halo");
- //debugPrint(mips, sharedSmallClusters, "MIP", "halo");
- //debugPrint(clumps, sharedSmallClusters, "Clump", "halo");
- //debugPrint(photonClusters, sharedSmallClusters, "Photon", "halo");
- //debugPrint(largeClustersWithoutSkeletons, sharedSmallClusters, "Large", "halo");
+ if (m_debug) {
+ System.out.println("Checking for small cluster sharing... "+smallClustersExcludingSeedsEtc.size()+" clusters -> "+sharedSmallClusters.listAllSharedClusters().size()+" SharedCluster objects");
+ //debugPrint(smallClusterSeeds, sharedSmallClusters, "Small", "halo");
+ //debugPrint(mips, sharedSmallClusters, "MIP", "halo");
+ //debugPrint(clumps, sharedSmallClusters, "Clump", "halo");
+ //debugPrint(photonClusters, sharedSmallClusters, "Photon", "halo");
+ //debugPrint(largeClustersWithoutSkeletons, sharedSmallClusters, "Large", "halo");
+ }
// Cache potential links (warning, there may be a lot!)
resetPotentialLinks();
@@ -540,9 +558,7 @@
double thresholdForProximity = 50.0; // 5cm. FIXME: Don't hard-code.
double thresholdForProximityClump = 75.0; // FIXME: Don't hard-code.
- boolean cheatOnScoring = false;
-
- if (cheatOnScoring) {
+ if (m_cheatOnScoring) {
initPotentialLinks_cheating(linkableClusters);
} else {
// Links between main components (clumps, mips) within a large cluster
@@ -781,7 +797,7 @@
}
- printStatus("FINAL STATUS:", tracksSortedByMomentum, allSharedClusters, newMapTrackToShowerComponents, newMapShowerComponentToTrack, newMapTrackToThreshold, newMapTrackToTolerance, photonClusters, mips, clumps, largeClustersWithoutSkeletons, smallClusterSeeds, newMapTrackToVetoedAdditions);
+ if (m_debug) { printStatus("FINAL STATUS:", tracksSortedByMomentum, allSharedClusters, newMapTrackToShowerComponents, newMapShowerComponentToTrack, newMapTrackToThreshold, newMapTrackToTolerance, photonClusters, mips, clumps, largeClustersWithoutSkeletons, smallClusterSeeds, newMapTrackToVetoedAdditions); }
long timeAfterFinalStatusPrintout = Calendar.getInstance().getTimeInMillis(); // DEBUG
if (m_debugTiming) { System.out.println("DEBUG: Time to print out final status took "+(timeAfterFinalStatusPrintout-timeAtEndOfLastEventIteration)+" ms"); }
////////////////////////////////////////////////////////////////////////
@@ -855,16 +871,18 @@
// Add to the output list
outputParticleList.add(part);
outputPhotonParticleList.add(part);
- // Debug printout
- MCParticle dominantParticle = quoteDominantParticle(clus);
- double truthMomentum = dominantParticle.getMomentum().magnitude();
- int truthPDG = dominantParticle.getPDGID();
- double effic = quoteEfficiency_P(dominantParticle, clus, allSharedClusters);
- double purity = quotePurity_P(dominantParticle, clus, allSharedClusters);
- double energyNoShare = m_photonCalib.getEnergy(clus);
- double efficNoShare = quoteEfficiency_P(dominantParticle, clus);
- double purityNoShare = quotePurity_P(dominantParticle, clus);
- System.out.println("DEBUG: Looking at photon cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity+". Excluding shared hits, E="+energyNoShare+", effic="+efficNoShare+", purity="+purityNoShare);
+ if (m_debug) {
+ // Debug printout
+ MCParticle dominantParticle = quoteDominantParticle(clus);
+ double truthMomentum = dominantParticle.getMomentum().magnitude();
+ int truthPDG = dominantParticle.getPDGID();
+ double effic = quoteEfficiency_P(dominantParticle, clus, allSharedClusters);
+ double purity = quotePurity_P(dominantParticle, clus, allSharedClusters);
+ double energyNoShare = m_photonCalib.getEnergy(clus);
+ double efficNoShare = quoteEfficiency_P(dominantParticle, clus);
+ double purityNoShare = quotePurity_P(dominantParticle, clus);
+ System.out.println("DEBUG: Looking at photon cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity+". Excluding shared hits, E="+energyNoShare+", effic="+efficNoShare+", purity="+purityNoShare);
+ }
}
}
@@ -951,8 +969,8 @@
// Add to the output list
outputParticleList.add(part);
outputNeutralHadronParticleList.add(part);
- // Debug printout
- {
+ if (m_debug) {
+ // Debug printout
MCParticle dominantParticle = quoteDominantParticle(clus);
double truthEnergy = dominantParticle.getEnergy();
double truthMomentum = dominantParticle.getMomentum().magnitude();
@@ -966,7 +984,7 @@
makePlotsOfEoverP(outputChargedParticleList, newMapTrackToShowerComponents, allSharedClusters);
- System.out.println("RECLUSTER: ALL DONE");
+ if (m_debug) { System.out.println("RECLUSTER: ALL DONE"); }
m_event = null;
}
@@ -1170,6 +1188,76 @@
totalChiSquared += trackChiSquared;
}
System.out.println("Total CHI^2 = "+totalChiSquared+" with NDF ~ "+tracksSortedByMomentum.size());
+ System.out.println("Here is a summary of all the energy deposits for each track:");
+ for (Track tr : tracksSortedByMomentum) {
+ double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
+ HitMap hitsEcal = ((HitMap)(m_event.get("inputHitMapEcal")));
+ HitMap hitsHcal = ((HitMap)(m_event.get("inputHitMapHcal")));
+ MCParticle trackTruth = ((BaseTrackMC)(tr)).getMCParticle();
+ List<CalorimeterHit> allTruthHitsInEcal = findHitsFromTruth_P(trackTruth, hitsEcal.values());
+ List<CalorimeterHit> allTruthHitsInHcal = findHitsFromTruth_P(trackTruth, hitsHcal.values());
+ int countCoreHitsEcal = 0;
+ int countCoreHitsHcal = 0;
+ int countSharedHitsEcal = 0;
+ int countSharedHitsHcal = 0;
+ double countWeightedSharedHitsEcal = 0.0;
+ double countWeightedSharedHitsHcal = 0.0;
+ int countUnmatchedSharedHitsEcal = 0;
+ int countUnmatchedSharedHitsHcal = 0;
+ Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
+ Set<Cluster> subClusters = new HashSet<Cluster>();
+ for (Cluster component : showerComponents) {
+ Set<Cluster> subClustersOfThisComponent = recursivelyFindSubClusters(component);
+ subClusters.addAll(subClustersOfThisComponent);
+ }
+ for (Cluster clus : showerComponents) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ if (allTruthHitsInEcal.contains(hit)) {
+ countCoreHitsEcal++;
+ } else if (allTruthHitsInHcal.contains(hit)) {
+ countCoreHitsHcal++;
+ }
+ }
+ }
+ // Investigate shared clusters
+ for (SharedClusterGroup shares : allSharedClusters) {
+ List<SharedCluster> sharedClusters = shares.listAllSharedClusters();
+ for (SharedCluster sharedCluster : sharedClusters) {
+ // First, check whether this shared cluster has any true hits...
+ int countTruthHitsEcalInThisSharedCluster = 0;
+ int countTruthHitsHcalInThisSharedCluster = 0;
+ for (CalorimeterHit hit : sharedCluster.getCluster().getCalorimeterHits()) {
+ if (allTruthHitsInEcal.contains(hit)) {
+ countTruthHitsEcalInThisSharedCluster++;
+ } else if (allTruthHitsInHcal.contains(hit)) {
+ countTruthHitsHcalInThisSharedCluster++;
+ }
+ }
+ if (countTruthHitsEcalInThisSharedCluster + countTruthHitsHcalInThisSharedCluster > 0) {
+ // There were some truth hits. Check the weight.
+ double sumWeightForThisSharedCluster = 0.0;
+ Set<Cluster> targetsInSharedCluster = sharedCluster.getTargetClusters();
+ for (Cluster targetInSharedCluster : targetsInSharedCluster) {
+ if (subClusters.contains(targetInSharedCluster)) {
+ sumWeightForThisSharedCluster += sharedCluster.getNormalizedWeight(targetInSharedCluster);
+ }
+ }
+ if (sumWeightForThisSharedCluster > 0.0) {
+ countSharedHitsEcal += countTruthHitsEcalInThisSharedCluster;
+ countSharedHitsHcal += countTruthHitsHcalInThisSharedCluster;
+ countWeightedSharedHitsEcal += ((double)(countTruthHitsEcalInThisSharedCluster)) * sumWeightForThisSharedCluster;
+ countWeightedSharedHitsHcal += ((double)(countTruthHitsHcalInThisSharedCluster)) * sumWeightForThisSharedCluster;
+ } else {
+ if (targetsInSharedCluster == null || targetsInSharedCluster.size()==0) {
+ countUnmatchedSharedHitsEcal += countTruthHitsEcalInThisSharedCluster;
+ countUnmatchedSharedHitsHcal += countTruthHitsHcalInThisSharedCluster;
+ }
+ }
+ }
+ }
+ }
+ System.out.println(" Track with p="+trackMomentum+" has "+allTruthHitsInEcal.size()+" ECAL + "+allTruthHitsInHcal.size()+" HCAL hits total. Of these, core has "+countCoreHitsEcal+" ECAL + "+countCoreHitsHcal+" HCAL hits and halo has "+countSharedHitsEcal+" ECAL + "+countSharedHitsHcal+" HCAL (weighted to "+countWeightedSharedHitsEcal+" + "+countWeightedSharedHitsHcal+"). Unmatched shared hits: "+countUnmatchedSharedHitsEcal+" ECAL + "+countUnmatchedSharedHitsHcal+" HCAL.");
+ }
for (Cluster clus : photonClusters) {
debugPrint("Photon", clus,tracksSortedByMomentum,newMapTrackToShowerComponents,newMapShowerComponentToTrack,allSharedClusters,newMapTrackToVetoedAdditions);
debugPrint(clus, newMapShowerComponentToTrack, photonClusters, mips, clumps,largeClustersWithoutSkeletons);
@@ -1636,7 +1724,8 @@
for (Cluster clus : output) {
if (clus.getCalorimeterHits().size() < 6 && m_mips.contains(clus)) {
// VETO
- System.out.println("DEBUG: Skipping daughters of cluster for implied-cluster-addition because it is a mip with only "+clus.getCalorimeterHits().size()+" hits.");
+ if (m_debug) { System.out.println("DEBUG: Skipping daughters of cluster for implied-cluster-addition because it is a mip with only "+clus.getCalorimeterHits().size()+" hits."); }
+ continue;
}
List<ScoredLink> potentialLinks = m_potentialLinks.get(clus);
@@ -1912,8 +2001,7 @@
// Start building cluster
double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
Cluster seed = tracksMatchedToClusters.get(tr);
- {
- // DEBUG stuff
+ if (m_debug) {
System.out.println("Rerecluster: starting with a track of p="+trackMomentum+", connected to a seed with "+seed.getCalorimeterHits().size()+" hits...");
}
Set<Cluster> allShowerComponents = new HashSet<Cluster>();
@@ -1925,12 +2013,13 @@
for (int iTier=0; iTier<25; iTier++) { // 25 is arbitrary cap in case we get stuck.
long timeAtStartOfThisTierLoop = Calendar.getInstance().getTimeInMillis(); // DEBUG
List<Cluster> componentsInPreviousTier = tieredShowerComponents.get(iTier);
- System.out.println("DEBUG: Shower of track with p="+trackMomentum+" has "+componentsInPreviousTier.size()+" clusters in tier "+iTier+". Expanding...");
+ if (m_debug) {System.out.println("DEBUG: Shower of track with p="+trackMomentum+" has "+componentsInPreviousTier.size()+" clusters in tier "+iTier+". Expanding..."); }
List<Cluster> componentsInNextTier = new Vector<Cluster>();
for (Cluster miniSeed : componentsInPreviousTier) {
if (miniSeed != seed && miniSeed.getCalorimeterHits().size() < 6 && m_mips.contains(miniSeed)) {
// VETO
- System.out.println("DEBUG: Skipping daughters of cluster for cluster-addition because it is a mip with only "+miniSeed.getCalorimeterHits().size()+" hits.");
+ if (m_debug) { System.out.println("DEBUG: Skipping daughters of cluster for cluster-addition because it is a mip with only "+miniSeed.getCalorimeterHits().size()+" hits."); }
+ continue;
}
long timeAtStartOfThisMiniSeedLoop = Calendar.getInstance().getTimeInMillis(); // DEBUG
// Search for links above threshold