lcsim/src/org/lcsim/contrib/uiowa/structural
diff -u -r1.1 -r1.2
--- CheatFragmentIdentifier.java 16 Dec 2005 21:11:39 -0000 1.1
+++ CheatFragmentIdentifier.java 26 Dec 2005 21:26:50 -0000 1.2
@@ -8,34 +8,100 @@
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.SimCalorimeterHit;
-
import util.HitCountAssociator;
+/**
+ * A cheating class to determine whether a given cluster is a fragment.
+ */
+
public class CheatFragmentIdentifier implements FragmentIdentifier
{
+ /**
+ * Constructor.
+ *
+ * @parameter clusterListName The name of the list of Cluster objects to which we'll compare potential fragments.
+ */
public CheatFragmentIdentifier(String clusterListName)
{
m_clusterListName = clusterListName;
}
+ /**
+ * Determine whether a particular cluster is a fragment.
+ * The algorithm is to identify the dominant MCParticle
+ * contributing to [clus] and then determine whether
+ * there is another cluster in the previously specified
+ * list of Cluster objects which has a larger contribution
+ * from the same MCParticle. If so, returns true (fragment).
+ *
+ * The cluster [clus] MUST be a member of the list of clusters.
+ *
+ * @parameter clus Cluster to test
+ * @parameter event EventHeader, from which we pull the list of clusters
+ */
public boolean isFragment(Cluster clus, EventHeader event)
{
List<Cluster> clusterList = event.get(Cluster.class, m_clusterListName);
- if (! clusterList.contains(clus)) {
- throw new AssertionError("Cluster not in list");
- }
+ if (! clusterList.contains(clus)) { throw new AssertionError("Cluster not in list"); }
- //System.out.println("DEBUG: Cluster list of size "+clusterList.size()+" contains cluster ["+clus+"]");
- //for (Cluster tmpClus : clusterList) {
- //System.out.println("DEBUG: Cluster list entry: ["+tmpClus+"] with "+tmpClus.getCalorimeterHits().size()+" hits");
- //}
+ // Which is the dominant MC particle?
+ MCParticle alphaParticle = findAlphaParticle(clus, event);
+ // OK, now is there another cluster to which it contributes more hits?
+ Cluster dominantCluster = findDominantCluster(alphaParticle, clusterList);
+
+ // It's a fragment if it's not the dominant cluster:
+ return (dominantCluster != clus);
+ }
+
+ /**
+ * A small utility routine. If [clus] is a fragment,
+ * returns its parent cluster. Otherwise, returns [clus].
+ */
+ public Cluster findParent(Cluster clus, EventHeader event)
+ {
+ if (isFragment(clus, event)) {
+ // Which is the dominant MC particle?
+ MCParticle alphaParticle = findAlphaParticle(clus, event);
+
+ // OK, now is there another cluster to which it contributes more hits?
+ List<Cluster> clusterList = event.get(Cluster.class, m_clusterListName);
+ Cluster dominantCluster = findDominantCluster(alphaParticle, clusterList);
+
+ // This is supposed to be a fragment, so the dominant cluster
+ // had better not be the same as clus.
+ if (clus == dominantCluster) { throw new AssertionError("BUG: Fragment, but is dominant cluster"); }
+
+ return dominantCluster; // the parent
+ } else {
+ return clus;
+ }
+ }
+
+
+ protected Cluster findDominantCluster(MCParticle alphaParticle, List<Cluster> clusterList)
+ {
+ int maxHitCount = 0;
+ Cluster bestCluster = null;
+ for (Cluster otherClus : clusterList) {
+ int hitCount = countHits(otherClus, alphaParticle);
+ if (hitCount > maxHitCount || bestCluster == null) {
+ maxHitCount = hitCount;
+ bestCluster = otherClus;
+ }
+ }
+ return bestCluster;
+ }
+
+ protected MCParticle findAlphaParticle(Cluster clus, EventHeader event)
+ {
+ MCParticle alphaParticle = null;
+ List<Cluster> clusterList = event.get(Cluster.class, m_clusterListName);
+ if (! clusterList.contains(clus)) { throw new AssertionError("Cluster not in list"); }
HitCountAssociator assoc = new HitCountAssociator(event);
assoc.setClusters(clusterList);
- // Which is the dominant MC particle?
MCParticle testAlphaParticle = assoc.associateClusterToMCParticles(clus).iterator().next();
- MCParticle alphaParticle = null;
int alphaHits = 0;
for (MCParticle part : assoc.associateClusterToMCParticles(clus)) {
int countHits = countHits(clus, part);
@@ -51,24 +117,8 @@
System.out.println("BUG: ["+tmpPart+"]: "+tmpPart.getType().getName()+" with energy "+tmpPart.getEnergy()+" and relevant hits = "+countHits(clus, tmpPart));
}
throw new AssertionError("BUG! Compare alphaParticle ["+alphaParticle+"] which is "+alphaParticle.getType().getName()+" with energy "+alphaParticle.getEnergy()+" vs testAlphaParticle ["+testAlphaParticle+"] which is "+testAlphaParticle.getType().getName()+" with energy "+testAlphaParticle.getEnergy()+". alphaHits="+alphaHits+", but hits from testAlphaParticle are "+countHits(clus, testAlphaParticle)) ;
- }
-
- // OK, now is there another cluster to which it contributes more hits?
- for (Cluster otherClus : clusterList) {
- int hitCount = countHits(otherClus, alphaParticle);
- if (hitCount > alphaHits && otherClus != clus) {
- //System.out.println("DEBUG: CheatFragmentIdentifier: cluster ["+clus+"] with "+alphaHits+" hits is a fragment (alpha is "+alphaParticle.getType().getName()+" with energy "+alphaParticle.getEnergy()+"; beaten by a cluster ["+otherClus+"] with "+hitCount+" hits)");
- return true; // fragment
- } else if (hitCount <= alphaHits && otherClus != clus) {
- if (hitCount>0) {
- //System.out.println("DEBUG: CheatFragmentIdentifier: found another cluster with same alpha particle, but fewer (or equal) hits ("+hitCount+" vs "+alphaHits+")");
- }
- }
- }
-
- // ... no, there isn't.
- //System.out.println("CheatFragmentIdentifier: cluster is not a fragment -- alpha particle is "+alphaParticle.getType().getName()+" with energy "+alphaParticle.getEnergy()+". No other cluster in a list of "+clusterList.size()+" matched better.");
- return false; // not a fragment.
+ }
+ return alphaParticle;
}
protected int countHits (Cluster clus, MCParticle part)
lcsim/src/org/lcsim/contrib/uiowa/structural
diff -u -r1.1 -r1.2
--- TestFragmentIdentifier.java 24 Dec 2005 19:22:58 -0000 1.1
+++ TestFragmentIdentifier.java 26 Dec 2005 21:26:50 -0000 1.2
@@ -19,7 +19,7 @@
* Test the performance of a FragmentIdentifier by
* comparing it to another, cheating FragmentIdentifier.
*
- * @version $Id: TestFragmentIdentifier.java,v 1.1 2005/12/24 19:22:58 mcharles Exp $
+ * @version $Id: TestFragmentIdentifier.java,v 1.2 2005/12/26 21:26:50 mcharles Exp $
*/
public class TestFragmentIdentifier implements FragmentIdentifier
@@ -45,6 +45,11 @@
m_hnumHitsVsDocaWhenFragNon = m_histoFactory.createCloud2D("numHitsVsDocaWhenWrong (test=frag, cheat=nonfrag)");
m_hnumHitsVsDocaWhenNonFrag = m_histoFactory.createCloud2D("numHitsVsDocaWhenWrong (test=nonfrag, cheat=frag)");
m_hnumHitsVsDocaWhenNonNon = m_histoFactory.createCloud2D("numHitsVsDocaWhenCorrect (test=nonfrag, cheat=nonfrag)");
+ m_hdocaToParentWhenFragFrag = m_histoFactory.createCloud1D("docaToParentWhenCorrect (test=frag, cheat=frag)");
+ m_hdocaToParentWhenNonFrag = m_histoFactory.createCloud1D("docaToParentWhenWrong (test=nonfrag, cheat=frag)");
+ m_hdistToParentWhenFragFrag = m_histoFactory.createCloud1D("distToParentWhenCorrect (test=frag, cheat=frag)");
+ m_hdistToParentWhenNonFrag = m_histoFactory.createCloud1D("distToParentWhenWrong (test=nonfrag, cheat=frag)");
+
} catch (IOException ioe1) {
ioe1.printStackTrace();
}
@@ -97,6 +102,40 @@
if (numHits>=4) { m_hdocaWhenNonNon.fill(doca); m_hnumHitsVsDocaWhenNonNon.fill(numHits,doca); }
}
+ if (cheat_isFragment) {
+ if (m_cheatID instanceof CheatFragmentIdentifier) {
+ CheatFragmentIdentifier cheatID = (CheatFragmentIdentifier)(m_cheatID) ;
+ Cluster parent = cheatID.findParent(clus, event);
+ if (parent == clus) { throw new AssertionError("Fragment but no parent!"); }
+
+ {
+ Hep3Vector[] parentPosAndDir = MiscUtilities.getPositionAndDirection(parent);
+ Hep3Vector[] linePosAndDir = MiscUtilities.getPositionAndDirection(clus);
+ if (parentPosAndDir != null && linePosAndDir != null && parentPosAndDir[0] != null && linePosAndDir[0] != null) {
+ Hep3Vector displacement = VecOp.sub(parentPosAndDir[0], linePosAndDir[0]);
+ double distance = displacement.magnitude();
+ if (test_isFragment) {
+ m_hdistToParentWhenFragFrag.fill(distance);
+ } else {
+ m_hdistToParentWhenNonFrag.fill(distance);
+ }
+ }
+ }
+
+ int numHitsParent = parent.getCalorimeterHits().size();
+ if (numHits>=4 && numHitsParent>=4) {
+ Hep3Vector[] parentPosAndDir = MiscUtilities.getPositionAndDirection(parent);
+ Hep3Vector[] linePosAndDir = MiscUtilities.getPositionAndDirection(clus);
+ double docaToParent = MiscUtilities.findDOCAToPoint(linePosAndDir[0], linePosAndDir[1], parentPosAndDir[0]);
+ if (test_isFragment) {
+ m_hdocaToParentWhenFragFrag.fill(docaToParent);
+ } else {
+ m_hdocaToParentWhenNonFrag.fill(docaToParent);
+ }
+ }
+ }
+ }
+
try {
m_tree.commit();
} catch(IOException ioe1) {
@@ -126,5 +165,10 @@
ICloud2D m_hnumHitsVsDocaWhenNonFrag ;
ICloud2D m_hnumHitsVsDocaWhenNonNon ;
+ ICloud1D m_hdocaToParentWhenFragFrag;
+ ICloud1D m_hdocaToParentWhenNonFrag;
+ ICloud1D m_hdistToParentWhenFragFrag;
+ ICloud1D m_hdistToParentWhenNonFrag;
+
}
