lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.22 -r1.23
--- ReclusterDriver.java 13 Jun 2008 22:46:38 -0000 1.22
+++ ReclusterDriver.java 17 Jun 2008 16:43:32 -0000 1.23
@@ -37,7 +37,7 @@
*
* This version is PRELIMINARY.
*
- * @version $Id: ReclusterDriver.java,v 1.22 2008/06/13 22:46:38 mcharles Exp $
+ * @version $Id: ReclusterDriver.java,v 1.23 2008/06/17 16:43:32 mcharles Exp $
* @author Mat Charles
*/
@@ -554,7 +554,7 @@
long timeAtEndOfThisTrackResidualCheck_EnergyCalculation = Calendar.getInstance().getTimeInMillis(); // DEBUG
if (m_debugTiming) { System.out.println("DEBUG: Computing residuals for this track: Energy calculation took "+(timeAtEndOfThisTrackResidualCheck_EnergyCalculation-timeAtStartOfThisTrackResidualCheck_EnergyCalculation)+" ms"); }
double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
- double sigma = 0.7 * Math.sqrt(trackMomentum);
+ double sigma = estimatedEnergyUncertainty(tr);
double normalizedResidual = (clusterEnergy-trackMomentum)/sigma;
if (normalizedResidual < -1.0) {
// FIXME: This cut-off shouldn't be hard-coded
@@ -744,11 +744,7 @@
outputChargedParticleList.add(part);
// Now consider E/p veto
double clusterEnergy = energy(showerComponents, allSharedClusters);
- double sigma = 0.7 * Math.sqrt(trackMomentumMag);
- if (trackMomentumMag < 1.0) {
- // Don't lowball uncertainty for low-momentum tracks
- sigma = 0.7;
- }
+ double sigma = estimatedEnergyUncertainty(tr);
double normalizedResidual = (clusterEnergy-trackMomentumMag)/sigma;
if (normalizedResidual > -windowEoverPveto && normalizedResidual < windowEoverPveto) {
// Passes E/p check
@@ -1203,13 +1199,23 @@
}
}
+ protected void printStatus(String desc, List<Track> tracksSortedByMomentum, List<SharedClusterGroup> allSharedClusters, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, Map<Cluster, Track> newMapShowerComponentToTrack, Map<Track,Double> newMapTrackToThreshold, Map<Track,Double> newMapTrackToTolerance, List<Cluster> photonClusters, List<Cluster> mips, List<Cluster> clumps, List<Cluster> largeClustersWithoutSkeletons, List<Cluster> smallClusterSeeds, Map<Track, Set<Cluster>> newMapTrackToVetoedAdditions) {
+ printStatus(desc, tracksSortedByMomentum, allSharedClusters, newMapTrackToShowerComponents, newMapShowerComponentToTrack, newMapTrackToThreshold, newMapTrackToTolerance, null, null, null, photonClusters, mips, clumps, largeClustersWithoutSkeletons, smallClusterSeeds, newMapTrackToVetoedAdditions);
+ }
+
protected void printStatus(String desc,
List<Track> tracksSortedByMomentum,
List<SharedClusterGroup> allSharedClusters,
+ // Track stuff
Map<Track, Set<Cluster>> newMapTrackToShowerComponents,
Map<Cluster, Track> newMapShowerComponentToTrack,
Map<Track,Double> newMapTrackToThreshold,
Map<Track,Double> newMapTrackToTolerance,
+ // Jet stuff
+ Map<Set<Track>, Set<Cluster>> newMapJetToShowerComponents,
+ Map<Cluster, Set<Track>> newMapShowerComponentToJet,
+ Map<Track, Set<Track>> mapTrackToJet,
+ // Cluster stuff
List<Cluster> photonClusters,
List<Cluster> mips,
List<Cluster> clumps,
@@ -1219,27 +1225,62 @@
)
{
System.out.println(desc);
+ Set<Set<Track>> jets = null;
+ if (newMapJetToShowerComponents != null) {
+ jets = newMapJetToShowerComponents.keySet();
+ } else {
+ jets = new HashSet<Set<Track>>(); // leave empty
+ }
double totalChiSquared = 0.0;
+ int totalNDF = 0;
for (Track tr : tracksSortedByMomentum) {
- Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
- double clusterEnergy = energy(showerComponents, allSharedClusters);
- double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
- double sigma = 0.7 * Math.sqrt(trackMomentum);
- double normalizedResidual = (clusterEnergy-trackMomentum)/sigma;
- List<Track> tmpTrackList = new Vector<Track>();
- tmpTrackList.add(tr);
- double realEffic = quoteEfficiency_T(tmpTrackList, showerComponents, allSharedClusters);
- double realPurity = quotePurity_T(tmpTrackList, showerComponents, allSharedClusters);
- double coreEffic = quoteEfficiency_T(tmpTrackList, showerComponents);
- double corePurity = quotePurity_T(tmpTrackList, showerComponents);
- System.out.println(" Track: threshold="+newMapTrackToThreshold.get(tr)
- +", tolerance="+newMapTrackToTolerance.get(tr)
- +", E/p is "+clusterEnergy+" / "+trackMomentum+" => NormResid = "+normalizedResidual
- +". Core: effic="+coreEffic+", purity="+corePurity
- +". Real: effic="+realEffic+", purity="+realPurity);
+ boolean fromJet = (mapTrackToJet != null && mapTrackToJet.get(tr) != null);
+ Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
+ double clusterEnergy = energy(showerComponents, allSharedClusters);
+ double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
+ double sigma = estimatedEnergyUncertainty(tr);
+ double normalizedResidual = (clusterEnergy-trackMomentum)/sigma;
+ List<Track> tmpTrackList = new Vector<Track>();
+ tmpTrackList.add(tr);
+ double realEffic = quoteEfficiency_T(tmpTrackList, showerComponents, allSharedClusters);
+ double realPurity = quotePurity_T(tmpTrackList, showerComponents, allSharedClusters);
+ double coreEffic = quoteEfficiency_T(tmpTrackList, showerComponents);
+ double corePurity = quotePurity_T(tmpTrackList, showerComponents);
+ String printme = new String();
+ if (fromJet) {
+ printme += "JET";
+ } else {
+ printme += " ";
double trackChiSquared = normalizedResidual * normalizedResidual;
totalChiSquared += trackChiSquared;
+ totalNDF++;
+ }
+ printme += "Track: threshold="+newMapTrackToThreshold.get(tr).floatValue()
+ +", tolerance="+newMapTrackToTolerance.get(tr).floatValue()
+ +", E/p is "+((float)(clusterEnergy))+" / "+((float)(trackMomentum))+" => NormResid = "+((float)(normalizedResidual))
+ +". Core: effic="+((float)(coreEffic))+", purity="+((float)(corePurity))
+ +". Real: effic="+((float)(realEffic))+", purity="+((float)(realPurity));
+ System.out.println(printme);
+ }
+ for (Set<Track> jet : jets) {
+ Set<Cluster> showerComponents = newMapJetToShowerComponents.get(jet);
+ double clusterEnergy = energy(showerComponents, allSharedClusters);
+ double jetMomentum = jetScalarMomentum(jet);
+ double sigma = estimatedEnergyUncertainty(jet);
+ double normalizedResidual = (clusterEnergy-jetMomentum)/sigma;
+ double realEffic = quoteEfficiency_T(jet, showerComponents, allSharedClusters);
+ double realPurity = quotePurity_T(jet, showerComponents, allSharedClusters);
+ double coreEffic = quoteEfficiency_T(jet, showerComponents);
+ double corePurity = quotePurity_T(jet, showerComponents);
+ System.out.println(" Jet: "
+ +" E/p is "+((float)(clusterEnergy))+" / "+((float)(jetMomentum))+" => NormResid = "+((float)(normalizedResidual))
+ +". Core: effic="+((float)(coreEffic))+", purity="+((float)(corePurity))
+ +". Real: effic="+((float)(realEffic))+", purity="+((float)(realPurity)));
+ double jetChiSquared = normalizedResidual * normalizedResidual;
+ totalChiSquared += jetChiSquared;
+ totalNDF++;
}
+
System.out.println("Total CHI^2 = "+totalChiSquared+" with NDF ~ "+tracksSortedByMomentum.size());
if (m_debug) {
System.out.println("Here is a summary of all the energy deposits for each track:");
@@ -1341,6 +1382,29 @@
}
}
+ protected double estimatedEnergyUncertainty(Track tr) {
+ double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
+ double sigma = 0.7 * Math.sqrt(trackMomentum);
+ if (trackMomentum < 1.0) {
+ // Don't lowball uncertainty for low-momentum tracks
+ sigma = 0.7;
+ }
+ return sigma;
+ }
+ protected double estimatedEnergyUncertainty(Collection<Track> tracks) {
+ double sumSquaredSigma = 0.0;
+ for (Track tr : tracks) {
+ double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
+ double sigmaSquared = 0.7 * 0.7 * trackMomentum;
+ if (trackMomentum < 1.0) {
+ // Don't lowball uncertainty for low-momentum tracks
+ sigmaSquared = 0.7 * 0.7;
+ }
+ sumSquaredSigma += sigmaSquared;
+ }
+ return Math.sqrt(sumSquaredSigma);
+ }
+
protected boolean testEoverP_oneSided(double clusterEnergy, Track tr, double tolerance) {
double[] trackThreeMomentum = tr.getMomentum();
double trackMomentumSq = trackThreeMomentum[0]*trackThreeMomentum[0] + trackThreeMomentum[1]*trackThreeMomentum[1] + trackThreeMomentum[2]*trackThreeMomentum[2];
@@ -1358,7 +1422,7 @@
double trackMomentumSq = trackMomentum * trackMomentum;
double trackEnergySq = trackMomentumSq + trackMassSq;
double trackEnergy = Math.sqrt(trackEnergySq);
- double trackEnergySigma = 0.7 * Math.sqrt(trackEnergy);
+ double trackEnergySigma = 0.7 * Math.sqrt(trackEnergy); // hmm...
double normalisedResidual = (clusterEnergy - trackEnergy)/trackEnergySigma;
return (normalisedResidual < tolerance); // not abs
}
@@ -1367,7 +1431,7 @@
double trackMomentumSq = trackMomentum * trackMomentum;
double trackEnergySq = trackMomentumSq + trackMassSq;
double trackEnergy = Math.sqrt(trackEnergySq);
- double trackEnergySigma = 0.7 * Math.sqrt(trackEnergy);
+ double trackEnergySigma = 0.7 * Math.sqrt(trackEnergy); // hmm...
double normalisedResidual = (clusterEnergy - trackEnergy)/trackEnergySigma;
return (Math.abs(normalisedResidual) < tolerance); // abs
}
@@ -1475,7 +1539,7 @@
trackList.add(tr);
return quoteEfficiency_T(trackList, clusterList);
}
- protected double quoteEfficiency_T(List<Track> trackList, Collection<Cluster> clusterList) {
+ protected double quoteEfficiency_T(Collection<Track> trackList, Collection<Cluster> clusterList) {
Cluster combined = makeCombinedCluster(clusterList);
return quoteEfficiency_T(trackList, combined);
}
@@ -1484,7 +1548,7 @@
trackList.add(tr);
return quotePurity_T(trackList, clusterList);
}
- protected double quotePurity_T(List<Track> trackList, Collection<Cluster> clusterList) {
+ protected double quotePurity_T(Collection<Track> trackList, Collection<Cluster> clusterList) {
Cluster combined = makeCombinedCluster(clusterList);
return quotePurity_T(trackList, combined);
}
@@ -1494,7 +1558,7 @@
trackList.add(tr);
return quoteEfficiency_T(trackList, clusterList, allSharedClusters);
}
- protected double quoteEfficiency_T(List<Track> trackList, Collection<Cluster> clusterList, List<SharedClusterGroup> allSharedClusters) {
+ protected double quoteEfficiency_T(Collection<Track> trackList, Collection<Cluster> clusterList, List<SharedClusterGroup> allSharedClusters) {
Cluster combined = makeCombinedCluster(clusterList);
return quoteEfficiency_T(trackList, combined, allSharedClusters);
}
@@ -1503,7 +1567,7 @@
trackList.add(tr);
return quotePurity_T(trackList, clusterList, allSharedClusters);
}
- protected double quotePurity_T(List<Track> trackList, Collection<Cluster> clusterList, List<SharedClusterGroup> allSharedClusters) {
+ protected double quotePurity_T(Collection<Track> trackList, Collection<Cluster> clusterList, List<SharedClusterGroup> allSharedClusters) {
Cluster combined = makeCombinedCluster(clusterList);
return quotePurity_T(trackList, combined, allSharedClusters);
}
@@ -1513,7 +1577,7 @@
partList.add(part);
return quotePurity_P(partList, clus, allSharedClusters);
}
- protected double quotePurity_T(List<Track> trackList, Cluster clus, List<SharedClusterGroup> allSharedClusters) {
+ protected double quotePurity_T(Collection<Track> trackList, Cluster clus, List<SharedClusterGroup> allSharedClusters) {
Vector<MCParticle> partList = new Vector<MCParticle>();
for (Track tr : trackList) {
List<MCParticle> particlesOfThisTrack = getMCParticle(tr);
@@ -1571,7 +1635,7 @@
partList.add(part);
return quoteEfficiency_P(partList, clus, allSharedClusters);
}
- protected double quoteEfficiency_T(List<Track> trackList, Cluster clus, List<SharedClusterGroup> allSharedClusters) {
+ protected double quoteEfficiency_T(Collection<Track> trackList, Cluster clus, List<SharedClusterGroup> allSharedClusters) {
Set<MCParticle> partSet = new HashSet<MCParticle>();
for (Track tr : trackList) {
List<MCParticle> particlesOfThisTrack = getMCParticle(tr);
@@ -1620,7 +1684,7 @@
double num = truthHitsInClusterCore.size() + countWeightedSharedTruthHits;
return num/denom;
}
- protected double quoteEfficiency_T(List<Track> trackList, Cluster clus) {
+ protected double quoteEfficiency_T(Collection<Track> trackList, Cluster clus) {
HitMap hitsEcal = ((HitMap)(m_event.get("inputHitMapEcal")));
HitMap hitsHcal = ((HitMap)(m_event.get("inputHitMapHcal")));
List<CalorimeterHit> allHits = new Vector<CalorimeterHit>();
@@ -1632,7 +1696,7 @@
double num = truthHitsInClusterCore.size();
return num/denom;
}
- protected double quotePurity_T(List<Track> trackList, Cluster clus) {
+ protected double quotePurity_T(Collection<Track> trackList, Cluster clus) {
Set<CalorimeterHit> truthHitsInCluster = findHitsFromTruth_T(trackList, clus.getCalorimeterHits());
double num = truthHitsInCluster.size();
double denom = clus.getCalorimeterHits().size();
@@ -1675,7 +1739,7 @@
}
return outputHits;
}
- protected Set<CalorimeterHit> findHitsFromTruth_T(List<Track> trackList, Collection<CalorimeterHit> hits) {
+ protected Set<CalorimeterHit> findHitsFromTruth_T(Collection<Track> trackList, Collection<CalorimeterHit> hits) {
Set<CalorimeterHit> outputHits = new HashSet<CalorimeterHit>();
for (Track tr : trackList) {
List<CalorimeterHit> foundHits = findHitsFromTruth_T(tr, hits);
@@ -3064,6 +3128,24 @@
}
return truthMap;
}
+ protected Map<MCParticle, List<SimCalorimeterHit>> truthInListFromHitList(List<CalorimeterHit> hits, List<MCParticle> mcList) {
+ Map<MCParticle, List<SimCalorimeterHit>> truthMap = truthFromHitList (hits);
+ Map<MCParticle, List<SimCalorimeterHit>> outputMap = new HashMap<MCParticle, List<SimCalorimeterHit>>();
+ for (MCParticle part : truthMap.keySet()) {
+ List<SimCalorimeterHit> truthHits = truthMap.get(part);
+ List<MCParticle> parents = findParentsInList(part, mcList);
+ for (MCParticle parent : parents) {
+ List<SimCalorimeterHit> parentHits = outputMap.get(parent);
+ if (parentHits == null) {
+ parentHits = new Vector<SimCalorimeterHit>();
+ outputMap.put(parent, parentHits);
+ }
+ parentHits.addAll(truthHits);
+ }
+ }
+ return outputMap;
+ }
+
protected int getLayer(CalorimeterHit hit) {
org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
@@ -3071,4 +3153,42 @@
int layer = id.getLayer();
return layer;
}
+
+ protected Hep3Vector jetThreeMomentum(Set<Track> jet) {
+ Hep3Vector totalJetMomentum = new BasicHep3Vector(0,0,0);
+ for (Track tr : jet) {
+ Hep3Vector trackMomentum = momentum(tr);
+ totalJetMomentum = VecOp.add(totalJetMomentum, trackMomentum);
+ }
+ return totalJetMomentum;
+ }
+ protected double jetScalarMomentum(Set<Track> jet) {
+ double totalJetEnergy = 0.0;
+ for (Track tr : jet) {
+ Hep3Vector trackMomentum = momentum(tr);
+ totalJetEnergy += trackMomentum.magnitude();
+ }
+ return totalJetEnergy;
+ }
+
+ protected Hep3Vector momentum(Track tr) {
+ Hep3Vector totalMomentum = null;
+ if (tr instanceof BaseTrackMC) {
+ totalMomentum = ((BaseTrackMC)(tr)).getMCParticle().getMomentum();
+ } else if (tr instanceof MultipleTrackTrack) {
+ totalMomentum = new BasicHep3Vector(0,0,0);
+ for (Track subtr : tr.getTracks()) {
+ Hep3Vector subTrackMomentum = null;
+ if (subtr instanceof BaseTrackMC) {
+ subTrackMomentum = ((BaseTrackMC)(subtr)).getMCParticle().getMomentum();
+ } else {
+ subTrackMomentum = new BasicHep3Vector(tr.getMomentum());
+ }
+ totalMomentum = VecOp.add(totalMomentum, subTrackMomentum);
+ }
+ } else {
+ totalMomentum = new BasicHep3Vector(tr.getMomentum());
+ }
+ return totalMomentum;
+ }
}