lcsim/src/org/lcsim/recon/cheater
diff -u -r1.5 -r1.6
--- ReconCheater.java 18 Oct 2005 17:09:50 -0000 1.5
+++ ReconCheater.java 26 Oct 2005 19:04:03 -0000 1.6
@@ -67,6 +67,7 @@
boolean acceptNuclearInteractionNeutrals;
boolean allowRadiation;
+ double DecayDistance;
double ECalResolution, HCalResolution;
double ECalSampling, HCalSampling, HCalDigital; // Hits/GeV
double pTrackMin, EClusterMin, NDigitalMin;
@@ -123,6 +124,7 @@
acceptNuclearInteractionNeutrals = cheatingTable.acceptNuclearInteractionNeutrals();
allowRadiation = cheatingTable.allowRadiation();
+ DecayDistance = cheatingTable.getDecayDistance();
ECalResolution = cheatingTable.getECalResolution();
ECalSampling = cheatingTable.getECalSampling();
HCalResolution = cheatingTable.getHCalResolution();
@@ -271,20 +273,12 @@
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.
+ // 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) ||
+ (!allowDecays && decayProduct && isLongLived(primaryParent)) ||
(!allowNuclearInteractions && interactionProduct)) {
skippedTracks.add(track);
continue;
@@ -327,19 +321,35 @@
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:";
+ text = " Total of "+extraTrackList.size()+" tracks from predecays, 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);
}
+ // Remove extra primary (PRE-DECAY) tracks.
+ List<CheatTrack> primaryTracks = new ArrayList<CheatTrack>();
+ for (CheatTrack track : extraTrackList) {
+ MCParticle mcp = track.getMCParticle();
+ if (getStatusCode(mcp).equals("GEN_PREDECAY") && !isLongLived(mcp)) {
+ primaryTracks.add(track); skippedTracks.add(track);
+ }
+ }
+ for (CheatTrack track : primaryTracks) extraTrackList.remove(track);
+ if (firstEvents && primaryTracks.size()>0) {
+ text = " Removed "+primaryTracks.size()+" predecay tracks:";
+ for (CheatTrack t : primaryTracks) { 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);
}
+ */
// Look for clusters in ClusterCheater list from neutral MCParticles.
if (firstEvents) {
@@ -391,12 +401,12 @@
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 (!found && !checkForParent(getParent(mcp))) unAssociated = true;
}
if (!unAssociated) continue;
// Skip decay and interaction products if switched off.
- if ((!allowDecays && decayProduct) ||
+ if ((!allowDecays && decayProduct && isLongLived(primaryParent)) ||
(!allowNuclearInteractions && interactionProduct)) {
skippedClusters.add(cluster);
continue;
@@ -412,10 +422,9 @@
}
// or add as an extra (unexpected) cluster.
else
- if ((!decayProduct || acceptDecayNeutrals) &&
- (!interactionProduct || acceptNuclearInteractionNeutrals)) {
+ //if ((!decayProduct || acceptDecayNeutrals) &&
+ //(!interactionProduct || acceptNuclearInteractionNeutrals)) {
extraClusterList.add(cluster);
- }
}
}
}
@@ -433,6 +442,7 @@
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);
@@ -440,6 +450,8 @@
else lostParticles.add(mcp);
}
if (firstEvents) {
+ text = " Record interacting, decaying or lost particles.";
+ System.out.println(text); //if (first) System.err.println(text);
if (decayedParticles.size()>0) {
text = " Total of "+decayedParticles.size()+" decayed particles:";
for (MCParticle p : decayedParticles) text+=" "+map.get(p);
@@ -456,83 +468,75 @@
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>();
+ List<MCParticle> missingDecayParents = new ArrayList<MCParticle>();
+ List<MCParticle> missingInteractionParents = 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);
- }
- }
+ // Check for missing decay and interation parents.
+ for (CheatTrack track : skippedTracks) {
+ MCParticle mcp = track.getMCParticle();
+ if (getStatusCode(mcp).equals("GEN_PREDECAY")) continue;
+ MCParticle parent = primaryMCParents.get(mcp);
+ if (checkForParent(parent)) continue;
+ if (checkDecay(parent) && !missingDecayParents.contains(parent)) missingDecayParents.add(parent);
+ else if (checkInteraction(parent) && !missingInteractionParents.contains(parent)) missingInteractionParents.add(parent);
+ }
+ for (CheatCluster cluster : skippedClusters) {
+ MCParticle mcp = cluster.getMCParticle();
+ MCParticle parent = primaryMCParents.get(mcp);
+ if (checkForParent(parent)) continue;
+ if (checkDecay(parent) && !missingDecayParents.contains(parent)) missingDecayParents.add(parent);
+ else if (checkInteraction(parent) && !missingInteractionParents.contains(parent)) missingInteractionParents.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;
+ if (firstEvents) {
+ text = " Total of "+missingDecayParents.size()+" missing decay parents:";
+ for (MCParticle p : missingDecayParents) text+=" "+map.get(p);
+ System.out.println(text); if (missingDecayParents.size()>0) System.err.println(text);
+ 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);
+ }
+
+ // Pick up missing decay and interaction parents.
+ for (MCParticle parent : missingDecayParents) {
+ MCParticle primaryParent = findPrimaryParent(parent); primaryMCParents.put(parent,primaryParent);
+ if (nEvt==89) System.err.println(" Missing parent #"+map.get(parent)+" primary parent #"+map.get(primaryParent));
+ if (missingDecayParents.contains(primaryParent)) continue;
+ if (!decayedParticles.contains(parent)) decayedParticles.add(parent);
+ if (!allowDecays) {
+ if (!missingParents.contains(parent)) missingParents.add(parent);
+ }
+ }
+ // Add missing interacted particle only if its parent is not already included.
+ for (MCParticle parent : missingInteractionParents) {
+ MCParticle primaryParent = findPrimaryParent(parent); primaryMCParents.put(parent,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);
+ List<MCParticle> grandParents = parent.getParents();
+ MCParticle grandParent = grandParents.get(0);
+ int i=0; boolean found = false;
+ while (grandParent!=primaryParent && i<10 && !found) {
+ if (missingInteractionParents.contains(grandParent) || missingParents.contains(grandParent)) found = true;
+ grandParents = grandParent.getParents(); grandParent = grandParents.get(0); i++;
+ }
+ if (found) continue;
+
+ if (!interactedParticles.contains(parent)) interactedParticles.add(parent);
+ if (!allowNuclearInteractions) {
+ if (!missingParents.contains(parent)) missingParents.add(parent);
}
}
-
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;
@@ -576,22 +580,30 @@
else {
// Add interacting or decaying parents.
for (MCParticle mcp : missingParents) {
- MCParticle primaryParent = findPrimaryParent(mcp);
- primaryMCParents.put(mcp,primaryParent);
+ MCParticle primaryParent = primaryMCParents.get(mcp);
+ boolean decayed = checkDecay(mcp);
+ boolean interacted = false;
+ if (!decayed) interacted = checkInteraction(mcp);
boolean decayProduct = checkDecay(primaryParent);
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) {
- if (allowDecays && mcp.getCharge()==0.) continue; // Skip decaying neutrals, e.g. K0_s.
+ if (decayed) {
+ //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 (interacted) {
+ //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))+".");
} else if (decayProduct) {
- if (!allowDecays) continue;
+ //if (!allowDecays) continue;
if (firstEvents)
System.out.println(" Adding decay product MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+".");
} else {
- if (allowNuclearInteractions) continue;
+ //if (allowNuclearInteractions) continue;
if (firstEvents)
- System.out.println(" Adding interacted MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+
+ System.out.println(" Adding MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+
" whose primary parent is #"+map.get(primaryMCParents.get(mcp))+".");
+ System.err.println(" Adding MCParticle #"+map.get(mcp)+" "+mcp.getPDGID()+".");
}
particles.add(new CheatReconstructedParticle(mcp));
finalEnergy += mcp.getEnergy();
@@ -665,24 +677,28 @@
double e = mcp.getEnergy();
int PDGID = mcp.getPDGID();
CheatCluster cluster = neutrals.get(mcp);
+ double eRatio = cluster.getRawEnergy()/mcp.getEnergy();
int nhits = 0;
for (CalorimeterHit hit : cluster.getCalorimeterHits())
if (hit.getTime()<100.) nhits++;
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,e);
+ aida.cloud1D("photon ratio").fill(eRatio);
+ aida.cloud2D("photon Energy vs. ratio").fill(eRatio,e);
}
else if (PDGID==130) {
- aida.cloud1D("Klong ratio").fill(ratio);
- aida.cloud2D("Klong Energy vs. ratio").fill(ratio,e);
+ aida.cloud1D("Klong E ratio").fill(eRatio);
+ aida.cloud2D("Klong Energy vs. E ratio").fill(eRatio,e);
+ aida.cloud1D("Klong Hit ratio").fill(ratio);
+ aida.cloud2D("Klong Energy vs. Hit 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,e);
+ aida.cloud1D("neutron E ratio").fill(eRatio);
+ aida.cloud2D("neutron Energy vs. E ratio").fill(eRatio,e);
+ aida.cloud1D("neutron Hit ratio").fill(ratio);
+ aida.cloud2D("neutron Energy vs. Hit ratio").fill(ratio,e);
if (e<1.0) aida.cloud1D("neutron low Energy").fill(e);
}
}
@@ -692,10 +708,15 @@
for (ReconstructedParticle particle : particles) {
CheatReconstructedParticle crp = (CheatReconstructedParticle)particle;
MCParticle mcp = crp.getMCParticle();
+ double e = mcp.getEnergy();
+ if (e<EClusterMin) continue;
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 (checkDecay(mcp)) {
+ aida.cloud1D("Final: Decay").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());
@@ -711,11 +732,18 @@
}
// 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;
+ //MCParticle primaryParent = findPrimaryParent(mcp); primaryMCParents.put(mcp,primaryParent);
+ //if (primaryParent==null) System.err.println(" Null primary parent for particle #"+map.get(mcp));
+ //else if (checkDecay(mcp) || checkDecay(primaryParent)) continue;
+ if (checkDecay(mcp)) {
+ if (!decayedParticles.contains(mcp)) decayedParticles.add(mcp);
+ continue;
+ }
+ boolean interacted = checkInteraction(mcp);
+ if (interacted) { if (!interactedParticles.contains(mcp)) interactedParticles.add(mcp); }
+ else lostParticles.add(mcp);
+
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;
@@ -723,10 +751,11 @@
boolean found = false;
for (ReconstructedParticle particle : particles) {
MCParticle crp = ((CheatReconstructedParticle)particle).getMCParticle();
- if (crp==mcp || crp==primaryMCParents.get(mcp)) found = true;
+ //if (crp==mcp || crp==primaryMCParents.get(mcp)) found = true;
+ if (crp==mcp) found = true;
}
if (!found) {
- if (!checkInteraction(mcp)) {
+ if (!interacted) {
tree.cd("Lost");
double cosTh = mcp.getPZ()/p;
aida.cloud1D("Charged: cosTh").fill(cosTh);
@@ -745,7 +774,6 @@
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();
@@ -756,19 +784,28 @@
// 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;
+ //MCParticle primaryParent = findPrimaryParent(mcp); primaryMCParents.put(mcp,primaryParent);
+ //if (primaryParent==null) System.err.println(" Null primary parent for particle #"+map.get(mcp));
+ //else if (checkDecay(mcp) || checkDecay(primaryParent)) continue;
+ if (checkDecay(mcp)) {
+ if (!decayedParticles.contains(mcp)) decayedParticles.add(mcp);
+ continue;
+ }
+ boolean interacted = checkInteraction(mcp);
+ if (interacted) { if (!interactedParticles.contains(mcp)) interactedParticles.add(mcp); }
+ else lostParticles.add(mcp);
+
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 (crp==mcp || crp==primaryMCParents.get(mcp)) found = true;
+ if (crp==mcp) found = true;
}
if (!found) {
- if (!checkInteraction(mcp)) {
+ if (!interacted) {
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());
@@ -802,28 +839,44 @@
}
}
}
- 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.);
+
+ // Decaying particles.
+ for (MCParticle mcp : decayedParticles) {
+ double e = mcp.getEnergy();
+ if (e<EClusterMin) continue;
+ int PDGID = mcp.getPDGID();
+
+ Hep3Vector endPoint = mcp.getEndPoint();
+ double distance = Math.sqrt(endPoint.x()*endPoint.x()+endPoint.y()*endPoint.y()+endPoint.z()*endPoint.z());
+ aida.cloud1D("Decayed: distance").fill(distance);
+ if (distance>=DecayDistance) aida.cloud1D("Decayed: LongLived").fill(PDGID);
}
// Interacting particles.
- int nPhotonConversions = 0, nNuclearInteractions = 0;
+ int nInteracted = 0, 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++;
+ String statusCode = getStatusCode(mcp); int PDGID = mcp.getPDGID();
+ if (statusCode.equals("GEN_FINAL_STATE")) {
+ if (Math.abs(PDGID)==211) nInteracted++;
+ // Photon conversions.
+ if (PDGID==22) {
+ nPhotonConversions++;
+ }
+ else {
+ nNuclearInteractions++;
+ }
}
}
+ if (allowNuclearInteractions) nProduced += nInteracted;
+ aida.cloud1D("# pi Produced").fill(nProduced);
+ aida.cloud1D("# pi Interacted").fill(nInteracted);
+ if (nProduced>0) {
+ double ratio = (double) nInteracted * 100./ nProduced;
+ if (ratio>100.) System.err.println(" Evt. #"+nEvt+" Int. ratio = "+df.format(ratio));
+ aida.cloud1D("int. ratio %").fill(ratio);
+ }
aida.cloud1D("# PhotonConversions").fill(nPhotonConversions);
aida.cloud1D("# NuclearInteractions").fill(nNuclearInteractions);
tree.cd("..");
@@ -850,12 +903,13 @@
/** Checks if a given particle decayed. */
boolean checkDecay(MCParticle mcp)
{
- if (getStatusCode(mcp).equals("GEN_PREDECAY")) return true;;
boolean decayed = false;
int PDGID = mcp.getPDGID();
boolean pi0 = PDGID==111;
if (PDGID==0 || pi0 || PDGID==92) return false;
+ if (getStatusCode(mcp).equals("GEN_PREDECAY")) return true;
+
double mass2 = 0.;
List<MCParticle> daughters = mcp.getDaughters();
if (daughters.size()<2) return false;
@@ -900,6 +954,21 @@
}
return hasNeutralDaughters;
}
+ boolean isLongLived(MCParticle mcp)
+ {
+ boolean longLived = false;
+ double distance = 0.;
+ try {
+ Hep3Vector endPoint = mcp.getEndPoint();
+ Hep3Vector origin = mcp.getOrigin();
+ double dx = endPoint.x() - origin.x(), dy = endPoint.y() - origin.y(), dz = endPoint.z() - origin.z();
+ distance = Math.sqrt(dx*dx+dy*dy+dz*dz);
+ }
+ catch (java.lang.RuntimeException ex) { }
+ if (distance>=DecayDistance) longLived = true;
+ return longLived;
+ }
+
/** Checks if a given particle interacted. */
boolean checkInteraction(MCParticle mcp)
{
@@ -964,46 +1033,55 @@
{
boolean inList = false;
if (charged.containsKey(parent) || neutrals.containsKey(parent)) inList = true;
- for (CheatTrack t : extraTrackList) { MCParticle p = t.getMCParticle(); if (map.get(p)==map.get(parent)) inList = true; }
- for (CheatCluster c : extraClusterList) { MCParticle p = c.getMCParticle(); if (map.get(p)==map.get(parent)) inList = true; }
+ for (CheatTrack t : extraTrackList) { MCParticle p = t.getMCParticle(); if (p==parent) inList = true; }
+ for (CheatCluster c : extraClusterList) { MCParticle p = c.getMCParticle(); if (p==parent) inList = true; }
+ if (getStatusCode(parent).equals("GEN_FINAL_STATE")) return inList;
- MCParticle grandParent = null;
+ MCParticle grandParent = null, primaryParent = primaryMCParents.get(parent);
List<MCParticle> grandParents = parent.getParents();
while (!grandParents.isEmpty() && !inList) {
grandParent = grandParents.get(0);
if (charged.containsKey(grandParent) || neutrals.containsKey(grandParent)) inList = true;
- for (CheatTrack t : extraTrackList) { MCParticle p = t.getMCParticle(); if (map.get(p)==map.get(grandParent)) inList = true; }
- for (CheatCluster c : extraClusterList) { MCParticle p = c.getMCParticle(); if (map.get(p)==map.get(grandParent)) inList = true; }
+ for (CheatTrack t : extraTrackList) { MCParticle p = t.getMCParticle(); if (p==grandParent) inList = true; }
+ for (CheatCluster c : extraClusterList) { MCParticle p = c.getMCParticle(); if (p==grandParent) inList = true; }
+ if (grandParent==primaryParent || getStatusCode(grandParent).equals("GEN_FINAL_STATE")) return inList;
grandParents = grandParent.getParents();
}
return inList;
}
- MCParticle findPrimaryParent(MCParticle mcp)
+ MCParticle getParent(MCParticle mcp)
{
- String statusCode = getStatusCode(mcp);
- boolean finalState = statusCode.equals("GEN_FINAL_STATE");
-
MCParticle parent = null;
List<MCParticle> parents = mcp.getParents();
if (!parents.isEmpty()) parent = parents.get(0);
- if (parent==null) {
- if (firstEvents) System.err.println(" Particle #"+map.get(mcp)+" ("+mcp.getPDGID()+") has no parent!");
- return parent;
- }
- boolean parentFinalState = getStatusCode(parent).equals("GEN_FINAL_STATE");
- int parentPDGID = parent.getPDGID();
-
- // Check for final state particles.
- if (!finalState) {
- // Get final state parent particle.
- List<MCParticle> grandParents = parent.getParents();
- while (!parentFinalState && !grandParents.isEmpty() && !getStatusCode(parent).equals("GEN_INITIAL")) {
- parent = grandParents.get(0); grandParents = parent.getParents();
- parentFinalState = getStatusCode(parent).equals("GEN_FINAL_STATE");
- parentPDGID = parent.getPDGID();
+ if (parent==null) System.err.println(" Evt.#"+nEvt+" Particle #"+map.get(mcp)+" ("+mcp.getPDGID()+") has no parent!");
+ return parent;
+ }
+ MCParticle findPrimaryParent(MCParticle mcp)
+ {
+ MCParticle parent = getParent(mcp);
+ if (parent==null) return parent;
+
+ String statusCode = getStatusCode(parent);
+ if (statusCode.equals("GEN_INITIAL")) return parent;
+
+ // Get initial state parent particle.
+ List<MCParticle> grandParents = parent.getParents();
+ MCParticle grandParent = grandParents.get(0);
+
+ MCParticle primaryParent = null;
+ while (!statusCode.equals("GEN_INITIAL")) {
+ if (statusCode.equals("GEN_FINAL_STATE")) {
+ if (allowDecays || !getStatusCode(grandParent).equals("GEN_PREDECAY")) return parent;
+ primaryParent = parent;
}
+ if (checkDecay(parent) && isLongLived(parent)) primaryParent = parent;
+ if (statusCode.equals("GEN_PREDECAY") && getStatusCode(grandParent).equals("GEN_INITIAL")) break;
+ parent = grandParent; grandParents = parent.getParents(); grandParent = grandParents.get(0);
+ statusCode = getStatusCode(parent);
}
+ if (primaryParent!=null) parent = primaryParent;
return parent;
}
@@ -1156,29 +1234,30 @@
for (ReconstructedParticle particle : particles) {
MCParticle mcp = ((CheatReconstructedParticle)particle).getMCParticle();
String statusCode = getStatusCode(mcp);
- if (!statusCode.equals("GEN_FINAL_STATE"))
+ if (!statusCode.equals("GEN_FINAL_STATE")) {
System.out.print(" ReconstructedParticle #"+map.get(mcp)+
" ("+mcp.getType().getName()+") status = "+statusCode);
- if (statusCode.equals("GEN_PREDECAY")) {
- double e = 0.;
- String list = "";
- List<MCParticle> daughters = mcp.getDaughters();
- for (MCParticle daughter : daughters) {
- list+=" "+map.get(daughter);
- if (getStatusCode(daughter).equals("GEN_FINAL_STATE")) e+= daughter.getEnergy();
- else {
- list+=" (";
- List<MCParticle> grandDaughters = daughter.getDaughters();
- for (MCParticle grandDaughter : grandDaughters) {
- list+=" "+map.get(grandDaughter);
- e+= grandDaughter.getEnergy();
+ if (statusCode.equals("GEN_PREDECAY")) {
+ double e = 0.;
+ String list = "";
+ List<MCParticle> daughters = mcp.getDaughters();
+ for (MCParticle daughter : daughters) {
+ list+=" "+map.get(daughter);
+ if (getStatusCode(daughter).equals("GEN_FINAL_STATE")) e+= daughter.getEnergy();
+ else {
+ list+=" (";
+ List<MCParticle> grandDaughters = daughter.getDaughters();
+ for (MCParticle grandDaughter : grandDaughters) {
+ list+=" "+map.get(grandDaughter);
+ e+= grandDaughter.getEnergy();
+ }
+ list+=")";
}
- list+=")";
}
+ System.out.print("\t E = "+df.format(mcp.getEnergy())+" / final state:"+list+" E = "+df.format(e));
}
- System.out.print("\t E = "+df.format(mcp.getEnergy())+" / final state:"+list+" E = "+df.format(e));
+ System.out.println("");
}
- System.out.println("");
text+=" "+map.get(mcp);
}
System.out.println(text); //if (first) System.err.println(text);
@@ -1275,11 +1354,12 @@
//if (firstEvents || (finalEnergy-totalEventEnergy)>10.) {
if (firstEvents) {
//String text = " Total event energy = "+df.format(finalEnergy);
- String text = " Evt #"+nEvt+": Total event energy = "+df.format(finalEnergy)+"/"+df.format(totalEventEnergy);
+ String text = " Evt #"+nEvt+": Total event energy = "+df.format(finalEnergy)+"/"+df.format(totalEventEnergy);
System.out.println(text+".\n"); System.err.println(text);
}
}
+
/**
* Used for standalone reconstruction
*/