lcsim/src/org/lcsim/recon/cheater
diff -u -r1.4 -r1.5
--- ReconCheater.java 28 Sep 2005 00:35:51 -0000 1.4
+++ ReconCheater.java 18 Oct 2005 17:09:50 -0000 1.5
@@ -32,9 +32,12 @@
import org.lcsim.util.loop.LCSimLoop;
import org.lcsim.util.aida.AIDA;
+// HEP analysis and framework packages
+import hep.aida.*; // Abstract analysis interfaces
import hep.physics.particle.*; // Particle,...
import hep.physics.vec.*; // Hep3Vector
+// Java packages
import java.io.File;
import java.net.URL;
import java.text.*; // DecimalFormat
@@ -50,16 +53,23 @@
public class ReconCheater extends Driver
implements ConditionsListener
{
- CheatingTable Cheating;
+ CheatingTable cheatingTable;
- boolean useFullTruth = false; // Use primary MCParticles.
- boolean useTruth = false; // Use MCParticles energies.
- boolean useECalParameterization = false, useHCalParameterization = true;
- boolean useReconPhotons = true, useReconNeutralHadrons = false;
- boolean allowNuclearInteractions = true;
-
- double ECalResolution = 0.18, HCalResolution = 0.60;
- double ECalSampling = 0.012, HCalSampling = 0.008, HCalDigital = 11.7; // Hits/GeV
+ boolean useFullTruth; // Use primary MCParticles.
+ boolean useTruth; // Use MCParticles energies.
+ boolean useECalParameterization, useHCalParameterization;
+ boolean useReconPhotons, useReconNeutralHadrons;
+ boolean allowDecays;
+ boolean acceptDecayProducts; // Note: Charged tracks affected by TrackingCheater.
+ boolean acceptDecayNeutrals;
+ boolean allowNuclearInteractions;
+ boolean acceptNuclearInteractionProducts;
+ boolean acceptNuclearInteractionNeutrals;
+ boolean allowRadiation;
+
+ double ECalResolution, HCalResolution;
+ double ECalSampling, HCalSampling, HCalDigital; // Hits/GeV
+ double pTrackMin, EClusterMin, NDigitalMin;
static boolean Hist = false;
@@ -73,11 +83,12 @@
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();
- String text = "ReconCheater version 0.7"+"\n";
+ String text = "ReconCheater version 0.8"+"\n";
System.out.println(text); System.err.println(text); // if (!XReconstruction.getBatchMode()) out.println(text);
if (!TrackingCheater.INITIALIZED) {
@@ -95,35 +106,69 @@
void init()
{
- useFullTruth = Cheating.useFullTruth();
- useTruth = Cheating.useTruth();
+ useFullTruth = cheatingTable.useFullTruth();
+ useTruth = cheatingTable.useTruth();
if (useFullTruth) useTruth = true;
- useECalParameterization = Cheating.useECalParameterization();
- useHCalParameterization = Cheating.useHCalParameterization();
- useReconPhotons = Cheating.useReconPhotons();
- useReconNeutralHadrons = Cheating.useReconNeutralHadrons();
- allowNuclearInteractions = Cheating.allowNuclearInteractions();
+ useECalParameterization = cheatingTable.useECalParameterization();
+ useHCalParameterization = cheatingTable.useHCalParameterization();
+ useReconPhotons = cheatingTable.useReconPhotons();
+ useReconNeutralHadrons = cheatingTable.useReconNeutralHadrons();
+
+ allowDecays = cheatingTable.allowDecays();
+ acceptDecayProducts = cheatingTable.acceptDecayProducts();
+ acceptDecayNeutrals = cheatingTable.acceptDecayNeutrals();
+ allowNuclearInteractions = cheatingTable.allowNuclearInteractions();
+ acceptNuclearInteractionProducts = cheatingTable.acceptNuclearInteractionProducts();
+ acceptNuclearInteractionNeutrals = cheatingTable.acceptNuclearInteractionNeutrals();
+ allowRadiation = cheatingTable.allowRadiation();
+
+ ECalResolution = cheatingTable.getECalResolution();
+ ECalSampling = cheatingTable.getECalSampling();
+ HCalResolution = cheatingTable.getHCalResolution();
+ HCalSampling = cheatingTable.getHCalSampling();
+ HCalDigital = cheatingTable.getHCalDigital();
+
+ pTrackMin = cheatingTable.getPTrackMin();
+ EClusterMin = cheatingTable.getEClusterMin();
+ NDigitalMin = cheatingTable.getNDigitalMin();
+
+ if (hist) {
+ AIDA aida = AIDA.defaultInstance();
+ ITree tree = aida.tree();
+ try {
+ tree.mkdir("ReconCheater");
+ tree.cd("ReconCheater");
+ tree.mkdir("Lost"); tree.mkdir("Interacted"); tree.mkdir("Conversion");
+ tree.cd("..");
+ }
+ catch (hep.aida.ref.tree.TreeObjectAlreadyExistException ex) {
+ }
+ }
}
protected int nEvt = 0;
boolean first = true, firstEvents = true;
+ boolean debug = false;
boolean hist;
+ /** Reconstruct an event by using the Monte Carlo truth to cheat. */
protected void process(EventHeader event)
{
super.process(event);
nEvt++;
if (firstEvents) System.out.println(" Event #"+nEvt+" (ReconCheater)");
+ String text = " Using reconstruction cheater.";
+ if (firstEvents) System.out.println(text); if (first) System.err.println(text);
- hist = getHistogramLevel() > 0;
- if (hist) System.err.println(" ReconCheater: hist = "+hist);
+ hist = hist || getHistogramLevel() > 0;
+ if (first && hist) System.err.println(" ReconCheater: hist = "+hist);
- if (Cheating == null) {
+ if (cheatingTable == null) {
ConditionsSet conditions = getConditionsManager().getConditions("Cheating");
conditions.addConditionsListener(this);
- Cheating = new CheatingTable(conditions);
+ cheatingTable = new CheatingTable(conditions);
init();
}
@@ -137,19 +182,25 @@
event.put("ReconCheater", particles, ReconstructedParticle.class, 0);
first = false;
- if (nEvt>=3) firstEvents = false;
+ if (nEvt>=3) firstEvents = debug = false;
}
public void conditionsChanged(ConditionsEvent event)
{
ConditionsSet conditions = getConditionsManager().getConditions("Cheating");
- Cheating = new CheatingTable(conditions);
+ cheatingTable = new CheatingTable(conditions);
init();
}
Map<MCParticle,CheatTrack> charged = null; // Tracks from MCParticles
Map<MCParticle,CheatCluster> neutrals = null; // and clusters.
+ Map<MCParticle,MCParticle> primaryMCParents = null;
+ List<MCParticle> decayedParticles = null;
+ List<MCParticle> interactedParticles = null;
+ List<MCParticle> lostParticles = null;
+ List<CheatTrack> extraTrackList = null;
+ List<CheatCluster> extraClusterList = null;
double finalEnergy;
/** Get reconstructed particles from Cheaters that match MCParticles,
@@ -157,333 +208,806 @@
*/
List<ReconstructedParticle> getReconstructedParticles(EventHeader event)
{
- String text = " Get selected particles.";
+ String text = " Get selected particles from reconstructed tracks and clusters, and MCParticles.";
if (firstEvents) System.out.println(text); if (first) System.err.println(text);
- // Get lists of charged and neutral MC particles.
+ // Get lists of charged and neutral final state MC particles.
List<MCParticle> chargedList = new ArrayList<MCParticle>();
List<MCParticle> neutralList = new ArrayList<MCParticle>();
getPythiaMCParticles(event,chargedList,neutralList);
+ // Reconstruct particles.
List<ReconstructedParticle> particles = new ArrayList();
+ // Keep records.
charged = new HashMap<MCParticle,CheatTrack>();
neutrals = new HashMap<MCParticle,CheatCluster>();
- Map<MCParticle,MCParticle> primaryMCParents = new HashMap<MCParticle,MCParticle>();
+ primaryMCParents = new HashMap<MCParticle,MCParticle>();
+ decayedParticles = new ArrayList<MCParticle>();
+ interactedParticles = new ArrayList<MCParticle>();
+ lostParticles = new ArrayList<MCParticle>();
// Look for tracks in TrackingCheater list from charged MCParticles.
- text = " Look for tracks from "+chargedList.size()+" charged MCParticles:";
- for (MCParticle p : chargedList) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
-
- List<CheatTrack> extraTrackList = new ArrayList<CheatTrack>();
+ if (firstEvents) {
+ text = " Look for tracks from "+chargedList.size()+" charged MCParticles:";
+ for (MCParticle p : chargedList) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ extraTrackList = new ArrayList<CheatTrack>();
List<CheatTrack> skippedTracks = new ArrayList<CheatTrack>();
List<List<CheatTrack>> trackCollections = event.get(CheatTrack.class);
- //if (firstEvents) System.out.println(" Size of track lists = "+trackCollections.size());
- finalEnergy = 0;
- for (List<CheatTrack> collection : trackCollections ) {
- String name = event.getMetaData(collection).getName();
- //if (first) System.out.println(" "+name+" has "+collection.size()+" tracks.");
+ for (List<CheatTrack> trackList : trackCollections ) {
+ String name = event.getMetaData(trackList).getName();
if (name.equals("CombinedTracks")) {
- for (CheatTrack track : collection) {
+ // 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());
- String statusCode = getStatusCode(mcp);
- boolean finalState = statusCode.equals("GEN_FINAL_STATE");
+ if (p<pTrackMin) continue;
- List<MCParticle> parents = mcp.getParents();
- MCParticle parent = null; if (!parents.isEmpty()) parent = parents.get(0);
- boolean parentFinalState = false; int parentPDGID = 0;
- if (parent!=null) {
- parentFinalState = parent.getGeneratorStatus()==Particle.FINAL_STATE;
- parentPDGID = parent.getPDGID();
-
- // Check for final state particles.
- if (!finalState) {
- // Get final state parent particle.
- List<MCParticle> grandParents = parent.getParents();
- while (!parentFinalState && !grandParents.isEmpty()) {
- parent = grandParents.get(0); grandParents = parent.getParents();
- parentFinalState = parent.getGeneratorStatus()==Particle.FINAL_STATE;
- parentPDGID = parent.getPDGID();
- }
- }
- primaryMCParents.put(mcp,parent);
+ // 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);
+
+ if (firstEvents) {
+ text = " "+mcp.getType().getName()+" #"+map.get(mcp);
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double r = Math.sqrt(endPoint.x()*endPoint.x()+endPoint.y()*endPoint.y());
+ if (decayed) System.out.println(" Track from"+text+" decaying at r="+df.format(r)+" mm");
+ else if (interacted)
+ System.out.println(" Track from"+text+" interacting at r="+df.format(r)+" mm");
+ }
+ // Skip tracks from decaying or interaction 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) ||
+ (!allowNuclearInteractions && interactionProduct)) {
+ skippedTracks.add(track);
+ continue;
}
+ // Record charged particle tracks that were found
if (chargedList.contains(mcp)) {
chargedList.remove(mcp);
- //if (!mcp.getSimulatorStatus().isBackscatter() &&
- // !mcp.getSimulatorStatus().vertexIsNotEndpointOfParent() &&
- //!mcp.getSimulatorStatus().isDecayedInTracker() &&
- // !mcp.getSimulatorStatus().isDecayedInCalorimeter()) {
- if (charged.get(mcp)==null) charged.put(mcp,track);
- else System.err.println(" track is already in map!");
- finalEnergy += mcp.getEnergy();
- }
- // Add as an extra track.
- else {
- // Skip if ...
- //if (finalState &&
- if (!(parentFinalState && (Math.abs(parentPDGID)==211 || Math.abs(parentPDGID)==321)) &&
- !(parentFinalState && (Math.abs(parentPDGID)==2212 || Math.abs(parentPDGID)==2112))) {
-
- if (p>0.100 &&
- !mcp.getSimulatorStatus().isBackscatter() &&
- !mcp.getSimulatorStatus().vertexIsNotEndpointOfParent() &&
- !mcp.getSimulatorStatus().isDecayedInCalorimeter()) {
- extraTrackList.add(track);
- }
- }
- else
- if (p>0.100) skippedTracks.add(track);
+
+ // 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);
}
}
}
- text = " Total of "+charged.keySet().size()+" charged particle tracks found:";
- for (MCParticle p : charged.keySet()) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
- text = " Total of "+chargedList.size()+" interacting charged particles:";
- for (MCParticle p : chargedList) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
- text = " Total of "+extraTrackList.size()+" tracks from interactions or decays:";
- for (CheatTrack t : extraTrackList) { MCParticle p = t.getMCParticle(); text+=" "+map.get(p); }
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
- text = " Total of "+skippedTracks.size()+" tracks skipped:";
- for (CheatTrack t : skippedTracks) { MCParticle p = t.getMCParticle(); text+=" "+map.get(p); }
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
+ // Remove secondaries if primary track has been reconstructed.
+ List<CheatTrack> secondaryTracks = new ArrayList<CheatTrack>();
+ for (CheatTrack track : extraTrackList) {
+ MCParticle mcp = track.getMCParticle();
+ MCParticle primaryParent = primaryMCParents.get(mcp);
+ List<MCParticle> parents = mcp.getParents();
+ MCParticle parent = parents.get(0);
+ boolean found = charged.containsKey(parent) || charged.containsKey(primaryParent);
+ for (CheatTrack t : extraTrackList)
+ if (t.getMCParticle()==parent || t.getMCParticle()==primaryParent) found = true;
+ if (found) {
+ secondaryTracks.add(track); skippedTracks.add(track);
+ }
+ }
+ for (CheatTrack track : secondaryTracks) extraTrackList.remove(track);
- // Look for clusters in ClusterCheater list from neutral MCParticles.
- text = " Look for clusters from "+neutralList.size()+" neutral MCParticles:";
- for (MCParticle p : neutralList) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
+ if (firstEvents) {
+ text = " Total of "+charged.keySet().size()+" charged particle tracks found:";
+ for (MCParticle p : charged.keySet()) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ text = " Total of "+chargedList.size()+" interacting or lost charged particles:";
+ for (MCParticle p : chargedList) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ text = " Total of "+extraTrackList.size()+" tracks from interactions or decays:";
+ for (CheatTrack t : extraTrackList) { MCParticle p = t.getMCParticle(); text+=" "+map.get(p); }
+ System.out.println(text); //if (first) System.err.println(text);
+ text = " Total of "+skippedTracks.size()+" tracks skipped:";
+ for (CheatTrack t : skippedTracks) { MCParticle p = t.getMCParticle(); text+=" "+map.get(p); }
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ // Record interacting, decaying or lost particles.
+ for (MCParticle mcp : chargedList) {
+ if (checkDecay(mcp) && !decayedParticles.contains(mcp)) decayedParticles.add(mcp);
+ else if (checkInteraction(mcp) && !interactedParticles.contains(mcp)) interactedParticles.add(mcp);
+ else lostParticles.add(mcp);
+ }
- List<CheatCluster> extraClusterList = new ArrayList<CheatCluster>();
+ // Look for clusters in ClusterCheater list from neutral MCParticles.
+ if (firstEvents) {
+ text = " Look for clusters from "+neutralList.size()+" neutral MCParticles:";
+ for (MCParticle p : neutralList) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ extraClusterList = new ArrayList<CheatCluster>();
List<CheatCluster> skippedClusters = new ArrayList<CheatCluster>();
- List<List<CheatCluster>> collections = event.get(CheatCluster.class);
- //if (firstEvents) System.out.println(" Size of cluster lists = "+collections.size());
- for (List<CheatCluster> collection : collections ) {
- String name = event.getMetaData(collection).getName();
- //if (first) System.out.println(" "+name+" has "+collection.size()+" clusters.");
+ 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 : collection) {
+ 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());
- String statusCode = getStatusCode(mcp);
- boolean finalState = statusCode.equals("GEN_FINAL_STATE");
+ double eMeasured = getMeasuredEnergy(cluster);
+ if (eMeasured<EClusterMin) continue;
- List<MCParticle> parents = mcp.getParents();
- MCParticle parent = null; if (!parents.isEmpty()) parent = parents.get(0);
- boolean parentFinalState = false; int parentPDGID = 0;
- if (parent!=null) {
- parentFinalState = parent.getGeneratorStatus()==Particle.FINAL_STATE;
- parentPDGID = parent.getPDGID();
-
- // Check for final state particles.
- if (!finalState) {
- // Get final state parent particle.
- List<MCParticle> grandParents = parent.getParents();
- while (!parentFinalState && !grandParents.isEmpty()) {
- parent = grandParents.get(0); grandParents = parent.getParents();
- parentFinalState = parent.getGeneratorStatus()==Particle.FINAL_STATE;
- parentPDGID = parent.getPDGID();
- }
- }
- primaryMCParents.put(mcp,parent);
+ // 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) unAssociated = true;
}
+ if (!unAssociated) continue;
- // Check if generator neutral list contains cluster MCParticle.
+ // Skip decay and interaction products if switched off.
+ if ((!allowDecays && decayProduct) ||
+ (!allowNuclearInteractions && interactionProduct)) {
+ skippedClusters.add(cluster);
+ continue;
+ }
+
+ // Record neutral particle clusters that were found
if (neutralList.contains(mcp)) {
+ neutralList.remove(mcp);
- // Check if neutral interacted or decayed.
- boolean interacted = false;
- for (CheatTrack t : extraTrackList) { if (t.getMCParticle()==mcp) interacted = true; }
- if (!interacted) {
- // Neutral is removed from list only if it didn't interact.
- neutralList.remove(mcp);
-
- // Add cluster to neutrals list.
- if (neutrals.get(mcp)==null) neutrals.put(mcp,cluster);
- else System.err.println(" cluster is already in map!");
- finalEnergy += mcp.getEnergy();
- }
+ // 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!");
}
- // Add as an extra (unexpected) cluster.
- else {
- if (!(parentFinalState && (Math.abs(parentPDGID)==211 || Math.abs(parentPDGID)==321)) &&
- !(parentFinalState && (Math.abs(parentPDGID)==2212 || Math.abs(parentPDGID)==2112))) {
- if (charge==0 && p>0.100 &&
- !mcp.getSimulatorStatus().isBackscatter() &&
- !mcp.getSimulatorStatus().vertexIsNotEndpointOfParent() &&
- !mcp.getSimulatorStatus().isDecayedInCalorimeter()) {
- extraClusterList.add(cluster);
- }
- }
- else
- if (p>0.100) skippedClusters.add(cluster);
+ // or add as an extra (unexpected) cluster.
+ else
+ if ((!decayProduct || acceptDecayNeutrals) &&
+ (!interactionProduct || acceptNuclearInteractionNeutrals)) {
+ extraClusterList.add(cluster);
}
}
}
}
- text = " Total of "+neutrals.keySet().size()+" neutral clusters found:";
- for (MCParticle p : neutrals.keySet()) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
- text = " Total of "+neutralList.size()+" interacting or lost neutrals:";
- for (MCParticle p : neutralList) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
- text = " Total of "+extraClusterList.size()+" clusters from interactions or decays:";
- for (CheatCluster t : extraClusterList) { MCParticle p = t.getMCParticle(); text+=" "+map.get(p); }
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
- text = " Total of "+skippedClusters.size()+" clusters skipped:";
- for (CheatCluster t : skippedClusters) { MCParticle p = t.getMCParticle(); text+=" "+map.get(p); }
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
-
- // Find parents of interaction tracks and clusters.
- if (firstEvents) System.out.println(" Find parents of interaction tracks and clusters.");
- List<MCParticle> parents = new ArrayList<MCParticle>();
+ if (firstEvents) {
+ text = " Total of "+neutrals.keySet().size()+" neutral clusters found:";
+ for (MCParticle p : neutrals.keySet()) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ text = " Total of "+neutralList.size()+" interacting or lost neutrals:";
+ for (MCParticle p : neutralList) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ text = " Total of "+extraClusterList.size()+" clusters from interactions or decays:";
+ for (CheatCluster c : extraClusterList) { MCParticle p = c.getMCParticle(); text+=" "+map.get(p); }
+ System.out.println(text); //if (first) System.err.println(text);
+ text = " Total of "+skippedClusters.size()+" clusters skipped:";
+ for (CheatCluster c : skippedClusters) { MCParticle p = c.getMCParticle(); text+=" "+map.get(p); }
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ // Record interacting, decaying or lost particles.
+ for (MCParticle mcp : neutralList) {
+ if (checkDecay(mcp)) decayedParticles.add(mcp);
+ else if (checkInteraction(mcp)) interactedParticles.add(mcp);
+ else lostParticles.add(mcp);
+ }
+ if (firstEvents) {
+ if (decayedParticles.size()>0) {
+ text = " Total of "+decayedParticles.size()+" decayed particles:";
+ for (MCParticle p : decayedParticles) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ if (interactedParticles.size()>0) {
+ text = " Total of "+interactedParticles.size()+" interacted particles:";
+ for (MCParticle p : interactedParticles) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ if (lostParticles.size()>0) {
+ text = " Total of "+lostParticles.size()+" lost particles:";
+ for (MCParticle p : lostParticles) text+=" "+map.get(p);
+ System.out.println(text); //if (first) System.err.println(text);
+ }
+ }
- for (CheatTrack track : skippedTracks) {
- MCParticle mcp = track.getMCParticle();
- MCParticle parent = primaryMCParents.get(mcp);
- if (parent==null) continue;
- if (!charged.containsKey(parent) && !neutrals.containsKey(parent) &&
- !parents.contains(parent)) {
- parents.add(parent);
- if (chargedList.contains(parent)) chargedList.remove(parent);
- if (neutralList.contains(parent)) neutralList.remove(parent);
- }
- }
- for (CheatCluster cluster : skippedClusters) {
- MCParticle mcp = cluster.getMCParticle();
- MCParticle parent = primaryMCParents.get(mcp);
- if (parent==null) continue;
- if (!charged.containsKey(parent) && !neutrals.containsKey(parent) &&
- !parents.contains(parent)) {
- parents.add(parent);
- if (chargedList.contains(parent)) chargedList.remove(parent);
- if (neutralList.contains(parent)) neutralList.remove(parent);
- }
- }
- //if (firstEvents) System.out.println(" Total of "+parents.size()+" interacting parents.");
- text = " Total of "+parents.size()+" interacting parents:";
- for (MCParticle p : parents) text+=" "+map.get(p);
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
+ // Find missing parents of decay or interaction tracks and clusters.
+ if (firstEvents) System.out.println(" Find missing parents of decay or interaction tracks and clusters.");
+ List<MCParticle> missingParents = new ArrayList<MCParticle>();
+
+ if (!allowDecays) {
+ // Add missing final state parents.
+ for (MCParticle mcp : decayedParticles) {
+ if (!checkForParent(mcp) && !missingParents.contains(mcp)) {
+ missingParents.add(mcp);
+ }
+ }
+ List<MCParticle> missingDecayParents = new ArrayList<MCParticle>();
+ // Add remaining missing decay parents.
+ for (CheatTrack track : skippedTracks) {
+ MCParticle mcp = track.getMCParticle();
+ MCParticle parent = primaryMCParents.get(mcp);
+ if (checkDecay(parent)) missingDecayParents.add(parent);
+ }
+ for (CheatCluster cluster : skippedClusters) {
+ MCParticle mcp = cluster.getMCParticle();
+ MCParticle parent = primaryMCParents.get(mcp);
+ if (checkDecay(parent)) missingDecayParents.add(parent);
+ }
+ if (firstEvents) {
+ text = " Total of "+missingDecayParents.size()+" missing decay parents:";
+ for (MCParticle p : missingDecayParents) text+=" "+map.get(p);
+ System.out.println(text); System.err.println(text);
+ }
+ for (MCParticle parent : missingDecayParents) {
+ // Check if parent is already recorded.
+ if (!checkForParent(parent) && !missingParents.contains(parent)) {
+ missingParents.add(parent);
+ }
+ }
+ }
+ if (!allowNuclearInteractions) {
+ List<MCParticle> missingInteractionParents = new ArrayList<MCParticle>();
+ for (MCParticle mcp : interactedParticles) {
+ MCParticle parent = primaryMCParents.get(mcp);
+ if (parent==null) { parent = findPrimaryParent(mcp); primaryMCParents.put(mcp,parent); }
+ if (!checkForParent(mcp) && !missingInteractionParents.contains(mcp)) {
+ missingInteractionParents.add(mcp);
+ }
+ }
+ if (firstEvents) {
+ text = " Total of "+missingInteractionParents.size()+" missing interaction parent(s):";
+ for (MCParticle p : missingInteractionParents) text+=" "+map.get(p);
+ System.out.println(text); if (missingInteractionParents.size()>0) System.err.println(text);
+ }
+ // Add missing interacted particle only if its parent(s) are not already included.
+ for (MCParticle mcp : missingInteractionParents) {
+ if (firstEvents) System.err.println(" MCParticle = "+map.get(mcp));
+ MCParticle primaryParent = primaryMCParents.get(mcp);
+ if (firstEvents) System.err.println(" PrimaryParent = "+map.get(primaryParent));
+ if (missingInteractionParents.contains(primaryParent) || missingParents.contains(primaryParent)) continue;
+
+ List<MCParticle> parents = mcp.getParents();
+ MCParticle parent = parents.get(0);
+ int i=0; boolean found = false;
+ while (parent!=primaryParent && i<10 && !found) {
+ if (firstEvents) System.err.println(" Parent = "+map.get(parent)); i++;
+ if (missingInteractionParents.contains(parent) || missingParents.contains(parent)) found = true;
+ parents = parent.getParents(); parent = parents.get(0);
+ }
+ if (found) continue;
+ if (!missingParents.contains(mcp))
+ missingParents.add(mcp);
+ }
+ }
+ if (firstEvents) {
+ text = " Total of "+missingParents.size()+" unaccounted for decay or interacting parent(s):";
+ for (MCParticle p : missingParents) text+=" "+map.get(p);
+ System.out.println(text); if (missingParents.size()>0) System.err.println(text);
+ }
// Select particles:
if (firstEvents) System.out.println(" Select particles to be reconstructed.");
+ finalEnergy = 0;
- // Add charged particles that were found in TrackingCheater list.
+ // Add charged particles that were found in Track list.
for (MCParticle mcp : charged.keySet()) {
CheatTrack track = charged.get(mcp);
+ for (ReconstructedParticle particle : particles) {
+ MCParticle p = ((CheatReconstructedParticle)particle).getMCParticle();
+ if (map.get(p)==map.get(mcp)) System.err.println(" Particle #"+map.get(mcp)+" already in list");
+ }
if (useFullTruth) particles.add(new CheatReconstructedParticle(mcp));
else particles.add(new CheatReconstructedParticle(track,(Particle)mcp));
+ finalEnergy += mcp.getEnergy();
}
- // Add neutral particles that were found in ClusterCheater list.
+ // Add neutral particles that were found in Cluster list.
for (MCParticle mcp : neutrals.keySet()) {
CheatCluster cluster = neutrals.get(mcp);
+ for (ReconstructedParticle particle : particles) {
+ MCParticle p = ((CheatReconstructedParticle)particle).getMCParticle();
+ if (map.get(p)==map.get(mcp)) System.err.println(" Particle #"+map.get(mcp)+" already in list");
+ }
if (useFullTruth) particles.add(new CheatReconstructedParticle(mcp));
else particles.add(new CheatReconstructedParticle(cluster,getMeasuredEnergy(cluster),(Particle)mcp));
+ finalEnergy += mcp.getEnergy();
}
// Charged and neutral lists contain MCParticles that interacted or decayed.
if (useFullTruth) {
// Add interacting charged particles.
for (MCParticle mcp : chargedList) {
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) continue;
particles.add(new CheatReconstructedParticle(mcp));
finalEnergy += mcp.getEnergy();
}
// Add interacting neutral particles.
for (MCParticle mcp : neutralList) {
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) continue;
particles.add(new CheatReconstructedParticle(mcp));
finalEnergy += mcp.getEnergy();
}
}
else {
- // Add interacting parents.
- for (MCParticle mcp : parents) {
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) continue;
+ // Add interacting or decaying parents.
+ for (MCParticle mcp : missingParents) {
+ MCParticle primaryParent = findPrimaryParent(mcp);
+ primaryMCParents.put(mcp,primaryParent);
+ boolean decayProduct = checkDecay(primaryParent);
+ if (getStatusCode(mcp).equals("GEN_PREDECAY")) {
+ if (allowDecays && mcp.getCharge()==0.) continue; // Skip decaying neutrals, e.g. K0_s.
+ if (firstEvents)
+ System.out.println(" Adding decayed MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+".");
+ } else if (decayProduct) {
+ if (!allowDecays) continue;
+ if (firstEvents)
+ System.out.println(" Adding decay product MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+".");
+ } else {
+ if (allowNuclearInteractions) continue;
+ if (firstEvents)
+ System.out.println(" Adding interacted MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+
+ " whose primary parent is #"+map.get(primaryMCParents.get(mcp))+".");
+ }
particles.add(new CheatReconstructedParticle(mcp));
finalEnergy += mcp.getEnergy();
- if (firstEvents)
- System.out.println(" Adding interacting MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+".");
}
// Add extra tracks.
for (CheatTrack track : extraTrackList) {
MCParticle mcp = track.getMCParticle();
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) continue;
+ boolean found = false;
+ for (ReconstructedParticle particle : particles) {
+ MCParticle crp = ((CheatReconstructedParticle)particle).getMCParticle();
+ if (firstEvents) {
+ if (crp==mcp) System.err.println(" Particle #"+map.get(mcp)+" already in list");
+ if (!acceptNuclearInteractionProducts) {
+ //if (crp==primaryMCParents.get(mcp))
+ //System.err.println(" Particle #"+map.get(mcp)+"'s parent is already in list");
+ if (primaryMCParents.get(crp)==mcp)
+ System.err.println(" Particle #"+map.get(mcp)+"'s daughter is already in list");
+ }
+ }
+ if (crp==mcp || (!acceptNuclearInteractionProducts && crp==primaryMCParents.get(mcp))) found = true;
+ }
+ if (found) continue;
particles.add(new CheatReconstructedParticle(track,(Particle)mcp));
finalEnergy += mcp.getEnergy();
}
// Add extra clusters.
for (CheatCluster cluster : extraClusterList) {
MCParticle mcp = cluster.getMCParticle();
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) continue;
+ boolean found = false;
+ for (ReconstructedParticle particle : particles) {
+ MCParticle crp = ((CheatReconstructedParticle)particle).getMCParticle();
+ if (firstEvents) {
+ if (crp==mcp) System.err.println(" Particle #"+map.get(mcp)+" already in list");
+ if (!acceptNuclearInteractionNeutrals) {
+ //if (crp==primaryMCParents.get(mcp))
+ //System.err.println(" Particle #"+map.get(mcp)+"'s parent is already in list");
+ if (primaryMCParents.get(crp)==mcp)
+ System.err.println(" Particle #"+map.get(mcp)+"'s daughter is already in list");
+ }
+ }
+ if (crp==mcp || (!acceptNuclearInteractionNeutrals && crp==primaryMCParents.get(mcp))) found = true;
+ }
+ if (found) continue;
particles.add(new CheatReconstructedParticle(cluster,getMeasuredEnergy(cluster),(Particle)mcp));
finalEnergy += mcp.getEnergy();
}
}
- // Check particles.
- text = " Total of "+particles.size()+" reconstructed particles:";
- for (ReconstructedParticle particle : particles) {
- MCParticle mcp = ((CheatReconstructedParticle)particle).getMCParticle();
- text+=" "+map.get(mcp);
- }
- if (firstEvents) System.out.println(text); //if (first) System.err.println(text);
+ // Check reconstructed particles.
+ checkReconstructedParticles(particles,chargedList,neutralList);
- // Total energy
- if (firstEvents) {
- text = " Total event energy = "+df.format(finalEnergy);
- System.out.println(text+".\n"); //System.err.println(text);
- }
-
- if (hist) analyzeReconstructedParticles(particles);
+ if (hist) analyzeReconstructedParticles(particles,chargedList,neutralList);
return particles;
}
- void analyzeReconstructedParticles(List<ReconstructedParticle> particles)
+ void analyzeReconstructedParticles(List<ReconstructedParticle> particles, List<MCParticle> chargedList,
+ List<MCParticle> neutralList)
{
AIDA aida = AIDA.defaultInstance();
+ ITree tree = aida.tree();
+ tree.cd("ReconCheater");
aida.cloud1D("# particles").fill(particles.size());
aida.cloud1D("Final Energy").fill(finalEnergy);
+ aida.cloud1D("Energy diff").fill(finalEnergy-totalEventEnergy);
+ // Reconstructed neutrals.
for (MCParticle mcp : neutrals.keySet()) {
+ double e = mcp.getEnergy();
int PDGID = mcp.getPDGID();
CheatCluster cluster = neutrals.get(mcp);
int nhits = 0;
for (CalorimeterHit hit : cluster.getCalorimeterHits())
if (hit.getTime()<100.) nhits++;
- double ratio = nhits/mcp.getEnergy();
+ double ratio = nhits/e;
if (PDGID==22) {
ratio = cluster.getRawEnergy()/mcp.getEnergy();
aida.cloud1D("photon ratio").fill(ratio);
- aida.cloud2D("photon Energy vs. ratio").fill(ratio,mcp.getEnergy());
+ aida.cloud2D("photon Energy vs. ratio").fill(ratio,e);
}
else if (PDGID==130) {
aida.cloud1D("Klong ratio").fill(ratio);
- aida.cloud2D("Klong Energy vs. ratio").fill(ratio,mcp.getEnergy());
+ aida.cloud2D("Klong Energy vs. ratio").fill(ratio,e);
+ if (e<1.0) aida.cloud1D("Klong low Energy").fill(e);
}
else if (Math.abs(PDGID)==2112) {
aida.cloud1D("neutron ratio").fill(ratio);
- aida.cloud2D("neutron Energy vs. ratio").fill(ratio,mcp.getEnergy());
+ aida.cloud2D("neutron Energy vs. ratio").fill(ratio,e);
+ if (e<1.0) aida.cloud1D("neutron low Energy").fill(e);
+ }
+ }
+
+ // Reconstructed particles.
+ int nProduced = 0;
+ for (ReconstructedParticle particle : particles) {
+ CheatReconstructedParticle crp = (CheatReconstructedParticle)particle;
+ MCParticle mcp = crp.getMCParticle();
+ String statusCode = getStatusCode(mcp); int PDGID = mcp.getPDGID();
+ if (statusCode.equals("GEN_FINAL_STATE")) {
+ if (Math.abs(PDGID)==211) nProduced++;
+ aida.cloud1D("Final: PDGID").fill(PDGID);
+ //if (firstEvents && Math.abs(PDGID)>2000)
+ // System.err.println(" ReconstructedParticle #"+map.get(mcp)+" ("+mcp.getType().getName()+
+ // ") status = "+statusCode+" parent = "+primaryMCParents.get(mcp).getType().getName());
+ if (Math.abs(PDGID)>5000)
+ System.err.println(" ReconstructedParticle #"+map.get(mcp)+" ("+mcp.getType().getName()+
+ ") status = "+statusCode+" parent = "+primaryMCParents.get(mcp).getType().getName());
+
+ }
+ else
+ if (statusCode.equals("GEN_PREDECAY")) {
+ aida.cloud1D("Decay: PDGID").fill(PDGID);
+ }
+ }
+
+ // Interacting or lost charged particles.
+ int nInteracted = 0;
+ for (MCParticle mcp : chargedList) {
+ MCParticle primaryParent = findPrimaryParent(mcp);
+ if (primaryParent==null) System.err.println(" Null primary parent for particle #"+map.get(mcp));
+ else if (checkDecay(mcp) || checkDecay(primaryParent)) continue;
+ double pt = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY());
+ double p = Math.sqrt(pt*pt+mcp.getPZ()*mcp.getPZ());
+ if (p<pTrackMin) continue;
+ int PDGID = mcp.getPDGID();
+ boolean found = false;
+ for (ReconstructedParticle particle : particles) {
+ MCParticle crp = ((CheatReconstructedParticle)particle).getMCParticle();
+ if (crp==mcp || crp==primaryMCParents.get(mcp)) found = true;
+ }
+ if (!found) {
+ if (!checkInteraction(mcp)) {
+ tree.cd("Lost");
+ double cosTh = mcp.getPZ()/p;
+ aida.cloud1D("Charged: cosTh").fill(cosTh);
+ tree.cd("..");
+ } else {
+ tree.cd("Interacted");
+ double e = mcp.getEnergy();
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double r = Math.sqrt(endPoint.x()*endPoint.x()+endPoint.y()*endPoint.y());
+ aida.cloud1D("Charged: radius").fill(r);
+ List<MCParticle> daughters = mcp.getDaughters();
+ double eInt = getDaughterEnergy(daughters);
+ int nCharged = 0;
+ for (MCParticle daughter : daughters) {
+ if (daughter.getPDGID()==0) continue;
+ if (daughter.getCharge()!=0) nCharged++;
+ }
+ aida.cloud1D("Charged: PDGID").fill(PDGID);
+ if (Math.abs(PDGID)==211) nInteracted++;
+ aida.cloud1D("Charged: E ratio").fill(eInt/e);
+ aida.cloud1D("Charged: # charged").fill(nCharged);
+ //aida.cloud1D("Charged: ").fill();
+ tree.cd("..");
+ }
+ }
+ }
+
+ // Interacting or lost neutral particles.
+ for (MCParticle mcp : neutralList) {
+ MCParticle primaryParent = findPrimaryParent(mcp);
+ if (primaryParent==null) System.err.println(" Null primary parent for particle #"+map.get(mcp));
+ else if (checkDecay(mcp) || checkDecay(primaryParent)) continue;
+ double e = mcp.getEnergy();
+ if (e<EClusterMin) continue;
+ int PDGID = mcp.getPDGID();
+ boolean found = false;
+ for (ReconstructedParticle particle : particles) {
+ MCParticle crp = ((CheatReconstructedParticle)particle).getMCParticle();
+ if (crp==mcp || crp==primaryMCParents.get(mcp)) found = true;
+ }
+ if (!found) {
+ if (!checkInteraction(mcp)) {
+ tree.cd("Lost");
+ double pt = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY());
+ double p = Math.sqrt(pt*pt+mcp.getPZ()*mcp.getPZ());
+ double cosTh = mcp.getPZ()/p;
+ aida.cloud1D("Neutral: cosTh").fill(cosTh);
+ tree.cd("..");
+ } else
+ // Photon conversions.
+ if (PDGID==22) {
+ tree.cd("Conversion");
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double r = Math.sqrt(endPoint.x()*endPoint.x()+endPoint.y()*endPoint.y());
+ aida.cloud1D("radius").fill(r);
+ tree.cd("..");
+ } else {
+ tree.cd("Interacted");
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double r = Math.sqrt(endPoint.x()*endPoint.x()+endPoint.y()*endPoint.y());
+ aida.cloud1D("Neutral: radius").fill(r);
+ List<MCParticle> daughters = mcp.getDaughters();
+ double eInt = getDaughterEnergy(daughters);
+ int nCharged = 0;
+ for (MCParticle daughter : daughters) {
+ if (daughter.getPDGID()==0) continue;
+ if (daughter.getCharge()!=0) nCharged++;
+ }
+ aida.cloud1D("Neutral: E ratio").fill(eInt/e);
+ aida.cloud1D("Neutral: # charged").fill(nCharged);
+ //aida.cloud1D("Neutral: ").fill();
+ tree.cd("..");
+ }
}
}
+ nProduced += nInteracted;
+ aida.cloud1D("# pi produced").fill(nProduced);
+ aida.cloud1D("# pi interacted").fill(nInteracted);
+ if (nProduced>0) {
+ double ratio = nInteracted / nProduced;
+ aida.cloud1D("int. ratio %").fill(ratio*100.);
+ }
+
+ // Interacting particles.
+ int nPhotonConversions = 0, nNuclearInteractions = 0;
+ for (MCParticle mcp : interactedParticles) {
+ double e = mcp.getEnergy();
+ if (e<EClusterMin) continue;
+ int PDGID = mcp.getPDGID();
+ // Photon conversions.
+ if (PDGID==22) {
+ nPhotonConversions++;
+ }
+ else {
+ nNuclearInteractions++;
+ }
+ }
+ aida.cloud1D("# PhotonConversions").fill(nPhotonConversions);
+ aida.cloud1D("# NuclearInteractions").fill(nNuclearInteractions);
+ tree.cd("..");
+ }
+ /** Gets the total energy of a particle's daughters. */
+ double getDaughterEnergy(MCParticle mcp)
+ {
+ return getDaughterEnergy(mcp.getDaughters());
+ }
[truncated at 1000 lines; 440 more skipped]