lcsim/src/org/lcsim/contrib/uiowa
diff -N CheatTrackClusterMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CheatTrackClusterMatcher.java 21 Sep 2007 02:29:05 -0000 1.1
@@ -0,0 +1,389 @@
+//package org.lcsim.recon.pfa.identifier;
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import hep.physics.vec.*;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.Track;
+import org.lcsim.util.*;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.recon.ztracking.cheater.CheatTrack;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import hep.physics.particle.Particle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.pfa.identifier.TrackClusterMatcher;
+
+public class CheatTrackClusterMatcher extends Driver implements TrackClusterMatcher
+{
+ String m_mcListName;
+
+ public CheatTrackClusterMatcher(String mcList) {
+ m_mcListName = mcList;
+ }
+
+ protected EventHeader m_event;
+ public void process(EventHeader event) {
+ m_event = event;
+ initGeometry(event);
+ }
+
+ // Main interface
+ public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters)
+ {
+ // Need truth info for this track
+ Particle truth = null;
+ if (tr instanceof ReconTrack) {
+ ReconTrack truthTrack = (ReconTrack) (tr);
+ truth = truthTrack.getMCParticle(); // returns a Particle
+ } else if (tr instanceof CheatTrack) {
+ CheatTrack truthTrack = (CheatTrack) (tr);
+ truth = truthTrack.getMCParticle(); // returns a MCParticle
+ } else {
+ throw new AssertionError("Don't know how to extract cheat information from this track of class "+tr.getClass().getName());
+ }
+ if (truth == null) {
+ // Track has no truth information
+ throw new AssertionError("CheatTrackClusterMatcher called for Track with no truth information");
+ }
+ MCParticle mctruth = (MCParticle) (truth);
+
+ if (m_debug) { System.out.println("DEBUG: Attempting to match a track with p="+new BasicHep3Vector(tr.getMomentum()).magnitude()+" (from an MCParticle of type "+mctruth.getType()+" and p="+mctruth.getMomentum().magnitude()+") to a list of "+clusters.size()+" clusters."); }
+
+ // Need truth info for whole event
+ List<MCParticle> mcList = m_event.get(MCParticle.class, m_mcListName);
+
+ // Look up truth parents. If the tracks and MC list use the same truth
+ // list then this should be a 1:1 mapping.
+ List<MCParticle> truthParentsInList = findParentsInList(mctruth, mcList);
+
+ // Check which clusters contributions from the MC parents:
+ Map<Cluster, List<CalorimeterHit>> clustersWithContributions = new HashMap<Cluster,List<CalorimeterHit>>();
+ for (Cluster clus : clusters) {
+ List<CalorimeterHit> hits = findContributedHits(clus, truthParentsInList, mcList);
+ if (hits.size()>0) {
+ clustersWithContributions.put(clus,hits);
+ }
+ }
+
+ // How many clusters did we match up with?
+ int nMatchedClusters = clustersWithContributions.keySet().size();
+ if (m_debug) { System.out.println("DEBUG: Based on truth information, matched track to "+nMatchedClusters+" clusters."); }
+ if (nMatchedClusters == 0) {
+ // No match -- for example, if track didn't reach the
+ // calorimeters.
+ if (m_debug) { System.out.println("DEBUG: No matched clusters => return NULL."); }
+ return null;
+ } else if (nMatchedClusters == 1) {
+ // Unique/exact match
+ Cluster match = clustersWithContributions.keySet().iterator().next();
+ if (m_debug) { System.out.println("DEBUG: Unique match => return cluster with "+match.getCalorimeterHits().size()+" hits."); }
+ return match;
+ } else {
+ // Non-unique match. Now we need to be cleverer...
+ Set<Cluster> matchedClusters = clustersWithContributions.keySet();
+ if (m_debug) { System.out.println("DEBUG: Looking at "+matchedClusters.size()+" matched clusters to see which ones work best..."); }
+ // Look at innermost ECAL hits.
+ // * If there are NO clusters with a hit in either/both
+ // of the first two layers, return whichever cluster has
+ // the innermost hit (not being too fussy if >1 cluster
+ // has the same innermost layer since it's a bad match
+ // in any case).
+ // * If there is ONE cluster with a hit in either/both
+ // of the first two layers, return that cluster.
+ // * If there are TWO OR MORE clusters with a hit in either/both
+ // of the first two layers, think some more.
+ Cluster clusterWithInnermostHitLayer = null;
+ CalorimeterHit globalInnermostHitInECAL = null;
+ List<Cluster> clustersWithHitsInFirstTwoLayers = new Vector<Cluster>();
+ for (Cluster clus : matchedClusters) {
+ List<CalorimeterHit> hits = clustersWithContributions.get(clus);
+ CalorimeterHit innermostHitInECAL = findInnermostHitInECAL(hits);
+ if (innermostHitInECAL != null) {
+ // There is an ECAL innermost hit...
+ int layer = getLayer(innermostHitInECAL);
+ // Check if in first two layers
+ if (layer < 2) {
+ clustersWithHitsInFirstTwoLayers.add(clus);
+ }
+ // Check if is current innermost cluster
+ if (clusterWithInnermostHitLayer == null || layer < getLayer(globalInnermostHitInECAL)) {
+ globalInnermostHitInECAL = innermostHitInECAL;
+ clusterWithInnermostHitLayer = clus;
+ }
+ }
+ if (m_debug) {
+ String printme = new String("DEBUG: Scanned a cluster with "+hits.size()+" hits: ");
+ if (innermostHitInECAL != null) {
+ printme += " Has ECAL hit(s). Innermost hit in layer ";
+ printme += getLayer(innermostHitInECAL);
+ if (getLayer(innermostHitInECAL) < 2) {
+ printme += " which is in first two layers";
+ }
+ if (clusterWithInnermostHitLayer == clus) {
+ printme += " so this is the innermost cluster (global innermost hit layer="+getLayer(globalInnermostHitInECAL)+")";
+ }
+ } else {
+ printme += " No ECAL hits.";
+ }
+ System.out.println(printme);
+ }
+ }
+ if (clustersWithHitsInFirstTwoLayers.size()==0) {
+ // No clusters with hits in first two layers => return innermost
+ // (This can be null if no clusters had any ECAL hits -- that's OK.)
+ if (m_debug) { System.out.println("DEBUG: Didn't find any clusters with hits in first two layers => returning innermost with "+clusterWithInnermostHitLayer.getCalorimeterHits().size()+" hits"); }
+ return clusterWithInnermostHitLayer;
+ } else if (clustersWithHitsInFirstTwoLayers.size()==1) {
+ // Exactly one sensible match. It should be both the innermost
+ // cluster and the only cluster with hit(s) in first two ECAL layers.
+ if (clustersWithHitsInFirstTwoLayers.get(0) != clusterWithInnermostHitLayer) { throw new AssertionError("BUG"); }
+ if (clusterWithInnermostHitLayer == null) { throw new AssertionError("BUG"); }
+ if (m_debug) { System.out.println("DEBUG: Found only one cluster with a hit in either/both of first two layers => returning that one with "+clusterWithInnermostHitLayer.getCalorimeterHits().size()+" hits"); }
+ return clusterWithInnermostHitLayer;
+ } else {
+ // Ugly case: >1 sensible match
+ if (m_debug) {
+ System.out.println("DEBUG: For track with p="+Math.sqrt(tr.getPX()*tr.getPX()+tr.getPY()*tr.getPY()+tr.getPZ()*tr.getPZ())+", there are "+clustersWithHitsInFirstTwoLayers.size()+" matched clusters with hits in first 2 ECAL layers.");
+ System.out.println("DEBUG: MCParticle is a "+mctruth.getPDGID()+" with momentum "+mctruth.getMomentum().magnitude());
+ System.out.println("DEBUG: Here are the truthParentsInList:");
+ for (MCParticle part : truthParentsInList) {
+ System.out.println("DEBUG: A "+part.getPDGID()+" with momentum "+part.getMomentum().magnitude());
+ }
+ }
+ // Make a HelixSwimmer to propagate the track
+ HelixSwimmer swimmer = new HelixSwimmer(m_fieldStrength[2]);
+ swimmer.setTrack(tr);
+ // Try swimming to the barrel:
+ double alphaBarrel = swimToBarrel(swimmer);
+ boolean validBarrel = false;
+ // Try swimming to the endcap:
+ double alphaEndcap = swimToEndcap(swimmer);
+ boolean validEndcap = false;
+ // Find intercept point
+ double alpha = Double.NaN;
+ if (isValidBarrelIntercept(swimmer, alphaBarrel)) {
+ alpha = alphaBarrel;
+ validBarrel = true;
+ } else if (isValidEndcapIntercept(swimmer, alphaEndcap)) {
+ alpha = alphaEndcap;
+ validEndcap = true;
+ }
+ Hep3Vector trackPoint = null;
+ if (Double.isNaN(alpha)) {
+ // Track extrapolation failed
+ // Not much hope here... just pick the one with the most matched hits.
+ if (m_debug) { System.out.println("DEBUG: Track extrapolation to ECAL failed."); }
+ Cluster bestMatch = null;
+ for (Cluster clus : clustersWithHitsInFirstTwoLayers) {
+ List<CalorimeterHit> hits = clustersWithContributions.get(clus);
+ int nMatchedHits = hits.size();
+ if (bestMatch==null || nMatchedHits > clustersWithContributions.get(bestMatch).size()) {
+ bestMatch = clus;
+ }
+ }
+ if (m_debug) { System.out.println("DEBUG: Returning a cluster with "+clustersWithContributions.get(bestMatch).size()+" / "+bestMatch.getCalorimeterHits().size()+" hits matched as my best guess."); }
+ return bestMatch;
+ } else {
+ trackPoint = swimmer.getPointAtDistance(alpha);
+ if (m_debug) { System.out.println("DEBUG: Extrapolated track; Track intercept point at ("+trackPoint.x()+","+trackPoint.y()+","+trackPoint.z()+")"); }
+ }
+ // Now find the nearest one...
+ Cluster nearestClusterToInterceptPoint = null;
+ CalorimeterHit nearestHitToInterceptPoint = null;
+ for (Cluster clus : clustersWithHitsInFirstTwoLayers) {
+ List<CalorimeterHit> hits = clustersWithContributions.get(clus);
+ CalorimeterHit innermostHitInECAL = findInnermostHitInECAL(hits);
+ if (m_debug) { System.out.println("DEBUG: Cluster with "+hits.size()+" / "+clus.getCalorimeterHits().size()+" hits matched, and first ECAL hit in layer "+getLayer(innermostHitInECAL)+" of "+innermostHitInECAL.getSubdetector().getName()); }
+ CalorimeterHit nearestHitInClusterToInterceptPoint = null;
+ for (CalorimeterHit hit : hits) {
+ int layer = getLayer(hit);
+ if (layer==0 || layer==1) {
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ String name = subdet.getName();
+ if (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0) {
+ // EM -- OK
+ double distance = proximity(trackPoint, hit);
+ if (m_debug) { System.out.println("DEBUG: Hit in ECAL layer "+layer+" distance from intercept point: "+distance); }
+ if (nearestHitToInterceptPoint==null || distance < proximity(trackPoint, nearestHitToInterceptPoint)) {
+ if (m_debug) { System.out.println("DEBUG: This is the new nearest hit"); }
+ nearestHitToInterceptPoint = hit;
+ nearestClusterToInterceptPoint = clus;
+ }
+ }
+ }
+ }
+ }
+ if (nearestClusterToInterceptPoint==null) { throw new AssertionError("Inconsistency"); }
+ if (m_debug) { System.out.println("DEBUG: Nearest cluster has proximity of "+proximity(trackPoint, nearestHitToInterceptPoint)+" and "+nearestClusterToInterceptPoint.getCalorimeterHits().size()+" hits. Returning it."); }
+ return nearestClusterToInterceptPoint;
+ }
+ }
+
+ }
+
+
+ // Utility routines
+ // ----------------
+
+ protected CalorimeterHit findInnermostHitInECAL(Cluster clus) {
+ return findInnermostHitInECAL(clus.getCalorimeterHits());
+ }
+
+ protected CalorimeterHit findInnermostHitInECAL(List<CalorimeterHit> hits) {
+ CalorimeterHit innermostHit = null;
+ for (CalorimeterHit hit : hits) {
+ int layer = getLayer(hit);
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0) {
+ // EM -- OK
+ if (innermostHit==null || getLayer(innermostHit)>layer) {
+ innermostHit = hit;
+ }
+ }
+ }
+ return innermostHit;
+ }
+
+ protected int getLayer(CalorimeterHit hit) {
+ if (hit==null) { throw new AssertionError("hit is null"); }
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+ protected List<CalorimeterHit> findContributedHits(Cluster clus, List<MCParticle> relevantParticles, List<MCParticle> allParticles) {
+ Vector<CalorimeterHit> output = new Vector<CalorimeterHit>();
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean relevantContributionToHit = false;
+ Set<MCParticle> particles = findMCParticles(hit, allParticles);
+ for (MCParticle part : particles) {
+ if (relevantParticles.contains(part)) {
+ relevantContributionToHit = true;
+ break;
+ }
+ }
+ if (relevantContributionToHit) {
+ output.add(hit);
+ }
+ }
+ return output;
+ }
+
+ 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;
+ }
+
+ protected double proximity(Hep3Vector point, CalorimeterHit hit) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ return distance;
+ }
+
+ protected double swimToBarrel(HelixSwimmer swimmer) {
+ // Look for a hit in the first layer of the ECAL barrel
+ return swimmer.getDistanceToRadius(m_ECAL_barrel_r);
+ }
+ protected double swimToEndcap(HelixSwimmer swimmer) {
+ // Look for a hit in the first layer of the ECAL endcap
+ double distanceToEndcap1 = swimmer.getDistanceToZ(m_ECAL_endcap_z);
+ double distanceToEndcap2 = swimmer.getDistanceToZ(-m_ECAL_endcap_z);
+ if (distanceToEndcap1>0) {
+ return distanceToEndcap1;
+ } else {
+ return distanceToEndcap2;
+ }
+ }
+ protected boolean isValidBarrelIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have -m_ECAL_barrel_z <= z <= +m_ECAL_barrel_z (within errors)
+ double uncertainty = 0.0;
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double z = intercept.z();
+ boolean zInRange = (z >= m_ECAL_barrel_zmin-uncertainty && z <= m_ECAL_barrel_zmax+uncertainty);
+ return zInRange;
+ }
+ protected boolean isValidEndcapIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have m_ECAL_endcap_rmin <= r <= m_ECAL_endcap_rmax (within errors)
+ double uncertainty = 0.0;
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double r = Math.sqrt(intercept.x()*intercept.x() + intercept.y()*intercept.y());
+ boolean rInRange = (r >= m_ECAL_endcap_rmin-uncertainty && r <= m_ECAL_endcap_rmax+uncertainty);
+ return rInRange;
+ }
+
+ public void initGeometry(EventHeader event)
+ {
+ Detector det = event.getDetector();
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ CylindricalCalorimeter eme = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ m_ECAL_barrel_zmin = emb.getZMin();
+ m_ECAL_barrel_zmax = emb.getZMax();
+ m_ECAL_barrel_r = emb.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_z = eme.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_rmin = eme.getInnerRadius();
+ m_ECAL_endcap_rmax = eme.getOuterRadius();
+ double[] zero = {0, 0, 0};
+ m_fieldStrength = det.getFieldMap().getField(zero);
+ m_init = true;
+ if (m_debug) {
+ System.out.println(this.getClass().getName()+": Init: ECAL barrel zmin="+m_ECAL_barrel_zmin);
+ System.out.println(this.getClass().getName()+": Init: ECAL barrel zmax="+m_ECAL_barrel_zmax);
+ System.out.println(this.getClass().getName()+": Init: ECAL barrel r="+m_ECAL_barrel_r);
+ System.out.println(this.getClass().getName()+": Init: ECAL endcap z="+m_ECAL_endcap_z);
+ System.out.println(this.getClass().getName()+": Init: ECAL endcap rmin="+m_ECAL_endcap_rmin);
+ System.out.println(this.getClass().getName()+": Init: ECAL endcap rmax="+m_ECAL_endcap_rmax);
+ }
+ }
+
+ public void setDebug(boolean debug) { m_debug = debug; }
+
+ protected boolean m_init = false;
+ protected double m_ECAL_barrel_zmin;
+ protected double m_ECAL_barrel_zmax;
+ protected double m_ECAL_barrel_r;
+ protected double m_ECAL_endcap_z;
+ protected double m_ECAL_endcap_rmin;
+ protected double m_ECAL_endcap_rmax;
+ protected boolean m_debug = false;
+ protected double[] m_fieldStrength;
+}
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.19 -r1.20
--- NonTrivialPFA.java 20 Sep 2007 21:11:03 -0000 1.19
+++ NonTrivialPFA.java 21 Sep 2007 02:29:05 -0000 1.20
@@ -346,7 +346,16 @@
// Check track matches are sensible; don't match multiple
// tracks to the same skeleton.
String eventSplitSkeletonClusters = "splitSkeletons";
- add(new CheckSkeletonsForMultipleTracks(evalWrapper, trackList, eventSkeletonClusters, eventSplitSkeletonClusters, eventMips, eventClumps));
+ {
+ LocalHelixExtrapolationTrackClusterMatcher extrapolate = new LocalHelixExtrapolationTrackClusterMatcher();
+ extrapolate.setCutSeparation(14.0); // about two cells
+ CheatTrackClusterMatcher cheater = new CheatTrackClusterMatcher("PerfectRecoMCParticles");
+ CheckSkeletonsForMultipleTracks separate = new CheckSkeletonsForMultipleTracks(evalWrapper, trackList, eventSkeletonClusters, eventSplitSkeletonClusters, eventMips, eventClumps, extrapolate);
+ //CheckSkeletonsForMultipleTracks separate = new CheckSkeletonsForMultipleTracks(evalWrapper, trackList, eventSkeletonClusters, eventSplitSkeletonClusters, eventMips, eventClumps, cheater);
+ add(cheater);
+ add(extrapolate);
+ add(separate);
+ }
// Add halo of nearby hits
prefix = "halo__";
@@ -720,6 +729,35 @@
add(new HaloAssigner(skeletonClusterList, inputHitMap, haloClusterList, outputHitMap));
}
+ protected void addCheatingTrackMatcher(String prefix, String trackList, String inputMIPList, String inputHaloClusterList, String inputSmallClusterList, String outputParticleList, boolean checkEoverP, ClusterEnergyCalculator calibration, boolean debug)
+ {
+ System.out.println("WARNING: Cheating on track matching");
+ // Merge input cluster lists
+ ListAddDriver<Cluster> mergeClusters = new ListAddDriver<Cluster>(Cluster.class);
+ mergeClusters.addInputList(inputHaloClusterList);
+ if (inputSmallClusterList != null) {
+ mergeClusters.addInputList(inputSmallClusterList);
+ }
+ mergeClusters.setOutputList(prefix+"inputClusterList");
+ add(mergeClusters);
+
+ // Try the clusters generically:
+ SimpleChargedParticleMaker hadID = new SimpleChargedParticleMaker();
+ //List<String> hitCollections = new Vector<String>();
+ //hitCollections.add("EcalBarrDigiHits");
+ //hitCollections.add("EcalEndcapDigiHits");
+ CheatTrackClusterMatcher clusMatch = new CheatTrackClusterMatcher("PerfectRecoMCParticles");
+ clusMatch.setDebug(false);
+ add(clusMatch);
+ hadID.setTrackMatcher(clusMatch);
+ hadID.setInputTrackList(trackList);
+ hadID.setOutputTrackList(prefix+"leftoverTracks");
+ hadID.setInputClusterList(prefix+"inputClusterList");
+ hadID.setOutputParticleList(outputParticleList);
+ hadID.setDebug(debug);
+ add(hadID);
+ }
+
protected void addTrackMatcher(String prefix, String trackList, String inputMIPList, String inputHaloClusterList, String inputSmallClusterList, String outputParticleList, boolean checkEoverP, ClusterEnergyCalculator calibration, boolean debug)
{
// First try the MIPs...
@@ -893,6 +931,7 @@
String tempChargedParticleList = prefix+"chargedHadronParticlesAfterFragmentHandling";
boolean applyEoverPcut = true;
addTrackMatcher(prefix+"FindCharged__", trackList, inputMIPList, clusterList, null, tempChargedParticleList, applyEoverPcut, calibrationForEoverP, trackDebug);
+ //addCheatingTrackMatcher(prefix+"FindCharged__", trackList, inputMIPList, clusterList, null, tempChargedParticleList, applyEoverPcut, calibrationForEoverP, trackDebug);
if (trackDebug) {
System.out.println(prefix+": Will write out charged particles as '"+tempChargedParticleList+"'");
}
