lcsim/src/org/lcsim/recon/cheater
diff -u -r1.6 -r1.7
--- ReconCheater.java 26 Oct 2005 19:04:03 -0000 1.6
+++ ReconCheater.java 1 Nov 2005 16:44:37 -0000 1.7
@@ -70,22 +70,19 @@
double DecayDistance;
double ECalResolution, HCalResolution;
double ECalSampling, HCalSampling, HCalDigital; // Hits/GeV
- double pTrackMin, EClusterMin, NDigitalMin;
+ double pTrackMin, EClusterMin;
+ double ECalEnergyMin, HCalEnergyMin, NDigitalMin;
- static boolean Hist = false;
+ boolean usePerfectEnergyFlow;
+
+ double Distance2XCluster, Distance4XCluster;
+
+ static boolean hist = false;
/** */
public ReconCheater()
{
- this(Hist);
- }
-
- public ReconCheater(boolean hist)
- {
- // super(hist);
- this.hist = hist;
-
System.out.println("");
//String text = "ReconCheater version 1.0"+"\n"+
// " Running hep.lcd version "+hep.lcd.Version.getVersion();
@@ -104,6 +101,8 @@
df.setMaximumFractionDigits(3);
}
+ public void setHist(boolean hist) { this.hist = hist; initHist(); }
+
void init()
{
@@ -133,14 +132,27 @@
pTrackMin = cheatingTable.getPTrackMin();
EClusterMin = cheatingTable.getEClusterMin();
+ ECalEnergyMin = cheatingTable.getECalEnergyMin();
+ HCalEnergyMin = cheatingTable.getHCalEnergyMin();
NDigitalMin = cheatingTable.getNDigitalMin();
+ usePerfectEnergyFlow = cheatingTable.usePerfectEnergyFlow();
+
+ Distance2XCluster = cheatingTable.getDistance2XCluster();
+ Distance4XCluster = cheatingTable.getDistance4XCluster();
+
+ initHist();
+ }
+
+ void initHist()
+ {
if (hist) {
AIDA aida = AIDA.defaultInstance();
ITree tree = aida.tree();
try {
tree.mkdir("ReconCheater");
tree.cd("ReconCheater");
+ tree.mkdir("Clusters");
tree.mkdir("Lost"); tree.mkdir("Interacted"); tree.mkdir("Conversion");
tree.cd("..");
}
@@ -151,8 +163,8 @@
protected int nEvt = 0;
boolean first = true, firstEvents = true;
- boolean debug = false;
- boolean hist;
+ boolean debug = true;
+ int imcp = -1;
/** Reconstruct an event by using the Monte Carlo truth to cheat. */
protected void process(EventHeader event)
@@ -194,7 +206,6 @@
init();
}
-
Map<MCParticle,CheatTrack> charged = null; // Tracks from MCParticles
Map<MCParticle,CheatCluster> neutrals = null; // and clusters.
Map<MCParticle,MCParticle> primaryMCParents = null;
@@ -239,64 +250,65 @@
extraTrackList = new ArrayList<CheatTrack>();
List<CheatTrack> skippedTracks = new ArrayList<CheatTrack>();
+ List<CheatTrack> trackList = null;
List<List<CheatTrack>> trackCollections = event.get(CheatTrack.class);
- for (List<CheatTrack> trackList : trackCollections ) {
- String name = event.getMetaData(trackList).getName();
+ for (List<CheatTrack> collection : trackCollections ) {
+ String name = event.getMetaData(collection).getName();
+ if (name.equals("CombinedTracks")) trackList = collection;
+ }
+ if (trackList==null) return particles;
+
+ // Work through track list.
+ for (CheatTrack track : trackList) {
+ MCParticle mcp = track.getMCParticle(); int PDGID = mcp.getPDGID();
+ double p = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()+mcp.getPZ()*mcp.getPZ());
+ if (p<pTrackMin) continue;
- if (name.equals("CombinedTracks")) {
- // Work through track list.
- for (CheatTrack track : trackList) {
- MCParticle mcp = track.getMCParticle(); int PDGID = mcp.getPDGID();
- double p = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()+mcp.getPZ()*mcp.getPZ());
- if (p<pTrackMin) continue;
-
- // Keep track of primary parentage.
- MCParticle primaryParent = findPrimaryParent(mcp);
- primaryMCParents.put(mcp,primaryParent);
-
- // Ignore radiation, compton scattering, ...
- if (!allowRadiation && mcp.getSimulatorStatus().vertexIsNotEndpointOfParent()) continue;
- if (!allowRadiation && checkRadiationParent(mcp)) continue;
- if (mcp.getSimulatorStatus().isBackscatter()) continue;
- // Ignore recoil nucleus.
- if (PDGID==0) continue;
-
- // Record particles that decay or interact.
- boolean decayed = checkDecay(mcp);
- boolean interacted = false;
- if (decayed) decayedParticles.add(mcp);
- else {
- interacted = checkInteraction(mcp);
- if (interacted) interactedParticles.add(mcp);
- }
- // Check for decay or interaction products.
- boolean decayProduct = checkDecay(primaryParent);
- boolean interactionProduct = checkInteractionParent(mcp);
-
- // Skip tracks from decaying or interacting particles.
- if ((allowDecays && decayed && hasNeutralDaughters(mcp)) || // But keep pi's from pi->mu decays, etc.
- (allowNuclearInteractions && interacted) || // Assume interacting tracks will not be used.
-
- // Skip decay and interaction products if switched off.
- (!allowDecays && decayProduct && isLongLived(primaryParent)) ||
- (!allowNuclearInteractions && interactionProduct)) {
- skippedTracks.add(track);
- continue;
- }
+ // Keep track of primary parentage.
+ MCParticle primaryParent = findPrimaryParent(mcp);
+ primaryMCParents.put(mcp,primaryParent);
+
+ // Ignore radiation, compton scattering, ...
+ if (!allowRadiation && mcp.getSimulatorStatus().vertexIsNotEndpointOfParent()) continue;
+ if (!allowRadiation && checkRadiationParent(mcp)) continue;
+ if (mcp.getSimulatorStatus().isBackscatter()) continue;
+ // Ignore recoil nucleus.
+ if (PDGID==0) continue;
+
+ // Record particles that decay or interact.
+ boolean decayed = checkDecay(mcp);
+ boolean interacted = false;
+ if (decayed) decayedParticles.add(mcp);
+ else {
+ interacted = checkInteraction(mcp);
+ if (interacted) interactedParticles.add(mcp);
+ }
+ // Check for decay or interaction products.
+ boolean decayProduct = checkDecay(primaryParent);
+ boolean interactionProduct = checkInteractionParent(mcp);
+
+ // Skip tracks from decaying or interacting particles.
+ if ((allowDecays && decayed && hasNeutralDaughters(mcp)) || // But keep pi's from pi->mu decays, etc.
+ (allowNuclearInteractions && interacted) || // Assume interacting tracks will not be used.
+
+ // Skip decay and interaction products if switched off.
+ (!allowDecays && decayProduct && isLongLived(primaryParent)) ||
+ (!allowNuclearInteractions && interactionProduct)) {
+ skippedTracks.add(track);
+ continue;
+ }
- // Record charged particle tracks that were found
- if (chargedList.contains(mcp)) {
- chargedList.remove(mcp);
-
- // Add reconstructed track to charged list.
- if (charged.get(mcp)==null) charged.put(mcp,track);
- else System.err.println(ClassName+" track is already in map!");
- }
- // or add as an extra track.
- else
- extraTrackList.add(track);
- }
+ // Record charged particle tracks that were found
+ if (chargedList.contains(mcp)) {
+ chargedList.remove(mcp);
+
+ // Add reconstructed track to charged list.
+ if (charged.get(mcp)==null) charged.put(mcp,track);
+ else System.err.println(ClassName+" track is already in map!");
}
+ // or add as an extra track.
+ else
+ extraTrackList.add(track);
}
// Remove secondaries if primary track has been reconstructed.
List<CheatTrack> secondaryTracks = new ArrayList<CheatTrack>();
@@ -332,7 +344,7 @@
List<CheatTrack> primaryTracks = new ArrayList<CheatTrack>();
for (CheatTrack track : extraTrackList) {
MCParticle mcp = track.getMCParticle();
- if (getStatusCode(mcp).equals("GEN_PREDECAY") && !isLongLived(mcp)) {
+ if (getStatusCode(mcp).equals("GEN_PREDECAY") && (allowDecays || !isLongLived(mcp))) {
primaryTracks.add(track); skippedTracks.add(track);
}
}
@@ -360,73 +372,81 @@
extraClusterList = new ArrayList<CheatCluster>();
List<CheatCluster> skippedClusters = new ArrayList<CheatCluster>();
+ List<CheatCluster> clusterList = null;
List<List<CheatCluster>> clusterCollections = event.get(CheatCluster.class);
- for (List<CheatCluster> clusterList : clusterCollections ) {
- String name = event.getMetaData(clusterList).getName();
-
- if (name.equals("RefinedClusters")) {
- for (CheatCluster cluster : clusterList) {
- MCParticle mcp = cluster.getMCParticle(); int PDGID = mcp.getPDGID();
- double charge = mcp.getCharge();
- double p = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()+mcp.getPZ()*mcp.getPZ());
- double eMeasured = getMeasuredEnergy(cluster);
- if (eMeasured<EClusterMin) continue;
-
- // Keep track of primary parentage.
- MCParticle primaryParent = findPrimaryParent(mcp);
- primaryMCParents.put(mcp,primaryParent);
-
- // Ignore radiation, compton scattering, ...
- if (!allowRadiation && mcp.getSimulatorStatus().vertexIsNotEndpointOfParent()) continue;
- if (!allowRadiation && checkRadiationParent(mcp)) continue;
- if (mcp.getSimulatorStatus().isBackscatter()) continue;
- // Ignore recoil nucleus.
- if (PDGID==0) continue;
-
- // Record particles that decay or interact.
- boolean decayed = checkDecay(mcp);
- boolean interacted = false;
- if (decayed && !decayedParticles.contains(mcp)) decayedParticles.add(mcp);
- else {
- interacted = checkInteraction(mcp);
- if (interacted && !interactedParticles.contains(mcp)) interactedParticles.add(mcp);
- }
- // Check for decay or interaction products.
- boolean decayProduct = checkDecay(primaryParent);
- boolean interactionProduct = checkInteractionParent(mcp);
-
- // Check if cluster is associated with a reconstructed track.
- boolean unAssociated = charge==0;
- if (!unAssociated) {
- boolean found = charged.containsKey(mcp);
- if (!found) for (CheatTrack t : extraTrackList) if (t.getMCParticle()==mcp) found = true;
- if (!found) for (CheatTrack t : skippedTracks) if (t.getMCParticle()==mcp) found = true;
- if (!found && !checkForParent(getParent(mcp))) unAssociated = true;
- }
- if (!unAssociated) continue;
-
- // Skip decay and interaction products if switched off.
- if ((!allowDecays && decayProduct && isLongLived(primaryParent)) ||
- (!allowNuclearInteractions && interactionProduct)) {
- skippedClusters.add(cluster);
- continue;
- }
+ for (List<CheatCluster> collection : clusterCollections ) {
+ String name = event.getMetaData(collection).getName();
+ if (name.equals("RefinedClusters")) clusterList = collection;
+ }
+ if (clusterList==null) return particles;
+
+ // Analyze clusters.
+ if (hist) analyzeReconstructedClusters(clusterList);
+
+ // Remove nearby clusters.
+ if (!usePerfectEnergyFlow) clusterList = getIsolatedClusters(clusterList);
+
+ // Work through cluster list.
+ for (CheatCluster cluster : clusterList) {
+ MCParticle mcp = cluster.getMCParticle(); int PDGID = mcp.getPDGID();
+ double charge = mcp.getCharge();
+ double p = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()+mcp.getPZ()*mcp.getPZ());
+ double eMeasured = getMeasuredEnergy(cluster);
+ if (eMeasured<EClusterMin) continue;
+
+ // Keep track of primary parentage.
+ MCParticle primaryParent = findPrimaryParent(mcp);
+ primaryMCParents.put(mcp,primaryParent);
+
+ // Ignore radiation, compton scattering, ...
+ if (!allowRadiation && mcp.getSimulatorStatus().vertexIsNotEndpointOfParent()) continue;
+ if (!allowRadiation && checkRadiationParent(mcp)) continue;
+ if (mcp.getSimulatorStatus().isBackscatter()) continue;
+ // Ignore recoil nucleus.
+ if (PDGID==0) continue;
+
+ // Record particles that decay or interact.
+ boolean decayed = checkDecay(mcp);
+ boolean interacted = false;
+ if (decayed && !decayedParticles.contains(mcp)) decayedParticles.add(mcp);
+ else {
+ interacted = checkInteraction(mcp);
+ if (interacted && !interactedParticles.contains(mcp)) interactedParticles.add(mcp);
+ }
+ // Check for decay or interaction products.
+ boolean decayProduct = checkDecay(primaryParent);
+ boolean interactionProduct = checkInteractionParent(mcp);
+
+ // Check if cluster is associated with a reconstructed track.
+ boolean unAssociated = charge==0;
+ if (!unAssociated) {
+ boolean found = charged.containsKey(mcp);
+ if (!found) for (CheatTrack t : extraTrackList) if (t.getMCParticle()==mcp) found = true;
+ if (!found) for (CheatTrack t : skippedTracks) if (t.getMCParticle()==mcp) found = true;
+ if (!found && !checkForParent(getParent(mcp))) unAssociated = true;
+ }
+ if (!unAssociated) continue;
+
+ // Skip decay and interaction products if switched off.
+ if ((!allowDecays && decayProduct && isLongLived(primaryParent)) ||
+ (!allowNuclearInteractions && interactionProduct)) {
+ skippedClusters.add(cluster);
+ continue;
+ }
- // Record neutral particle clusters that were found
- if (neutralList.contains(mcp)) {
- neutralList.remove(mcp);
-
- // Add reconstructed cluster to neutrals list.
- if (neutrals.get(mcp)==null) neutrals.put(mcp,cluster);
- else System.err.println(ClassName+" cluster is already in map!");
- }
- // or add as an extra (unexpected) cluster.
- else
- //if ((!decayProduct || acceptDecayNeutrals) &&
- //(!interactionProduct || acceptNuclearInteractionNeutrals)) {
- extraClusterList.add(cluster);
- }
+ // Record neutral particle clusters that were found
+ if (neutralList.contains(mcp)) {
+ neutralList.remove(mcp);
+
+ // Add reconstructed cluster to neutrals list.
+ if (neutrals.get(mcp)==null) neutrals.put(mcp,cluster);
+ else System.err.println(ClassName+" cluster is already in map!");
}
+ // or add as an extra (unexpected) cluster.
+ else
+ //if ((!decayProduct || acceptDecayNeutrals) &&
+ //(!interactionProduct || acceptNuclearInteractionNeutrals)) {
+ extraClusterList.add(cluster);
}
if (firstEvents) {
text = " Total of "+neutrals.keySet().size()+" neutral clusters found:";
@@ -661,6 +681,112 @@
return particles;
}
+ List<CheatCluster> getIsolatedClusters(List<CheatCluster> clusterList)
+ {
+ List<CheatCluster> result = new ArrayList<CheatCluster>();
+
+ // Work through cluster list.
+ for (CheatCluster cluster : clusterList) {
+ MCParticle mcp = cluster.getMCParticle(); int PDGID = mcp.getPDGID();
+ boolean Klong = PDGID==130;
+ boolean neutron = Math.abs(PDGID)==2112;
+ if (Klong || neutron) {
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double energy = getMeasuredEnergy(cluster);
+ if (energy<EClusterMin) continue;
+
+ // Check for decay or interaction products.
+ MCParticle primaryParent = findPrimaryParent(mcp); primaryMCParents.put(mcp,primaryParent);
+ boolean decayProduct = checkDecay(primaryParent);
+ boolean interactionProduct = checkInteractionParent(mcp);
+
+ // Skip decay and interaction products if switched off.
+ if ((!allowDecays && decayProduct && isLongLived(primaryParent)) ||
+ (!allowNuclearInteractions && interactionProduct))
+ continue;
+
+ // Determine distance to nearest large cluster.
+ double distanceLarge2XCluster = 1000., distanceLarge4XCluster = 1000.;
+ for (CheatCluster cluster2 : clusterList) {
+ if (cluster2==cluster) continue;
+ MCParticle mcp2 = cluster2.getMCParticle(); int PDGID2 = mcp2.getPDGID();
+ Hep3Vector endPoint2 = mcp2.getEndPoint();
+ double energy2 = getMeasuredEnergy(cluster2);
+ if (energy2<EClusterMin) continue;
+
+ double dx = endPoint2.x()-endPoint.x(), dy = endPoint2.y()-endPoint.y(), dz = endPoint2.z()-endPoint.z();
+ double d = Math.sqrt(dx*dx + dy*dy + dz*dz);
+ if ((energy2>(2.*energy)) && d<distanceLarge2XCluster) distanceLarge2XCluster = d;
+ if ((energy2>(4.*energy)) && d<distanceLarge4XCluster) distanceLarge4XCluster = d;
+ }
+ if (distanceLarge2XCluster<Distance2XCluster) continue;
+ if (distanceLarge4XCluster<Distance4XCluster) continue;
+ }
+ result.add(cluster);
+ }
+ return result;
+ }
+
+
+ void analyzeReconstructedClusters(List<CheatCluster> clusterList)
+ {
+ AIDA aida = AIDA.defaultInstance();
+ ITree tree = aida.tree();
+ tree.cd("ReconCheater"); tree.cd("Clusters");
+
+ // Work through cluster list.
+ int nNeutrals = 0, nRemainingNeutrals = 0;
+ for (CheatCluster cluster : clusterList) {
+ MCParticle mcp = cluster.getMCParticle(); int PDGID = mcp.getPDGID();
+ boolean Klong = PDGID==130;
+ boolean neutron = Math.abs(PDGID)==2112;
+ if (!Klong && !neutron) continue;
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double energy = getMeasuredEnergy(cluster);
+ if (energy<EClusterMin) continue;
+
+ // Check for decay or interaction products.
+ MCParticle primaryParent = findPrimaryParent(mcp); primaryMCParents.put(mcp,primaryParent);
+ boolean decayProduct = checkDecay(primaryParent);
+ boolean interactionProduct = checkInteractionParent(mcp);
+
+ // Skip decay and interaction products if switched off.
+ if ((!allowDecays && decayProduct && isLongLived(primaryParent)) ||
+ (!allowNuclearInteractions && interactionProduct))
+ continue;
+ nNeutrals++;
+
+ // Plot distance to nearest large cluster.
+ double distance = 1000., distanceLarge2XCluster = 1000., distanceLarge4XCluster = 1000.;
+ for (CheatCluster cluster2 : clusterList) {
+ if (cluster2==cluster) continue;
+ MCParticle mcp2 = cluster2.getMCParticle(); int PDGID2 = mcp2.getPDGID();
+ Hep3Vector endPoint2 = mcp2.getEndPoint();
+ double energy2 = getMeasuredEnergy(cluster2);
+ if (energy2<EClusterMin) continue;
+
+ double dx = endPoint2.x()-endPoint.x(), dy = endPoint2.y()-endPoint.y(), dz = endPoint2.z()-endPoint.z();
+ double d = Math.sqrt(dx*dx + dy*dy + dz*dz);
+ if (d<100.) {
+ aida.cloud1D("Nearby ratio").fill(energy2/energy);
+ }
+ if ((energy2>(energy/2.)) && d<distance) distance = d;
+ if ((energy2>(2.*energy)) && d<distanceLarge2XCluster) distanceLarge2XCluster = d;
+ if ((energy2>(4.*energy)) && d<distanceLarge4XCluster) distanceLarge4XCluster = d;
+ }
+ if (distance<1000.) aida.cloud1D("Nearest 0.5X cluster").fill(distance);
+ if (distanceLarge2XCluster<1000.) aida.cloud1D("Nearest Large2X cluster").fill(distanceLarge2XCluster);
+ if (distanceLarge4XCluster<1000.) aida.cloud1D("Nearest Large4X cluster").fill(distanceLarge4XCluster);
+
+ // Apply isolation cuts.
+ if (distanceLarge2XCluster<Distance2XCluster) continue;
+ if (distanceLarge4XCluster<Distance4XCluster) continue;
+ nRemainingNeutrals++;
+ }
+ aida.cloud1D("# neutral hadrons").fill(nNeutrals);
+ aida.cloud1D("# remaining neutrals").fill(nRemainingNeutrals);
+ tree.cd(".."); tree.cd("..");
+ }
void analyzeReconstructedParticles(List<ReconstructedParticle> particles, List<MCParticle> chargedList,
List<MCParticle> neutralList)
{
@@ -1106,15 +1232,19 @@
boolean photon = mcp.getPDGID()==22;
if (photon) {
energy = energy/ECalSampling;
+ if (energy<ECalEnergyMin) energy = 0.;
}
else {
- if (HCalDigital==0) energy = energy/HCalSampling;
+ if (HCalDigital==0) {
+ energy = energy/HCalSampling;
+ if (energy<HCalEnergyMin) energy = 0.;
+ }
else {
int nhits = 0;
for (CalorimeterHit hit : cluster.getCalorimeterHits())
if (hit.getTime()<100.) nhits++;
energy = nhits/HCalDigital;
- if (nhits<NDigitalMin) energy = 0.;
+ if (nhits<NDigitalMin || energy<HCalEnergyMin) energy = 0.;
}
}
Hep3Vector pV = mcp.getMomentum();
@@ -1124,16 +1254,18 @@
if (photon) {
if (useECalParameterization) {
r = 1. - random.nextGaussian()*ECalResolution/Math.sqrt(mcp.getEnergy());
- if (r<0.1) r = 0.1;
+ if (r<0.1) r = 0.1; energy = r * pV.magnitude();
+ if (energy<ECalEnergyMin) energy = 0.;
}
}
else
if (useHCalParameterization) {
r = 1. - random.nextGaussian()*HCalResolution/Math.sqrt(mcp.getEnergy());
- if (r<0.1) r = 0.1;
+ if (r<0.1) r = 0.1; energy = r * pV.magnitude();
+ if (energy<HCalEnergyMin) energy = 0.;
}
}
- return r * pV.magnitude();
+ return energy;
}
Map<MCParticle,Integer> map = null;
@@ -1158,7 +1290,7 @@
double totalEventEnergy;
/** Get charged and neutral lists for MCParticles generated by the Pythia Monte Carlo. */
- public void getPythiaMCParticles(EventHeader event, List<MCParticle> chargedList, List<MCParticle> neutralList)
+ void getPythiaMCParticles(EventHeader event, List<MCParticle> chargedList, List<MCParticle> neutralList)
{
// Get Monte Carlo particles from LCSIM event.
String text = " Get Pythia Monte Carlo charged and neutral lists.";
