lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.7 -r1.8
--- ReclusterDTreeDriver.java 4 Jan 2008 19:25:25 -0000 1.7
+++ ReclusterDTreeDriver.java 9 Jan 2008 18:34:21 -0000 1.8
@@ -23,6 +23,8 @@
public class ReclusterDTreeDriver extends ReclusterDriver {
+ protected boolean m_cheatOnPhotonsMisidentifiedAsNeutralHadrons = false;
+
protected String m_dTreeClusterListName;
public ReclusterDTreeDriver(String dTreeClusterList, String trackList, String mcList) {
@@ -621,35 +623,7 @@
}
}
- // Photons (...)
- for (Cluster clus : photons) {
- Track tr = newMapShowerComponentToTrack.get(clus);
- if (tr == null) {
- // write out photon
- BaseReconstructedParticle part = new BaseReconstructedParticle();
- part.addCluster(clus);
- // Add fuzzy hits
- Cluster sharedHitClus = makeClusterOfSharedHits(clus, allSharedClusters);
- if (sharedHitClus.getCalorimeterHits().size()>0) { part.addCluster(sharedHitClus); }
- double clusterEnergy = energy(clus, allSharedClusters, m_photonCalib);
- Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
- Hep3Vector threeMomentum = VecOp.mult(clusterEnergy, VecOp.unit(pos)); // assume it comes from the IP
- HepLorentzVector fourMomentum = new BasicHepLorentzVector(clusterEnergy, threeMomentum);
- part.set4Vector(fourMomentum);
- part.setReferencePoint(0,0,0);
- part.setCharge(0);
- // Set the PID and mass
- part.addParticleID(pid_photon);
- part.setParticleIdUsed(pid_photon);
- part.setMass(0.0);
- // Add to the output list
- outputParticleList.add(part);
- outputPhotonParticleList.add(part);
- }
- }
-
- // NHad
-
+ // Find neutral hadron pieces
List<Cluster> unmatchedClusterPieces = new Vector<Cluster>();
for (Cluster clus : linkableClusters) {
Track tr = newMapShowerComponentToTrack.get(clus);
@@ -699,10 +673,64 @@
neutralClusterCores.add(newShower);
}
}
-
+
+ // Now, optionally, we can cheat for a specific case.
+ // Sometimes photons (or photon fragments) are misidentified as neutral
+ // hadrons. In those cases, the energy is overestimated. Optionally,
+ // we can go back and change the ID to photon based on truth info.
+ Set<Cluster> extraClustersToTreatAsPhotons = new HashSet<Cluster>();
+ if (m_cheatOnPhotonsMisidentifiedAsNeutralHadrons) {
+ System.out.println("WARNING: Cheating on neutral hadron -> photon ID");
+ for (Cluster clus : neutralClusterCores) {
+ // Require exactly one subcluster (since real photons don't have structure)
+ if (clus.getClusters().size() == 1) {
+ // It's a candidate
+ MCParticle domPartOfClus = quoteDominantParticle(clus);
+ int domPDG = domPartOfClus.getPDGID();
+ if (domPDG == 22) {
+ extraClustersToTreatAsPhotons.add(clus);
+ }
+ }
+ }
+ }
+
+ // Photons (...)
+ List<Cluster> photonsToUse = new Vector<Cluster>(photons);
+ photonsToUse.addAll(extraClustersToTreatAsPhotons);
+ for (Cluster clus : photonsToUse) {
+ Track tr = newMapShowerComponentToTrack.get(clus);
+ if (tr == null) {
+ // write out photon
+ BaseReconstructedParticle part = new BaseReconstructedParticle();
+ part.addCluster(clus);
+ // Add fuzzy hits
+ Cluster sharedHitClus = makeClusterOfSharedHits(clus, allSharedClusters);
+ if (sharedHitClus.getCalorimeterHits().size()>0) { part.addCluster(sharedHitClus); }
+ double clusterEnergy = energy(clus, allSharedClusters, m_photonCalib);
+ Hep3Vector pos = new BasicHep3Vector(clus.getPosition());
+ Hep3Vector threeMomentum = VecOp.mult(clusterEnergy, VecOp.unit(pos)); // assume it comes from the IP
+ HepLorentzVector fourMomentum = new BasicHepLorentzVector(clusterEnergy, threeMomentum);
+ part.set4Vector(fourMomentum);
+ part.setReferencePoint(0,0,0);
+ part.setCharge(0);
+ // Set the PID and mass
+ part.addParticleID(pid_photon);
+ part.setParticleIdUsed(pid_photon);
+ part.setMass(0.0);
+ // Add to the output list
+ outputParticleList.add(part);
+ outputPhotonParticleList.add(part);
+ if (m_debug) { debugPrintPhotonParticleInfo(clus, allSharedClusters); }
+ }
+ }
+
// Make reconstructed particles from the neutral hadron showers:
for (Cluster clus : neutralClusterCores) {
+ if (extraClustersToTreatAsPhotons.contains(clus)) {
+ // Skip this (already used as a photon)
+ continue;
+ }
// write out neutral hadron
BaseReconstructedParticle part = new BaseReconstructedParticle();
part.addCluster(clus);
@@ -725,6 +753,7 @@
// Add to the output list
outputParticleList.add(part);
outputNeutralHadronParticleList.add(part);
+ if (m_debug) { debugPrintNeutralHadronParticleInfo(clus, allSharedClusters); }
}
List<ReconstructedParticle> outputParticleListWithEoverPveto = new Vector<ReconstructedParticle>();
@@ -732,6 +761,10 @@
outputParticleListWithEoverPveto.addAll(outputNeutralHadronParticleList);
outputParticleListWithEoverPveto.addAll(outputChargedParticleListWithEoverPveto);
+ if (m_debug) {
+ debugPrintHitStatus(allSharedClusters, tracksSortedByMomentum, newMapTrackToShowerComponents, newMapShowerComponentToTrack, photons, neutralClusterCores);
+ }
+
// Write out
m_event.put(m_outputParticleListName, outputParticleList);
m_event.put(m_outputParticleListName+"_withEoverPveto", outputParticleListWithEoverPveto);
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.14 -r1.15
--- ReclusterDriver.java 4 Jan 2008 19:20:27 -0000 1.14
+++ ReclusterDriver.java 9 Jan 2008 18:34:21 -0000 1.15
@@ -32,7 +32,7 @@
*
* This version is PRELIMINARY.
*
- * @version $Id: ReclusterDriver.java,v 1.14 2008/01/04 19:20:27 mcharles Exp $
+ * @version $Id: ReclusterDriver.java,v 1.15 2008/01/09 18:34:21 mcharles Exp $
* @author Mat Charles
*/
@@ -780,16 +780,7 @@
outputParticleList.add(part);
outputPhotonParticleList.add(part);
if (m_debug) {
- // Debug printout
- MCParticle dominantParticle = quoteDominantParticle(clus);
- double truthMomentum = dominantParticle.getMomentum().magnitude();
- int truthPDG = dominantParticle.getPDGID();
- double effic = quoteEfficiency_P(dominantParticle, clus, allSharedClusters);
- double purity = quotePurity_P(dominantParticle, clus, allSharedClusters);
- double energyNoShare = m_photonCalib.getEnergy(clus);
- double efficNoShare = quoteEfficiency_P(dominantParticle, clus);
- double purityNoShare = quotePurity_P(dominantParticle, clus);
- System.out.println("DEBUG: Looking at photon cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity+". Excluding shared hits, E="+energyNoShare+", effic="+efficNoShare+", purity="+purityNoShare);
+ debugPrintPhotonParticleInfo(clus, allSharedClusters);
}
}
}
@@ -880,14 +871,7 @@
outputParticleList.add(part);
outputNeutralHadronParticleList.add(part);
if (m_debug) {
- // Debug printout
- MCParticle dominantParticle = quoteDominantParticle(clus);
- double truthEnergy = dominantParticle.getEnergy();
- double truthMomentum = dominantParticle.getMomentum().magnitude();
- int truthPDG = dominantParticle.getPDGID();
- double purity = quotePurity_P(dominantParticle, clus, allSharedClusters);
- double effic = quoteEfficiency_P(dominantParticle, clus, allSharedClusters);
- System.out.println("DEBUG: Looking at neutral cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity);
+ debugPrintNeutralHadronParticleInfo(clus, allSharedClusters);
}
}
@@ -1223,6 +1207,8 @@
}
}
// Investigate shared clusters
+ Set<CalorimeterHit> allTruthHitsInEcalShared = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> allTruthHitsInHcalShared = new HashSet<CalorimeterHit>();
for (SharedClusterGroup shares : allSharedClusters) {
List<SharedCluster> sharedClusters = shares.listAllSharedClusters();
for (SharedCluster sharedCluster : sharedClusters) {
@@ -1232,8 +1218,10 @@
for (CalorimeterHit hit : sharedCluster.getCluster().getCalorimeterHits()) {
if (allTruthHitsInEcal.contains(hit)) {
countTruthHitsEcalInThisSharedCluster++;
+ allTruthHitsInEcalShared.add(hit);
} else if (allTruthHitsInHcal.contains(hit)) {
countTruthHitsHcalInThisSharedCluster++;
+ allTruthHitsInHcalShared.add(hit);
}
}
if (countTruthHitsEcalInThisSharedCluster + countTruthHitsHcalInThisSharedCluster > 0) {
@@ -1259,32 +1247,10 @@
}
}
}
- System.out.println(" Track with p="+trackMomentum+" has "+allTruthHitsInEcal.size()+" ECAL + "+allTruthHitsInHcal.size()+" HCAL hits total. Of these, core has "+countCoreHitsEcal+" ECAL + "+countCoreHitsHcal+" HCAL hits and halo has "+countSharedHitsEcal+" ECAL + "+countSharedHitsHcal+" HCAL (weighted to "+countWeightedSharedHitsEcal+" + "+countWeightedSharedHitsHcal+"). Unmatched shared hits: "+countUnmatchedSharedHitsEcal+" ECAL + "+countUnmatchedSharedHitsHcal+" HCAL.");
+ System.out.println(" Track with p="+trackMomentum+" has "+allTruthHitsInEcal.size()+" ECAL + "+allTruthHitsInHcal.size()+" HCAL hits total. Of these, core has "+countCoreHitsEcal+" ECAL + "+countCoreHitsHcal+" HCAL hits and halo has "+countSharedHitsEcal+" ECAL + "+countSharedHitsHcal+" HCAL (weighted to "+countWeightedSharedHitsEcal+" + "+countWeightedSharedHitsHcal+"). Unmatched shared hits: "+countUnmatchedSharedHitsEcal+" ECAL + "+countUnmatchedSharedHitsHcal+" HCAL. Total shared hit count: "+allTruthHitsInEcalShared.size()+" ECAL + "+allTruthHitsInHcalShared.size()+" HCAL");
// MORE DEBUG
- System.out.println("Core clusters:");
- for (Cluster clus : subClusters) {
- MCParticle domPart = quoteDominantParticle(clus);
- System.out.println(" -> "+clus.getCalorimeterHits().size()+" (dominant: "+domPart.getPDGID()+" with p="+domPart.getMomentum().magnitude());
- }
- // MORE DEBUG
- System.out.println("Core contributions:");
- Map<MCParticle,Set<CalorimeterHit>> mapTruthToHitSets = new HashMap<MCParticle,Set<CalorimeterHit>>();
- for (Cluster clus : subClusters) {
- for (CalorimeterHit hit : clus.getCalorimeterHits()) {
- SimCalorimeterHit simhit = (SimCalorimeterHit)(hit);
- int nMCParticles = simhit.getMCParticleCount();
- for (int iMC=0; iMC<nMCParticles; iMC++) {
- MCParticle part = simhit.getMCParticle(iMC);
- Set<CalorimeterHit> hitSet = mapTruthToHitSets.get(part);
- if (hitSet == null) { hitSet = new HashSet<CalorimeterHit>(); mapTruthToHitSets.put(part,hitSet); }
- hitSet.add(hit);
- }
- }
- }
- for (MCParticle part : mapTruthToHitSets.keySet()) {
- Set<CalorimeterHit> hitSet = mapTruthToHitSets.get(part);
- System.out.println(" -> "+part.getPDGID()+" with p="+part.getMomentum().magnitude()+" has "+hitSet.size()+" core hits.");
- }
+ debugPrintCoreClusterInfo(subClusters);
+ debugPrintCoreContributions(subClusters);
}
for (Cluster clus : photonClusters) {
debugPrint("Photon", clus,tracksSortedByMomentum,newMapTrackToShowerComponents,newMapShowerComponentToTrack,allSharedClusters,newMapTrackToVetoedAdditions);
@@ -1655,6 +1621,14 @@
List<CalorimeterHit> hitsOfThisParticle = findHitsFromTruth_P(part, hits);
outputHits.addAll(hitsOfThisParticle);
}
+ // Verify no duplication
+ Set<Long> uniqueIDs = new HashSet<Long>();
+ for (CalorimeterHit hit : outputHits) {
+ uniqueIDs.add(hit.getCellID());
+ }
+ if (uniqueIDs.size() != outputHits.size()) {
+ throw new AssertionError("When looking for truth hits of track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+", I found "+outputHits.size()+" hits in a collection of "+hits.size()+", but only "+uniqueIDs.size()+" were unique.");
+ }
return outputHits;
}
protected List<CalorimeterHit> findHitsFromTruth_P(MCParticle part, Collection<CalorimeterHit> hitsToSearch) {
@@ -1738,7 +1712,11 @@
continue;
} else {
// Ambiguity
- throw new AssertionError("Ambiguity in subcluster match");
+ if (m_debug) {
+ throw new AssertionError("Ambiguity in subcluster match");
+ } else {
+ System.out.println("ERROR: Ambiguity in subcluster match");
+ }
}
}
return foundSubClusters;
@@ -1846,7 +1824,6 @@
if (m_mips != null) {
if (clus.getCalorimeterHits().size() < 6 && m_mips.contains(clus)) {
// VETO
- if (m_debug) { System.out.println("DEBUG: Skipping daughters of cluster for implied-cluster-addition because it is a mip with only "+clus.getCalorimeterHits().size()+" hits."); }
continue;
}
}
@@ -1862,7 +1839,7 @@
if (!allShowerComponents.contains(potentialNewCluster) && !output.contains(potentialNewCluster) && !clustersToAddThisPass.contains(potentialNewCluster)) {
// We haven't used this one yet -- should add it if it's available
if (unavailableClusters.contains(potentialNewCluster)) {
- if (m_debug) { System.out.println("DEBUG: Skipping cluster because it's already been used"); }
+ //if (m_debug) { System.out.println("DEBUG: Skipping cluster because it's already been used"); }
} else {
clustersToAddThisPass.add(potentialNewCluster);
}
@@ -1990,7 +1967,7 @@
return m_sharedClusters;
}
}
- private class SharedCluster {
+ protected class SharedCluster {
Cluster m_rawCluster = null;
Map<Cluster, Double> m_sharedBetween = null;
double m_sumOfWeights = 0.0;
@@ -2458,6 +2435,68 @@
super.suspend();
}
+ protected void debugPrintPhotonParticleInfo(Cluster clus, List<SharedClusterGroup> allSharedClusters) {
+ double clusterEnergy = energy(clus, allSharedClusters, m_photonCalib);
+ MCParticle dominantParticle = quoteDominantParticle(clus);
+ double truthMomentum = dominantParticle.getMomentum().magnitude();
+ int truthPDG = dominantParticle.getPDGID();
+ double effic = quoteEfficiency_P(dominantParticle, clus, allSharedClusters);
+ double purity = quotePurity_P(dominantParticle, clus, allSharedClusters);
+ double energyNoShare = m_photonCalib.getEnergy(clus);
+ double efficNoShare = quoteEfficiency_P(dominantParticle, clus);
+ double purityNoShare = quotePurity_P(dominantParticle, clus);
+ System.out.println("DEBUG: Looking at photon cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity+". Excluding shared hits, E="+energyNoShare+", effic="+efficNoShare+", purity="+purityNoShare);
+ // More debug info:
+ Set<Cluster> subClusters = new HashSet<Cluster>();
+ subClusters.add(clus);
+ debugPrintCoreClusterInfo(subClusters);
+ debugPrintCoreContributions(subClusters);
+ }
+
+ protected void debugPrintCoreClusterInfo(Set<Cluster> subClusters) {
+ System.out.println("Core clusters:");
+ for (Cluster clus : subClusters) {
+ MCParticle domPart = quoteDominantParticle(clus);
+ System.out.println(" -> "+clus.getCalorimeterHits().size()+" (dominant: "+domPart.getPDGID()+" with p="+domPart.getMomentum().magnitude());
+ }
+ }
+ protected void debugPrintCoreContributions(Set<Cluster> subClusters) {
+ System.out.println("Core contributions:");
+ Map<MCParticle,Set<CalorimeterHit>> mapTruthToHitSets = new HashMap<MCParticle,Set<CalorimeterHit>>();
+ for (Cluster clus : subClusters) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ SimCalorimeterHit simhit = (SimCalorimeterHit)(hit);
+ int nMCParticles = simhit.getMCParticleCount();
+ for (int iMC=0; iMC<nMCParticles; iMC++) {
+ MCParticle part = simhit.getMCParticle(iMC);
+ Set<CalorimeterHit> hitSet = mapTruthToHitSets.get(part);
+ if (hitSet == null) { hitSet = new HashSet<CalorimeterHit>(); mapTruthToHitSets.put(part,hitSet); }
+ hitSet.add(hit);
+ }
+ }
+ }
+ for (MCParticle part : mapTruthToHitSets.keySet()) {
+ Set<CalorimeterHit> hitSet = mapTruthToHitSets.get(part);
+ System.out.println(" -> "+part.getPDGID()+" with p="+part.getMomentum().magnitude()+" has "+hitSet.size()+" core hits.");
+ }
+ }
+
+ protected void debugPrintNeutralHadronParticleInfo(Cluster clus, List<SharedClusterGroup> allSharedClusters) {
+ double clusterEnergy = energy(clus, allSharedClusters, m_neutralCalib);
+ MCParticle dominantParticle = quoteDominantParticle(clus);
+ double truthEnergy = dominantParticle.getEnergy();
+ double truthMomentum = dominantParticle.getMomentum().magnitude();
+ int truthPDG = dominantParticle.getPDGID();
+ double purity = quotePurity_P(dominantParticle, clus, allSharedClusters);
+ double effic = quoteEfficiency_P(dominantParticle, clus, allSharedClusters);
+ System.out.println("DEBUG: Looking at neutral cluster with measured energy = "+clusterEnergy+". Dominant particle is "+truthPDG+" with p="+truthMomentum+". Cluster has effic="+effic+", purity="+purity);
+ // More debug info:
+ Set<Cluster> subClusters = recursivelyFindSubClusters(clus);
+ subClusters.remove(clus); // Ignore the wrapper cluster which isn't really a piece of the core -- it's the whole thing
+ debugPrintCoreClusterInfo(subClusters);
+ debugPrintCoreContributions(subClusters);
+ }
+
private void debugPrint(String typeName,
Cluster clus,
List<Track> tracksSortedByMomentum,
@@ -2747,5 +2786,134 @@
System.out.println(printme);
}
}
+
+ protected void debugPrintHitStatus(List<SharedClusterGroup> allSharedClusters, List<Track> tracksSortedByMomentum, Map<Track, Set<Cluster>> newMapTrackToShowerComponents, Map<Cluster, Track> newMapShowerComponentToTrack, List<Cluster> photons, List<Cluster> neutralClusterCores) {
+ // DEBUG: Check the status of all hits in the event.
+ HitMap hitMapEcal = ((HitMap)(m_event.get("inputHitMapEcal")));
+ HitMap hitMapHcal = ((HitMap)(m_event.get("inputHitMapHcal")));
+ Set<CalorimeterHit> hitsEcal = new HashSet<CalorimeterHit>(hitMapEcal.values());
+ Set<CalorimeterHit> hitsHcal = new HashSet<CalorimeterHit>(hitMapHcal.values());
+ Set<CalorimeterHit> hitsEcalSharedWithTargets = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> hitsHcalSharedWithTargets = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> hitsEcalSharedWithoutTargets = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> hitsHcalSharedWithoutTargets = new HashSet<CalorimeterHit>();
+
+ for (SharedClusterGroup shares : allSharedClusters) {
+ List<SharedCluster> sharedClusters = shares.listAllSharedClusters();
+ for (SharedCluster sharedCluster : sharedClusters) {
+ Set<Cluster> targets = sharedCluster.getTargetClusters();
+ for (CalorimeterHit hit : sharedCluster.getCluster().getCalorimeterHits()) {
+ boolean hitInEcal = hitsEcal.contains(hit);
+ boolean hitInHcal = hitsHcal.contains(hit);
+ if (hitInEcal && hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (!hitInEcal && !hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (targets.size()==0) {
+ if (hitInEcal) {
+ hitsEcalSharedWithoutTargets.add(hit);
+ } else {
+ hitsHcalSharedWithoutTargets.add(hit);
+ }
+ } else {
+ if (hitInEcal) {
+ hitsEcalSharedWithTargets.add(hit);
+ } else {
+ hitsHcalSharedWithTargets.add(hit);
+ }
+ }
+ }
+ }
+ }
+
+ Set<CalorimeterHit> hitsEcalInChargedCore = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> hitsHcalInChargedCore = new HashSet<CalorimeterHit>();
+ for (Track tr : tracksSortedByMomentum) {
+ Set<Cluster> showerComponents = newMapTrackToShowerComponents.get(tr);
+ for (Cluster clus : showerComponents) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean hitInEcal = hitsEcal.contains(hit);
+ boolean hitInHcal = hitsHcal.contains(hit);
+ if (hitInEcal && hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (!hitInEcal && !hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (hitInEcal) {
+ hitsEcalInChargedCore.add(hit);
+ } else {
+ hitsHcalInChargedCore.add(hit);
+ }
+ }
+ }
+ }
+
+ Set<CalorimeterHit> hitsEcalInPhotonCore = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> hitsHcalInPhotonCore = new HashSet<CalorimeterHit>();
+ for (Cluster clus : photons) {
+ Track tr = newMapShowerComponentToTrack.get(clus);
+ if (tr == null) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean hitInEcal = hitsEcal.contains(hit);
+ boolean hitInHcal = hitsHcal.contains(hit);
+ if (hitInEcal && hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (!hitInEcal && !hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (hitInEcal) {
+ hitsEcalInPhotonCore.add(hit);
+ } else {
+ hitsHcalInPhotonCore.add(hit);
+ }
+ }
+ }
+ }
+
+ Set<CalorimeterHit> hitsEcalInNeutralCore = new HashSet<CalorimeterHit>();
+ Set<CalorimeterHit> hitsHcalInNeutralCore = new HashSet<CalorimeterHit>();
+ for (Cluster clus : neutralClusterCores) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean hitInEcal = hitsEcal.contains(hit);
+ boolean hitInHcal = hitsHcal.contains(hit);
+ if (hitInEcal && hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (!hitInEcal && !hitInHcal) { throw new AssertionError("Book-keeping error"); }
+ if (hitInEcal) {
+ hitsEcalInNeutralCore.add(hit);
+ } else {
+ hitsHcalInNeutralCore.add(hit);
+ }
+ }
+ }
+
+ int countIdentifiedHitsEcal = hitsEcalSharedWithTargets.size() + hitsEcalSharedWithoutTargets.size() + hitsEcalInChargedCore.size() + hitsEcalInPhotonCore.size() + hitsEcalInNeutralCore.size();
+ int countIdentifiedHitsHcal = hitsHcalSharedWithTargets.size() + hitsHcalSharedWithoutTargets.size() + hitsHcalInChargedCore.size() + hitsHcalInPhotonCore.size() + hitsHcalInNeutralCore.size();
+
+ System.out.println("Of the "+hitsEcal.size()+" ECAL hits, I identified "+countIdentifiedHitsEcal+":");
+ System.out.println(" "+hitsEcalSharedWithTargets.size()+" are shared and have a target");
+ System.out.println(" "+hitsEcalSharedWithoutTargets.size()+" are shared and have no target");
+ System.out.println(" "+hitsEcalInChargedCore.size()+" are in the core of a charged shower");
+ System.out.println(" "+hitsEcalInPhotonCore.size()+" are in the core of a photon shower");
+ System.out.println(" "+hitsEcalInNeutralCore.size()+" are in the core of a neutral hadron shower");
+ System.out.println("Of the "+hitsHcal.size()+" HCAL hits, I identified "+countIdentifiedHitsHcal+":");
+ System.out.println(" "+hitsHcalSharedWithTargets.size()+" are shared and have a target");
+ System.out.println(" "+hitsHcalSharedWithoutTargets.size()+" are shared and have no target");
+ System.out.println(" "+hitsHcalInChargedCore.size()+" are in the core of a charged shower");
+ System.out.println(" "+hitsHcalInPhotonCore.size()+" are in the core of a photon shower");
+ System.out.println(" "+hitsHcalInNeutralCore.size()+" are in the core of a neutral hadron shower");
+
+ for (Track tr : tracksSortedByMomentum) {
+ List<CalorimeterHit> trackTruthHitsInEcal = findHitsFromTruth_T(tr, hitsEcal);
+ List<CalorimeterHit> trackTruthHitsInHcal = findHitsFromTruth_T(tr, hitsHcal);
+ System.out.println("For track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" with "+trackTruthHitsInEcal.size()+" ECAL + "+trackTruthHitsInHcal.size()+" HCAL hits:");
+ int count;
+ String printme = new String(" ECAL hits: ");
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInEcal) { if (hitsEcalSharedWithTargets.contains(hit)) { count++; } } ; printme += " sharedWithTarget="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInEcal) { if (hitsEcalSharedWithoutTargets.contains(hit)) { count++; } } ; printme += " sharedWithoutTarget="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInEcal) { if (hitsEcalInChargedCore.contains(hit)) { count++; } } ; printme += " chargedCore="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInEcal) { if (hitsEcalInPhotonCore.contains(hit)) { count++; } } ; printme += " photonCore="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInEcal) { if (hitsEcalInNeutralCore.contains(hit)) { count++; } } ; printme += " neutralCore="+count;
+ System.out.println(printme);
+ printme = new String(" HCAL hits: ");
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInHcal) { if (hitsHcalSharedWithTargets.contains(hit)) { count++; } } ; printme += " sharedWithTarget="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInHcal) { if (hitsHcalSharedWithoutTargets.contains(hit)) { count++; } } ; printme += " sharedWithoutTarget="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInHcal) { if (hitsHcalInChargedCore.contains(hit)) { count++; } } ; printme += " chargedCore="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInHcal) { if (hitsHcalInPhotonCore.contains(hit)) { count++; } } ; printme += " photonCore="+count;
+ count = 0; for (CalorimeterHit hit : trackTruthHitsInHcal) { if (hitsHcalInNeutralCore.contains(hit)) { count++; } } ; printme += " neutralCore="+count;
+ System.out.println(printme);
+ }
+ }
}
