lcsim/src/org/lcsim/recon/cheater
diff -N ReconCheater.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReconCheater.java 26 Sep 2005 22:47:06 -0000 1.1
@@ -0,0 +1,629 @@
+/* ReconCheater.java
+ *
+ * ...
+ *
+ * Written by Mike Ronan, LBNL (Berkeley) Fall 2005.
+ * Modified by:
+ *
+ */
+
+package org.lcsim.recon.cheater;
+
+// LCSIM framework packages
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.util.MCParticleClassifier;
+import org.lcsim.conditions.ConditionsEvent;
+import org.lcsim.conditions.ConditionsListener;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
+import org.lcsim.recon.cluster.cheat.ClusterCheater;
+import org.lcsim.recon.cluster.util.RefinedCluster;
+import org.lcsim.recon.cluster.cheat.CheatClusterDriver;
+import org.lcsim.recon.ztracking.FoundTrack;
+import org.lcsim.recon.ztracking.cheater.CheatTrack;
+import org.lcsim.recon.ztracking.cheater.TrackingCheater;
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+
+import org.lcsim.util.loop.LCSimLoop;
+import org.lcsim.util.aida.AIDA;
+
+import hep.physics.particle.*; // Particle,...
+import hep.physics.vec.*; // Hep3Vector
+
+import java.io.File;
+import java.net.URL;
+import java.text.*; // DecimalFormat
+import java.util.List;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Random;
+
+
+/** ... */
+public class ReconCheater extends Driver
+{
+ 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
+
+ static boolean Hist = false;
+
+
+ /** */
+ public ReconCheater()
+ {
+ this(Hist);
+ }
+
+ public ReconCheater(boolean hist)
+ {
+ // super(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";
+ System.out.println(text); System.err.println(text); // if (!XReconstruction.getBatchMode()) out.println(text);
+
+ if (!TrackingCheater.INITIALIZED) {
+ System.out.println(" Charged tracking using cheater - ");
+ TrackingCheater trackFinder = new TrackingCheater();
+ add(trackFinder);
+ }
+ System.out.println(" EM and Had clustering using cheater - ");
+ CheatClusterDriver clusterFinder = new CheatClusterDriver();
+ add(clusterFinder);
+
+ if (useFullTruth) useTruth = true;
+ df.setMaximumFractionDigits(3);
+ }
+
+ protected int nEvt = 0;
+ boolean first = true, firstEvents = true;
+ boolean hist;
+
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ nEvt++;
+ if (firstEvents) System.out.println(" Event #"+nEvt+" (ReconCheater)");
+
+ hist = getHistogramLevel() > 0;
+ if (hist) System.err.println(" ReconCheater: hist = "+hist);
+ /*
+ if (IDEff == null)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("IDEfficiency");
+ conditions.addConditionsListener(this);
+ IDEff = new IDResolutionTables(conditions);
+ }
+ */
+ // Get MCParticle mapping from event.
+ map = getMCParticleMapping(event);
+
+ // Get reconstructed particles from Track and Cluster Cheaters.
+ List<ReconstructedParticle> particles = getReconstructedParticles(event);
+
+ // Add the reconstructed particles to the event.
+ event.put("ReconCheater", particles, ReconstructedParticle.class, 0);
+
+ first = false;
+ if (nEvt>=3) firstEvents = false;
+ }
+ /*
+ public void conditionsChanged(ConditionsEvent event)
+ {
+ ConditionsSet conditions = getConditionsManager().getConditions("IDEfficiency");
+ IDEff = new IDResolutionTables(conditions);
+ }
+ */
+
+ Map<MCParticle,CheatTrack> charged = null; // Tracks from MCParticles
+ Map<MCParticle,CheatCluster> neutrals = null; // and clusters.
+ double finalEnergy;
+
+ /** Get reconstructed particles from Cheaters that match MCParticles,
+ * then select particles to add to them.
+ */
+ List<ReconstructedParticle> getReconstructedParticles(EventHeader event)
+ {
+ String text = " Get selected particles.";
+ if (firstEvents) System.out.println(text); if (first) System.err.println(text);
+
+ // Get lists of charged and neutral MC particles.
+ List<MCParticle> chargedList = new ArrayList<MCParticle>();
+ List<MCParticle> neutralList = new ArrayList<MCParticle>();
+
+ getPythiaMCParticles(event,chargedList,neutralList);
+
+ List<ReconstructedParticle> particles = new ArrayList();
+
+ charged = new HashMap<MCParticle,CheatTrack>();
+ neutrals = new HashMap<MCParticle,CheatCluster>();
+ Map<MCParticle,MCParticle> primaryMCParents = new HashMap<MCParticle,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>();
+ 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.");
+
+ if (name.equals("CombinedTracks")) {
+ for (CheatTrack track : collection) {
+ 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");
+
+ 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);
+ }
+
+ 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);
+ }
+ }
+ }
+ }
+ 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);
+
+ // 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);
+
+ List<CheatCluster> 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.");
+
+ if (name.equals("RefinedClusters")) {
+ for (CheatCluster cluster : collection) {
+ 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");
+
+ 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);
+ }
+
+ // Check if generator neutral list contains cluster MCParticle.
+ if (neutralList.contains(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 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);
+ }
+ }
+ }
+ }
+ 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>();
+
+ 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);
+
+
+ // Select particles:
+ if (firstEvents) System.out.println(" Select particles to be reconstructed.");
+
+ // Add charged particles that were found in TrackingCheater list.
+ for (MCParticle mcp : charged.keySet()) {
+ CheatTrack track = charged.get(mcp);
+ if (useFullTruth) particles.add(new CheatReconstructedParticle(mcp));
+ else particles.add(new CheatReconstructedParticle(track,(Particle)mcp));
+ }
+ // Add neutral particles that were found in ClusterCheater list.
+ for (MCParticle mcp : neutrals.keySet()) {
+ CheatCluster cluster = neutrals.get(mcp);
+ double e = mcp.getEnergy(); if (e<=mcp.getType().getMass()) continue;
+ if (useFullTruth) particles.add(new CheatReconstructedParticle(mcp));
+ else particles.add(new CheatReconstructedParticle(cluster,getMeasuredEnergy(cluster),(Particle)mcp));
+ }
+
+ // Charged and neutral lists contain MCParticles that interacted or decayed.
+ if (useFullTruth) {
+ // Add interacting charged particles.
+ for (MCParticle mcp : chargedList) {
+ particles.add(new CheatReconstructedParticle(mcp));
+ finalEnergy += mcp.getEnergy();
+ }
+ // Add interacting neutral particles.
+ for (MCParticle mcp : neutralList) {
+ double e = mcp.getEnergy(); if (e<=mcp.getType().getMass()) continue;
+ particles.add(new CheatReconstructedParticle(mcp));
+ finalEnergy += mcp.getEnergy();
+ }
+ }
+ else {
+ // Add interacting parents.
+ for (MCParticle mcp : parents) {
+ 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();
+ particles.add(new CheatReconstructedParticle(track,(Particle)mcp));
+ finalEnergy += mcp.getEnergy();
+ }
+ // Add extra clusters.
+ for (CheatCluster cluster : extraClusterList) {
+ MCParticle mcp = cluster.getMCParticle();
+ double e = mcp.getEnergy(); if (e<=mcp.getType().getMass()) 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);
+
+ // Total energy
+ if (firstEvents) System.out.println(" Total event energy = "+df.format(finalEnergy)+".\n");
+
+ if (hist) analyzeReconstructedParticles(particles);
+
+ return particles;
+ }
+
+ void analyzeReconstructedParticles(List<ReconstructedParticle> particles)
+ {
+ AIDA aida = AIDA.defaultInstance();
+
+ aida.cloud1D("# particles").fill(particles.size());
+ aida.cloud1D("Final Energy").fill(finalEnergy);
+
+ for (MCParticle mcp : neutrals.keySet()) {
+ 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();
+
+ if (PDGID==22) {
+ ratio = cluster.getRawEnergy()/mcp.getEnergy();
+ aida.cloud1D("photon ratio").fill(ratio);
+ aida.cloud2D("photon Energy vs. ratio").fill(ratio,mcp.getEnergy());
+ }
+ else if (PDGID==130) {
+ aida.cloud1D("Klong ratio").fill(ratio);
+ aida.cloud2D("Klong Energy vs. ratio").fill(ratio,mcp.getEnergy());
+ }
+ else if (Math.abs(PDGID)==2112) {
+ aida.cloud1D("neutron ratio").fill(ratio);
+ aida.cloud2D("neutron Energy vs. ratio").fill(ratio,mcp.getEnergy());
+ }
+ }
+ }
+
+ MCParticleClassifier mcpClassifier = new MCParticleClassifier();
+
+ String getStatusCode(MCParticle mcp)
+ {
+ String code = ""+mcpClassifier.getClassification(mcp);
+ return code;
+ }
+
+ Random random = new Random();
+
+ double getMeasuredEnergy(CheatCluster cluster)
+ {
+ MCParticle mcp = cluster.getMCParticle();
+ double energy = cluster.getRawEnergy();
+ if (mcp.getPDGID()==22) {
+ energy = energy/ECalSampling;
+ }
+ else {
+ //energy = energy/HCalSampling;
+ int nhits = 0;
+ for (CalorimeterHit hit : cluster.getCalorimeterHits())
+ if (hit.getTime()<100.) nhits++;
+ energy = nhits/HCalDigital;
+ }
+ Hep3Vector pV = mcp.getMomentum();
+ double r = energy/pV.magnitude();
+ if (useTruth) r = 1.;
+ else {
+ if (mcp.getPDGID()==22) {
+ if (useECalParameterization) {
+ r = 1. - random.nextGaussian()*ECalResolution/Math.sqrt(mcp.getEnergy());
+ if (r<0.1) r = 0.1;
+ }
+ }
+ else
+ if (useHCalParameterization) {
+ r = 1. - random.nextGaussian()*HCalResolution/Math.sqrt(mcp.getEnergy());
+ if (r<0.1) r = 0.1;
+ }
+ }
+ return r * pV.magnitude();
+ }
+
+ Map<MCParticle,Integer> map = null;
+
+ /** Get ... */
+ Map<MCParticle,Integer> getMCParticleMapping(EventHeader event)
+ {
+ Map<MCParticle,Integer> map = new HashMap<MCParticle,Integer>();
+
+ // Get Monte Carlo particles from LCSIM event
+ List<MCParticle> particleList = event.getMCParticles();
+ int n=-1;
+ for (MCParticle mcp : particleList) {
+ Particle particle = (Particle) mcp; n++;
+ map.put(mcp,n);
+ }
+ return map;
+ }
+
+ /** Get Pythia Monte Carlo charged and neutral lists. */
+ public 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.";
+ if (firstEvents) System.out.println(text); if (first) System.err.println(text);
+
+ List<MCParticle> particleList = event.getMCParticles();
+
+ int NMCParticles = 0;
+
+ int n=-1,ipn=0; boolean initial = true; double finalEnergy=0.;
+ for (MCParticle mcp : particleList) {
+ Particle particle = (Particle) mcp; n++;
+
+ String statusCode = ""+getStatusCode(mcp);
+ boolean finalState = statusCode.equals("GEN_FINAL_STATE");
+
+ //if (finalState) {
+ int PDGID = particle.getPDGID();
+ String Name = particle.getType().getName();
+
+ int parentPDGID = 0; String parentNumber = "";
+ List<MCParticle> parents = particle.getParents();
+ Particle parent = null; if (!parents.isEmpty()) parent = parents.get(0);
+ if (parent!=null) { parentPDGID = parent.getPDGID(); parentNumber = " (#"+map.get(parent)+")"; }
+
+ boolean neutrino = Math.abs(PDGID)==NuEID || Math.abs(PDGID)==NuMuID || Math.abs(PDGID)==NuTauID;
+ if (neutrino) {
+ parentPDGID = parent.getPDGID();
+ while (parentPDGID == PDGID) {
+ List<MCParticle> grandParents = parent.getParents();
+ Particle grandParent = null; if (!grandParents.isEmpty()) grandParent = grandParents.get(0);
+ parent = grandParent; if (parent!=null) parentPDGID = parent.getPDGID();
+ }
+ //continue;
+ }
+ double p = Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()+mcp.getPZ()*mcp.getPZ());
+ text = " Particle #"+n+" (PDGID "+PDGID+") "+Name+" (p="+df.format(p)+") \tstatus = "+
+ statusCode+" parent = "+parentNumber+parentPDGID;
+ if (firstEvents) System.out.println(text);
+
+ boolean quark = Math.abs(PDGID)<=6;
+ boolean gluon = PDGID==21;
+ boolean unknown = PDGID==0 || PDGID==92;
+
+ if (finalState && !neutrino && !(quark||gluon||unknown)) {
+ NMCParticles++;
+ Hep3Vector PV = particle.getMomentum();
+ double Pt = Math.sqrt(PV.x()*PV.x()+PV.y()*PV.y());
+ double P = Math.sqrt(Pt*Pt+PV.z()*PV.z());
+ double CosTh = -1; if (P!=0.) CosTh = PV.z()/P;
+ finalEnergy += particle.getEnergy();
+
+ if (Math.abs(CosTh)<InitialCosThMin) {
+ ipn++;
+ if (particle.getCharge()==0) neutralList.add(mcp);
+ else chargedList.add(mcp);
+ }
+ }
+ }
+ if (firstEvents) System.out.println(" Found "+NMCParticles+" final state Monte Carlo particles"+
+ " with total energy = "+df.format(finalEnergy)+".");
+ }
+
+
+ /**
+ * Used for standalone reconstruction
+ */
+ public static void main(String[] argv)
+ {
+ String JobName = "ReconCheater "+argv[0];
+ try {
+ /*
+ URL url = new URL("http://www.lcsim.org/datasamples/slic_sdjan03_K0S_Theta90_1-10GeV.slcio");
+ FileCache cache = new FileCache();
+ File file = cache.getCachedFile(url);
+ */
+ File file = new File(argv[1]);
+ //LCIODriver writer = new LCIODriver("output.slcio");
+
+ LCSimLoop loop = new LCSimLoop();
+ loop.setLCIORecordSource(file);
+ loop.add(new ReconCheater());
+ //Test LoopProcessor processor = new LoopProcessor(argv,loop);
+ // LoopProcessor processor = new LoopProcessor(argv,nEvtFile,loop,parameters);
+ loop.loop(-1);
+ loop.dispose();
+
+ AIDA.defaultInstance().saveAs("new.aida");
+ }
+ catch (Exception x) { }
+ {
+ }
+ }
+
+ protected DecimalFormat df = new DecimalFormat();
+
+ protected final static double InitialCosThMin = 0.9999;
+
+ protected final static int NuEID = 12;
+ protected final static int NuMuID = 14;
+ protected final static int NuTauID = 16;
+}