lcsim/src/org/lcsim/recon/pfa/identifier
diff -N TrackClusterMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackClusterMatcher.java 10 Aug 2006 22:57:02 -0000 1.1
@@ -0,0 +1,28 @@
+package org.lcsim.recon.pfa.identifier;
+
+import java.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.Track;
+
+/**
+ * Interface for matching Tracks to Clusters
+ *
+ * @version $Id: TrackClusterMatcher.java,v 1.1 2006/08/10 22:57:02 mcharles Exp $
+ */
+
+public interface TrackClusterMatcher
+{
+ /**
+ * Attempt to match a Track <code>tr</code> to a Cluster
+ * from the list <code>clusters</code>. The return value
+ * is the matched Cluster, or null if there is no acceptable
+ * match found.
+ *
+ * Implementations typically return a cluster from the list supplied.
+ * They may also return a new cluster which contains clusters
+ * (via Cluster.getClusters()) from this list (and no other hits/clusters).
+ *
+ * @version $Id: TrackClusterMatcher.java,v 1.1 2006/08/10 22:57:02 mcharles Exp $
+ */
+ public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters);
+}
lcsim/src/org/lcsim/recon/pfa/identifier
diff -N SimpleTrackClusterMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SimpleTrackClusterMatcher.java 10 Aug 2006 22:57:02 -0000 1.1
@@ -0,0 +1,273 @@
+package org.lcsim.recon.pfa.identifier;
+
+import java.util.*;
+import hep.physics.vec.*;
+
+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;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.util.Driver;
+
+/**
+ * Attempt to match a Track to a Cluster, based on the intercept point
+ * on the ECAL inner surface.
+ *
+ * @version $Id: SimpleTrackClusterMatcher.java,v 1.1 2006/08/10 22:57:02 mcharles Exp $
+ */
+
+public class SimpleTrackClusterMatcher extends Driver implements TrackClusterMatcher
+{
+
+ /**
+ * Match this track to a cluster from the list supplied.
+ */
+ public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters)
+ {
+ // 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;
+
+ // Fixme: Here we should check that the track really does go all the
+ // way to the ECAL instead of stopping/decaying/interacting earlier.
+ double alpha = Double.NaN;
+ if (isValidBarrelIntercept(swimmer, alphaBarrel)) {
+ alpha = alphaBarrel;
+ validBarrel = true;
+ } else if (isValidEndcapIntercept(swimmer, alphaEndcap)) {
+ alpha = alphaEndcap;
+ validEndcap = true;
+ }
+
+ // Did we make a successful extrapolation?
+ if (Double.isNaN(alpha)) {
+ // No -- failed
+ return null;
+ }
+ if ( !(validEndcap || validBarrel) ) {
+ // Invalid state
+ throw new AssertionError("Invalid state: alpha is not NaN, but not a valid barrel or endcap intercept");
+ }
+ if ( validEndcap && validBarrel ) {
+ throw new AssertionError("Invalid state: barrel="+validBarrel+", endcap="+validEndcap);
+ }
+
+ // If we reach here, we extrapolated to the barrel or
+ // to the endcap successfully.
+
+ Cluster matchedCluster = findMatchedCluster(tr, swimmer, alpha, clusters);
+ if (matchedCluster != null) {
+ // Matched OK
+ if (m_debug) {
+ Hep3Vector momentumVec = new BasicHep3Vector(tr.getMomentum());
+ System.out.println("DEBUG: Extrapolated track to cluster (momentum = "+momentumVec.magnitude()+")");
+ }
+ return matchedCluster;
+ } else {
+ // No match found
+ if (m_debug) {
+ Hep3Vector momentumVec = new BasicHep3Vector(tr.getMomentum());
+ System.out.println("DEBUG: Failed to extrapolate track (momentum = "+momentumVec.magnitude()+")");
+ }
+ return null;
+ }
+ }
+
+ 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 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) {
+ // Geometrically, it looks good.
+ // Is it sensible in terms of energy?
+
+ // // FIXME: Farm the calibration out to an external class properly
+ // double estimatedClusterEnergy = SimpleNeutralParticleMaker.estimateClusterEnergy(clus);
+ // // We don't expect an exact match due to resolution, energy lost
+ // // to fragments etc., but a good portion of the energy should be
+ // // in the cluster
+ // double[] trackMomentum = tr.getMomentum();
+ // double trackMomentumMag = Math.sqrt(trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
+ // boolean fractionEnergyOK = estimatedClusterEnergy > 0.5*trackMomentumMag;
+ // boolean absoluteEnergyOK = (trackMomentumMag - estimatedClusterEnergy < 0.2); // deliberately one-sided
+ // boolean energyOK = fractionEnergyOK || absoluteEnergyOK;
+ boolean energyOK = true; // For now, always pass
+
+ if (energyOK) {
+ // Matches OK
+ return nearbyCluster;
+ } else {
+ // This cluster isn't a good match -- ignore it.
+ }
+ }
+ }
+ // 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;
+ }
+
+ public void process(EventHeader event) {
+ initGeometry(event);
+ }
+
+ 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);
+ }
+ }
+
+ 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;
+
+ 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;
+ }
+}
lcsim/src/org/lcsim/recon/pfa/identifier
diff -N SimpleTrackMIPClusterMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SimpleTrackMIPClusterMatcher.java 10 Aug 2006 22:57:02 -0000 1.1
@@ -0,0 +1,59 @@
+package org.lcsim.recon.pfa.identifier;
+
+import java.util.*;
+import hep.physics.vec.*;
+
+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;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+
+/**
+ * Attempt to match a Track to a MIP-like Cluster, based on the intercept point
+ * on the ECAL inner surface and on the direction of the track at the
+ * intercept point.
+ *
+ * @version $Id: SimpleTrackMIPClusterMatcher.java,v 1.1 2006/08/10 22:57:02 mcharles Exp $
+ */
+
+public class SimpleTrackMIPClusterMatcher extends SimpleTrackClusterMatcher
+{
+ protected Cluster findMatchedCluster(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;
+ }
+}
