4 added files
lcsim/src/org/lcsim/recon/cluster/mst
diff -N MSTPhotonFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MSTPhotonFinder.java 11 Aug 2006 23:31:25 -0000 1.1
@@ -0,0 +1,170 @@
+package org.lcsim.recon.cluster.mst;
+
+import java.util.*;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.recon.cluster.util.Clusterer;
+import org.lcsim.util.decision.*;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.recon.cluster.util.ClusterSizeDecision;
+import org.lcsim.recon.cluster.util.ClusterFirstLayerDecision;
+import org.lcsim.recon.cluster.util.ClusterLayerSeparationDecision;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+/**
+ * Class to reconstruct and select photon clusters in the ECAL.
+ * Conforms to the Clusterer interface.
+ *
+ * @version $Id: MSTPhotonFinder.java,v 1.1 2006/08/11 23:31:25 mcharles Exp $
+ */
+
+public class MSTPhotonFinder implements Clusterer
+{
+ /** Simple constructor */
+ public MSTPhotonFinder() {
+ m_shapeDec.setMinimum(null,null,null);
+ m_shapeDec.setMaximum(null,null,null);
+ }
+
+ // Config stuff
+ protected List<Cluster> m_MIPVetoList = null;
+ protected double m_coreThreshold = 7.5; // default (mm)
+ protected double m_fragThreshold = 75.0; // default (mm)
+ protected int m_layerProximityThreshold = 2; // default
+ protected int m_coreSizeMinimum = 10;
+ protected int m_fragSizeMaximum = 6;
+ protected int m_coreFirstLayerRange = 6;
+ protected DecisionMakerSingle<Cluster> m_inputDec = null;
+ protected PhotonShapeDecision m_shapeDec = new PhotonShapeDecision();
+
+ /** Set minimum 3D distance between hits in photon core */
+ public void setCoreThreshold(double newThreshold) { m_coreThreshold = newThreshold; }
+
+ /** Set minimum transverse distance between photon core center and fragment center */
+ public void setFragmentThreshold(double newThreshold) { m_fragThreshold = newThreshold; }
+
+ /**
+ * For a fragment to be added, require that there is at
+ * least one fragment hit within must be within [layerThreshold]
+ * layers of a core hit. Default is 2.
+ */
+ public void setLayerProximityThreshold(int layerThreshold) { m_layerProximityThreshold = layerThreshold; }
+
+ /**
+ * Set input hit decision (e.g. energy threshold) if desired.
+ * As per MSTClusterDriver interface, the hits are wrapped into
+ * one-hit clusters and a DecisionMakerSingle<Cluster> is applied.
+ * This applies to both cores and fragments.
+ */
+ public void setInputHitDecision(DecisionMakerSingle<Cluster> dec) { m_inputDec = dec; }
+
+ /** Cores must have at least [min] hits. */
+ public void setCoreSizeMinimum(int min) { m_coreSizeMinimum = min; }
+
+ /** Fragments must have at most [min] hits. */
+ public void setFragmentSizeMaximum(int max) { m_fragSizeMaximum = max; }
+
+ /** Cores must have a hit in the first [nLayers] layers. */
+ public void setCoreFirstLayerRange(int nLayers) { m_coreFirstLayerRange = nLayers; }
+
+ /** Ron's Clusterer interface */
+ public List<Cluster> createClusters(Map<Long,CalorimeterHit> map) {
+ List<CalorimeterHit> list = new Vector<CalorimeterHit>();
+ list.addAll(map.values());
+ return createClusters(list);
+ }
+
+ /**
+ * Find photon candidate clusters, checking for consistency
+ * with photon shape. The caller supplies a list of hits in
+ * the ECAL from which the MIPs have already been removed.
+ *
+ */
+ public List<Cluster> createClusters(List<CalorimeterHit> hits)
+ {
+ // First step: Find clusters with MST
+ // ----------------------------------
+
+ // Values are hard-coded for now
+ MSTClusterDriver mstClusterer = new MSTClusterDriver();
+ mstClusterer.setThreshold(m_coreThreshold);
+ mstClusterer.registerMetrics(new GeometricalDistance());
+ if (m_inputDec != null) {
+ mstClusterer.setInputDecision(m_inputDec);
+ }
+ List<Cluster> photonClusters = mstClusterer.createClusters(hits);
+
+ // Second step: Separate photon cores from fragments
+ // -------------------------------------------------
+
+ // Set up the decision-maker to identify valid photon cores
+ AndDecisionMakerSingle<Cluster> coreInputDecision = new AndDecisionMakerSingle<Cluster>();
+ coreInputDecision.addDecisionMaker(new ClusterSizeDecision(m_coreSizeMinimum)); // Require >= n hits
+ coreInputDecision.addDecisionMaker(new ClusterFirstLayerDecision(m_coreFirstLayerRange)); // Require first hit in first n layers
+ if (m_MIPVetoList != null) {
+ // NOT YET IMPLEMENTED
+ //List<Cluster> mipList = event.get(Cluster.class, m_MIPVetoList);
+ //coreInputDecision.addDecisionMaker(new ProximityDecision(mipList, 0.75)); // Require isolation w.r.t. MIPs
+ }
+ coreInputDecision.addDecisionMaker(m_shapeDec); // Require photon shape
+
+ // Set up the decision-maker to identify fragments.
+ // We should require that they are
+ // (a) inconsistent with being a core
+ // (b) consistent with being a fragment
+ AndDecisionMakerSingle<Cluster> fragInputDecision = new AndDecisionMakerSingle<Cluster>();
+ fragInputDecision.addDecisionMaker(new NotDecisionMakerSingle<Cluster> (coreInputDecision)); // Inconsistent with being a core
+ fragInputDecision.addDecisionMaker(new NotDecisionMakerSingle<Cluster> (new ClusterSizeDecision(m_fragSizeMaximum+1))); // Require < n+1 hits
+
+ // OK, now split into cores and non-cores
+ ListFilter findCores = new ListFilter(coreInputDecision);
+ ListFilter findFrags = new ListFilter(fragInputDecision);
+ List<Cluster> cores = findCores.filterList(photonClusters);
+ List<Cluster> frags = findFrags.filterList(photonClusters);
+ if (m_debug) { System.out.println("DEBUG: Found "+cores.size()+" cores and "+frags.size()+" fragments."); }
+
+ // Third step: Find fragments surrounding each photon
+ // --------------------------------------------------
+
+ DecisionMakerPair<Cluster,Cluster> pairDecision = new ClusterLayerSeparationDecision(m_layerProximityThreshold);
+ Metrics metric = new TransverseDistanceMetric();
+
+ Map<Cluster,Cluster> closestCore = new HashMap<Cluster,Cluster>();
+
+ // For each fragment, find the best-match core
+ for (Cluster frag : frags) {
+ // Loop over cores to find the best match
+ double dMin = -1;
+ Cluster bestCore = null;
+ for (Cluster core : cores) {
+ if (pairDecision.valid(core, frag)) {
+ double d = metric.getDistance(core, frag);
+ if (bestCore==null || d < dMin) {
+ // This is the best match so far
+ bestCore = core;
+ dMin = d;
+ }
+ }
+ }
+ if (bestCore != null && dMin < m_fragThreshold) {
+ // Found a valid match
+ closestCore.put(frag,bestCore);
+ }
+ }
+
+ // Fourth step: Add nearby fragments to core
+ // -----------------------------------------
+ for (Cluster frag : closestCore.keySet()) {
+ // Cores were produced by MST earlier, so should be BasicCluster
+ Cluster matchedCore = closestCore.get(frag);
+ BasicCluster basicMatchedCore = (BasicCluster) (matchedCore);
+ basicMatchedCore.addCluster(frag);
+ }
+
+ return cores;
+ }
+
+ protected boolean m_debug = false;
+ public void setDebug(boolean debug) { m_debug = debug; }
+}
lcsim/src/org/lcsim/recon/cluster/mst
diff -N MSTPhotonFinderDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MSTPhotonFinderDriver.java 11 Aug 2006 23:31:25 -0000 1.1
@@ -0,0 +1,101 @@
+package org.lcsim.recon.cluster.mst;
+
+import java.util.*;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.util.Driver;
+import org.lcsim.util.decision.DecisionMakerSingle;
+
+/**
+ * Driver wrapper around the MSTPhotonFinder class.
+ *
+ * @version $Id: MSTPhotonFinderDriver.java,v 1.1 2006/08/11 23:31:25 mcharles Exp $
+ */
+
+public class MSTPhotonFinderDriver extends Driver
+{
+ /** Simple constructor. */
+ public MSTPhotonFinderDriver() {
+ m_finder = new MSTPhotonFinder();
+ }
+
+ /** Processs one event, looking for photon candidates. */
+ public void process(EventHeader event) {
+ super.process(event);
+ HitMap hits = (HitMap) (event.get(m_inputHitMapName));
+ List<Cluster> photons = m_finder.createClusters(hits);
+ event.put(m_outputClusterListName, photons);
+ // Build output HitMap
+ HitMap outputHitMap = new HitMap(hits);
+ // Remove photon hits
+ for (Cluster clus : photons) {
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Long cellID = new Long(hit.getCellID());
+ outputHitMap.remove(cellID);
+ }
+ }
+ if (m_debug) { System.out.println("DEBUG: In this event, found "+photons.size()+" photons and wrote them out as '"+m_outputClusterListName+"'."); }
+ event.put(m_outputHitMapName, outputHitMap);
+ }
+
+ /** Read this HitMap in from the event and look for photons in those hits. */
+ public void setInputHitMap(String name) { m_inputHitMapName = name; }
+ /** Write out unused hits to the event under this name. */
+ public void setOutputHitMap(String name) { m_outputHitMapName = name; }
+ /** Write out photon clusters to the event under this name. */
+ public void setOutputClusterList(String name) { m_outputClusterListName = name ; }
+
+ protected MSTPhotonFinder m_finder;
+ protected String m_outputHitMapName;
+ protected String m_inputHitMapName;
+ protected String m_outputClusterListName;
+
+ // Config stuff
+
+ /** Set minimum 3D distance between hits in photon core */
+ public void setCoreThreshold(double newThreshold) {
+ m_finder.setCoreThreshold(newThreshold);
+ }
+ /** Set minimum transverse distance between photon core center and fragment center */
+ public void setFragmentThreshold(double newThreshold) {
+ m_finder.setFragmentThreshold(newThreshold);
+ }
+ /**
+ * For a fragment to be added, require that there is at
+ * least one fragment hit within must be within [layerThreshold]
+ * layers of a core hit. Default is 2.
+ */
+ public void setLayerProximityThreshold(int layerThreshold) {
+ m_finder.setLayerProximityThreshold(layerThreshold);
+ }
+ /**
+ * Set input hit decision (e.g. energy threshold) if desired.
+ * As per MSTClusterDriver interface, the hits are wrapped into
+ * one-hit clusters and a DecisionMakerSingle<Cluster> is applied.
+ * This applies to both cores and fragments.
+ */
+ public void setInputHitDecision(DecisionMakerSingle<Cluster> dec) {
+ m_finder.setInputHitDecision(dec);
+ }
+ /** Cores must have at least [min] hits. */
+ public void setCoreSizeMinimum(int min) {
+ m_finder.setCoreSizeMinimum(min);
+ }
+ /** Fragments must have at most [min] hits. */
+ public void setFragmentSizeMaximum(int max) {
+ m_finder.setFragmentSizeMaximum(max);
+ }
+ /** Cores must have a hit in the first [nLayers] layers. */
+ public void setCoreFirstLayerRange(int nLayers) {
+ m_finder.setCoreFirstLayerRange(nLayers);
+ }
+
+ protected boolean m_debug = false;
+ public void setDebug(boolean debug) {
+ m_debug = debug;
+ m_finder.setDebug(debug);
+ }
+}
lcsim/src/org/lcsim/recon/cluster/mst
diff -N TransverseDistanceMetric.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TransverseDistanceMetric.java 11 Aug 2006 23:31:25 -0000 1.1
@@ -0,0 +1,78 @@
+package org.lcsim.recon.cluster.mst;
+
+import java.io.*;
+import java.util.*;
+import hep.physics.vec.*;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+
+
+/**
+ * A distance metric that computes the transverse distance
+ * between a large cluster with well-defined direction
+ * (e.g. a MIP or a photon) and a fragment.
+ *
+ * The first cluster is the large one.
+ * Let its primary eigenvector be d1.
+ * Let its center-of-energy be at p1.
+ * The second cluster is the fragment.
+ * Let its center-of-energy be at p2.
+ *
+ * Let the angle between (d1) and (p2-p1) be theta:
+ * (d1).(p2-p1) = |d1| |p2-p1| cos(theta)
+ * [Note that there is an ambiguity between (theta) and (pi-theta).]
+ * Then the metric M is given by
+ * M^2 = |p2-p1|^2 sin^2(theta)
+ * = |p2-p1|^2 [1-cos^2(theta)]
+ * = |p2-p1|^2 - |p2-p1|^2 cos^2(theta)
+ * = |p2-p1|^2 - [(d1).(p2-p1)/|d1|]^2
+ */
+public class TransverseDistanceMetric implements Metrics {
+ /**
+ * Simple constructor
+ */
+ public TransverseDistanceMetric(){}
+
+ /**
+ * Compute the transverse distance from the large cluster to the
+ * other, as described above.
+ */
+ public double getDistance (Cluster cluster1, Cluster cluster2)
+ {
+ // Find the centers and direction:
+ Hep3Vector p1 = new BasicHep3Vector(cluster1.getPosition());
+ Hep3Vector p2 = new BasicHep3Vector(cluster2.getPosition());
+ Hep3Vector d1 = findDirection(cluster1);
+
+ // Vector arithmetic:
+ Hep3Vector p2_minus_p1 = VecOp.sub(p2,p1);
+ double p2_minus_p1_squared = VecOp.dot(p2_minus_p1, p2_minus_p1);
+ double d1_dot_p2_minus_p1 = VecOp.dot(d1, p2_minus_p1);
+ double d1_squared = VecOp.dot(d1,d1);
+
+ // Compute the metric with the formula above:
+ double m_squared = p2_minus_p1_squared - d1_dot_p2_minus_p1*d1_dot_p2_minus_p1/d1_squared ;
+ return Math.sqrt(m_squared);
+ }
+
+ private Hep3Vector findDirection(Cluster clus)
+ {
+ if (clus.getCalorimeterHits().size()<4) {
+ throw new AssertionError("Too few hits to calculate direction of cluster: there were "+clus.getCalorimeterHits().size()+", but I needed 4.");
+ }
+
+ 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 direction = new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]);
+ return direction;
+ }
+}
lcsim/src/org/lcsim/recon/cluster/mst
diff -N PhotonShapeDecision.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PhotonShapeDecision.java 11 Aug 2006 23:31:25 -0000 1.1
@@ -0,0 +1,89 @@
+package org.lcsim.recon.cluster.mst;
+
+import java.util.*;
+import hep.aida.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+import org.lcsim.util.decision.DecisionMakerSingle;
+import org.lcsim.event.CalorimeterHit;
+
+/**
+ * Decision-maker class which attempts to identify photons based
+ * on their normalized eigenvalues. Default is no cuts.
+ */
+
+public class PhotonShapeDecision implements DecisionMakerSingle<Cluster>
+{
+ /** Simple constructor. */
+ public PhotonShapeDecision() {
+ }
+
+ /** Main interface. */
+ public boolean valid(Cluster clus) {
+ calculateValues(clus);
+ for (int i=0; i<3; i++){
+ if (m_max[i] != null) {
+ if (m_value[i] > m_max[i].doubleValue()) {
+ return false;
+ }
+ }
+ if (m_min[i] != null) {
+ if (m_value[i] < m_min[i].doubleValue()) {
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+ /** Set maximum values for cuts on normalized eigenvectors. */
+ public void setMaximum(double max0, double max1, double max2) {
+ m_max[0] = new Double(max0);
+ m_max[1] = new Double(max1);
+ m_max[2] = new Double(max2);
+ }
+ /** Set maximum values for cuts on normalized eigenvectors. A null value means no cut. */
+ public void setMaximum(Double max0, Double max1, Double max2) {
+ m_max[0] = max0;
+ m_max[1] = max1;
+ m_max[2] = max2;
+ }
+
+ /** Set minimum values for cuts on normalized eigenvectors. */
+ public void setMinimum(double min0, double min1, double min2) {
+ m_min[0] = new Double(min0);
+ m_min[1] = new Double(min1);
+ m_min[2] = new Double(min2);
+ }
+ /** Set minimum values for cuts on normalized eigenvectors. A null value means no cut. */
+ public void setMinimum(Double min0, Double min1, Double min2) {
+ m_min[0] = min0;
+ m_min[1] = min1;
+ m_min[2] = min2;
+ }
+
+ /** Internal utility routine. */
+ protected void calculateValues(Cluster clus) {
+ if (clus.getCalorimeterHits().size() < 4) {
+ // Too few hits to calculate the tensor
+ throw new AssertionError(this.getClass().getName()+".calculateValues() called on a cluster with only "+clus.getCalorimeterHits().size()+" hits; need at least 4.");
+ }
+ BasicCluster copy = new BasicCluster();
+ copy.addCluster(clus);
+ TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+ copy.setPropertyCalculator(calc);
+ copy.calculateProperties();
+ double[] nev = calc.getNormalizedEigenvalues();
+ if (nev == null) {
+ throw new AssertionError("Normalized eigenvalues not calculated");
+ }
+ m_value[0] = nev[0];
+ m_value[1] = (nev[2]+nev[1])*0.5;
+ m_value[2] = (nev[2]-nev[1])*0.5;
+ }
+
+ protected Double[] m_max = {null,null,null};
+ protected Double[] m_min = {null,null,null};
+ protected double[] m_value = {0,0,0};
+}
CVSspam 0.2.8
