lcsim/src/org/lcsim/recon/pfa/output
diff -N ConfusionPlotter.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ConfusionPlotter.java 15 Jul 2006 20:19:47 -0000 1.1
@@ -0,0 +1,439 @@
+package org.lcsim.recon.pfa.output;
+
+import java.util.List;
+import java.util.Vector;
+import java.util.Set;
+import java.util.HashSet;
+import java.io.IOException;
+
+import hep.aida.ITree;
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.ICloud1D;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
+
+import hep.physics.particle.Particle;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.Track;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.CalorimeterHit;
+
+/**
+ * Plot the confusion distributions.
+ *
+ * There are two categories of reconstructed particle:
+ * particles with a track (charged), and
+ * particles without a track (neutral).
+ *
+ * There are three categories of truth particle:
+ * charged with an associated track,
+ * charged with no associated track, and
+ * neutral.
+ *
+ * For each truth category, we determine what fraction of the
+ * particle assignments are right or wrong.
+ *
+ * @version $Id: ConfusionPlotter.java,v 1.1 2006/07/15 20:19:47 mcharles Exp $
+ */
+
+public class ConfusionPlotter extends Driver
+{
+ boolean m_debug = false;
+ String m_inputRecoParticleListName;
+ String m_inputTrueParticleListName;
+ public ConfusionPlotter(String recoParticleList, String trueParticleList) {
+ m_inputRecoParticleListName = recoParticleList;
+ m_inputTrueParticleListName = trueParticleList;
+ initPlots();
+ }
+
+ public void setDebug(boolean debug) { m_debug = debug; }
+
+ public void process(EventHeader event)
+ {
+ // Here is the list of reconstructed particles:
+ List<ReconstructedParticle> inputRecoList = event.get(ReconstructedParticle.class, m_inputRecoParticleListName);
+
+ // Here is the list of truth particles:
+ List<MCParticle> inputTrueList = event.get(MCParticle.class, m_inputTrueParticleListName);
+
+ // We need to figure out which truth particles have tracks
+ // that were matched to clusters (i.e. used in PFA).
+ List<Particle> truthParticlesWithMatchedTracks = new Vector<Particle>();
+ for (ReconstructedParticle inputParticle : inputRecoList) {
+ List<Track> tracks = inputParticle.getTracks();
+ if (tracks != null) {
+ for (Track track : tracks) {
+ if (track instanceof ReconTrack) {
+ ReconTrack rt = (ReconTrack) (track);
+ Particle mc = rt.getMCParticle();
+ truthParticlesWithMatchedTracks.add(mc);
+ } else {
+ throw new AssertionError("ERROR in "+this.getClass()+": I can't get truth information from this track: "+track);
+ }
+ }
+ }
+ }
+
+ // Book-keeping:
+ flushCache();
+
+ // Loop over all the reconstructed particles and
+ // identify contributions from different truth categories.
+ for (ReconstructedParticle inputParticle : inputRecoList) {
+ // Did we reconstruct this particle as charged or neutral?
+ boolean hasTracks = (inputParticle.getTracks() != null && inputParticle.getTracks().size() > 0);
+ boolean hasCharge = (Math.abs(inputParticle.getCharge()) > 0.5);
+ boolean recoIsCharged = hasTracks || hasCharge;
+ // What contributed to the calorimeter cluster?
+ List<Cluster> clusters = inputParticle.getClusters();
+ for (Cluster clus : clusters) {
+ // For each cluster, find the contributing truth particles:
+ Set<MCParticle> contributingParticles = findMCParticles(clus, inputTrueList);
+ for (MCParticle truthPart : contributingParticles) {
+ // Figure out whether the truth particle has charge and/or a track:
+ boolean truthIsCharged = (Math.abs(truthPart.getCharge())>0.5);
+ boolean truthHasTrack = (truthParticlesWithMatchedTracks.contains(truthPart));
+ if (!truthIsCharged && truthHasTrack) { throw new AssertionError("ERROR: MCParticle "+truthPart+" is neutral but has a track"); }
+ // Count hits and energy:
+ int hitsHCAL = countHitsInClusterHCAL(truthPart, clus, inputTrueList);
+ double energyECAL = correctedEnergyInClusterECAL(truthPart, clus, inputTrueList);
+ // Record:
+ if (truthIsCharged && truthHasTrack) {
+ if (recoIsCharged) {
+ energyECAL_ChargedParticleWithTrack_ClusterWithTrack += energyECAL;
+ hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack += hitsHCAL;
+ } else {
+ energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack += energyECAL;
+ hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack += hitsHCAL;
+ }
+ } else if (truthIsCharged && !truthHasTrack) {
+ if (recoIsCharged) {
+ energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack += energyECAL;
+ hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack += hitsHCAL;
+ } else {
+ energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack += energyECAL;
+ hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack += hitsHCAL;
+ }
+ } else if (!truthIsCharged && !truthHasTrack) {
+ if (recoIsCharged) {
+ energyECAL_NeutralParticle_ClusterWithTrack += energyECAL;
+ hitsHCAL_NeutralParticle_ClusterWithTrack += hitsHCAL;
+ } else {
+ energyECAL_NeutralParticle_ClusterWithNoTrack += energyECAL;
+ hitsHCAL_NeutralParticle_ClusterWithNoTrack += hitsHCAL;
+ }
+ } else {
+ throw new AssertionError("I don't know how to handle this case: truthIsCharged="+truthIsCharged+", truthHasTrack="+truthHasTrack);
+ }
+ }
+ }
+ }
+
+ if (m_debug) {
+ debugPrintout();
+ }
+ fillPlots();
+ }
+
+ public void suspend() {
+ try {
+ m_tree.commit();
+ } catch(IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+ super.suspend();
+ }
+
+ protected void initPlots()
+ {
+ IAnalysisFactory af = IAnalysisFactory.create();
+ try {
+ m_tree = af.createTreeFactory().create("confusion.aida","xml",false,true);
+ m_histoFactory = af.createHistogramFactory(m_tree);
+ m_hECAL_charged_track_to_track = m_histoFactory.createCloud1D("hECAL_charged_track_to_track");
+ m_hECAL_charged_track_to_notrack = m_histoFactory.createCloud1D("hECAL_charged_track_to_notrack");
+ m_hECAL_charged_notrack_to_track = m_histoFactory.createCloud1D("hECAL_charged_notrack_to_track");
+ m_hECAL_charged_notrack_to_notrack = m_histoFactory.createCloud1D("hECAL_charged_notrack_to_notrack");
+ m_hHCAL_charged_track_to_track = m_histoFactory.createCloud1D("hHCAL_charged_track_to_track");
+ m_hHCAL_charged_track_to_notrack = m_histoFactory.createCloud1D("hHCAL_charged_track_to_notrack");
+ m_hHCAL_charged_notrack_to_track = m_histoFactory.createCloud1D("hHCAL_charged_notrack_to_track");
+ m_hHCAL_charged_notrack_to_notrack = m_histoFactory.createCloud1D("hHCAL_charged_notrack_to_notrack");
+ m_hECAL_neutral_to_track = m_histoFactory.createCloud1D("hECAL_neutral_to_track");
+ m_hECAL_neutral_to_notrack = m_histoFactory.createCloud1D("hECAL_neutral_to_notrack");
+ m_hHCAL_neutral_to_track = m_histoFactory.createCloud1D("hHCAL_neutral_to_track");
+ m_hHCAL_neutral_to_notrack = m_histoFactory.createCloud1D("hHCAL_neutral_to_notrack");
+ m_hECAL_charged_track_fractionbad = m_histoFactory.createCloud1D("hECAL_charged_track_fractionbad");
+ m_hECAL_charged_notrack_fractionbad = m_histoFactory.createCloud1D("hECAL_charged_notrack_fractionbad");
+ m_hECAL_neutral_fractionbad = m_histoFactory.createCloud1D("hECAL_neutral_fractionbad");
+ m_hHCAL_charged_track_fractionbad = m_histoFactory.createCloud1D("hHCAL_charged_track_fractionbad");
+ m_hHCAL_charged_notrack_fractionbad = m_histoFactory.createCloud1D("hHCAL_charged_notrack_fractionbad");
+ m_hHCAL_neutral_fractionbad = m_histoFactory.createCloud1D("hHCAL_neutral_fractionbad");
+ } catch (IOException ioe1) {
+ ioe1.printStackTrace();
+ }
+ }
+
+ ITree m_tree = null;
+ IHistogramFactory m_histoFactory = null;
+ ICloud1D m_hECAL_charged_track_to_track ;
+ ICloud1D m_hECAL_charged_track_to_notrack ;
+ ICloud1D m_hECAL_charged_notrack_to_track ;
+ ICloud1D m_hECAL_charged_notrack_to_notrack ;
+ ICloud1D m_hHCAL_charged_track_to_track ;
+ ICloud1D m_hHCAL_charged_track_to_notrack ;
+ ICloud1D m_hHCAL_charged_notrack_to_track ;
+ ICloud1D m_hHCAL_charged_notrack_to_notrack ;
+ ICloud1D m_hECAL_neutral_to_track ;
+ ICloud1D m_hECAL_neutral_to_notrack ;
+ ICloud1D m_hHCAL_neutral_to_track ;
+ ICloud1D m_hHCAL_neutral_to_notrack ;
+
+ ICloud1D m_hECAL_charged_track_fractionbad ;
+ ICloud1D m_hECAL_charged_notrack_fractionbad ;
+ ICloud1D m_hECAL_neutral_fractionbad ;
+ ICloud1D m_hHCAL_charged_track_fractionbad ;
+ ICloud1D m_hHCAL_charged_notrack_fractionbad ;
+ ICloud1D m_hHCAL_neutral_fractionbad ;
+
+
+ double energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack = 0.0;
+ double energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack = 0.0;
+ double energyECAL_ChargedParticleWithTrack_ClusterWithTrack = 0.0;
+ double energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack = 0.0;
+ double energyECAL_NeutralParticle_ClusterWithTrack = 0.0;
+ double energyECAL_NeutralParticle_ClusterWithNoTrack = 0.0;
+ int hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack = 0;
+ int hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack = 0;
+ int hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack = 0;
+ int hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack = 0;
+ int hitsHCAL_NeutralParticle_ClusterWithTrack = 0;
+ int hitsHCAL_NeutralParticle_ClusterWithNoTrack = 0;
+
+ protected void flushCache() {
+ energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack = 0.0;
+ energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack = 0.0;
+ energyECAL_ChargedParticleWithTrack_ClusterWithTrack = 0.0;
+ energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack = 0.0;
+ energyECAL_NeutralParticle_ClusterWithTrack = 0.0;
+ energyECAL_NeutralParticle_ClusterWithNoTrack = 0.0;
+ hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack = 0;
+ hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack = 0;
+ hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack = 0;
+ hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack = 0;
+ hitsHCAL_NeutralParticle_ClusterWithTrack = 0;
+ hitsHCAL_NeutralParticle_ClusterWithNoTrack = 0;
+ }
+
+ protected void debugPrintout() {
+ // Some debug printout:
+ System.out.println("In this event, ECAL energy for charged particles with tracks: "
+ +energyECAL_ChargedParticleWithTrack_ClusterWithTrack+" (particle with track -> cluster with track) "
+ +energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack+" (particle with track -> cluster with no track) ");
+ System.out.println("In this event, ECAL energy for charged particles with no track: "
+ +energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack+" (particle with no track -> cluster with track) "
+ +energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack+" (particle with no track -> cluster with no track) ");
+ System.out.println("In this event, ECAL energy for neutral particles: "
+ +energyECAL_NeutralParticle_ClusterWithTrack+" (particle -> cluster with track) "
+ +energyECAL_NeutralParticle_ClusterWithNoTrack+" (particle -> cluster with no track) ");
+ System.out.println("In this event, HCAL hits for charged particles with tracks: "
+ +hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack+" (particle with track -> cluster with track) "
+ +hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack+" (particle with track -> cluster with no track) ");
+ System.out.println("In this event, HCAL hits for charged particles with no track: "
+ +hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack+" (particle with no track -> cluster with track) "
+ +hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack+" (particle with no track -> cluster with no track) ");
+ System.out.println("In this event, HCAL hits for neutral particles: "
+ +hitsHCAL_NeutralParticle_ClusterWithTrack+" (particle -> cluster with track) "
+ +hitsHCAL_NeutralParticle_ClusterWithNoTrack+" (particle -> cluster with no track) ");
+ }
+
+ protected void fillPlots()
+ {
+ // More book-keeping:
+ double badECAL_track_to_notrack = energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack;
+ double goodECAL_track_to_track = energyECAL_ChargedParticleWithTrack_ClusterWithTrack;
+ double badECAL_notrack_to_track = energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack + energyECAL_NeutralParticle_ClusterWithTrack;
+ double goodECAL_notrack_to_notrack = energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack + energyECAL_NeutralParticle_ClusterWithNoTrack;
+ int badHCAL_track_to_notrack = hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack;
+ int goodHCAL_track_to_track = hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack;
+ int badHCAL_notrack_to_track = hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack + hitsHCAL_NeutralParticle_ClusterWithTrack;
+ int goodHCAL_notrack_to_notrack = hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack + hitsHCAL_NeutralParticle_ClusterWithNoTrack;
+
+ // Fill plots:
+ m_hECAL_charged_track_to_track .fill(energyECAL_ChargedParticleWithTrack_ClusterWithTrack);
+ m_hECAL_charged_track_to_notrack .fill(energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack);
+ m_hECAL_charged_notrack_to_track .fill(energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack);
+ m_hECAL_charged_notrack_to_notrack .fill(energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack);
+ m_hHCAL_charged_track_to_track .fill(hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack);
+ m_hHCAL_charged_track_to_notrack .fill(hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack);
+ m_hHCAL_charged_notrack_to_track .fill(hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack);
+ m_hHCAL_charged_notrack_to_notrack .fill(hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack);
+ m_hECAL_neutral_to_track .fill(energyECAL_NeutralParticle_ClusterWithTrack);
+ m_hECAL_neutral_to_notrack .fill(energyECAL_NeutralParticle_ClusterWithNoTrack);
+ m_hHCAL_neutral_to_track .fill(hitsHCAL_NeutralParticle_ClusterWithTrack);
+ m_hHCAL_neutral_to_notrack .fill(hitsHCAL_NeutralParticle_ClusterWithNoTrack);
+
+ double fraction_ECAL_chargedtrack_bad = 0.0;
+ double fraction_ECAL_chargednotrack_bad = 0.0;
+ double fraction_ECAL_neutral_bad = 0.0;
+ double fraction_HCAL_chargedtrack_bad = 0.0;
+ double fraction_HCAL_chargednotrack_bad = 0.0;
+ double fraction_HCAL_neutral_bad = 0.0;
+
+ if (energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack>0.0) {
+ fraction_ECAL_chargedtrack_bad = (energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack / (energyECAL_ChargedParticleWithTrack_ClusterWithTrack+energyECAL_ChargedParticleWithTrack_ClusterWithNoTrack));
+ }
+ if (energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack>0.0) {
+ fraction_ECAL_chargednotrack_bad =(energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack / (energyECAL_ChargedParticleWithNoTrack_ClusterWithTrack+energyECAL_ChargedParticleWithNoTrack_ClusterWithNoTrack));
+ }
+ if (energyECAL_NeutralParticle_ClusterWithTrack>0.0) {
+ fraction_ECAL_neutral_bad = (energyECAL_NeutralParticle_ClusterWithTrack / (energyECAL_NeutralParticle_ClusterWithTrack+energyECAL_NeutralParticle_ClusterWithNoTrack));
+ }
+ if (hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack > 0) {
+ fraction_HCAL_chargedtrack_bad = ( ((double)(hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack)) / ((double)(hitsHCAL_ChargedParticleWithTrack_ClusterWithNoTrack+hitsHCAL_ChargedParticleWithTrack_ClusterWithTrack)) );
+ }
+ if (hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack > 0) {
+ fraction_HCAL_chargednotrack_bad = ( ((double)(hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack)) / ((double)(hitsHCAL_ChargedParticleWithNoTrack_ClusterWithTrack + hitsHCAL_ChargedParticleWithNoTrack_ClusterWithNoTrack)) );
+ }
+ if (hitsHCAL_NeutralParticle_ClusterWithTrack > 0) {
+ fraction_HCAL_neutral_bad = ( ((double)(hitsHCAL_NeutralParticle_ClusterWithTrack)) / ((double)(hitsHCAL_NeutralParticle_ClusterWithTrack + hitsHCAL_NeutralParticle_ClusterWithNoTrack)) );
+ }
+
+ m_hECAL_charged_track_fractionbad .fill(fraction_ECAL_chargedtrack_bad);
+ m_hECAL_charged_notrack_fractionbad.fill(fraction_ECAL_chargednotrack_bad);
+ m_hECAL_neutral_fractionbad .fill(fraction_ECAL_neutral_bad);
+ m_hHCAL_charged_track_fractionbad .fill(fraction_HCAL_chargedtrack_bad);
+ m_hHCAL_charged_notrack_fractionbad.fill(fraction_HCAL_chargednotrack_bad);
+ m_hHCAL_neutral_fractionbad .fill(fraction_HCAL_neutral_bad);
+ }
+
+ // FIXME: This belongs in an external class
+ protected int countHitsInClusterHCAL(MCParticle truthPart, Cluster clus, List<MCParticle> inputTrueList) {
+ int countHitsHCAL = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean useThisHit = hitMatch(truthPart, (SimCalorimeterHit)(hit), inputTrueList);
+ if (useThisHit) {
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("HADBarrel") == 0 || name.compareTo("HADEndcap") == 0) {
+ // HCAL
+ if (hit.getTime() < 100) {
+ // Within first 100 ns => OK
+ countHitsHCAL++;
+ }
+ }
+ }
+ }
+ return countHitsHCAL;
+ }
+
+ // FIXME: This belongs in an external class
+ protected double correctedEnergyInClusterECAL(MCParticle truthPart, Cluster clus, List<MCParticle> inputTrueList) {
+ double energyECAL = 0.0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ boolean useThisHit = hitMatch(truthPart, (SimCalorimeterHit)(hit), inputTrueList);
+ if (useThisHit) {
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0) {
+ // ECAL
+ SimCalorimeterHit simHit = (SimCalorimeterHit) hit;
+ double rawEnergy = hitMatchEnergy(truthPart, simHit, inputTrueList);
+ double rawEnergyTotal = hit.getRawEnergy();
+ double frac = (rawEnergy / rawEnergyTotal);
+ energyECAL += frac * hit.getCorrectedEnergy();
+ }
+ }
+ }
+ return energyECAL;
+ }
+
+ // FIXME: This belongs in an external class
+ protected boolean hitMatch(MCParticle part, SimCalorimeterHit hit, List<MCParticle> inputTrueList) {
+ Set<MCParticle> matchingParticles = findMCParticles(hit, inputTrueList);
+ for (MCParticle hitPart : matchingParticles) {
+ if (part == hitPart) {
+ return true;
+ }
+ }
+ return false;
+ }
+
+ // FIXME: This belongs in an external class
+ protected double hitMatchEnergy(MCParticle part, SimCalorimeterHit hit, List<MCParticle> inputTrueList) {
+ int nContributingParticles = hit.getMCParticleCount();
+ double energySum = 0.0;
+ for (int i=0; i<nContributingParticles; i++) {
+ MCParticle hitPart = hit.getMCParticle(i);
+ List<MCParticle> parentsInList = findParentsInList(hitPart, inputTrueList);
+ if (parentsInList.contains(part) || part==null) {
+ energySum += hit.getContributedEnergy(i);
+ }
+ }
+ return energySum;
+ }
+
+ /**
+ * Internal utility routine, belongs in another class
+ */
+ protected Set<MCParticle> findMCParticles(Cluster clus, List<MCParticle> mcList)
+ {
+ Set clusterParticles = new HashSet<MCParticle>();
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Set<MCParticle> hitParticles = findMCParticles(hit, mcList);
+ clusterParticles.addAll(hitParticles);
+ }
+ return clusterParticles;
+ }
+
+ /**
+ * Internal utility routine, belongs in another class
+ */
+ protected Set<MCParticle> findMCParticles(CalorimeterHit hit, List<MCParticle> mcList)
+ {
+ if ( ! (hit instanceof SimCalorimeterHit) ) {
+ throw new AssertionError("Non-simulated hit!");
+ } else {
+ SimCalorimeterHit simHit = (SimCalorimeterHit) (hit);
+ Set<MCParticle> contributingParticlesFromList = new HashSet<MCParticle>();
+ int nContributingParticles = simHit.getMCParticleCount();
+ for (int i=0; i<nContributingParticles; i++) {
+ MCParticle part = simHit.getMCParticle(i);
+ List<MCParticle> parentsInList = findParentsInList(part, mcList);
+ contributingParticlesFromList.addAll(parentsInList);
+ }
+ return contributingParticlesFromList;
+ }
+ }
+
+ /**
+ * Internal utility routine, belongs in another class
+ */
+ protected List<MCParticle> findParentsInList(MCParticle part, List<MCParticle> mcList)
+ {
+ List<MCParticle> outputList = new Vector<MCParticle>();
+ if (mcList.contains(part)) {
+ // Already in there
+ outputList.add(part);
+ } else {
+ // Not in there -- recurse up through parents
+ List<MCParticle> parents = part.getParents();
+ if (parents.size()==0) {
+ // Ran out of options -- add nothing and return below
+ } else {
+ for (MCParticle parent : parents) {
+ List<MCParticle> ancestorsInList = findParentsInList(parent, mcList);
+ outputList.addAll(ancestorsInList);
+ }
+ }
+ }
+ return outputList;
+ }
+}
