lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.5 -r1.6
--- ReclusterDriver.java 19 Nov 2007 22:19:02 -0000 1.5
+++ ReclusterDriver.java 21 Nov 2007 19:48:37 -0000 1.6
@@ -30,7 +30,7 @@
*
* This version is PRELIMINARY.
*
- * @version $Id: ReclusterDriver.java,v 1.5 2007/11/19 22:19:02 mcharles Exp $
+ * @version $Id: ReclusterDriver.java,v 1.6 2007/11/21 19:48:37 mcharles Exp $
* @author Mat Charles
*/
@@ -504,238 +504,68 @@
//debugPrint(largeClustersWithoutSkeletons, sharedSmallClusters, "Large", "halo");
// Cache potential links (warning, there may be a lot!)
resetPotentialLinks();
- // First look for mip-clump links, only working within a large, contiguous cluster:
+
+ // Scoring constants
double scaleFactorTrackToTrack = 1.0;
double scaleFactorTrackToClump = 1.0;
double scaleFactorTrackToSmallSeed = 1.0;
+ double thresholdForProximity = 50.0; // 5cm. FIXME: Don't hard-code.
+ double thresholdForProximityClump = 75.0; // FIXME: Don't hard-code.
- // Links between main components (clumps, mips) within a large cluster
- // -------------------------------------------------------------------
- System.out.println("DEBUG: There are "+largeClusters.size()+" large clusters in the event which can contain skeletons.");
- for (Cluster largeClus : largeClusters) {
- List<Cluster> subClustersMips = new Vector<Cluster>();
- List<Cluster> subClustersClumps = new Vector<Cluster>();
- for (Cluster subClus : largeClus.getClusters()) {
- for (Cluster subsubClus : subClus.getClusters()) {
- if (mips.contains(subsubClus)) { subClustersMips.add(subsubClus); }
- if (clumps.contains(subsubClus)) { subClustersClumps.add(subsubClus); }
- }
- }
- if (m_debugLinkScores) { System.out.println("DEBUG: Large cluster with "+largeClus.getCalorimeterHits().size()+" hits contains "+subClustersMips.size()+" MIPs and "+subClustersClumps.size()+" clumps."); }
- // Potential mip-mip links [careful -- don't double-count]
- for (int i=0; i<subClustersMips.size(); i++) {
- Cluster mip1 = subClustersMips.get(i);
- for (int j=i+1; j<subClustersMips.size(); j++) {
- Cluster mip2 = subClustersMips.get(j);
- double likelihood = m_eval.getLinkLikelihoodTrackToTrack(mip1,mip2);
- double score = likelihood / scaleFactorTrackToTrack;
- if (score > 1.0) { score = 1.0; } else if (score < 0.0) { score = 0.0; }
- addPotentialLink(mip1, mip2, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip1.getCalorimeterHits().size()+"] -- MIP["+mip2.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(mip1, mip2, "MIP", "MIP", score); }
- }
- }
- // Potential mip-clump links
- for (Cluster mip : subClustersMips) {
- for (Cluster clump : subClustersClumps) {
- double likelihood = m_eval.getLinkLikelihoodTrackToClump(mip,clump);
- double score = likelihood / scaleFactorTrackToClump;
- if (score > 1.0) { score = 1.0; } else if (score < 0.0) { score = 0.0; }
- addPotentialLink(mip, clump, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip.getCalorimeterHits().size()+"] -- Clump["+clump.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(mip, clump, "MIP", "Clump", score); }
- }
- }
- // Some debug info
- if (m_debugLinkScores) {
- if (subClustersMips.size()+subClustersClumps.size()==1 || subClustersMips.size()==0) {
- for (Cluster clus : subClustersMips) { System.out.println("DEBUG: Unlinkable MIP["+clus.getCalorimeterHits().size()+"]"); }
- for (Cluster clus : subClustersClumps) { System.out.println("DEBUG: Unlinkable Clump["+clus.getCalorimeterHits().size()+"]"); }
- }
- }
- }
-
- // Potential links to other things (harder to pin down...)
- // -------------------------------------------------------
+ boolean cheatOnScoring = false;
- // Link from a MIP to another MIP/clump which is *not* in the same large cluster.
- // This can happen if:
- // 1) One is part of a secondary shower
- // 2) Linkage across the ECAL/HCAL boundary failed
- // MIP <--> MIP
- for (int i=0; i<mips.size(); i++) {
- Cluster mip1 = mips.get(i);
- for (int j=i+1; j<mips.size(); j++) {
- Cluster mip2 = mips.get(j);
- boolean linkAlreadyExists = checkForLink(mip1, mip2);
- if (!linkAlreadyExists) {
- double likelihood = m_eval.getLinkLikelihoodTrackToTrack(mip1,mip2);
- double score = likelihood / scaleFactorTrackToTrack;
- // Penalty factor for going outside main cluster, based on angular separation
- Hep3Vector clus1Position = new BasicHep3Vector(mip1.getPosition());
- Hep3Vector clus2Position = new BasicHep3Vector(mip2.getPosition());
- double penaltyFactor = VecOp.dot(VecOp.unit(clus1Position), VecOp.unit(clus2Position));
- if (penaltyFactor > 0.8) { penaltyFactor = 0.8; }
- score *= penaltyFactor;
- if (score > 1.0) { score = 1.0; } else if (score < 0.0) { score = 0.0; }
- addPotentialLink(mip1, mip2, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip1.getCalorimeterHits().size()+"] -- MIP["+mip2.getCalorimeterHits().size()+"] = "+score+" (not in same large cluster)"); }
- if (m_debugLinkScores) { debugPrintLink(mip1, mip2, "MIP", "MIP", score); }
- }
- }
- }
- // MIP <--> Clump
- for (Cluster mip : mips) {
- for (Cluster clump : clumps) {
- boolean linkAlreadyExists = checkForLink(mip, clump);
- if (!linkAlreadyExists) {
- double likelihood = m_eval.getLinkLikelihoodTrackToClump(mip,clump);
- double score = likelihood / scaleFactorTrackToClump;
- // Penalty factor for going outside main cluster, based on angular separation
- Hep3Vector clus1Position = new BasicHep3Vector(mip.getPosition());
- Hep3Vector clus2Position = new BasicHep3Vector(clump.getPosition());
- double penaltyFactor = VecOp.dot(VecOp.unit(clus1Position), VecOp.unit(clus2Position));
- if (penaltyFactor > 0.8) { penaltyFactor = 0.8; }
- score *= penaltyFactor;
- if (score > 1.0) { score = 1.0; } else if (score < 0.0) { score = 0.0; }
- addPotentialLink(mip, clump, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip.getCalorimeterHits().size()+"] -- Clump["+clump.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(mip, clump, "MIP", "Clump", score); }
- }
- }
- }
+ if (cheatOnScoring) {
+ initPotentialLinks_cheating(linkableClusters);
+ } else {
+ // Links between main components (clumps, mips) within a large cluster
+ // -------------------------------------------------------------------
+ System.out.println("DEBUG: There are "+largeClusters.size()+" large clusters in the event which can contain skeletons.");
+ initPotentialLinks_SkeletonsWithinLargeClus(largeClusters, mips, clumps, scaleFactorTrackToTrack, scaleFactorTrackToClump);
+
+ // Potential links to other things (harder to pin down...)
+ // -------------------------------------------------------
+
+ // Link from a MIP to another MIP/clump which is *not* in the same large cluster.
+ // This can happen if:
+ // 1) One is part of a secondary shower
+ // 2) Linkage across the ECAL/HCAL boundary failed
+ // MIP <--> MIP
+ initPotentialLinks_MipMip(mips, scaleFactorTrackToTrack, true);
+ // MIP <--> Clump
+ initPotentialLinks_MipClump(mips, clumps, scaleFactorTrackToClump, true);
- // Link from MIP to a nearby cluster (use pointing)
- for (Cluster mip : mips) {
- double thresholdForProximity = 50.0; // 5cm. FIXME: Don't hard-code.
+ // Link from MIP to a nearby cluster (use pointing)
// MIP <--> small cluster seed
- for (Cluster clus : smallClusterSeeds) {
- double score = scoreOnProximityAndPointing(mip, clus, thresholdForProximity);
- addPotentialLink(mip, clus, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip.getCalorimeterHits().size()+"] -- SmallSeed["+clus.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(mip, clus, "MIP", "SmallSeed", score); }
- if (m_debugLinkScores) {
- double likelihood = m_eval.getLinkLikelihoodTrackToClump(mip,clus);
- double dist = proximity(mip, clus);
- double scaledDistance = dist / thresholdForProximity;
- double penaltyFactor = 1.0 / (scaledDistance * scaledDistance);
- double recalculatedScore = likelihood;
- if (penaltyFactor < 1.0) {
- recalculatedScore *= penaltyFactor;
- }
- System.out.println("DEBUG: proximity = "+dist+" and threshold = "+thresholdForProximity+" => penaltyFactor = "+penaltyFactor);
- System.out.println("DEBUG: likelihood = "+likelihood+" => combined score = "+recalculatedScore);
- }
- }
+ initPotentialLinks_MipMisc(mips, smallClusterSeeds, thresholdForProximity, "SmallSeed");
// MIP <--> photon
- for (Cluster clus : photonClusters) {
- boolean photonIsSeed = seeds.contains(clus);
- double score = scoreOnProximityAndPointing(mip, clus, thresholdForProximity);
- if (photonIsSeed) {
- addPotentialLink(mip, clus, score);
- if (m_debugLinkScores) { debugPrintLink(mip, clus, "MIP", "Photon", score); }
- } else {
- String printme = new String("DEBUG: Skipped making link from MIP["+mip.getCalorimeterHits().size()+"] -- Photon["+clus.getCalorimeterHits().size()+"] because photon is not a seed.");
- MCParticle domPart1 = quoteDominantParticle(mip);
- MCParticle domPart2 = quoteDominantParticle(clus);
- if (domPart1 == domPart2) {
- int domPDG = domPart1.getPDGID();
- double domMom = domPart1.getMomentum().magnitude();
- printme += " MISTAKE! They both come from "+domPDG+" with p="+domMom;
- double likelihood = m_eval.getLinkLikelihoodTrackToClump(mip,clus);
- double distance = proximity(mip, clus);
- printme += ". Score would have been "+score+" since likelihood="+likelihood+" and proximity="+distance;
- }
- System.out.println(printme);
- }
+ List<Cluster> seedPhotonClusters = new Vector<Cluster>();
+ for (Cluster photon : photonClusters) {
+ if (seeds.contains(photon)) {
+ seedPhotonClusters.add(photon);
+ }
}
+ initPotentialLinks_MipMisc(mips, seedPhotonClusters, thresholdForProximity, "Photon");
// MIP <--> large cluster
- for (Cluster clus : largeClustersWithoutSkeletons) {
- double score = scoreOnProximityAndPointing(mip, clus, thresholdForProximity);
- addPotentialLink(mip, clus, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip.getCalorimeterHits().size()+"] -- LargeClus["+clus.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(mip, clus, "MIP", "LargeClus", score); }
- }
- }
- // Link certain things based on proximity and angle when there is
- // little other information. Use these quantities:
- // 1) Proximity
- // 2) Dot product of these vectors: (IP to inner clus) and (inner clus to outer clus)
- // For genuinely large clusters (e.g. clumps with 50 hits) we also have pointing information
- for (Cluster smallSeed : smallClusterSeeds) {
- Hep3Vector smallSeedPosition = new BasicHep3Vector(smallSeed.getPosition());
- double thresholdForProximity = 50.0; // 5cm. FIXME: Don't hard-code.
+ initPotentialLinks_MipMisc(mips, largeClustersWithoutSkeletons, thresholdForProximity, "LargeClus");
+
+ // Link certain things based on proximity and angle when there is
+ // little other information. Use these quantities:
+ // 1) Proximity
+ // 2) Dot product of these vectors: (IP to inner clus) and (inner clus to outer clus)
+ // For genuinely large clusters (e.g. clumps with 50 hits) we also have pointing information
+
// Small seeds <--> clumps
- for (Cluster clus : clumps) {
- double score = scoreOnProximityAndAngle(smallSeed, clus, thresholdForProximity);
- addPotentialLink(smallSeed, clus, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: SmallSeed["+smallSeed.getCalorimeterHits().size()+"] -- Clump["+clus.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(smallSeed, clus, "SmallSeed", "Clump", score); }
- }
+ initPotentialLinks_MiscMisc(smallClusterSeeds, clumps, thresholdForProximity, "SmallSeed", "Clump");
// Small seeds <--> photons
- for (Cluster clus : photonClusters) {
- double score = scoreOnProximityAndAngle(smallSeed, clus, thresholdForProximity);
- addPotentialLink(smallSeed, clus, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: SmallSeed["+smallSeed.getCalorimeterHits().size()+"] -- Photon["+clus.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(smallSeed, clus, "SmallSeed", "Photon", score); }
- if (m_debugLinkScores) {
- double dist = proximity(smallSeed, clus);
- double scaledDist = dist / thresholdForProximity;
- double proximityScore = 1.0 / (scaledDist * scaledDist);
- Hep3Vector clus1Position = new BasicHep3Vector(smallSeed.getPosition());
- Hep3Vector clus2Position = new BasicHep3Vector(clus.getPosition());
- Hep3Vector positionOfInnerCluster = null;
- Hep3Vector positionOfOuterCluster = null;
- if (clus2Position.magnitude() > clus1Position.magnitude()) {
- positionOfInnerCluster = clus1Position;
- positionOfOuterCluster = clus2Position;
- } else {
- positionOfInnerCluster = clus2Position;
- positionOfOuterCluster = clus1Position;
- }
- Hep3Vector displacementVector = VecOp.sub(positionOfOuterCluster, positionOfInnerCluster);
- double dotProduct = VecOp.dot(VecOp.unit(displacementVector), VecOp.unit(positionOfInnerCluster));
- double dotProductScore = dotProduct;
- if (dotProduct < 0.0) { dotProductScore = 0.0; }
- double recalculatedScore = proximityScore * dotProductScore;
- System.out.println("DEBUG: proximity="+dist+" and threshold="+thresholdForProximity+" => proximityScore="+proximityScore);
- System.out.println("DEBUG: dotProduct="+dotProduct+" => score="+recalculatedScore);
- }
- }
+ initPotentialLinks_MiscMisc(smallClusterSeeds, photonClusters, thresholdForProximity, "SmallSeed", "Photon");
// Small seeds <--> large clusters
- for (Cluster clus : largeClustersWithoutSkeletons) {
- double score = scoreOnProximityAndAngle(smallSeed, clus, thresholdForProximity);
- addPotentialLink(smallSeed, clus, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: SmallSeed["+smallSeed.getCalorimeterHits().size()+"] -- LargeClus["+clus.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(smallSeed, clus, "SmallSeed", "LargeClus", score); }
- }
- }
- // Clump <--> clump
- for (int i=0; i<clumps.size(); i++) {
- double thresholdForProximity = 75.0; // FIXME: Don't hard-code.
- Cluster clump1 = clumps.get(i);
- for (int j=i+1; j<clumps.size(); j++) {
- Cluster clump2 = clumps.get(j);
- double score = scoreOnProximityAndAngle(clump1, clump2, thresholdForProximity);
- score *= 0.8; // Arbitrary penalty [FIXME]
- addPotentialLink(clump1, clump2, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: Clump["+clump1.getCalorimeterHits().size()+"] -- Clump["+clump2.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(clump1, clump2, "Clump", "Clump", score); }
- }
- }
- // Clump <--> large cluster
- for (Cluster clump : clumps) {
- double thresholdForProximity = 75.0; // FIXME: Don't hard-code.
- for (Cluster clus : largeClustersWithoutSkeletons) {
- double score = scoreOnProximityAndAngle(clump, clus, thresholdForProximity);
- score *= 0.8; // Arbitrary penalty [FIXME]
- addPotentialLink(clump, clus, score);
- //if (m_debugLinkScores) { System.out.println("DEBUG: Clump["+clump.getCalorimeterHits().size()+"] -- LargeClus["+clus.getCalorimeterHits().size()+"] = "+score); }
- if (m_debugLinkScores) { debugPrintLink(clump, clus, "Clump", "LargeClus", score); }
- }
+ initPotentialLinks_MiscMisc(smallClusterSeeds, largeClustersWithoutSkeletons, thresholdForProximity, "SmallSeed", "LargeClus");
+ // Clump <--> clump
+ initPotentialLinks_MiscSelf(clumps, thresholdForProximityClump, "Clump", true);
+ // Clump <--> large cluster
+ initPotentialLinks_MiscMisc(clumps, largeClustersWithoutSkeletons, thresholdForProximityClump, "Clump", "LargeClus");
}
- // [other legitimate link types]
- // [...]
+
// OK. Now sort the potential links by score:
sortLinks();
// Try to build skeletons, working with tracks in momentum-sorted order...
@@ -955,6 +785,7 @@
// write out charged shower
BaseReconstructedParticle part = new BaseReconstructedParticle();
for (Cluster clus : showerComponents) { part.addCluster(clus); }
+ // [FIXME: ADD FUZZY SIMCALORIMETERHITS]
part.addTrack(tr);
part.setCharge(tr.getCharge());
// Get the particle's three-momentum.
@@ -981,6 +812,7 @@
// write out photon
BaseReconstructedParticle part = new BaseReconstructedParticle();
part.addCluster(clus);
+ // [FIXME: ADD FUZZY SIMCALORIMETERHITS]
double clusterEnergy = energy(clus, allSharedClusters, m_photonCalib);
Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
Hep3Vector threeMomentum = VecOp.mult(clusterEnergy, VecOp.unit(pos)); // assume it comes from the IP
@@ -1046,15 +878,17 @@
for (Cluster miniSeed : componentsInPreviousTier) {
// Search for links above threshold
List<ScoredLink> potentialLinks = m_potentialLinks.get(miniSeed);
- for (ScoredLink link : potentialLinks) {
- if (link.score() < threshold) { break ; }
- Cluster newClus = link.counterpart(miniSeed);
- if (newClus == null) { throw new AssertionError("Null link!"); }
- // Ensure that we haven't already added component
- if (unusedUnmatchedClusterPieces.contains(newClus)) {
- componentsInNextTier.add(newClus);
- unusedUnmatchedClusterPieces.remove(newClus);
- piecesForThisCluster.add(newClus);
+ if (potentialLinks != null) {
+ for (ScoredLink link : potentialLinks) {
+ if (link.score() < threshold) { break ; }
+ Cluster newClus = link.counterpart(miniSeed);
+ if (newClus == null) { throw new AssertionError("Null link!"); }
+ // Ensure that we haven't already added component
+ if (unusedUnmatchedClusterPieces.contains(newClus)) {
+ componentsInNextTier.add(newClus);
+ unusedUnmatchedClusterPieces.remove(newClus);
+ piecesForThisCluster.add(newClus);
+ }
}
}
}
@@ -1072,6 +906,7 @@
// write out neutral hadron
BaseReconstructedParticle part = new BaseReconstructedParticle();
part.addCluster(clus);
+ // [FIXME: ADD FUZZY SIMCALORIMETERHITS]
double clusterEnergy = energy(clus, allSharedClusters, m_neutralCalib);
double clusterMomentumMagSq = clusterEnergy*clusterEnergy - mass_K0*mass_K0;
if (clusterMomentumMagSq<0.0) { clusterMomentumMagSq = 0.0; }
@@ -1107,6 +942,168 @@
m_event = null;
}
+ protected void initPotentialLinks_cheating(Collection<Cluster> clusters) {
+ Cluster[] tmpArray = new Cluster[1];
+ Cluster[] clusterArray = clusters.toArray(tmpArray);
+ if (clusterArray.length != clusters.size()) { throw new AssertionError("Book-keeping error"); }
+ for (int i=0; i<clusterArray.length; i++) {
+ Cluster clus1 = clusterArray[i];
+ MCParticle domPart1 = quoteDominantParticle(clus1);
+ for (int j=i+1; j<clusterArray.length; j++) {
+ Cluster clus2 = clusterArray[j];
+ if (clus1 == clus2) { throw new AssertionError("Book-keeping error"); }
+ MCParticle domPart2 = quoteDominantParticle(clus2);
+ if (domPart1 == domPart2) {
+ addPotentialLink(clus1, clus2, 1.0);
+ }
+ }
+ }
+ }
+
+ protected void initPotentialLinks_MipMip(List<Cluster> mips, double scaleFactorTrackToTrack, boolean applyPenaltyFactor)
+ {
+ for (int i=0; i<mips.size(); i++) {
+ Cluster mip1 = mips.get(i);
+ for (int j=i+1; j<mips.size(); j++) {
+ Cluster mip2 = mips.get(j);
+ boolean linkAlreadyExists = checkForLink(mip1, mip2);
+ if (!linkAlreadyExists) {
+ double likelihood = m_eval.getLinkLikelihoodTrackToTrack(mip1,mip2);
+ double score = likelihood / scaleFactorTrackToTrack;
+ // Penalty factor for going outside main cluster, based on angular separation
+ Hep3Vector clus1Position = new BasicHep3Vector(mip1.getPosition());
+ Hep3Vector clus2Position = new BasicHep3Vector(mip2.getPosition());
+ if (applyPenaltyFactor) {
+ double penaltyFactor = VecOp.dot(VecOp.unit(clus1Position), VecOp.unit(clus2Position));
+ if (penaltyFactor > 0.8) { penaltyFactor = 0.8; }
+ score *= penaltyFactor;
+ }
+ if (score > 1.0) { score = 1.0; } else if (score < 0.0) { score = 0.0; }
+ addPotentialLink(mip1, mip2, score);
+ //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip1.getCalorimeterHits().size()+"] -- MIP["+mip2.getCalorimeterHits().size()+"] = "+score+" (not in same large cluster)"); }
+ if (m_debugLinkScores) { debugPrintLink(mip1, mip2, "MIP", "MIP", score); }
+ }
+ }
+ }
+ }
+
+ protected void initPotentialLinks_MipClump(List<Cluster> mips, List<Cluster> clumps, double scaleFactorTrackToClump, boolean applyPenaltyFactor)
+ {
+ for (Cluster mip : mips) {
+ for (Cluster clump : clumps) {
+ boolean linkAlreadyExists = checkForLink(mip, clump);
+ if (!linkAlreadyExists) {
+ double likelihood = m_eval.getLinkLikelihoodTrackToClump(mip,clump);
+ double score = likelihood / scaleFactorTrackToClump;
+ if (applyPenaltyFactor) {
+ // Penalty factor for going outside main cluster, based on angular separation
+ Hep3Vector clus1Position = new BasicHep3Vector(mip.getPosition());
+ Hep3Vector clus2Position = new BasicHep3Vector(clump.getPosition());
+ double penaltyFactor = VecOp.dot(VecOp.unit(clus1Position), VecOp.unit(clus2Position));
+ if (penaltyFactor > 0.8) { penaltyFactor = 0.8; }
+ score *= penaltyFactor;
+ }
+ if (score > 1.0) { score = 1.0; } else if (score < 0.0) { score = 0.0; }
+ addPotentialLink(mip, clump, score);
+ //if (m_debugLinkScores) { System.out.println("DEBUG: MIP["+mip.getCalorimeterHits().size()+"] -- Clump["+clump.getCalorimeterHits().size()+"] = "+score); }
+ if (m_debugLinkScores) { debugPrintLink(mip, clump, "MIP", "Clump", score); }
+ }
+ }
+ }
+ }
+
+ protected void initPotentialLinks_MipMisc(List<Cluster> mips, List<Cluster> miscClusters, double thresholdForProximity, String typeName)
+ {
+ for (Cluster mip : mips) {
+ for (Cluster clus : miscClusters) {
+ boolean linkAlreadyExists = checkForLink(mip, clus);
+ if (!linkAlreadyExists) {
+ double score = scoreOnProximityAndPointing(mip, clus, thresholdForProximity);
+ addPotentialLink(mip, clus, score);
+ if (m_debugLinkScores) { debugPrintLink(mip, clus, "MIP", typeName, score); }
+ if (m_debugLinkScores) {
+ double likelihood = m_eval.getLinkLikelihoodTrackToClump(mip,clus);
+ double dist = proximity(mip, clus);
+ double scaledDistance = dist / thresholdForProximity;
+ double penaltyFactor = 1.0 / (scaledDistance * scaledDistance);
+ double recalculatedScore = likelihood;
+ if (penaltyFactor < 1.0) {
+ recalculatedScore *= penaltyFactor;
+ }
+ System.out.println("DEBUG: proximity = "+dist+" and threshold = "+thresholdForProximity+" => penaltyFactor = "+penaltyFactor);
+ System.out.println("DEBUG: likelihood = "+likelihood+" => combined score = "+recalculatedScore);
+ }
+ }
+ }
+ }
+ }
+
+ protected void initPotentialLinks_MiscMisc(List<Cluster> list1, List<Cluster> list2, double thresholdForProximity, String type1, String type2)
+ {
+ for (Cluster clus1 : list1) {
+ if (list2.contains(clus1)) { throw new AssertionError("Book-keeping error"); }
+ for (Cluster clus2 : list2) {
+ boolean linkAlreadyExists = checkForLink(clus1, clus2);
+ if (!linkAlreadyExists) {
+ double score = scoreOnProximityAndAngle(clus1, clus2, thresholdForProximity);
+ addPotentialLink(clus1, clus2, score);
+ if (m_debugLinkScores) { debugPrintLink(clus1, clus2, type1, type2, score); }
+ }
+ }
+ }
+ }
+
+ protected void initPotentialLinks_MiscSelf(List<Cluster> clusterList, double thresholdForProximity, String typeName, boolean applyPenalty)
+ {
+ for (int i=0; i<clusterList.size(); i++) {
+ Cluster clus1 = clusterList.get(i);
+ for (int j=i+1; j<clusterList.size(); j++) {
+ Cluster clus2 = clusterList.get(j);
+ boolean linkAlreadyExists = checkForLink(clus1, clus2);
+ if (!linkAlreadyExists) {
+ double score = scoreOnProximityAndAngle(clus1, clus2, thresholdForProximity);
+ if (applyPenalty) {
+ // Penalty factor for going outside main cluster, based on angular separation
+ Hep3Vector clus1Position = new BasicHep3Vector(clus1.getPosition());
+ Hep3Vector clus2Position = new BasicHep3Vector(clus2.getPosition());
+ double penaltyFactor = VecOp.dot(VecOp.unit(clus1Position), VecOp.unit(clus2Position));
+ if (penaltyFactor > 0.8) { penaltyFactor = 0.8; }
+ score *= penaltyFactor;
+ }
+ addPotentialLink(clus1, clus2, score);
+ if (m_debugLinkScores) { debugPrintLink(clus1, clus2, typeName, typeName, score); }
+ }
+ }
+ }
+ }
+
+ protected void initPotentialLinks_SkeletonsWithinLargeClus(List<Cluster> largeClusters, List<Cluster> mips, List<Cluster> clumps, double scaleFactorTrackToTrack, double scaleFactorTrackToClump)
+ {
+ for (Cluster largeClus : largeClusters) {
+ List<Cluster> subClustersMips = new Vector<Cluster>();
+ List<Cluster> subClustersClumps = new Vector<Cluster>();
+ for (Cluster subClus : largeClus.getClusters()) {
+ for (Cluster subsubClus : subClus.getClusters()) {
+ if (mips.contains(subsubClus)) { subClustersMips.add(subsubClus); }
+ if (clumps.contains(subsubClus)) { subClustersClumps.add(subsubClus); }
+ }
+ }
+ if (m_debugLinkScores) { System.out.println("DEBUG: Large cluster with "+largeClus.getCalorimeterHits().size()+" hits contains "+subClustersMips.size()+" MIPs and "+subClustersClumps.size()+" clumps."); }
+ // Potential mip-mip links
+ initPotentialLinks_MipMip(subClustersMips, scaleFactorTrackToTrack, false);
+ // Potential mip-clump links
+ initPotentialLinks_MipClump(subClustersMips, subClustersClumps, scaleFactorTrackToClump, false);
+
+ // Some debug info
+ if (m_debugLinkScores) {
+ if (subClustersMips.size()+subClustersClumps.size()==1 || subClustersMips.size()==0) {
+ for (Cluster clus : subClustersMips) { System.out.println("DEBUG: Unlinkable MIP["+clus.getCalorimeterHits().size()+"]"); }
+ for (Cluster clus : subClustersClumps) { System.out.println("DEBUG: Unlinkable Clump["+clus.getCalorimeterHits().size()+"]"); }
+ }
+ }
+ }
+ }
+
protected void printStatus(String desc,
List<Track> tracksSortedByMomentum,
List<SharedClusterGroup> allSharedClusters,
@@ -1897,69 +1894,71 @@
List<ScoredLink> potentialLinks = m_potentialLinks.get(miniSeed);
long timeAfterLookingUpPotentialLinks = Calendar.getInstance().getTimeInMillis(); // DEBUG
if (m_debugTiming) { System.out.println("DEBUG: Getting potential links took "+(timeAfterLookingUpPotentialLinks-timeAtStartOfThisMiniSeedLoop)+" ms"); }
- for (ScoredLink link : potentialLinks) {
- if (link.score() < threshold) { break ; }
- long timeAtStartOfThisPotentialLink = Calendar.getInstance().getTimeInMillis(); // DEBUG
- Cluster newClus = link.counterpart(miniSeed);
- if (newClus == null) { throw new AssertionError("Null link!"); }
- // Ensure that we haven't already added component
- long timeBeforePreviousUseTests = Calendar.getInstance().getTimeInMillis(); // DEBUG
- boolean testNotAlreadyInBaseShower = !(allShowerComponents.contains(newClus));
- boolean testNotAlreadyInNextTier = !(componentsInNextTier.contains(newClus));
- boolean testNotAlreadyAssigned = (newMapShowerComponentToTrack.get(newClus)==null);
- boolean testNotSeed = !(allSeedsInEvent.contains(newClus));
- long timeAfterPreviousUseTests = Calendar.getInstance().getTimeInMillis(); // DEBUG
- if (m_debugTiming) { System.out.println("DEBUG: Testing if previously used took "+(timeAfterPreviousUseTests-timeBeforePreviousUseTests)+" ms (and prep took "+(timeBeforePreviousUseTests-timeAtStartOfThisPotentialLink)+" ms)"); }
- if (testNotAlreadyInBaseShower && testNotAlreadyInNextTier && testNotAlreadyAssigned && testNotSeed) {
- // E/p check with full propagation test
- long timeBeforeProbeFullPropagation = Calendar.getInstance().getTimeInMillis(); // DEBUG
- Set<Cluster> impliedAdditionalClusters = probeFullPropagation(allShowerComponents, link, miniSeed);
- // Check that full-propagation set doesn't contain any other seeds (that would be a sign that
- // something went badly wrong!)
- boolean impliedAdditionalClustersContainsIllegalSeed = false;
- for (Cluster impliedClus : impliedAdditionalClusters) {
- if (allSeedsInEvent.contains(impliedClus) && impliedClus != seed) {
- impliedAdditionalClustersContainsIllegalSeed = true;
- break;
+ if (potentialLinks != null) {
+ for (ScoredLink link : potentialLinks) {
+ if (link.score() < threshold) { break ; }
+ long timeAtStartOfThisPotentialLink = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ Cluster newClus = link.counterpart(miniSeed);
+ if (newClus == null) { throw new AssertionError("Null link!"); }
+ // Ensure that we haven't already added component
+ long timeBeforePreviousUseTests = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ boolean testNotAlreadyInBaseShower = !(allShowerComponents.contains(newClus));
+ boolean testNotAlreadyInNextTier = !(componentsInNextTier.contains(newClus));
+ boolean testNotAlreadyAssigned = (newMapShowerComponentToTrack.get(newClus)==null);
+ boolean testNotSeed = !(allSeedsInEvent.contains(newClus));
+ long timeAfterPreviousUseTests = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ if (m_debugTiming) { System.out.println("DEBUG: Testing if previously used took "+(timeAfterPreviousUseTests-timeBeforePreviousUseTests)+" ms (and prep took "+(timeBeforePreviousUseTests-timeAtStartOfThisPotentialLink)+" ms)"); }
+ if (testNotAlreadyInBaseShower && testNotAlreadyInNextTier && testNotAlreadyAssigned && testNotSeed) {
+ // E/p check with full propagation test
+ long timeBeforeProbeFullPropagation = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ Set<Cluster> impliedAdditionalClusters = probeFullPropagation(allShowerComponents, link, miniSeed);
+ // Check that full-propagation set doesn't contain any other seeds (that would be a sign that
+ // something went badly wrong!)
+ boolean impliedAdditionalClustersContainsIllegalSeed = false;
+ for (Cluster impliedClus : impliedAdditionalClusters) {
+ if (allSeedsInEvent.contains(impliedClus) && impliedClus != seed) {
+ impliedAdditionalClustersContainsIllegalSeed = true;
+ break;
+ }
}
- }
- long timeAfterProbeFullPropagation = Calendar.getInstance().getTimeInMillis(); // DEBUG
- if (m_debugTiming) { System.out.println("DEBUG: probeFullPropagation() took "+(timeAfterProbeFullPropagation-timeBeforeProbeFullPropagation)+" ms"); }
- if (impliedAdditionalClustersContainsIllegalSeed) {
- System.out.println("DEBUG: Link with score="+link.score()+" to a subcluster with "+newClus.getCalorimeterHits().size()+" hits implied adding a total of "+impliedAdditionalClusters.size()+" clusters -- rejected because this would have picked up another track's seed.");
- } else {
- // OK
- Set<Cluster> existingShowerPlusAdditionalClusters = new HashSet<Cluster>();
- existingShowerPlusAdditionalClusters.addAll(allShowerComponents);
- existingShowerPlusAdditionalClusters.addAll(impliedAdditionalClusters);
- existingShowerPlusAdditionalClusters.addAll(componentsInNextTier);
- long timeBeforeEnergyCalculation = Calendar.getInstance().getTimeInMillis(); // DEBUG
- double energyOfExistingShowerPlusAdditionalClusters = energy(existingShowerPlusAdditionalClusters, allSharedClusters);
- long timeAfterEnergyCalculation = Calendar.getInstance().getTimeInMillis(); // DEBUG
- if (m_debugTiming) { System.out.println("DEBUG: Energy calculation took "+(timeAfterEnergyCalculation-timeBeforeEnergyCalculation)+" ms (and prep took"+(timeBeforeEnergyCalculation-timeAfterProbeFullPropagation)+" ms)"); }
- boolean testValidEoverP = testEoverP(energyOfExistingShowerPlusAdditionalClusters, tr, tolerance);
- long timeAfterTestEoverP = Calendar.getInstance().getTimeInMillis(); // DEBUG
- if (m_debugTiming) { System.out.println("DEBUG: E/p test took "+(timeAfterTestEoverP-timeAfterEnergyCalculation)+" ms"); }
- if (testValidEoverP) {
- // OK -- add 'em
- componentsInNextTier.addAll(impliedAdditionalClusters);
- if (m_debugLinkScores) {
- String printme = "DEBUG: Link with score="+link.score()+" to a subcluster with "+newClus.getCalorimeterHits().size()+" hits implied adding a total of "+impliedAdditionalClusters.size()+" clusters for a new running total of "+energyOfExistingShowerPlusAdditionalClusters+". Implicitly added:";
- for (Cluster impliedClus : impliedAdditionalClusters) {
- printme += " ["+impliedClus.getCalorimeterHits().size()+"]";
- }
- printme += ".";
- System.out.println(printme);
- }
+ long timeAfterProbeFullPropagation = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ if (m_debugTiming) { System.out.println("DEBUG: probeFullPropagation() took "+(timeAfterProbeFullPropagation-timeBeforeProbeFullPropagation)+" ms"); }
+ if (impliedAdditionalClustersContainsIllegalSeed) {
+ System.out.println("DEBUG: Link with score="+link.score()+" to a subcluster with "+newClus.getCalorimeterHits().size()+" hits implied adding a total of "+impliedAdditionalClusters.size()+" clusters -- rejected because this would have picked up another track's seed.");
} else {
- // Nope -- failed
- if (m_debugLinkScores) { System.out.println("DEBUG: Link with score="+link.score()+" to a subcluster with "+newClus.getCalorimeterHits().size()+" hits implied adding a total of "+impliedAdditionalClusters.size()+" clusters -- rejected because new running total would be "+energyOfExistingShowerPlusAdditionalClusters); }
- vetoedAdditions.add(newClus);
+ // OK
+ Set<Cluster> existingShowerPlusAdditionalClusters = new HashSet<Cluster>();
+ existingShowerPlusAdditionalClusters.addAll(allShowerComponents);
+ existingShowerPlusAdditionalClusters.addAll(impliedAdditionalClusters);
+ existingShowerPlusAdditionalClusters.addAll(componentsInNextTier);
+ long timeBeforeEnergyCalculation = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ double energyOfExistingShowerPlusAdditionalClusters = energy(existingShowerPlusAdditionalClusters, allSharedClusters);
+ long timeAfterEnergyCalculation = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ if (m_debugTiming) { System.out.println("DEBUG: Energy calculation took "+(timeAfterEnergyCalculation-timeBeforeEnergyCalculation)+" ms (and prep took"+(timeBeforeEnergyCalculation-timeAfterProbeFullPropagation)+" ms)"); }
+ boolean testValidEoverP = testEoverP(energyOfExistingShowerPlusAdditionalClusters, tr, tolerance);
+ long timeAfterTestEoverP = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ if (m_debugTiming) { System.out.println("DEBUG: E/p test took "+(timeAfterTestEoverP-timeAfterEnergyCalculation)+" ms"); }
+ if (testValidEoverP) {
+ // OK -- add 'em
+ componentsInNextTier.addAll(impliedAdditionalClusters);
+ if (m_debugLinkScores) {
+ String printme = "DEBUG: Link with score="+link.score()+" to a subcluster with "+newClus.getCalorimeterHits().size()+" hits implied adding a total of "+impliedAdditionalClusters.size()+" clusters for a new running total of "+energyOfExistingShowerPlusAdditionalClusters+". Implicitly added:";
+ for (Cluster impliedClus : impliedAdditionalClusters) {
+ printme += " ["+impliedClus.getCalorimeterHits().size()+"]";
+ }
+ printme += ".";
+ System.out.println(printme);
+ }
+ } else {
+ // Nope -- failed
+ if (m_debugLinkScores) { System.out.println("DEBUG: Link with score="+link.score()+" to a subcluster with "+newClus.getCalorimeterHits().size()+" hits implied adding a total of "+impliedAdditionalClusters.size()+" clusters -- rejected because new running total would be "+energyOfExistingShowerPlusAdditionalClusters); }
+ vetoedAdditions.add(newClus);
+ }
}
}
+ long timeAtEndOfThisPotentialLink = Calendar.getInstance().getTimeInMillis(); // DEBUG
+ if (m_debugTiming) { System.out.println("DEBUG: Testing this potential link took "+(timeAtEndOfThisPotentialLink-timeAtStartOfThisPotentialLink)+" ms"); }
}
- long timeAtEndOfThisPotentialLink = Calendar.getInstance().getTimeInMillis(); // DEBUG
- if (m_debugTiming) { System.out.println("DEBUG: Testing this potential link took "+(timeAtEndOfThisPotentialLink-timeAtStartOfThisPotentialLink)+" ms"); }
}
long timeAtEndOfThisMiniSeedLoop = Calendar.getInstance().getTimeInMillis(); // DEBUG
if (m_debugTiming) { System.out.println("DEBUG: Testing this miniSeed took "+(timeAtEndOfThisMiniSeedLoop-timeAtStartOfThisMiniSeedLoop)+" ms"); }