lcsim/src/org/lcsim/recon/pfa/identifier
diff -u -r1.4 -r1.5
--- SimpleChargedParticleMaker.java 18 Jul 2006 00:44:02 -0000 1.4
+++ SimpleChargedParticleMaker.java 10 Aug 2006 22:57:34 -0000 1.5
@@ -24,123 +24,42 @@
*
* Currently, PID is done by cheating.
*
- * @version $Id: SimpleChargedParticleMaker.java,v 1.4 2006/07/18 00:44:02 mcharles Exp $
+ * @version $Id: SimpleChargedParticleMaker.java,v 1.5 2006/08/10 22:57:34 mcharles Exp $
*/
public class SimpleChargedParticleMaker extends Driver
{
/** Simple constructor. */
- public SimpleChargedParticleMaker() {}
+ public SimpleChargedParticleMaker() {
+ }
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 setDebug(boolean debug) { m_debug = debug; }
+ public void setTrackMatcher(TrackClusterMatcher matcher) { m_matcher = matcher; }
public void process(EventHeader event)
{
- initGeometry(event);
+ super.process(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:
Map<Track,Cluster> matchedTracks = new HashMap<Track,Cluster>();
- List<Track> outputTrackList = new Vector<Track>(); // initially empty
+ List<Track> unmatchedTracks = new Vector<Track>(inputTrackList);
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>();
- List<Track> unsuccessfullyExtrapolatedTracks = new Vector<Track>();
-
+ // Try to match each track to a cluster
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");
- }
- } else {
- unsuccessfullyExtrapolatedTracks.add(tr);
- }
- }
-
- // 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;
- if (m_debug) {
- Hep3Vector momentumVec = new BasicHep3Vector(tr.getMomentum());
- System.out.println("DEBUG: Extrapolated track to MIP (momentum = "+momentumVec.magnitude()+")");
- }
- break;
- }
- }
- } else {
- // If that didn't work, try any cluster...
- matchedCluster = findMatchedCluster(tr, swimmer, alpha.doubleValue(), inputClusterList);
- }
-
+ Cluster matchedCluster = m_matcher.matchTrackToCluster(tr, inputClusterList);
if (matchedCluster != null) {
- // We found a match
matchedTracks.put(tr, matchedCluster);
- if (m_debug) {
- Hep3Vector momentumVec = new BasicHep3Vector(tr.getMomentum());
- System.out.println("DEBUG: Extrapolated track to cluster (momentum = "+momentumVec.magnitude()+")");
- }
- } else {
- // No match -- put the track into the list of unused tracks
- outputTrackList.add(tr);
- if (m_debug) {
- Hep3Vector momentumVec = new BasicHep3Vector(tr.getMomentum());
- System.out.println("DEBUG: Failed to extrapolate track (momentum = "+momentumVec.magnitude()+")");
- }
+ unmatchedTracks.remove(tr);
+ // Should we check that the matched cluster is in the
+ // expected list?
}
}
@@ -154,277 +73,26 @@
// Already used this cluster in a particle => just add the track to that particle
LocalReconstructedParticle part = matchedClusters.get(clus);
part.addTrack(tr);
- double[] trackMomentum = tr.getMomentum();
- double trackMomentumMagSq = (trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
- double trackMass = 0.140; // Treat the second particle as a pion
- if (tr instanceof ReconTrack) {
- Particle truthParticle = ((ReconTrack)(tr)).getMCParticle();
- if (truthParticle != null) {
- trackMass = truthParticle.getMass();
- }
- }
- double trackEnergy = Math.sqrt(trackMomentumMagSq + trackMass*trackMass);
- part.setEnergy(part.getEnergy() + trackEnergy);
+ part.recomputeKinematics();
} else {
// This cluster hasn't been used yet -- initialize its particle
LocalReconstructedParticle part = new LocalReconstructedParticle();
part.addTrack(tr);
part.addCluster(clus);
- double[] trackMomentum = tr.getMomentum();
- double trackMomentumMagSq = (trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
- double mass = 0.140;
- if (tr instanceof ReconTrack) {
- Particle truthParticle = ((ReconTrack)(tr)).getMCParticle();
- if (truthParticle != null) {
- mass = truthParticle.getMass();
- }
- }
- double energy = Math.sqrt(trackMomentumMagSq + mass*mass);
- part.setEnergy(energy);
+ part.recomputeKinematics();
matchedClusters.put(clus, part);
}
- }
-
+ }
outputParticleList.addAll(matchedClusters.values());
// Write out
- event.put(m_outputTrackListName, outputTrackList);
+ event.put(m_outputTrackListName, unmatchedTracks);
event.put(m_outputParticleListName, outputParticleList);
-
- // Here is a hook for debug/plots
- checkEvent(successfullyExtrapolatedTracks, unsuccessfullyExtrapolatedTracks, outputTrackList, outputParticleList, inputClusterList, inputMIPList, swimmer);
}
- protected void checkEvent(Map<Track,Double> successfullyExtrapolatedTracks, List<Track> unsuccessfullyExtrapolatedTracks, List<Track> outputTrackList, List<ReconstructedParticle> outputParticleList, List<Cluster> inputClusterList, List<Cluster> inputMIPList, HelixSwimmer swimmer)
- {
- // nothing by default
- }
-
- 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) {
- // 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;
- }
-
- 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;
- protected boolean m_debug = false;
-
+ protected TrackClusterMatcher m_matcher;
String m_inputTrackListName;
String m_outputTrackListName;
- String m_inputMIPListName;
String m_inputClusterListName;
String m_outputParticleListName;
@@ -443,6 +111,23 @@
public void setCharge(double q) {
_charge = q;
}
+ public void recomputeKinematics() {
+ double energy = 0.0;
+ for (Track tr : this.getTracks()) {
+ double[] trackMomentum = tr.getMomentum();
+ double trackMomentumMagSq = (trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2]);
+ double mass = 0.140;
+ if (tr instanceof ReconTrack) {
+ Particle truthParticle = ((ReconTrack)(tr)).getMCParticle();
+ if (truthParticle != null) {
+ mass = truthParticle.getMass();
+ }
+ }
+ double trackEnergy = Math.sqrt(trackMomentumMagSq + mass*mass);
+ energy += trackEnergy;
+ }
+ this.setEnergy(energy);
+ }
}
}