lcsim/src/org/lcsim/recon/pfa/identifier
diff -u -r1.1 -r1.2
--- SimpleTrackClusterMatcher.java 10 Aug 2006 22:57:02 -0000 1.1
+++ SimpleTrackClusterMatcher.java 3 Nov 2006 19:01:03 -0000 1.2
@@ -13,22 +13,56 @@
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
import org.lcsim.util.Driver;
+import org.lcsim.event.MCParticle;
+import hep.physics.particle.Particle;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.cluster.util.ClusterEnergyCalculator;
/**
* 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 $
+ * Currently, the match criteria are hard-code. It would be better to
+ * supply them as DecisionMaker objects.
+ *
+ * @version $Id: SimpleTrackClusterMatcher.java,v 1.2 2006/11/03 19:01:03 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)
{
+ if (m_debug) {
+ System.out.println("DEBUG: SimpleTrackClusterMatcher.matchTrackToCluster invoked for a list of "+clusters.size()+" clusters.");
+ double[] trackMom = tr.getMomentum();
+ double trackMomMag = Math.sqrt(trackMom[0]*trackMom[0] + trackMom[1]*trackMom[1] + trackMom[2]*trackMom[2]);
+ System.out.println("DEBUG: Track has momentum "+trackMomMag+" --> ("+trackMom[0]+", "+trackMom[1]+", "+trackMom[2]+")");
+ System.out.println("DEBUG: Track has reference point ("+tr.getReferencePointX()+", "+tr.getReferencePointY()+", "+tr.getReferencePointZ()+")");
+ System.out.println("DEBUG: Track parameters are:"
+ +" d0="+tr.getTrackParameter(0)
+ +", phi0="+tr.getTrackParameter(1)
+ +", omega="+tr.getTrackParameter(2)
+ +", z0="+tr.getTrackParameter(3)
+ +", s="+tr.getTrackParameter(4)
+ );
+ if (tr instanceof org.lcsim.mc.fast.tracking.ReconTrack) {
+ org.lcsim.mc.fast.tracking.ReconTrack rectr = (org.lcsim.mc.fast.tracking.ReconTrack) (tr);
+ System.out.println("DEBUG: Track parameters (unsmeared from ReconTrack) are:"
+ +" d0="+rectr.getNotSmearedTrack().getTrackParameter(0)
+ +", phi0="+rectr.getNotSmearedTrack().getTrackParameter(1)
+ +", omega="+rectr.getNotSmearedTrack().getTrackParameter(2)
+ +", z0="+rectr.getNotSmearedTrack().getTrackParameter(3)
+ +", s="+rectr.getNotSmearedTrack().getTrackParameter(4)
+ );
+ } else {
+ System.out.println("Could not unsmear track because it is a "+tr.getClass().getName());
+ }
+ }
+
// Make a HelixSwimmer to propagate the track
HelixSwimmer swimmer = new HelixSwimmer(m_fieldStrength[2]);
swimmer.setTrack(tr);
@@ -54,6 +88,19 @@
// Did we make a successful extrapolation?
if (Double.isNaN(alpha)) {
// No -- failed
+ if (m_debug) {
+ System.out.println("DEBUG: "+this.getClass().getName()+": Failed to extrapolate: alphaBarrel="+alphaBarrel+" ("+validBarrel+"), alphaEndcap="+alphaEndcap+" ("+validEndcap+") => alpha="+alpha);
+ double[] trackMom = tr.getMomentum();
+ double trackMomMag = Math.sqrt(trackMom[0]*trackMom[0] + trackMom[1]*trackMom[1] + trackMom[2]*trackMom[2]);
+ if (tr instanceof org.lcsim.mc.fast.tracking.ReconTrack) {
+ Particle debugTruthMatch = ((ReconTrack)(tr)).getMCParticle();
+ if (debugTruthMatch != null) {
+ System.out.println("DEBUG: "+this.getClass().getName()+": No track match for track with p="+trackMomMag+" from truth "+debugTruthMatch.getType().getName()+" with E="+debugTruthMatch.getEnergy()+" because extrapolation failed.");
+ } else {
+ System.out.println("DEBUG: "+this.getClass().getName()+": No track match for track with p="+trackMomMag+" [no truth match] because extrapolation failed.");
+ }
+ }
+ }
return null;
}
if ( !(validEndcap || validBarrel) ) {
@@ -72,18 +119,45 @@
// Matched OK
if (m_debug) {
Hep3Vector momentumVec = new BasicHep3Vector(tr.getMomentum());
- System.out.println("DEBUG: Extrapolated track to cluster (momentum = "+momentumVec.magnitude()+")");
+ System.out.println("DEBUG: "+this.getClass().getName()+": Extrapolated track to cluster (momentum = "+momentumVec.magnitude()+")");
+ double[] trackMom = tr.getMomentum();
+ double trackMomMag = Math.sqrt(trackMom[0]*trackMom[0] + trackMom[1]*trackMom[1] + trackMom[2]*trackMom[2]);
+ if (tr instanceof org.lcsim.mc.fast.tracking.ReconTrack) {
+ Particle debugTruthMatch = ((ReconTrack)(tr)).getMCParticle();
+ if (debugTruthMatch != null) {
+ System.out.println("DEBUG: "+this.getClass().getName()+": Track match for track with p="+trackMomMag+" from truth "+debugTruthMatch.getType().getName()+" with E="+debugTruthMatch.getEnergy());
+ } else {
+ System.out.println("DEBUG: "+this.getClass().getName()+": Track match for track with p="+trackMomMag+" [no truth match]");
+ }
+ }
}
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()+")");
+ System.out.println("DEBUG: "+this.getClass().getName()+": Failed to extrapolate track (momentum = "+momentumVec.magnitude()+")");
+ double[] trackMom = tr.getMomentum();
+ double trackMomMag = Math.sqrt(trackMom[0]*trackMom[0] + trackMom[1]*trackMom[1] + trackMom[2]*trackMom[2]);
+ if (tr instanceof org.lcsim.mc.fast.tracking.ReconTrack) {
+ Particle debugTruthMatch = ((ReconTrack)(tr)).getMCParticle();
+ if (debugTruthMatch != null) {
+ System.out.println("DEBUG: "+this.getClass().getName()+": No track match for track with p="+trackMomMag+" from truth "+debugTruthMatch.getType().getName()+" with E="+debugTruthMatch.getEnergy());
+ } else {
+ System.out.println("DEBUG: "+this.getClass().getName()+": No track match for track with p="+trackMomMag+" [no truth match]");
+ }
+ }
}
+
return null;
}
}
+
+ protected boolean m_checkEoverP = false;
+ public void setCheckEoverP(boolean checkEoverP) {
+ m_checkEoverP = checkEoverP;
+ if (m_debug) { System.out.println("DEBUG: set m_checkEoverP to "+m_checkEoverP); }
+ }
protected double swimToBarrel(HelixSwimmer swimmer) {
// Look for a hit in the first layer of the ECAL barrel
@@ -110,6 +184,8 @@
protected Cluster findMatchedCluster(Track tr, HelixSwimmer swimmer, double alpha, List<Cluster> clusters)
{
+ if (m_debug) { System.out.println("DEBUG: SimpleTrackClusterMatched.findMatchedCluster() invoked for a list of "+clusters.size()+" clusters."); }
+
// Find the track intercept and direction
swimmer.setTrack(tr);
Hep3Vector trackPoint = swimmer.getPointAtDistance(alpha);
@@ -119,29 +195,59 @@
// Obtain geometrical info:
CalorimeterHit nearestHit = findNearestHit(trackPoint, nearbyCluster);
double separation = proximity(trackPoint, nearestHit);
- int firstLayerHit = getLayer(nearestHit);
+ int layerOfNearestHit = getLayer(nearestHit);
+ CalorimeterHit firstHitInECAL = findInnermostHitInECAL(nearbyCluster);
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 goodFirstLayer = (firstHitInECAL!=null && getLayer(firstHitInECAL) < 5);
boolean goodSeparation = (separation < 30.0);
boolean foundMatch = goodSubDet && goodFirstLayer && goodSeparation;
+ if (m_debug) {
+ String printme = new String();
+ printme += "Debug: Didn't match track to cluster since";
+ printme += " subdet="+subdet.getName()+" ["+goodSubDet+"] and";
+ if (firstHitInECAL!=null) {
+ printme += " firstlayer="+getLayer(firstHitInECAL)+" ["+goodFirstLayer+"] and";
+ } else {
+ printme += " firstlayer=null ["+goodFirstLayer+"] and";
+ }
+ printme += " separation="+separation+" ["+goodSeparation+"]";
+ System.out.println(printme);
+ String debugContributions = new String();
+ debugContributions += "DEBUG: Cluster contents:";
+ Map<MCParticle, List<CalorimeterHit>> tmpMap = new HashMap<MCParticle, List<CalorimeterHit>>();
+ for (CalorimeterHit hit : nearbyCluster.getCalorimeterHits()) {
+ SimCalorimeterHit simhit = (SimCalorimeterHit) (hit);
+ for (int i=0; i<simhit.getMCParticleCount(); i++) {
+ MCParticle hitPart = simhit.getMCParticle(i);
+ if ( ! (tmpMap.keySet().contains(hitPart)) ) {
+ tmpMap.put(hitPart, new Vector<CalorimeterHit>());
+ }
+ tmpMap.get(hitPart).add(hit);
+ }
+ }
+ for (MCParticle hitPart : tmpMap.keySet()) {
+ debugContributions += " ";
+ debugContributions += hitPart.getType().getName();
+ debugContributions += " (E=";
+ debugContributions += hitPart.getEnergy();
+ debugContributions += ", hits=";
+ debugContributions += tmpMap.get(hitPart).size();
+ debugContributions += ")";
+ }
+ System.out.println(debugContributions);
+ }
if (foundMatch) {
// Geometrically, it looks good.
// Is it sensible in terms of energy?
+
+ boolean energyOK = true; // By default, always pass
+ if (m_checkEoverP) {
+ boolean passesEoverPcut = checkEoverP(nearbyCluster, tr);
+ energyOK = energyOK && passesEoverPcut;
+ }
- // // 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;
@@ -154,6 +260,49 @@
return null;
}
+ protected boolean checkEoverP(Cluster nearbyCluster, Track tr) {
+ // 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
+
+ // Check energy and uncertainty from calorimeter:
+ double estimatedClusterEnergy = estimateClusterEnergy(nearbyCluster);
+ double estimatedClusterEnergyUncertainty = 0.7 * Math.sqrt(estimatedClusterEnergy); // 70%/sqrt(E) for hadrons
+
+ // Check energy from track
+ double[] trackMomentum = tr.getMomentum();
+ double trackMomentumMagSq = trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2];
+ double trackMomentumMag = Math.sqrt(trackMomentumMagSq);
+ double massPion = 0.14;
+ double massProton = 0.94;
+ double energyReleaseIfPion = Math.sqrt(trackMomentumMagSq + massPion*massPion);
+ double energyReleaseIfProton = trackMomentumMag;
+ double energyReleaseIfAntiproton = trackMomentumMag + massProton + massProton;
+
+ double allowedVariation = 3.0; // 3sigma
+ boolean energyDiffOK_pion = Math.abs((energyReleaseIfPion - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < allowedVariation;
+ boolean energyDiffOK_proton = Math.abs((energyReleaseIfProton - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < allowedVariation;
+ boolean energyDiffOK_antiproton = Math.abs((energyReleaseIfAntiproton - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < allowedVariation;
+ boolean energyDiffOK = energyDiffOK_pion || energyDiffOK_proton || energyDiffOK_antiproton;
+
+ //boolean fractionEnergyOK = estimatedClusterEnergy > 0.5*trackMomentumMag; // don't use
+ //boolean absoluteEnergyOK = (trackMomentumMag - estimatedClusterEnergy < 0.2); // deliberately one-sided (old... is this right?)
+
+ if (m_debug) {
+ System.out.println("DEBUG [gen]: Checked E/P; found cluster E="+estimatedClusterEnergy
+ +" +- "+estimatedClusterEnergyUncertainty
+ +" and track E="+energyReleaseIfPion+" ("+energyReleaseIfProton+"/"+energyReleaseIfAntiproton+")"
+ +" => "+energyDiffOK);
+ System.out.println("DEBUG: Comparing track with momentum "+trackMomentumMag+" to a cluster with estimated energy "+estimatedClusterEnergy+" +- "+estimatedClusterEnergyUncertainty);
+ System.out.println(" If pi+: |"+energyReleaseIfPion+" - "+estimatedClusterEnergy+"|/"+estimatedClusterEnergyUncertainty+" = "+Math.abs((energyReleaseIfPion - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty)+" => "+energyDiffOK_pion);
+ System.out.println(" If p: |"+energyReleaseIfProton+" - "+estimatedClusterEnergy+"|/"+estimatedClusterEnergyUncertainty+" = "+Math.abs((energyReleaseIfProton - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty)+" => "+energyDiffOK_proton);
+ System.out.println(" If pbar: |"+energyReleaseIfAntiproton+" - "+estimatedClusterEnergy+"|/"+estimatedClusterEnergyUncertainty+" = "+Math.abs((energyReleaseIfAntiproton - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty)+" => "+energyDiffOK_antiproton);
+ System.out.println(" => Overall: "+energyDiffOK);
+ }
+
+ return energyDiffOK;
+ }
+
protected List<Cluster> findNearestClusters(Hep3Vector point, List<Cluster> clusterList)
{
Map<Cluster,Double> mapClusterToDistance = new HashMap<Cluster, Double>();
@@ -176,11 +325,29 @@
double distance = VecOp.sub(hitPosition, point).magnitude();
if (distance<minDist || nearest==null) {
nearest = hit;
+ minDist = distance;
}
}
return nearest;
}
+ protected CalorimeterHit findInnermostHitInECAL(Cluster clus) {
+ CalorimeterHit innermostHit = null;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ 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 double proximity(Hep3Vector point, Cluster clus) {
CalorimeterHit nearestHit = findNearestHit(point, clus);
return proximity(point, nearestHit);
@@ -212,10 +379,16 @@
return layer;
}
+ protected EventHeader m_event;
public void process(EventHeader event) {
+ m_event = event;
initGeometry(event);
}
+ public void setDebug(boolean debug) {
+ m_debug = debug;
+ }
+
public void initGeometry(EventHeader event)
{
Detector det = event.getDetector();
@@ -270,4 +443,11 @@
}
Map<T,Double> m_map;
}
+
+ protected ClusterEnergyCalculator m_calib = null;
+ protected double estimateClusterEnergy(Cluster clus) {
+ return m_calib.getEnergy(clus);
+ }
+ /** Specify what energy calibration to use for E/P check. */
+ public void setCalibration(ClusterEnergyCalculator calib) { m_calib = calib; }
}
