lcsim/src/org/lcsim/contrib/uiowa
diff -N ReclusterDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReclusterDriver.java 14 Oct 2007 22:08:10 -0000 1.1
@@ -0,0 +1,830 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.util.*;
+import org.lcsim.util.hitmap.*;
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+import org.lcsim.recon.cluster.mipfinder.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.recon.cluster.clumpfinder.*;
+import org.lcsim.recon.pfa.identifier.*;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.base.*;
+import hep.physics.particle.Particle;
+
+public class ReclusterDriver extends Driver {
+
+ String m_mcList;
+ String m_inputSmallClusters;
+ String m_inputUnusedHits;
+ String m_inputSkeletonClusters;
+ String m_inputPhotonClusters;
+ String m_inputMuonParticles;
+ String m_inputTrackList;
+ String m_inputLargeClusters;
+ String m_inputMIPs;
+ String m_inputClumps;
+ String m_inputSplitSkeletonClusters;
+ EventHeader m_event;
+
+ public ReclusterDriver(String mcList, String trackList, String muonParticles, String photonClusters, String skeletonClusters, String smallClusters, String unusedHits, String largeClusters, String mips, String clumps, String splitSkeletonClusters) {
+ m_mcList = mcList;
+ m_inputTrackList = trackList;
+ m_inputMuonParticles = muonParticles;
+ m_inputPhotonClusters = photonClusters;
+ m_inputSkeletonClusters = skeletonClusters;
+ m_inputSmallClusters = smallClusters;
+ m_inputUnusedHits = unusedHits;
+ m_inputLargeClusters = largeClusters;
+ m_inputMIPs = mips;
+ m_inputClumps = clumps;
+ m_inputSplitSkeletonClusters = splitSkeletonClusters;
+
+ // For now, hard-coded
+ LocalHelixExtrapolationTrackClusterMatcher genMatch = new LocalHelixExtrapolationTrackClusterMatcher();
+ LocalHelixExtrapolationTrackMIPClusterMatcher mipMatch = new LocalHelixExtrapolationTrackMIPClusterMatcher();
+ add(mipMatch);
+ add(genMatch);
+ m_trackClusterMatcher = new DualActionTrackClusterMatcher(mipMatch, genMatch);
+ }
+
+ public void process(EventHeader event) {
+ super.process(event); // hmm...
+ m_event = event;
+
+ // Read in
+ List<Track> trackList = event.get(Track.class, m_inputTrackList);
+ List<ReconstructedParticle> muonParticles = event.get(ReconstructedParticle.class, m_inputMuonParticles);
+ List<Cluster> photonClusters = event.get(Cluster.class, m_inputPhotonClusters);
+ List<Cluster> skeletonClusters = event.get(Cluster.class, m_inputSkeletonClusters);
+ List<Cluster> smallClusters = event.get(Cluster.class, m_inputSmallClusters);
+ List<Cluster> largeClusters = event.get(Cluster.class, m_inputLargeClusters);
+ List<Cluster> mips = event.get(Cluster.class, m_inputMIPs);
+ List<Cluster> clumps = event.get(Cluster.class, m_inputClumps);
+ HitMap unusedHits = ((HitMap)(event.get(m_inputUnusedHits)));
+ List<Cluster> splitSkeletonClusters = event.get(Cluster.class, m_inputSplitSkeletonClusters);
+
+ 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>();
+ for (Cluster skel : skeletonClusters) {
+ for (Cluster clus : skel.getClusters()) {
+ if (mips.contains(clus)) {
+ mapMIPToSkeleton.put(clus, skel);
+ } else if (clumps.contains(clus)) {
+ mapClumpToSkeleton.put(clus, skel);
+ } else {
+ throw new AssertionError("Internal consistency failure");
+ }
+ }
+ }
+ Map<Cluster,Cluster> mapMIPToSplitSkeleton = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> mapClumpToSplitSkeleton = new HashMap<Cluster,Cluster>();
+ for (Cluster skel : splitSkeletonClusters) {
+ for (Cluster clus : skel.getClusters()) {
+ if (mips.contains(clus)) {
+ mapMIPToSplitSkeleton.put(clus, skel);
+ } else if (clumps.contains(clus)) {
+ mapClumpToSplitSkeleton.put(clus, skel);
+ } else {
+ throw new AssertionError("Internal consistency failure");
+ }
+ }
+ }
+
+
+ {
+ int countSmallClusterHits = 0;
+ for (Cluster clus : smallClusters) {
+ countSmallClusterHits += clus.getCalorimeterHits().size();
+ }
+ System.out.println("DEBUG: ReclusterDriver read in "+smallClusters.size()+" small clusters comprising "+countSmallClusterHits+" hits.");
+ }
+
+ // Many of those LargeClusters contain skeleton(s), but a few may not. Identify them.
+ List<Cluster> largeClustersWithoutSkeletons = new Vector<Cluster>();
+ for (Cluster largeCluster : largeClusters) {
+ List<CalorimeterHit> largeClusterHits = largeCluster.getCalorimeterHits();
+ boolean largeClusterContainsSkeleton = false;
+ for (Cluster skeletonCluster : skeletonClusters) {
+ for (CalorimeterHit hit : skeletonCluster.getCalorimeterHits()) {
+ if (largeClusterHits.contains(hit)) {
+ largeClusterContainsSkeleton = true;
+ break;
+ }
+ }
+ if (largeClusterContainsSkeleton) { break; }
+ }
+ if (!largeClusterContainsSkeleton) {
+ largeClustersWithoutSkeletons.add(largeCluster);
+ }
+ }
+
+ {
+ int countLargeHitsUnused = 0;
+ for (Cluster clus : largeClustersWithoutSkeletons) {
+ countLargeHitsUnused += clus.getCalorimeterHits().size();
+ }
+ System.out.println("DEBUG: ReclusterDriver found "+largeClustersWithoutSkeletons.size()+" large clusters without skeletons, comprising "+countLargeHitsUnused+" hits.");
+ }
+
+ // Pull those unused hits:
+ List<CalorimeterHit> unusedHitsIgnoringLargeClustersWithoutSkeletons = new Vector<CalorimeterHit>();
+ unusedHitsIgnoringLargeClustersWithoutSkeletons.addAll(unusedHits.values());
+ for (Cluster clus : largeClustersWithoutSkeletons) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ unusedHitsIgnoringLargeClustersWithoutSkeletons.remove(hit);
+ }
+ }
+ List<Cluster> wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons = new Vector<Cluster>();
+ for (CalorimeterHit hit : unusedHitsIgnoringLargeClustersWithoutSkeletons) {
+ BasicCluster clus = new BasicCluster();
+ clus.addHit(hit);
+ wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons.add(clus);
+ }
+
+ // Match tracks to clusters without regard for E/p.
+ // We allow them to match the following categories of cluster:
+ // mips & clumps (from skeletonClusters)
+ // largeClustersWithoutSkeletons
+ // photonClusters (could be electrons)
+ // smallClusters
+ List<Cluster> mixedClusters = new Vector<Cluster>();
+ mixedClusters.addAll(mips);
+ mixedClusters.addAll(clumps);
+ mixedClusters.addAll(largeClustersWithoutSkeletons);
+ mixedClusters.addAll(photonClusters); // could match [electrons]
+ mixedClusters.addAll(smallClusters);
+ Map<Track,Cluster> tracksMatchedToClusters = new HashMap<Track,Cluster>();
+ Map<Cluster, List<Track>> clustersMatchedToTracks = new HashMap<Cluster, List<Track>>();
+ for (Track tr : trackList) {
+ Cluster matchedCluster = m_trackClusterMatcher.matchTrackToCluster(tr, mixedClusters);
+ if (matchedCluster != null) {
+ tracksMatchedToClusters.put(tr, matchedCluster);
+ List<Track> clusTrList = clustersMatchedToTracks.get(matchedCluster);
+ if (clusTrList == null) { clusTrList = new Vector<Track>(); clustersMatchedToTracks.put(matchedCluster, clusTrList); }
+ clusTrList.add(tr);
+ }
+ }
+
+ // It can happen that some tracks don't match because they point to unused hits.
+ List<Track> unmatchedTracks = new Vector<Track>();
+ unmatchedTracks.addAll(trackList);
+ unmatchedTracks.removeAll(tracksMatchedToClusters.keySet());
+
+ // Some debug printout
+ System.out.println("DEBUG: "+unmatchedTracks.size()+" unmatched tracks remaining.");
+ 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 (smallClusters.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)";
+ System.out.println(printme);
+ }
+
+ // OK. Now pick out those clusters that have something connected to them...
+ Collection<Cluster> seeds = clustersMatchedToTracks.keySet();
+ Map<Cluster, Double> seedClusterTrackEnergy = new HashMap<Cluster,Double>();
+ for (Cluster seed : seeds) {
+ List<Track> seedTracks = clustersMatchedToTracks.get(seed);
+ double energySum = 0.0;
+ for (Track tr : seedTracks) {
+ double massSq = 0.140 * 0.140;
+ double[] threeMomentum = tr.getMomentum();
+ double momentumSq = threeMomentum[0]*threeMomentum[0] + threeMomentum[1]*threeMomentum[1] + threeMomentum[2]*threeMomentum[2];
+ double energySq = momentumSq + massSq;
+ energySum += Math.sqrt(energySq);
+ }
+ seedClusterTrackEnergy.put(seed, new Double(energySum));
+ }
+
+ // Set up initial thresholds
+ Map<Cluster,Double> mapSeedsToThresholds = new HashMap<Cluster,Double>();
+ for (Cluster seed : seeds) {
+ mapSeedsToThresholds.put(seed, new Double(0.9));
+ }
+
+ // Need to add other clusters to them somewhat intelligently...
+
+ Map<Cluster,Cluster> mapPhotonClustersToSeeds = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> mapSmallClustersToSeeds = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> mapLargeClustersToSeeds = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> mapMipsToSeeds = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> mapClumpsToSeeds = new HashMap<Cluster,Cluster>();
+ Map<Cluster,Cluster> mapHitsToSeeds = new HashMap<Cluster,Cluster>();
+
+ for (Cluster clus : photonClusters) {
+ if (seeds.contains(clus)) {
+ mapPhotonClustersToSeeds.put(clus,clus);
+ } else {
+ mapPhotonClustersToSeeds.put(clus,null);
+ }
+ }
+ for (Cluster clus : smallClusters) {
+ if (seeds.contains(clus)) {
+ mapSmallClustersToSeeds.put(clus,clus);
+ } else {
+ mapSmallClustersToSeeds.put(clus,null);
+ }
+ }
+ for (Cluster clus : largeClusters) {
+ if (seeds.contains(clus)) {
+ mapLargeClustersToSeeds.put(clus,clus);
+ } else {
+ mapLargeClustersToSeeds.put(clus,null);
+ }
+ }
+ for (Cluster clus : mips) {
+ if (seeds.contains(clus)) {
+ mapMipsToSeeds.put(clus,clus);
+ } else {
+ mapMipsToSeeds.put(clus,null);
+ }
+ }
+ for (Cluster clus : clumps) {
+ if (seeds.contains(clus)) {
+ mapClumpsToSeeds.put(clus,clus);
+ } else {
+ mapClumpsToSeeds.put(clus,null);
+ }
+ }
+ for (Cluster clus : wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons) {
+ if (seeds.contains(clus)) {
+ mapHitsToSeeds.put(clus,clus);
+ } else {
+ mapHitsToSeeds.put(clus,null);
+ }
+ }
+
+
+
+
+
+ printStatus("Initial state:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ // First pass: assign unambiguous split skeletons
+ for (Cluster skel : splitSkeletonClusters) {
+ List<Cluster> seedsInSkeleton = new Vector<Cluster>();
+ for (Cluster clus : skel.getClusters()) {
+ if (seeds.contains(clus)) {
+ seedsInSkeleton.add(clus);
+ }
+ }
+ if (seedsInSkeleton.size() == 1) {
+ // unique
+ Cluster seed = seedsInSkeleton.get(0);
+ for (Cluster clus : skel.getClusters()) {
+ if (mips.contains(clus)) {
+ mapMipsToSeeds.put(clus, seed);
+ } else if (clumps.contains(clus)) {
+ mapClumpsToSeeds.put(clus, seed);
+ } else {
+ throw new AssertionError("Internal consistency failure");
+ }
+ }
+ }
+ }
+
+ printStatus("After first pass:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ // Second pass: attach small clusters and hits to nearby seeds OR skeleton components OR large clusters OR photon clusters
+ // Create maps from small thing to big thing
+ double threshold = 50.0;
+ Map<Cluster, Cluster> mapAttachSmallClusters = new HashMap<Cluster,Cluster>();
+ Map<Cluster, Cluster> mapAttachHits = new HashMap<Cluster,Cluster>();
+ Map<List<Cluster>, Map<Cluster,Cluster>> targetMap = new HashMap<List<Cluster>, Map<Cluster,Cluster>>(); // e.g. from "mips" to "mapMipsToSeeds"
+ targetMap.put(mips, mapMipsToSeeds);
+ targetMap.put(clumps, mapClumpsToSeeds);
+ targetMap.put(largeClusters, mapLargeClustersToSeeds);
+ targetMap.put(photonClusters, mapPhotonClustersToSeeds);
+ List<Cluster> targets = new Vector<Cluster>();
+ for (List<Cluster> targetList : targetMap.keySet()) {
+ targets.addAll(targetList);
+ }
+
+ for (Cluster clus : smallClusters) {
+ double minDist = -1.0;
+ Cluster nearest = null;
+ for (Cluster target : targets) {
+ double dist = proximity(clus, target);
+ if (nearest == null || dist < minDist) {
+ minDist = dist;
+ nearest = target;
+ }
+ }
+ if (minDist < threshold) {
+ mapAttachSmallClusters.put(clus, nearest);
+ }
+ }
+ for (Cluster clus : wrappedUnusedHitsIgnoringLargeClustersWithoutSkeletons) {
+ double minDist = -1.0;
+ Cluster nearest = null;
+ for (Cluster target : targets) {
+ double dist = proximity(clus, target);
+ if (nearest == null || dist < minDist) {
+ minDist = dist;
+ nearest = target;
+ }
+ }
+ if (minDist < threshold) {
+ mapAttachHits.put(clus, nearest);
+ }
+ }
+
+ for (Cluster clus : mapAttachSmallClusters.keySet()) {
+ Cluster target = mapAttachSmallClusters.get(clus);
+ for (List<Cluster> targetList : targetMap.keySet()) {
+ if (targetList.contains(target)) {
+ mapSmallClustersToSeeds.put(clus, targetMap.get(targetList).get(target)); // maps clus to [seed of [target of [clus] ] ]
+ }
+ }
+ }
+
+ for (Cluster clus : mapAttachHits.keySet()) {
+ Cluster target = mapAttachHits.get(clus);
+ for (List<Cluster> targetList : targetMap.keySet()) {
+ if (targetList.contains(target)) {
+ mapHitsToSeeds.put(clus, targetMap.get(targetList).get(target)); // maps clus to [seed of [target of [clus] ] ]
+ }
+ }
+ }
+
+ printStatus("After second pass:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ // SCORING?
+ //
+ // Isolated hit to anything: Distance (1/r2), possibly relative position
+ // Small cluster to anything: Distance (1/r2), possibly relative position
+ // Skeleton component to skeleton component: tricky...
+ // Skeleton component to large cluster: tricky...
+ // Skeleton component to photon cluster: tricky...
+ // Large cluster to photon cluster or vice-versa: ?
+ //
+ // For small things, indirect proximity (small thing -> big thing -> seed) is OK.
+ // For big things, indirect proximity (big thing -> other big thing -> seed) is explosive.
+ //
+ // For skeleton components:
+ // * Assign a score to each pair (e.g. MIP <--> CLUMP)
+ // * For each valid route from seed, take product of scores,
+ // e.g. route 1: SEED <--> CLUMP = 0.4
+ // e.g. route 2: SEED <--> MIP <--> CLUMP = 0.9 * 1.0
+ // * Take maximum of scores.
+ // * For each pair, score caps at 1 -- e.g. likelihood in (1.0, 0.9) => score = 1, likelihood in (0.9, 0.0) => score = likelihood/0.9
+ // * This approach ensures that score for each indirectly linked cluster is
+ // always <= score for best path to get there. So if the distant cluster
+ // is above threshold, we are guaranteed a path to it is also above threshold.
+ //
+ // To keep things somewhat bounded, need to impose distance cut-offs.
+ // Scoring normalization: 1.0 indicates a link that you'd definitely want to make. Cap at 1.0
+
+ iterateOnce(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, targetMap, seeds);
+ printStatus("After third pass:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ adjustThresholds(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, seeds, seedClusterTrackEnergy);
+ iterateOnce(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, targetMap, seeds);
+ printStatus("After fourth pass:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ adjustThresholds(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, seeds, seedClusterTrackEnergy);
+ iterateOnce(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, targetMap, seeds);
+ printStatus("After fifth pass:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ adjustThresholds(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, seeds, seedClusterTrackEnergy);
+ iterateOnce(mapSeedsToThresholds, mapHitsToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapPhotonClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, targetMap, seeds);
+ printStatus("After sixth pass:", seeds, seedClusterTrackEnergy, mapPhotonClustersToSeeds, mapSmallClustersToSeeds, mapLargeClustersToSeeds, mapMipsToSeeds, mapClumpsToSeeds, mapHitsToSeeds, clustersMatchedToTracks, mapSeedsToThresholds);
+
+ /*
+ for (Cluster seed : seeds) {
+ List<Cluster> wholeCluster = new Vector<Cluster>();
+ wholeCluster.add(seed);
+
+ double trackEnergy = seedClusterTrackEnergy.get(seed);
+ seedMap.put(seed, wholeClsuter);
+ }
+ */
+ }
+
+ void adjustThresholds(Map<Cluster,Double> mapSeedsToThresholds,
+ Map<Cluster,Cluster> mapHitsToSeeds,
+ Map<Cluster,Cluster> mapSmallClustersToSeeds,
+ Map<Cluster,Cluster> mapLargeClustersToSeeds,
+ Map<Cluster,Cluster> mapPhotonClustersToSeeds,
+ Map<Cluster,Cluster> mapMipsToSeeds,
+ Map<Cluster,Cluster> mapClumpsToSeeds,
+ Collection<Cluster> seeds,
+ Map<Cluster,Double> seedClusterTrackEnergy)
+ {
+ for (Cluster seed : seeds) {
+ List<Cluster> matchedPhotons = reverseMap(seed, mapPhotonClustersToSeeds);
+ List<Cluster> matchedSmallClusters = reverseMap(seed, mapSmallClustersToSeeds);
+ List<Cluster> matchedLargeClusters = reverseMap(seed, mapLargeClustersToSeeds);
+ List<Cluster> matchedMips = reverseMap(seed, mapMipsToSeeds);
+ List<Cluster> matchedClumps = reverseMap(seed, mapClumpsToSeeds);
+ List<Cluster> matchedHits = reverseMap(seed, mapHitsToSeeds);
+
+ List<Cluster> combinedClusters = new Vector<Cluster>();
+ combinedClusters.addAll(matchedPhotons);
+ combinedClusters.addAll(matchedSmallClusters);
+ combinedClusters.addAll(matchedLargeClusters);
+ combinedClusters.addAll(matchedMips);
+ combinedClusters.addAll(matchedClumps);
+ combinedClusters.addAll(matchedHits);
+
+ double clusterEnergy = energy(combinedClusters);
+ double trackEnergy = seedClusterTrackEnergy.get(seed);
+ double trackEnergySigma = 0.7 * Math.sqrt(trackEnergy);
+ double normalisedResidual = (clusterEnergy - trackEnergy)/trackEnergySigma;
+
+ Double currentThreshold = mapSeedsToThresholds.get(seed);
+ double newThreshold = currentThreshold;
+
+ if ( normalisedResidual > 1.0) {
+ // Cluster is too big => increase theshold
+ double baseScale = 1.0 - currentThreshold;
+ if (baseScale > 0.1) { baseScale = 0.1; }
+ double fractionToMove = (normalisedResidual / 10.0);
+ if (fractionToMove > 0.5) { fractionToMove = 0.5; }
+ newThreshold += (baseScale * fractionToMove);
+ } else if (normalisedResidual < -1.0) {
+ // Too small => decrease threshold
+ double baseScale = 1.0 - currentThreshold;
+ if (baseScale > 0.1) { baseScale = 0.1; }
+ double fractionToMove = (-normalisedResidual / 10.0);
+ if (fractionToMove > 0.5) { fractionToMove = 0.5; }
+ newThreshold -= (baseScale * fractionToMove);
+ }
+
+ mapSeedsToThresholds.put(seed, newThreshold);
+ }
+ }
+
+ void iterateOnce(Map<Cluster,Double> mapSeedsToThresholds,
+ Map<Cluster,Cluster> mapHitsToSeeds,
+ Map<Cluster,Cluster> mapSmallClustersToSeeds,
+ Map<Cluster,Cluster> mapLargeClustersToSeeds,
+ Map<Cluster,Cluster> mapPhotonClustersToSeeds,
+ Map<Cluster,Cluster> mapMipsToSeeds,
+ Map<Cluster,Cluster> mapClumpsToSeeds,
+ Map<List<Cluster>, Map<Cluster,Cluster>> targetMap,
+ Collection<Cluster> seeds
+ ) {
+ // Compute scores
+
+ double const_hitScoreThreshold = 20.0;
+ double const_hitScoreScale = 40.0;
+ double const_hitMaxAllowedDistance = 200.0;
+
+ Map<Cluster, Map<Cluster,Double>> isolatedHitToSeedScoreMap = scoreOnProximity(const_hitScoreThreshold, const_hitScoreScale, const_hitMaxAllowedDistance, mapHitsToSeeds.keySet(), targetMap, seeds);
+ Map<Cluster, Map<Cluster,Double>> smallClusterToSeedScoreMap = scoreOnProximity(const_hitScoreThreshold, const_hitScoreScale, const_hitMaxAllowedDistance, mapSmallClustersToSeeds.keySet(), targetMap, seeds);
+ Map<Cluster, Map<Cluster,Double>> largeClusterToSeedScoreMap = scoreOnProximity(const_hitScoreThreshold, const_hitScoreScale, const_hitMaxAllowedDistance, mapLargeClustersToSeeds.keySet(), targetMap, seeds);
+ Map<Cluster, Map<Cluster,Double>> photonClusterToSeedScoreMap = scoreOnProximity(const_hitScoreThreshold, const_hitScoreScale, const_hitMaxAllowedDistance, mapPhotonClustersToSeeds.keySet(), targetMap, seeds);
+ Map<Cluster, Map<Cluster,Double>> mipToSeedScoreMap = scoreOnProximity(const_hitScoreThreshold, const_hitScoreScale, const_hitMaxAllowedDistance, mapMipsToSeeds.keySet(), targetMap, seeds);
+ Map<Cluster, Map<Cluster,Double>> clumpToSeedScoreMap = scoreOnProximity(const_hitScoreThreshold, const_hitScoreScale, const_hitMaxAllowedDistance, mapClumpsToSeeds.keySet(), targetMap, seeds);
+
+ // Update mappings based on scores
+ updateMappings(isolatedHitToSeedScoreMap, mapHitsToSeeds, mapSeedsToThresholds);
+ updateMappings(smallClusterToSeedScoreMap, mapSmallClustersToSeeds, mapSeedsToThresholds);
+ updateMappings(largeClusterToSeedScoreMap, mapLargeClustersToSeeds, mapSeedsToThresholds);
+ updateMappings(photonClusterToSeedScoreMap, mapPhotonClustersToSeeds, mapSeedsToThresholds);
+ updateMappings(mipToSeedScoreMap, mapMipsToSeeds, mapSeedsToThresholds);
+ updateMappings(clumpToSeedScoreMap, mapClumpsToSeeds, mapSeedsToThresholds);
+ }
+
+ protected void updateMappings(Map<Cluster, Map<Cluster,Double>> clusterToSeedScoreMap, Map<Cluster,Cluster> clusterToSeedMap, Map<Cluster,Double> mapSeedsToThresholds) {
+ Collection<Cluster> seeds = mapSeedsToThresholds.keySet();
+ for (Cluster clus : clusterToSeedMap.keySet()) {
+ if (seeds.contains(clus)) { continue; } // don't do anything silly
+ double maxReducedScore = 0.0;
+ Cluster bestSeedMatch = null;
+ Map<Cluster,Double> scoreMap = clusterToSeedScoreMap.get(clus);
+ for (Cluster seed : scoreMap.keySet()) {
+ Double seedScore = scoreMap.get(seed);
+ Double seedThreshold = mapSeedsToThresholds.get(seed);
+ double reducedScore = seedScore / seedThreshold;
+ if (reducedScore > 1.0) {
+ if (bestSeedMatch==null || reducedScore > maxReducedScore) {
+ maxReducedScore = reducedScore;
+ bestSeedMatch = seed;
+ }
+ }
+ }
+ clusterToSeedMap .put(clus, bestSeedMatch);
+ }
+ }
+
+ protected Map<Cluster, Map<Cluster,Double>> scoreOnProximity(double const_hitScoreThreshold, double const_hitScoreScale, double const_hitMaxAllowedDistance, Collection<Cluster> inputClusters, Map<List<Cluster>, Map<Cluster,Cluster>> targetMap, Collection<Cluster> seeds)
+ {
+ // Crummy inverse-square drop-off formula:
+ // f(x) = a / (x+b)^2
+ // f(x0) = 1
+ // f(x1) = 0.5
+ // => b = (sqrt(2)*x0 + x1) / (-1 -sqrt(2)) // taking negative solution
+ // a = (x0 + b)^2
+ double x0 = const_hitScoreThreshold;
+ double x1 = x0 + const_hitScoreScale;
+ double b = (1.41421*x0 - x1) / (-2.41421);
+ double a = (x0+b)*(x0+b);
+ Map<Cluster, Map<Cluster,Double>> scoreMap = new HashMap<Cluster, Map<Cluster,Double>>();
+ for (Cluster clus : inputClusters) {
+ if (seeds.contains(clus)) { continue; } // don't do anything silly
+ Map<Cluster, Double> scoresForThisCluster = new HashMap<Cluster,Double>();
+ for (List<Cluster> targetList : targetMap.keySet()) {
+ for (Cluster target : targetList) {
+ double dist = proximity(clus, target);
+ if (dist > const_hitMaxAllowedDistance) { continue; }
+ double score = a / ((dist + b)*(dist + b));
+ if (score > 1.0) {
+ if (dist > const_hitScoreThreshold*1.01) {
+ System.out.println("DEBUG: Dist="+dist+", a="+a+", b="+b+", so score="+score);
+ double score_x0 = a / ((x0+b)*(x0+b));
+ double score_x1 = a / ((x1+b)*(x1+b));
+ System.out.println("DEBUG: Score("+x0+") = "+score_x0+", Score("+x1+") = "+score_x1);
+ throw new AssertionError("Internal consistency failure");
+ }
+ score = 1.0;
+ }
+ Cluster seedOfTarget = targetMap.get(targetList).get(target);
+ if (seedOfTarget != null) {
+ Double previousScore = scoresForThisCluster.get(seedOfTarget);
+ if (previousScore == null || previousScore < score) {
+ scoresForThisCluster.put(seedOfTarget, new Double(score));
+ }
+ }
+ }
+ }
+ scoreMap.put(clus, scoresForThisCluster);
+ }
+ return scoreMap;
+ }
+
+ protected void printStatus(String desc, Collection<Cluster> seeds, Map<Cluster,Double> seedClusterTrackEnergy, Map<Cluster,Cluster> mapPhotonClustersToSeeds, Map<Cluster,Cluster> mapSmallClustersToSeeds, Map<Cluster,Cluster> mapLargeClustersToSeeds, Map<Cluster,Cluster> mapMipsToSeeds, Map<Cluster,Cluster> mapClumpsToSeeds, Map<Cluster,Cluster> mapHitsToSeeds, Map<Cluster,List<Track>> clustersMatchedToTracks, Map<Cluster,Double> mapSeedsToThresholds) {
+ System.out.println(desc);
+ for (Cluster seed : seeds) {
+ String printme = new String();
+ printme += " Seed cluster (track E="+seedClusterTrackEnergy.get(seed)+", threshold="+mapSeedsToThresholds.get(seed)+") --> ";
+ List<Cluster> matchedPhotons = reverseMap(seed, mapPhotonClustersToSeeds);
+ List<Cluster> matchedSmallClusters = reverseMap(seed, mapSmallClustersToSeeds);
+ List<Cluster> matchedLargeClusters = reverseMap(seed, mapLargeClustersToSeeds);
+ List<Cluster> matchedMips = reverseMap(seed, mapMipsToSeeds);
+ List<Cluster> matchedClumps = reverseMap(seed, mapClumpsToSeeds);
+ List<Cluster> matchedHits = reverseMap(seed, mapHitsToSeeds);
+ if (matchedPhotons.size() > 0) { printme += matchedPhotons.size()+" photons; "; }
+ if (matchedSmallClusters.size() > 0) { printme += matchedSmallClusters.size()+" smallClus; "; }
+ if (matchedLargeClusters.size() > 0) { printme += matchedLargeClusters.size()+" largeClus; "; }
+ if (matchedMips.size() > 0) { printme += matchedMips.size()+" mips; "; }
+ if (matchedClumps.size() > 0) { printme += matchedClumps.size()+" clumps; "; }
+ if (matchedHits.size() > 0) { printme += matchedHits.size()+" hits; "; }
+
+ List<Cluster> combinedClusters = new Vector<Cluster>();
+ combinedClusters.addAll(matchedPhotons);
+ combinedClusters.addAll(matchedSmallClusters);
+ combinedClusters.addAll(matchedLargeClusters);
+ combinedClusters.addAll(matchedMips);
+ combinedClusters.addAll(matchedClumps);
+ combinedClusters.addAll(matchedHits);
+ printme += " --> energy = "+energy(combinedClusters);
+ List<Track> trackList = clustersMatchedToTracks.get(seed);
+ if (trackList.size()>0) {
+ printme += " (effic="+quoteEfficiency(trackList, combinedClusters)+", purity="+quotePurity(trackList, combinedClusters)+")";
+ }
+
+ System.out.println(printme);
+ }
+
+ List<Cluster> unassignedPhotonClusters = new Vector<Cluster>();
+ List<Cluster> unassignedSmallClusters = new Vector<Cluster>();
+ List<Cluster> unassignedLargeClusters = new Vector<Cluster>();
+ List<Cluster> unassignedMips = new Vector<Cluster>();
+ List<Cluster> unassignedClumps = new Vector<Cluster>();
+ List<Cluster> unassignedHits = new Vector<Cluster>();
+ for (Cluster clus : mapPhotonClustersToSeeds.keySet()) {
+ if (mapPhotonClustersToSeeds.get(clus) == null) {
+ unassignedPhotonClusters.add(clus);
+ }
+ }
+ for (Cluster clus : mapSmallClustersToSeeds.keySet()) {
+ if (mapSmallClustersToSeeds.get(clus) == null) {
+ unassignedSmallClusters.add(clus);
+ }
+ }
+ for (Cluster clus : mapLargeClustersToSeeds.keySet()) {
+ if (mapLargeClustersToSeeds.get(clus) == null) {
+ unassignedLargeClusters.add(clus);
+ }
+ }
+ for (Cluster clus : mapMipsToSeeds.keySet()) {
+ if (mapMipsToSeeds.get(clus) == null) {
+ unassignedMips.add(clus);
+ }
+ }
+ for (Cluster clus : mapClumpsToSeeds.keySet()) {
+ if (mapClumpsToSeeds.get(clus) == null) {
+ unassignedClumps.add(clus);
+ }
+ }
+ for (Cluster clus : mapHitsToSeeds.keySet()) {
+ if (mapHitsToSeeds.get(clus) == null) {
+ unassignedHits.add(clus);
+ }
+ }
+ System.out.println("Unassigned: "+unassignedPhotonClusters.size()+" photons, "+unassignedSmallClusters.size()+" small clusters, "+unassignedLargeClusters.size()+" large clusters, "+unassignedMips.size()+" mips, "+unassignedClumps.size()+" clumps, "+unassignedHits.size()+" hits.");
+ }
+
+ protected List<Cluster> reverseMap(Cluster inputClus, Map<Cluster,Cluster> map) {
+ List<Cluster> output = new Vector<Cluster>();
+ for (Cluster clus : map.keySet()) {
+ if (map.get(clus) == inputClus) {
+ output.add(clus);
+ }
+ }
+ return output;
+ }
+
+ protected double energy(List<Cluster> inputClusters) {
+ Cluster magicCombinedCluster = makeCombinedCluster(inputClusters);
+ return m_calib.getEnergy(magicCombinedCluster);
+ }
+ protected Cluster makeCombinedCluster(List<Cluster> inputClusters) {
+ BasicCluster magicCombinedCluster = new BasicCluster();
+ for (Cluster clus : inputClusters) {
+ magicCombinedCluster.addCluster(clus);
+ }
+ return magicCombinedCluster;
+ }
+
+ 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"); }
+ double minDist = 0;
+ 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;
+ }
+ }
+ }
+ return minDist;
+ }
+
+ protected double quoteEfficiency(List<Track> trackList, List<Cluster> clusterList) {
+ Cluster combined = makeCombinedCluster(clusterList);
+ return quoteEfficiency(trackList, combined);
+ }
+ protected double quotePurity(List<Track> trackList, List<Cluster> clusterList) {
+ Cluster combined = makeCombinedCluster(clusterList);
+ return quotePurity(trackList, combined);
+ }
+ protected double quoteEfficiency(List<Track> trackList, Cluster clus) {
+ HitMap hitsEcal = ((HitMap)(m_event.get("inputHitMapEcal")));
+ HitMap hitsHcal = ((HitMap)(m_event.get("inputHitMapHcal")));
+ List<CalorimeterHit> allHits = new Vector<CalorimeterHit>();
+ allHits.addAll(hitsEcal.values());
+ allHits.addAll(hitsHcal.values());
+ Set<CalorimeterHit> truthHitsInEvent = findHitsFromTruth(trackList, allHits);
+ Set<CalorimeterHit> truthHitsInCluster = findHitsFromTruth(trackList, clus.getCalorimeterHits());
+ double denom = truthHitsInEvent.size();
+ double num = truthHitsInCluster.size();
+ return num/denom;
+ }
+ protected double quotePurity(List<Track> trackList, Cluster clus) {
+ Set<CalorimeterHit> truthHitsInCluster = findHitsFromTruth(trackList, clus.getCalorimeterHits());
+ double num = truthHitsInCluster.size();
+ double denom = clus.getCalorimeterHits().size();
+ return num/denom;
+ }
+
+ Set<CalorimeterHit> findHitsFromTruth(List<Track> trackList, List<CalorimeterHit> hits) {
+ Set<CalorimeterHit> outputHits = new HashSet<CalorimeterHit>();
+ for (Track tr : trackList) {
+ List<CalorimeterHit> foundHits = findHitsFromTruth(tr, hits);
+ outputHits.addAll(foundHits);
+ }
+ return outputHits;
+ }
+
+ List<CalorimeterHit> findHitsFromTruth(Track tr, List<CalorimeterHit> hits) {
+ MCParticle truthPart = null;
+ if (tr instanceof ReconTrack) {
+ truthPart = (MCParticle)(((ReconTrack)(tr)).getMCParticle());
+ } else if (tr instanceof BaseTrackMC) {
+ truthPart = ((BaseTrackMC)(tr)).getMCParticle();
+ }
+ return findHitsFromTruth(truthPart, hits);
+ }
+ List<CalorimeterHit> findHitsFromTruth(MCParticle part, List<CalorimeterHit> hitsToSearch) {
+ List<MCParticle> mcList = m_event.get(MCParticle.class, m_mcList);
+ List<MCParticle> truthParents = findParentsInList(part, mcList);
+ List<CalorimeterHit> truthHits = new Vector<CalorimeterHit>();
+ for (CalorimeterHit hit : hitsToSearch) {
+ Set<MCParticle> contributingParticles = findMCParticles(hit, mcList);
+ boolean match = false;
+ for (MCParticle particleInHit : contributingParticles) {
+ if (truthParents.contains(particleInHit)) {
+ match = true;
+ }
+ }
+ if (match) {
+ truthHits.add(hit);
+ }
+ }
+ return truthHits;
+ }
+
+ protected Set<MCParticle> findMCParticles(CalorimeterHit hit, List<MCParticle> mcList)
+ {
+ if ( ! (hit instanceof SimCalorimeterHit) ) {
+ throw new AssertionError("Non-simulated hit!");
+ } else {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) (hit);
+ Set<MCParticle> contributingParticlesFromList = new HashSet<MCParticle>();
+ int nContributingParticles = simHit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ MCParticle part = simHit.getMCParticle(i);
+ List<MCParticle> parentsInList = findParentsInList(part, mcList);
+ contributingParticlesFromList.addAll(parentsInList);
+ }
+ return contributingParticlesFromList;
+ }
+ }
+ protected List<MCParticle> findParentsInList(MCParticle part, List<MCParticle> mcList)
+ {
+ List<MCParticle> outputList = new Vector<MCParticle>();
+ if (mcList.contains(part)) {
+ // Already in there
+ outputList.add(part);
+ } else {
+ // Not in there -- recurse up through parents
+ List<MCParticle> parents = part.getParents();
+ if (parents.size()==0) {
+ // Ran out of options -- add nothing and return below
+ } else {
+ for (MCParticle parent : parents) {
+ List<MCParticle> ancestorsInList = findParentsInList(parent, mcList);
+ outputList.addAll(ancestorsInList);
+ }
+ }
+ }
+ return outputList;
+ }
+
+
+ TrackClusterMatcher m_trackClusterMatcher;
+
+ // Ugly...
+ DetailedNeutralHadronClusterEnergyCalculator m_calib = new DetailedNeutralHadronClusterEnergyCalculator();
+}
