Print

Print


Commit in lcsim/src/org/lcsim/recon/cluster on MAIN
cheat/CheatClusterDriver.java+51added 1.1
     /CheatClusterer.java+59added 1.1
fixedcone/FixedConeClusterDriver.java+320added 1.1
nn/NearestNeighborClusterDriver.java+109added 1.1
+539
4 added files
Refactored clustering code.

lcsim/src/org/lcsim/recon/cluster/cheat
CheatClusterDriver.java added at 1.1
diff -N CheatClusterDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ CheatClusterDriver.java	18 Jul 2005 03:11:09 -0000	1.1
@@ -0,0 +1,51 @@
+package org.lcsim.recon.cluster.cheat;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.util.Driver;
+
+
+/**
+ * Driver to create and write CheatClusters to the event
+ */
+
+public class CheatClusterDriver extends Driver
+{
+    CheatClusterer _clusterer;
+    
+    public CheatClusterDriver()
+    {
+        _clusterer = new CheatClusterer();
+    }
+    
+    public void process(EventHeader event)
+    {
+        //First look for clusters in individual collections
+        
+        List<List<SimCalorimeterHit>> collections = event.get(SimCalorimeterHit.class);
+        for (List<SimCalorimeterHit> collection : collections)
+        {
+            Map<MCParticle,CheatCluster> result = _clusterer.findClusters(collection);
+            String name = event.getMetaData(collection).getName();
+            if (result.size() > 0) event.put(name+"Clusters",new ArrayList(result.values()));
+        }
+        
+        // Then look for refined clusters combining hits from all collections
+        
+        List<List<CheatCluster>> clusters = event.get(CheatCluster.class);
+        Map<MCParticle, CheatCluster> refined = _clusterer.findRefinedClusters(clusters);
+        if (refined.size() > 0) event.put("RefinedClusters",new ArrayList(refined.values()));
+    }
+    
+    public String toString()
+    {
+        return "CheatClusterDriver";
+    }
+}

lcsim/src/org/lcsim/recon/cluster/cheat
CheatClusterer.java added at 1.1
diff -N CheatClusterer.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ CheatClusterer.java	18 Jul 2005 03:11:09 -0000	1.1
@@ -0,0 +1,59 @@
+package org.lcsim.recon.cluster.cheat;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+
+
+/**
+ * The Cluster cheater works by finding perfectly reconstructed clusters
+ * based by looking at the MC truth information associated with each hit.
+ *
+ * To simplify the coding for now, if more than one MC particle is associated with a
+ * single hit, the entire hit is associated with the cluster belonging to the MC
+ * particle that contributed the most energy to the hit. Note that there is currently
+ * no attempt to merge clusters that would in fact be impossible to resolve.
+ */
+
+public class CheatClusterer
+{
+    
+    public Map<MCParticle, CheatCluster> findClusters(List<SimCalorimeterHit> hits)
+    {
+        Map<MCParticle,CheatCluster> result = new HashMap<MCParticle,CheatCluster>();
+        for (SimCalorimeterHit hit : hits)
+        {
+            double eMax = hit.getContributedEnergy(0);
+            MCParticle pMax = hit.getMCParticle(0);
+            for (int i=1; i<hit.getMCParticleCount(); i++)
+            {
+                if (hit.getContributedEnergy(i) > eMax) pMax = hit.getMCParticle(i);
+            }
+            CheatCluster cc = result.get(pMax);
+            if (cc == null) result.put(pMax,cc = new CheatCluster(pMax));
+            cc.addHit(hit);
+        }
+        return result;
+    }
+    public Map<MCParticle, CheatCluster> findRefinedClusters(List<List<CheatCluster>> collections)
+    {
+        Map<MCParticle, CheatCluster> result = new HashMap<MCParticle,CheatCluster>();
+        for (List<CheatCluster> clusters : collections)
+        {
+            for (CheatCluster cluster : clusters)
+            {
+                MCParticle p = cluster.getMCParticle();
+                CheatCluster rc = result.get(p);
+                if (rc == null) result.put(p,rc = new CheatCluster(p));
+                rc.addCluster(cluster);
+            }
+        }
+        return result;
+    }
+    
+}

lcsim/src/org/lcsim/recon/cluster/fixedcone
FixedConeClusterDriver.java added at 1.1
diff -N FixedConeClusterDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ FixedConeClusterDriver.java	18 Jul 2005 03:11:09 -0000	1.1
@@ -0,0 +1,320 @@
+package org.lcsim.recon.cluster.fixedcone;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.event.util.CalorimeterHitEsort;
+import org.lcsim.geometry.CalorimeterIDDecoder;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.fourvec.Lorentz4Vector;
+import org.lcsim.util.fourvec.Momentum4Vector;
+
+/**
+ * FixedConeClusterer implements a
+ * <font face="symbol">q-f </font> cone clustering algorithm
+ * that assigns all neighboring hits to the same cluster if they fall
+ * within a radius R of the cluster axis. The axis is originally defined
+ * by a seed cell, and is iteratively updated as cells are added.
+ * This version of the ClusterBuilder splits overlapping clusters
+ * by assigning cells in the overlap region to the nearest cluster axis.
+ *
+ * @author Norman A. Graf
+ * @version 1.0
+ */
+
+public class FixedConeClusterDriver extends Driver
+{
+    private double _radius;
+    private double _seedEnergy;
+    private double _minEnergy;
+    private int _numLayers;
+    private double _samplingFraction;
+    private double[] _layerEnergy;
+    private FixedConeClusterer _clusterer;
+    
+    CalorimeterIDDecoder _decoder;
+    
+    private static final double PI=Math.PI;
+    private static final double TWOPI = 2.*PI;
+    
+    /**
+     * Constructor
+     *
+     * @param   radius The cone radius in <font face="symbol">q-f </font> space
+     * @param   seed   The minimum energy for a cone seed cell (in GeV)
+     * @param   minE   The minimum energy for a cluster (in GeV)
+     */
+    public FixedConeClusterDriver(double radius, double seed, double minE)
+    {
+        _radius = radius;
+        // overwrite later with sampling fraction correction
+        _seedEnergy = seed;
+        _minEnergy = minE;
+        _clusterer = new FixedConeClusterer(_radius, _seedEnergy, _minEnergy);
+    }
+    
+    
+    /**
+     * Processes an Event to find CalorimeterClusters
+     *
+     * @param   event  The Event to process
+     */
+    public void process(EventHeader event)
+    {
+        List<CalorimeterHit> collection = event.get(CalorimeterHit.class,"EcalBarrHits");
+        _decoder = (CalorimeterIDDecoder) event.getMetaData(collection).getIDDecoder();
+//        System.out.println(_decoder);
+        List<BasicCluster> clusters = _clusterer.makeClusters(collection,_decoder);
+        
+        String name = event.getMetaData(collection).getName();
+        
+//        System.out.println("found "+clusters.size()+" clusters");
+        if (clusters.size() > 0) event.put(name+"FixedConeClusters",clusters);    
+    }
+    
+    /**
+     * Takes a list of CalorimeterHits, clusters them, and returns a list of
+     * CalorimeterClusters
+     *
+     * @param   in  A collection of CalorimeterHits
+     * @return  List of found CalorimeterClusters
+     */
+ /*   
+    public List cluster(List<? extends CalorimeterHit> in, CalorimeterIDDecoder decoder)
+    {
+        List out = new ArrayList();
+        
+        double rsquared = _radius*_radius;
+        // sort the vector in descending energy for efficiency
+        // this starts with the highest energy seeds.
+        Collections.sort(in, new CalorimeterHitEsort());
+        
+        int nclus = 0;
+        int size = in.size();
+        boolean[] used = new boolean[size];
+        //   outer loop finds a seed
+        for(int i =0; i<size; ++i)
+        {
+            if (!used[i])
+            {
+                CalorimeterHit p = in.get(i);
+                if (p.getEnergy()>_seedEnergy)
+                {
+                    decoder.setID(p.getCellID());
+                    double cellE = p.getEnergy();
+                    double px = cellE*Math.cos(decoder.getPhi())*Math.sin(decoder.getTheta());
+                    double py = cellE*Math.sin(decoder.getPhi())*Math.sin(decoder.getTheta());
+                    double pz = cellE*Math.cos(decoder.getTheta());
+                    Lorentz4Vector sum = new Momentum4Vector(px,py,pz,cellE);
+                    double phiseed=sum.phi();
+                    double thetaseed=sum.theta();
+                    
+                    // constituent cells
+                    List<CalorimeterHit> members = new ArrayList<CalorimeterHit>();
+                    members.add(p);
+                    // inner loop adds neighboring cells to seed
+                    
+                    for (int j = i+1; j<size; ++j)
+                    {
+                        //                        if (in.get(j)!=null)
+                        if(!used[j])
+                        {
+                            CalorimeterHit p2 = in.get(j);
+                            decoder.setID(p2.getCellID());
+                            double phi = decoder.getPhi();
+                            double theta=decoder.getTheta();
+                            double dphi=phi-phiseed;
+                            
+                            if(dphi<-PI) dphi+=TWOPI;
+                            if(dphi>PI) dphi-=TWOPI;
+                            double dtheta = theta-thetaseed;
+                            double R2=dphi*dphi+dtheta*dtheta;
+                            if(R2< rsquared)
+                            {
+                                //  particle within cone
+                                cellE = p2.getEnergy();
+                                px = cellE*Math.cos(decoder.getPhi())*Math.sin(decoder.getTheta());
+                                py = cellE*Math.sin(decoder.getPhi())*Math.sin(decoder.getTheta());
+                                pz = cellE*Math.cos(decoder.getTheta());
+                                sum.plusEquals(new Momentum4Vector(px, py, pz, cellE));
+                                members.add(p2);
+                                // tag this element so we don't reuse it
+                                used[j]=true;
+                                
+                                //   recalculate cone center
+                                phiseed=sum.phi();
+                                thetaseed=sum.theta();
+                            }
+                        }
+                    }// end of inner loop
+                    
+                    
+                    // if energy of cluster is large enough add it to the list
+                    if(sum.E()>_minEnergy)
+                    {
+                        //                   HitsCluster clus = new HitsCluster(decoder,sum,members,_numLayers,_samplingFraction);
+                        CalorimeterCluster clus = new CalorimeterCluster(decoder,sum, members,30,1.);
+                        for(CalorimeterHit hit : members)
+                        {
+                            clus.addHit(hit);
+                        }
+                        out.add(clus);
+                        nclus++;
+                    }
+                    
+                }
+            }// end of outer loop
+        }
+        
+        if (nclus>1)
+        {
+            // sort the clusters in descending energy
+            Collections.sort(out, new HitsClusterESort());
+            // loop over the found clusters and look for overlaps
+            // i.e distance between clusters is less than 2*R
+            for(int i=0; i<out.size();++i  )
+            {
+                for(int j=i+1; j<out.size(); ++j)
+                {
+                    double dTheta = dTheta((CalorimeterCluster)out.get(i), (CalorimeterCluster)out.get(j));
+                    if (dTheta<2*_radius)
+                    {
+                        resolve((CalorimeterCluster)out.get(i), (CalorimeterCluster)out.get(j), decoder);
+                    }
+                }
+            }
+        }
+        
+//        System.out.println("found "+out.size()+ " clusters");
+        return out;
+    }
+*/    
+    
+    /**
+     * Calculate the angle between two CalorimeterClusters
+     *
+     * @param   c1  First CalorimeterCluster
+     * @param   c2  Second CalorimeterCluster
+     * @return     The angle between the two clusters
+     */
+/*    
+    public double dTheta(CalorimeterCluster c1, CalorimeterCluster c2)
+    {
+        Lorentz4Vector v1 = c1.vector();
+        Lorentz4Vector v2 = c2.vector();
+        double costheta = (v1.vec3dot(v2))/(v1.p()*v2.p());
+        return Math.acos(costheta);
+    }
+*/    
+    /**
+     * Given two overlapping clusters, assign cells to nearest axis.
+     * The cluster axes are <em> not </em> iterated.
+     * Cluster quantities are recalculated after the split.
+     *
+     * @param   c1 First CalorimeterCluster
+     * @param   c2 Second CalorimeterCluster
+     */
+  /*  
+    public void resolve(CalorimeterCluster c1, CalorimeterCluster c2, CalorimeterIDDecoder decoder)
+    {
+        // do not recalculate cluster axis until all reshuffling is done
+        // do not want the cones to shift
+        // this behavior may change in the future
+        
+        Lorentz4Vector v1 = c1.vector();
+        double phi1 = v1.phi();
+        double theta1 = v1.theta();
+        Lorentz4Vector v2 = c2.vector();
+        double phi2 = v2.phi();
+        double theta2 = v2.theta();
+        List<CalorimeterHit> cells1 = c1.cells();
+        List<CalorimeterHit> cells2 = c2.cells();
+        int size2 = cells2.size();
+        
+        List<CalorimeterHit> tmp = new ArrayList<CalorimeterHit>();
+        // loop over cells in first cluster...
+        for(int i=0; i<cells1.size(); ++i)
+        {
+            CalorimeterHit p2 = (CalorimeterHit) cells1.get(i);
+            decoder.setID(p2.getCellID());
+            double phi = decoder.getPhi();
+            double theta = decoder.getTheta();
+            //distance to cluster 1
+            double dphi1=phi-phi1;
+            if(dphi1<-PI) dphi1+=TWOPI;
+            if(dphi1>PI) dphi1-=TWOPI;
+            double dtheta1 = theta-theta1;
+            double R1=dphi1*dphi1+(dtheta1)*(dtheta1);
+            // distance to cluster 2
+            double dphi2=phi-phi2;
+            if(dphi2<-PI) dphi2+=TWOPI;
+            if(dphi2>PI) dphi2-=TWOPI;
+            double dtheta2 = theta-theta2;
+            double R2=dphi2*dphi2+(dtheta2)*(dtheta2);
+            
+            if (R2<R1)
+            {
+                //add this cell to the vector of cells to remove
+                tmp.add(p2);
+            }
+        }
+        for(int i=0; i<tmp.size(); ++i)
+        {
+            CalorimeterHit p2 = (CalorimeterHit) tmp.get(i);
+            // remove this cell from cluster1
+            cells1.remove(p2);
+            // add it to cluster2
+            cells2.add(p2);
+        }
+        tmp.clear();
+        // repeat for cluster 2
+        // only loop over cells in original cluster
+        for(int i=0; i<size2; ++i)
+        {
+            CalorimeterHit p2 = (CalorimeterHit) cells2.get(i);
+            decoder.setID(p2.getCellID());
+            double phi = decoder.getPhi();
+            double theta = decoder.getTheta();
+            //distance to cluster 1
+            double dphi1=phi-phi1;
+            if(dphi1<-PI) dphi1+=TWOPI;
+            if(dphi1>PI) dphi1-=TWOPI;
+            double dtheta1 = theta-theta1;
+            double R1=dphi1*dphi1+(dtheta1)*(dtheta1);
+            // distance to cluster 2
+            double dphi2=phi-phi2;
+            if(dphi2<-PI) dphi2+=TWOPI;
+            if(dphi2>PI) dphi2-=TWOPI;
+            double dtheta2 = theta-theta2;
+            double R2=dphi2*dphi2+(dtheta2)*(dtheta2);
+            
+            if (R1<R2)
+            {
+                //add this cell to the vector of cells to remove
+                tmp.add(p2);
+            }
+        }
+        for(int i=0; i<tmp.size(); ++i)
+        {
+            CalorimeterHit p2 = (CalorimeterHit) tmp.get(i);
+            // remove this cell from cluster2
+            cells2.remove(p2);
+            // add it to cluster1
+            cells1.add(p2);
+        }
+        c1.calculateVec();
+        c2.calculateVec();
+    }
+   */
+    
+    public String toString()
+    {
+        return "FixedConeClusterDriver with clusterer "+_clusterer;
+    }
+}

lcsim/src/org/lcsim/recon/cluster/nn
NearestNeighborClusterDriver.java added at 1.1
diff -N NearestNeighborClusterDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ NearestNeighborClusterDriver.java	18 Jul 2005 03:11:09 -0000	1.1
@@ -0,0 +1,109 @@
+package org.lcsim.recon.cluster.nn;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.geometry.CalorimeterIDDecoder;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.Driver;
+
+/**
+ * A simple driver which applies a nearest neighbor clustering algorithm
+ * to calorimeter cell collections contained in the event.
+ *
+ * @author Norman A. Graf with additions of threshold code by Eric J. Benavidez ([log in to unmask])
+ */
+public class NearestNeighborClusterDriver extends Driver
+{
+   // the minimum number of cells to qualify as a cluster
+   private int _minNcells;
+   
+   // the neighborhood in u to search
+   private int _dU;
+   // the neighborhood in v to search
+   private int _dV;
+   // the neighborood in layers to search
+   private int _dLayer;
+   // energy threshold
+   private double _thresh;
+   // Create clusters class
+   private NearestNeighborClusterer _clusterer;
+   
+   /**
+    * Creates a new instance of NearestNeighborClusterDriver with a domain of (1,1,1)
+    * and no threshold.
+    * @param ncells The minimum number of cells required to constitute a cluster.
+    */
+   public NearestNeighborClusterDriver(int ncells)
+   {
+      this(1, 1, 1, ncells, 0);
+   }
+   
+   /**
+    * Creates a new instance of NearestNeighborClusterDriver
+    * @param dU  The extent in the U coordinate which defines a neighborhood.
+    * @param dV  The extent in the V coordinate which defines a neighborhood.
+    * @param dLayer  The extent in the layer which defines a neighborhood.
+    * @param ncells The minimum number of cells required to constitute a cluster.
+    * @param threshold The energy threshold that must be met or exceeded in order
+    *   to be considered in the cluster.
+    */
+   public NearestNeighborClusterDriver(int dU, int dV, int dLayer, int ncells,
+            double threshold)
+   {
+      _dU = dU;
+      _dV = dV;
+      _dLayer = dLayer;
+      _minNcells = ncells;
+      _thresh = threshold;
+      _clusterer = new NearestNeighborClusterer(_dU,_dV,_dLayer,_minNcells,_thresh);
+   }
+   
+   /**
+    * Creates a new NearestNeighborClusterDriver with no energy threshold applied.
+    * @param dU  The extent in the U coordinate which defines a neighborhood.
+    * @param dV  The extent in the V coordinate which defines a neighborhood.
+    * @param dLayer  The extent in the layer which defines a neighborhood.
+    * @param ncells The minimum number of cells required to constitute a cluster.
+    */
+   
+   public NearestNeighborClusterDriver(int dU, int dV, int dLayer, int ncells)
+   {
+       this(dU, dV, dLayer, ncells, 0);
+   }
+   
+   /**
+    * Process an event. Extract calorimeter cells from collections contained in the
+    * event, cluster the cells using a nearest-neighbor algorithm, and add these
+    * cluster containers back to the event.
+    *
+    * @param event
+    */
+   public void process(EventHeader event)
+   {
+      //First look for clusters in individual collections
+      
+      List<List<CalorimeterHit>> collections = event.get(CalorimeterHit.class);
+      for (List<CalorimeterHit> collection: collections)
+      {
+         CalorimeterIDDecoder decoder = (CalorimeterIDDecoder) event.getMetaData(collection).getIDDecoder();
+         List<BasicCluster> clusters = _clusterer.findClusters(collection,decoder);
+         
+         String name = event.getMetaData(collection).getName();
+         if (clusters.size() > 0)
+             event.put(name+"NNClusters",clusters);
+      }
+   }
+   
+   public String toString()
+   {
+      return "NearestNeighborClusterDriver with clusterer "+_clusterer;
+   }
+   
+}
CVSspam 0.2.8