2 added + 1 modified, total 3 files
lcsim/src/org/lcsim/contrib/uiowa/template
diff -N ChargedHadronIdentifier.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ChargedHadronIdentifier.java 27 Jan 2006 23:54:01 -0000 1.1
@@ -0,0 +1,336 @@
+package template;
+
+import java.util.*;
+import hep.physics.vec.*;
+
+import org.lcsim.event.Track;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+
+public class ChargedHadronIdentifier extends Driver
+{
+ public ChargedHadronIdentifier() {
+ }
+
+ public void setInputTrackList(String name) { m_inputTrackListName = name; }
+ public void setOutputTrackList(String name){ m_outputTrackListName = name; }
+ public void setInputMIPList(String name){ m_inputMIPListName = name; }
+ public void setInputClusterList(String name){ m_inputClusterListName = name; }
+ public void setOutputParticleList(String name){ m_outputParticleListName = name; }
+
+ public void process(EventHeader event)
+ {
+ initGeometry(event);
+
+ // Inputs:
+ List<Track> inputTrackList = event.get(Track.class, m_inputTrackListName);
+ List<Cluster> inputMIPList = event.get(Cluster.class, m_inputMIPListName);
+ List<Cluster> inputClusterList = event.get(Cluster.class, m_inputClusterListName);
+
+ // Outputs:
+ List<Track> outputTrackList = new Vector<Track>(); // initially empty
+ List<ReconstructedParticle> outputParticleList = new Vector<ReconstructedParticle>();
+
+ // Get the field strength using Steve McGill's code:
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ double[] fieldStrength = det.getFieldMap().getField(zero);
+ HelixSwimmer swimmer = new HelixSwimmer(fieldStrength[2]);
+
+ // For each track, we'll try to extrapolate it to the ECAL surface.
+ // If we do, we store the appropriate value of alpha (distance to swim)
+ Map<Track, Double> successfullyExtrapolatedTracks = new HashMap<Track,Double>();
+
+ for (Track tr : inputTrackList) {
+ // Make a swimmer:
+ swimmer.setTrack(tr);
+ double alpha = Double.NaN;
+ // // Try swimming to the barrel:
+ double alphaBarrel = swimToBarrel(swimmer);
+ boolean validBarrel = false;
+ // Try swimming to the endcap:
+ double alphaEndcap = swimToEndcap(swimmer);
+ boolean validEndcap = false;
+ // Fixme: Here we should check that the track really does go all the
+ // way to the ECAL instead of stopping/decaying/interacting earlier.
+ // This used to be done with the checkDecayPoint() method in
+ // contrib.uiowa.structural.SwimToECAL
+ if (isValidBarrelIntercept(swimmer, alphaBarrel)) {
+ alpha = alphaBarrel;
+ validBarrel = true;
+ } else if (isValidEndcapIntercept(swimmer, alphaEndcap)) {
+ alpha = alphaEndcap;
+ validEndcap = true;
+ }
+
+ if ( ! Double.isNaN(alpha) ) {
+ // Found something.
+ if ( validEndcap || validBarrel) {
+ if ( !(validEndcap && validBarrel) ) {
+ successfullyExtrapolatedTracks.put(tr, new Double(alpha));
+ } else {
+ throw new AssertionError("Confusing intercept! barrel="+validBarrel+", endcap="+validEndcap);
+ }
+ } else {
+ throw new AssertionError("alpha is not NaN, but not a valid barrel or endcap intercept");
+ }
+ }
+ }
+
+ // See: structural.MatchHelixToCluster
+ for (Track tr : successfullyExtrapolatedTracks.keySet()) {
+ Double alpha = successfullyExtrapolatedTracks.get(tr);
+ Cluster matchedCluster = null; // This will be from inputClusterList
+ // First try to match it to a MIP...
+ Cluster matchedMIP = findMatchedMIP(tr, swimmer, alpha.doubleValue(), inputMIPList);
+ if (matchedMIP != null) {
+ // Matching MIP. Look up its parent cluster:
+ for (Cluster parent : inputClusterList) {
+ List<Cluster> clustersInsideParent = recursivelyFindSubClusters(parent);
+ if (clustersInsideParent.contains(matchedMIP)) {
+ // Found it
+ matchedCluster = parent;
+ break;
+ }
+ }
+ } else {
+ // If that didn't work, try any cluster...
+ matchedCluster = findMatchedCluster(tr, swimmer, alpha.doubleValue(), inputClusterList);
+ }
+
+ if (matchedCluster != null) {
+ // We found a match
+ // Make a new ReconstructedParticle
+ BasicReconstructedParticle part = new BasicReconstructedParticle();
+ part.addTrack(tr);
+ part.addCluster(matchedCluster);
+ outputParticleList.add(part);
+ } else {
+ // No match -- put the track into the list of unused tracks
+ outputTrackList.add(tr);
+ }
+ }
+
+ // Write out
+ event.put(m_outputTrackListName, outputTrackList);
+ event.put(m_outputParticleListName, outputParticleList);
+ }
+
+ 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
+ return swimmer.getDistanceToZ(m_ECAL_endcap_z);
+ }
+ protected boolean isValidBarrelIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have -m_ECAL_barrel_z <= z <= +m_ECAL_barrel_z
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double z = intercept.z();
+ boolean zInRange = (z >= m_ECAL_barrel_zmin && z <= m_ECAL_barrel_zmax);
+ return zInRange;
+ }
+ protected boolean isValidEndcapIntercept(HelixSwimmer swimmer, double alpha) {
+ // Must have m_ECAL_endcap_rmin <= r <= m_ECAL_endcap_rmax
+ Hep3Vector intercept = swimmer.getPointAtDistance(alpha);
+ double r = Math.sqrt(intercept.x()*intercept.x() + intercept.y()*intercept.y());
+ boolean rInRange = (r >= m_ECAL_endcap_rmin && r <= m_ECAL_endcap_rmax);
+ return rInRange;
+ }
+ protected List<Cluster> recursivelyFindSubClusters(Cluster clus)
+ {
+ List<Cluster> output = new Vector<Cluster>();
+ for (Cluster dau : clus.getClusters()) {
+ output.addAll(recursivelyFindSubClusters(dau));
+ }
+ output.add(clus);
+ return output;
+ }
+
+ protected void initGeometry(EventHeader event)
+ {
+ if (!m_init) {
+ 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();
+ m_init = true;
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_zmin="+m_ECAL_barrel_zmin);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_zmax="+m_ECAL_barrel_zmax);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_barrel_r="+m_ECAL_barrel_r);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_z="+m_ECAL_endcap_z);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_rmin="+m_ECAL_endcap_rmin);
+ System.out.println(this.getClass().getName()+": Init: m_ECAL_endcap_rmax="+m_ECAL_endcap_rmax);
+ }
+ }
+
+ protected Cluster findMatchedMIP(Track tr, HelixSwimmer swimmer, double alpha, List<Cluster> mips)
+ {
+ // Find the track intercept and direction
+ swimmer.setTrack(tr);
+ Hep3Vector trackPoint = swimmer.getPointAtDistance(alpha);
+ // Obtain the unit vector giving the tangent:
+ double delta = 0.1;
+ if (alpha < 0) { delta *= -1.0; }
+ Hep3Vector aLittleFurther = swimmer.getPointAtDistance(alpha+delta);
+ Hep3Vector tangent = VecOp.unit(VecOp.sub(aLittleFurther, trackPoint));
+
+ List<Cluster> nearestMIPs = findNearestClusters(trackPoint, mips);
+ for (Cluster nearbyMIP : nearestMIPs) {
+ // Obtain geometrical info:
+ CalorimeterHit nearestHit = findNearestHit(trackPoint, nearbyMIP);
+ double separation = proximity(trackPoint, nearestHit);
+ int firstLayerHit = getLayer(nearestHit);
+ double unitDotProduct = findUnitDotProduct(tangent, nearbyMIP);
+ org.lcsim.geometry.Subdetector subdet = nearestHit.getSubdetector();
+ // Make cuts:
+ boolean goodSubDet = (subdet.getName().compareTo("EMBarrel")==0) || (subdet.getName().compareTo("EMEndcap")==0);
+ boolean goodFirstLayer = (firstLayerHit < 5);
+ boolean goodDotProduct = (Math.abs(unitDotProduct) > 0.85);
+ boolean goodSeparation = (separation < 50.0);
+ boolean foundMatch = goodSubDet && goodFirstLayer && goodDotProduct && goodSeparation;
+ if (foundMatch) {
+ // OK, made a good match
+ return nearbyMIP;
+ }
+ }
+ // No match
+ return null;
+ }
+
+ protected Cluster findMatchedCluster(Track tr, HelixSwimmer swimmer, double alpha, List<Cluster> clusters)
+ {
+ // Find the track intercept and direction
+ swimmer.setTrack(tr);
+ Hep3Vector trackPoint = swimmer.getPointAtDistance(alpha);
+
+ List<Cluster> nearestClusters = findNearestClusters(trackPoint, clusters);
+ for (Cluster nearbyCluster : nearestClusters) {
+ // Obtain geometrical info:
+ CalorimeterHit nearestHit = findNearestHit(trackPoint, nearbyCluster);
+ double separation = proximity(trackPoint, nearestHit);
+ int firstLayerHit = getLayer(nearestHit);
+ org.lcsim.geometry.Subdetector subdet = nearestHit.getSubdetector();
+ // Make cuts:
+ boolean goodSubDet = (subdet.getName().compareTo("EMBarrel")==0) || (subdet.getName().compareTo("EMEndcap")==0);
+ boolean goodFirstLayer = (firstLayerHit < 5);
+ boolean goodSeparation = (separation < 30.0);
+ boolean foundMatch = goodSubDet && goodFirstLayer && goodSeparation;
+ if (foundMatch) {
+ // OK, made a good match
+ return nearbyCluster;
+ }
+ }
+ // No match
+ return null;
+ }
+
+ protected List<Cluster> findNearestClusters(Hep3Vector point, List<Cluster> clusterList)
+ {
+ Map<Cluster,Double> mapClusterToDistance = new HashMap<Cluster, Double>();
+ List<Cluster> sortedListOfClusters = new Vector<Cluster>();
+ for (Cluster clus : clusterList) {
+ double dist = proximity(point, clus);
+ mapClusterToDistance.put(clus, new Double(dist));
+ sortedListOfClusters.add(clus);
+ }
+ Comparator<Cluster> comp = new CompareMapping<Cluster>(mapClusterToDistance);
+ Collections.sort(sortedListOfClusters, comp);
+ return sortedListOfClusters;
+ }
+ protected CalorimeterHit findNearestHit(Hep3Vector point, Cluster clus)
+ {
+ CalorimeterHit nearest = null;
+ double minDist = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ if (distance<minDist || nearest==null) {
+ nearest = hit;
+ }
+ }
+ return nearest;
+ }
+
+ protected double proximity(Hep3Vector point, Cluster clus) {
+ CalorimeterHit nearestHit = findNearestHit(point, clus);
+ return proximity(point, nearestHit);
+ }
+ protected double proximity(Hep3Vector point, CalorimeterHit hit) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ return distance;
+ }
+
+ protected double findUnitDotProduct(Hep3Vector tangent, Cluster clus)
+ {
+ // Find the cluster direction
+ BasicCluster copy = new BasicCluster();
+ copy.addCluster(clus);
+ TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+ copy.setPropertyCalculator(calc);
+ copy.calculateProperties();
+ double[][]axes = calc.getPrincipleAxis();
+ Hep3Vector clusterDir = new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]);
+ // Get the dot product:
+ double unitDotProduct = VecOp.dot(tangent, clusterDir) / (tangent.magnitude() * clusterDir.magnitude());
+ return unitDotProduct;
+ }
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+ private class CompareMapping<T> implements Comparator<T> {
+ public CompareMapping(Map<T,Double> map) {
+ m_map = map;
+ }
+ public int compare(Object o1, Object o2) {
+ Cluster c1 = (Cluster) o1;
+ Cluster c2 = (Cluster) o2;
+ Double D1 = m_map.get(c1);
+ Double D2 = m_map.get(c2);
+ if (D1.equals(D2)) {
+ // Equal
+ return 0;
+ } else if (D1.doubleValue() < D2.doubleValue()) {
+ return -1;
+ } else {
+ return +1;
+ }
+ }
+ Map<T,Double> m_map;
+ }
+
+ 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;
+
+ String m_inputTrackListName;
+ String m_outputTrackListName;
+ String m_inputMIPListName;
+ String m_inputClusterListName;
+ String m_outputParticleListName;
+}
+
+
lcsim/src/org/lcsim/contrib/uiowa/template
diff -N HaloAssigner.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HaloAssigner.java 27 Jan 2006 23:54:01 -0000 1.1
@@ -0,0 +1,114 @@
+package template;
+
+import java.util.*;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+public class HaloAssigner extends Driver
+{
+ public HaloAssigner(String inputClusterListName, String inputHitMapName, String outputClusterListName, String outputHitMapName) {
+ m_inputClusterListName = inputClusterListName;
+ m_inputHitMapName = inputHitMapName;
+ m_outputClusterListName = outputClusterListName;
+ m_outputHitMapName = outputHitMapName;
+ }
+
+ public void process(EventHeader event)
+ {
+ // Get the input clusters and hits:
+ List<Cluster> inputSkeletons = event.get(Cluster.class, m_inputClusterListName);
+ Map<Long, CalorimeterHit> inputHitMap = (Map<Long, CalorimeterHit>) (event.get(m_inputHitMapName));
+
+ // Prepare the outputs:
+ Map<Long, CalorimeterHit> outputHitMap = new HashMap<Long,CalorimeterHit>(inputHitMap); // initially full
+ List<Cluster> outputClusterList = new Vector<Cluster>();
+
+ // For all unused hits, assign the hit to a cluster.
+ // We do this as a two-step process to ensure that the
+ // we're independent of the hit/cluster ordering.
+
+ // First pass: find assignments
+ Map<Cluster, List<CalorimeterHit>> foundAssignments = new HashMap<Cluster, List<CalorimeterHit>>();
+ for (CalorimeterHit hit : inputHitMap.values()) {
+ Cluster bestMatch = findBestCluster(hit, inputSkeletons);
+ List<CalorimeterHit> matchedHits = foundAssignments.get(bestMatch);
+ if (matchedHits == null) {
+ matchedHits = new Vector<CalorimeterHit>();
+ foundAssignments.put(bestMatch, matchedHits);
+ // Assigned => remove from output hitmap
+ outputHitMap.remove(hit.getCellID());
+ }
+ matchedHits.add(hit);
+ }
+
+ // Second pass: create expanded clusters
+ for (Cluster inputCluster : inputSkeletons) {
+ BasicCluster outputCluster = new BasicCluster();
+ outputCluster.addCluster(inputCluster);
+ List<CalorimeterHit> matchedHits = foundAssignments.get(inputCluster);
+ if (matchedHits != null) {
+ for (CalorimeterHit matchedHit : matchedHits) {
+ outputCluster.addHit(matchedHit);
+ }
+ }
+ outputClusterList.add(outputCluster);
+ }
+
+ event.put(m_outputClusterListName, outputClusterList);
+ event.put(m_outputHitMapName, outputHitMap);
+ }
+
+ protected Cluster findBestCluster(CalorimeterHit hit, List<Cluster> clusters)
+ {
+ Cluster bestMatch = null;
+ double minDistance = 0;
+ for (Cluster subCluster : clusters) {
+ // How far to this cluster?
+ double dist = distance(subCluster, hit);
+ if (bestMatch == null || dist < minDistance) {
+ minDistance = dist;
+ bestMatch = subCluster;
+ }
+ }
+ return bestMatch;
+ }
+
+ // This belongs outside this class.
+ private double distance(Cluster clus, CalorimeterHit hit)
+ {
+ // Loop over hits...
+ boolean firstCheck = true;
+ double minDistance = Double.NaN; // Will stay NaN if clus is empty
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+ for (CalorimeterHit hitInCluster : hits) {
+ double dist = distance(hit, hitInCluster);
+ if (firstCheck || dist<minDistance) {
+ minDistance = dist;
+ firstCheck = false;
+ }
+ }
+
+ return minDistance;
+ }
+
+ // This belongs outside this class.
+ private double distance(CalorimeterHit hit1, CalorimeterHit hit2)
+ {
+ Hep3Vector vect1 = new BasicHep3Vector(hit1.getPosition());
+ Hep3Vector vect2 = new BasicHep3Vector(hit2.getPosition());
+ Hep3Vector displacement = VecOp.sub(vect1, vect2);
+ return displacement.magnitude();
+ }
+ String m_inputClusterListName;
+ String m_inputHitMapName;
+ String m_outputClusterListName;
+ String m_outputHitMapName;
+}
lcsim/src/org/lcsim/contrib/uiowa/template
diff -u -r1.2 -r1.3
--- NonTrivialPFA.java 26 Jan 2006 01:28:00 -0000 1.2
+++ NonTrivialPFA.java 27 Jan 2006 23:54:01 -0000 1.3
@@ -133,7 +133,22 @@
// Outputs:
// * Fleshed-out skeletons (with halo added)
// * Modified hitmap with any remaining clustered hits
- // Step 6c: Extrapolate tracks to ECAL surface
+ add(new template.HaloAssigner("skeletons", "structural unused hits", "skeletons plus halo", "structural unused hits minus halo"));
+ // Step 6c: Extrapolate tracks to ECAL surface, form charged particles
+ // Inputs:
+ // * The tracks (which we'll extrapolate)
+ // * The MIPs (which are our best chance of a good association)
+ // * The fleshed-out skeletons (to match the tracks to -- or to their MIP components -- and form particles)
+ // Outputs:
+ // * Any unused tracks
+ // * The reconstructed particles
+ ChargedHadronIdentifier hadID = new ChargedHadronIdentifier();
+ hadID.setInputTrackList(EventHeader.TRACKS);
+ hadID.setOutputTrackList("leftover tracks");
+ hadID.setInputMIPList("mips");
+ hadID.setInputClusterList("skeletons plus halo");
+ hadID.setOutputParticleList("charged hadron particles");
+ add(hadID);
// Step 6d: Handle fragments
// Step 6e: Plots
}
@@ -146,6 +161,8 @@
add(new DebugInfoHitMap("hit map hcal without mips or clumps"));
add(new DebugInfoHitMap("ecal hit map after mst"));
add(new DebugInfoHitMap("hcal hit map after mst"));
+ add(new DebugInfoHitMap("structural unused hits"));
+ add(new DebugInfoHitMap("structural unused hits minus halo"));
add(new DebugInfoClusterList("mips ecal"));
add(new DebugInfoClusterList("mips hcal"));
@@ -156,10 +173,27 @@
add(new DebugInfoClusterList("mst clusters ecal"));
add(new DebugInfoClusterList("mst clusters hcal"));
add(new DebugInfoClusterList("mst clusters linked"));
+ add(new DebugInfoClusterList("skeletons"));
+ add(new DebugInfoClusterList("skeletons plus halo"));
add(new DebugInfoTrackList(EventHeader.TRACKS));
+ add(new DebugInfoTrackList("leftover tracks"));
+
+ add(new DebugInfoParticleList("charged hadron particles"));
}
+ public void process(EventHeader event) {
+ // Special handling of things that need per-event info.
+ // This probably belongs in its own driver.
+ if (m_perEventQuantities != null) {
+ for (StructuralLikelihoodQuantityWithEventInfo quant : m_perEventQuantities) {
+ quant.setEventInfo(event);
+ }
+ }
+ super.process(event);
+ }
+
+
protected void makeEventInfoList(LikelihoodEvaluator eval)
{
// Handle things that have per-event info:
CVSspam 0.2.8
