lcsim/src/org/lcsim/recon/cluster/structural
diff -N CheatFragmentIdentifier.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CheatFragmentIdentifier.java 2 Mar 2006 20:51:54 -0000 1.1
@@ -0,0 +1,86 @@
+package org.lcsim.recon.cluster.structural;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.event.*;
+import org.lcsim.recon.cluster.analysis.*;
+
+/**
+ * A cheating class to determine whether a given cluster is a fragment.
+ *
+ * @version $Id: CheatFragmentIdentifier.java,v 1.1 2006/03/02 20:51:54 mcharles Exp $
+ */
+public class CheatFragmentIdentifier implements FragmentIdentifier
+{
+ /**
+ * Constructor. The initialization strings are the same as used by
+ * <code>org.lcsim.recon.cluster.analysis.CreateClusterAnalysisLists</code>
+ *
+ * @param hitListNames String array of names of List<CalorimeterHit> in the event
+ * @param clusterListNames String array of names of List<Cluster> in the event.
+ * @param mcListName Name of the self-consistent List<MCParticle> in the event
+ * @param outputListMCParticleToCluster Name to write out the List<MCPClusterInfo> to the event
+ * @param outputListClusterToMCParticle Name to write out the List<ClusterMCPInfo> to the event
+ * @see org.lcsim.recon.cluster.analysis.CreateClusterAnalysisLists
+ */
+ public CheatFragmentIdentifier(String[] hitListNames, String[] clusterListNames, String mcListName, String outputListMCParticleToCluster, String outputListClusterToMCParticle)
+ {
+ m_clusterAssociatorOutputListMCParticleToCluster = outputListMCParticleToCluster;
+ m_clusterAssociatorOutputListClusterToMCParticle = outputListClusterToMCParticle;
+ m_clusterAssociator = new CreateClusterAnalysisLists(hitListNames, clusterListNames, mcListName, outputListMCParticleToCluster, outputListClusterToMCParticle);
+ }
+
+ /**
+ * Determine whether a particular cluster clus is a fragment
+ * using truth information. This is done using the tools in
+ * the org.lcsim.recon.cluster.analysis package.
+ *
+ * The cluster [clus] MUST be a member of the named list(s)
+ * of clusters supplied in the constructor.
+ *
+ * @parameter clus Cluster to test
+ * @parameter event EventHeader
+ */
+ public boolean isFragment(Cluster clus, EventHeader event)
+ {
+ // Check whether the event cache is out of date
+ if (m_clusterAssociatorLastEvent != event) {
+ m_clusterAssociator.CreateLists(event);
+ m_clusterAssociatorLastEvent = event;
+ }
+
+ // Look for the list of cluster info objects
+ List<ClusterMCPInfo> clusterInfoList = event.get(ClusterMCPInfo.class, m_clusterAssociatorOutputListClusterToMCParticle);
+
+ // Find the info:
+ ClusterMCPInfo info = null;
+ for (ClusterMCPInfo currentInfo : clusterInfoList) {
+ Cluster currentCluster = currentInfo.getCluster();
+ if (clus == currentCluster) {
+ info = currentInfo;
+ }
+ }
+ if (info == null) {
+ throw new java.lang.NullPointerException("ERROR: Info not found");
+ } else {
+ // Found it OK.
+ // We look up which the dominant contributing MCParticle is.
+ // Then we look at that MCParticle and see which cluster it
+ // puts the most energy in (the dominant cluster of that particle).
+ //
+ // The cluster we're testing is a primary if and only if it
+ // is the dominant cluster of its dominant MCParticle.
+
+ MCPClusterInfo maxParticleInfo = info.getMaxMCPC();
+ Cluster maxClusterOfMaxParticle = maxParticleInfo.getMaxCluster();
+ boolean isPrimary = (clus == maxClusterOfMaxParticle);
+ return (!isPrimary); // It's a fragment if it's not a primary
+ }
+ }
+
+ // Support for the association:
+ transient protected EventHeader m_clusterAssociatorLastEvent = null;
+ protected CreateClusterAnalysisLists m_clusterAssociator = null;
+ protected String m_clusterAssociatorOutputListMCParticleToCluster;
+ protected String m_clusterAssociatorOutputListClusterToMCParticle;
+}
lcsim/src/org/lcsim/recon/cluster/structural
diff -N FragmentIdentifier.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FragmentIdentifier.java 2 Mar 2006 20:51:54 -0000 1.1
@@ -0,0 +1,20 @@
+package org.lcsim.recon.cluster.structural;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+
+/**
+ * Determine whether a cluster is a fragment or not.
+ *
+ * @version $Id: FragmentIdentifier.java,v 1.1 2006/03/02 20:51:54 mcharles Exp $
+ */
+
+public interface FragmentIdentifier
+{
+ /**
+ * Attempt to determine whether the cluster clus is a fragment.
+ *
+ * @return true if we think it is a fragment, false if not.
+ */
+ public boolean isFragment(Cluster clus, EventHeader event);
+}
lcsim/src/org/lcsim/recon/cluster/structural
diff -N SimpleFragmentIdentifier.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SimpleFragmentIdentifier.java 2 Mar 2006 20:51:54 -0000 1.1
@@ -0,0 +1,164 @@
+package org.lcsim.recon.cluster.structural;
+
+import java.util.*;
+import hep.physics.vec.*;
+
+import org.lcsim.event.*;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+
+/**
+ * A simple cut-based FragmentIdentifier. The following algorithm is used:
+ *
+ * <OL>
+ * <LI> If there is a track associated to the cluster, it is a primary.
+ * <LI> If there is no track associated to the cluster:
+ * <OL> <LI> If nhits < 4, it is a fragment.
+ * <LI> If nhits > n, it is a primary.
+ * <LI> If 4 <= nhits <= n, then compute the distance of closest approach
+ * (DOCA) to the interaction point. It is a primary if the DOCA < d
+ * </OL>
+ * </OL>
+ *
+ * where n and d are supplied as parameters. The nhits>=4 requirement is
+ * hard-coded -- this is because computing the axes of the cluster
+ * (and hence the DOCA to the IP) requires at least four hits.
+ *
+ * To check whether the cluster has a track, the class loops over lists
+ * of ReconstructedParticle supplied by the user with the
+ * <code>addParticleList</code> method. If any charged particle contains
+ * the cluster being tested, it is considered to have a track and therefore
+ * to be a primary.
+ *
+ * @version $Id: SimpleFragmentIdentifier.java,v 1.1 2006/03/02 20:51:54 mcharles Exp $
+ */
+
+public class SimpleFragmentIdentifier implements FragmentIdentifier
+{
+ /**
+ * Constructor.
+ *
+ * @param n Parameter for cut on number of hits (see above)
+ * @param d Parameter for cut on DOCA to IP (see above)
+ */
+ public SimpleFragmentIdentifier(int n, double d) {
+ m_cutNumHits = n;
+ m_cutDOCA = d;
+ }
+
+ /**
+ * Try to determine whether this cluster is a fragment
+ *
+ * @param clus Cluster to test
+ * @param event Current event
+ */
+ public boolean isFragment(Cluster clus, EventHeader event)
+ {
+ // Look up the particles:
+ List<ReconstructedParticle> particleList = new Vector<ReconstructedParticle> ();
+ for (String name : m_particleListNames) {
+ List<ReconstructedParticle> currentList = event.get(ReconstructedParticle.class, name);
+ particleList.addAll(currentList);
+ }
+
+ // Is there a track associated with this cluster?
+ boolean hasTrack = false;
+ for (ReconstructedParticle part : particleList) {
+ if ( Math.abs(part.getCharge()) > 0.5 ) {
+ List<Cluster> clustersAssociatedWithTrack = recursivelyFindSubClusters(part);
+ if (clustersAssociatedWithTrack.contains(clus)) {
+ hasTrack = true;
+ break;
+ }
+ }
+ }
+
+ // Check # hits
+ int nHits = clus.getCalorimeterHits().size();
+ boolean nHitsGood = (nHits > m_cutNumHits);
+
+ if (nHits<4) {
+ // If too few hits (1, 2, maybe 3), DOCA is undefined.
+ // Treat these cases as fragments always (unless there's a track)
+ if (!hasTrack) {
+ return true; // fragment
+ } else {
+ return false; // has track => tiny but real cluster
+ }
+ } else {
+ // Enough hits to check DOCA.
+ BasicHep3Vector originPoint = new BasicHep3Vector(0,0,0);
+ Hep3Vector[] linePosAndDir = getPositionAndDirection(clus);
+ double DOCA = findDOCAToPoint(linePosAndDir[0], linePosAndDir[1], originPoint);
+ boolean docaGood = (DOCA < m_cutDOCA);
+ // It's a fragment if there's no track AND not enough hits AND doesn't point back to IP:
+ boolean isFragment = !hasTrack && !nHitsGood && !docaGood;
+ return isFragment;
+ }
+ }
+
+ /** Utility routine: Find all clusters contained in a ReconstructedParticle */
+ protected List<Cluster> recursivelyFindSubClusters(ReconstructedParticle part)
+ {
+ List<Cluster> allSubClusters = new Vector<Cluster>();
+ for (Cluster clus : part.getClusters()) {
+ List<Cluster> subClusters = recursivelyFindSubClusters(clus);
+ allSubClusters.addAll(subClusters);
+ }
+ return allSubClusters;
+ }
+ /** Utility routine: Find all clusters contained in a Cluster */
+ protected List<Cluster> recursivelyFindSubClusters(Cluster clus)
+ {
+ List<Cluster> output = new Vector<Cluster>();
+ for (Cluster dau : clus.getClusters()) {
+ output.addAll(recursivelyFindSubClusters(dau));
+ }
+ output.add(clus);
+ return output;
+ }
+
+ // Belongs in another class:
+ private double findDOCAToPoint(Hep3Vector linePoint, Hep3Vector lineDirection, Hep3Vector point)
+ {
+ Hep3Vector displacement = new BasicHep3Vector(point.x() - linePoint.x(), point.y() - linePoint.y(), point.z() - linePoint.z());
+ double dotProduct = VecOp.dot(displacement, lineDirection);
+ double doca = Math.sqrt(displacement.magnitudeSquared() - dotProduct*dotProduct);
+ return doca;
+ }
+ // Belongs in another class:
+ private Hep3Vector[] getPositionAndDirection(Cluster clus)
+ {
+ if (clus.getCalorimeterHits().size() < 4) {
+ // Too few hits to calculate the tensor
+ return null;
+ }
+
+ BasicCluster copy = new BasicCluster();
+ copy.addCluster(clus);
+ TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+ copy.setPropertyCalculator(calc);
+ copy.calculateProperties();
+ double[][]axes = calc.getPrincipleAxis();
+ if (axes == null) {
+ throw new AssertionError("Principal axes not calculated");
+ }
+ Hep3Vector[] posAndDir = new Hep3Vector[2];
+ posAndDir[0] = new BasicHep3Vector(copy.getPosition());
+ posAndDir[1] = new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]);
+ return posAndDir;
+ }
+
+ /**
+ * Add a list of particles (via its string name). Any clusters which belong
+ * to a charged particle in this list are considered primaries.
+ */
+ public void addParticleList(String name) {m_particleListNames.add(name);}
+
+ /** Wipe all information about particle lists. */
+ public void clearParticleLists() { m_particleListNames.clear(); }
+
+ List<String> m_particleListNames = new Vector<String>();
+ int m_cutNumHits;
+ double m_cutDOCA;
+}
