5 removed + 5 modified, total 10 files
lcsim/src/org/lcsim/recon/cluster/cheat
diff -N ClusterCheater.java
--- ClusterCheater.java 8 Jul 2005 17:15:05 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,45 +0,0 @@
-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 ClusterCheater extends Driver
-{
- CreateCheatClusters ccc;
- public ClusterCheater()
- {
- ccc = new CreateCheatClusters();
- }
-
- 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 = ccc.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 = ccc.findRefinedClusters(clusters);
- if (refined.size() > 0) event.put("RefinedClusters",new ArrayList(refined.values()));
- }
-}
lcsim/src/org/lcsim/recon/cluster/cheat
diff -N CreateCheatClusters.java
--- CreateCheatClusters.java 8 Jul 2005 17:15:05 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,59 +0,0 @@
-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 CreateCheatClusters
-{
-
- 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
diff -u -r1.4 -r1.5
--- FixedConeClusterer.java 8 Jul 2005 17:15:06 -0000 1.4
+++ FixedConeClusterer.java 18 Jul 2005 03:06:18 -0000 1.5
@@ -2,16 +2,13 @@
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.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.ClusterESort;
+import org.lcsim.recon.cluster.util.FixedConeClusterPropertyCalculator;
import org.lcsim.util.fourvec.Lorentz4Vector;
import org.lcsim.util.fourvec.Momentum4Vector;
@@ -28,58 +25,334 @@
* @version 1.0
*/
-public class FixedConeClusterer extends Driver
+public class FixedConeClusterer
{
private double _radius;
- private double _Eseed;
- private double _Emin;
+ private double _seedEnergy;
+ private double _minEnergy;
private int _numLayers;
private double _samplingFraction;
private double[] _layerEnergy;
- private CreateFCClusters clusterer;
- CalorimeterIDDecoder _decoder;
private static final double PI=Math.PI;
private static final double TWOPI = 2.*PI;
+ public FixedConeClusterPropertyCalculator _clusterPropertyCalculator;
+ private FixedConeDistanceMetric _dm;
+
+ public enum FixedConeDistanceMetric
+ { DOTPRODUCT, DPHIDCOSTHETA, DPHIDTHETA }
/**
* 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)
+ * @param distMetric The distance metric:
+ * 0 for dot-product method,
+ * 1 for R<sup>2</sup> = (d phi)<sup>2 + (d(cos theta))<sup>2</sup>
+ * (simplified version of dot-product method)
+ * 2 for R<sup>2</sup> = (d phi)<sup>2 + (d theta)<sup>2</sup>
+ * (backwards compatibility mode with old FixedConeClusterer implementation
+ * Anything else will be assumed to be 0.
*/
- public FixedConeClusterer(double radius, double seed, double minE)
+ public FixedConeClusterer(double radius, double seed, double minE, FixedConeDistanceMetric distMetric)
{
_radius = radius;
// overwrite later with sampling fraction correction
- _Eseed = seed;
- _Emin = minE;
- clusterer = new CreateFCClusters(_radius, _Eseed, _Emin);
- }
+ _seedEnergy = seed;
+ _minEnergy = minE;
+ _dm = distMetric;
+ }
+ /**
+ * Constructor with default distance metric (dot-product method)
+ *
+ * @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 FixedConeClusterer(double radius, double seed, double minE)
+ {
+ this(radius, seed, minE, FixedConeDistanceMetric.DOTPRODUCT);
+ }
+
+/*
+ public void setRadius(double radius)
+ {
+ _radius = radius;
+ }
+ public void setSeed(double seed)
+ {
+ _seedEnergy = seed;
+ }
+ public void setMinE(double minE)
+ {
+ _minEnergy = minE;
+ }
+ */
+
+ /**
+ * Make clusters from the input list
+ */
+ public List<BasicCluster> makeClusters(List<CalorimeterHit> in, CalorimeterIDDecoder decoder)
+ {
+ List<BasicCluster> out = new ArrayList<BasicCluster>();
+ _clusterPropertyCalculator = new FixedConeClusterPropertyCalculator(decoder);
+
+ 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 = 0;
+ boolean cond = false;
+
+ switch(_dm)
+ {
+ case DPHIDCOSTHETA: // R^2 = dphi^2 + d(cos theta)^2
+ double dcostheta = Math.cos(theta)-Math.cos(thetaseed);
+ R2 = dphi*dphi + dcostheta*dcostheta;
+ cond = (R2 < rsquared);
+ break;
+ case DPHIDTHETA: // R^2 = dphi^2 + dtheta^2
+ R2 = dphi*dphi + dtheta*dtheta;
+ cond = (R2 < rsquared);
+ break;
+ case DOTPRODUCT:
+ default: // dot product method
+ double dotp = Math.sin(theta)*Math.sin(thetaseed)*
+ (Math.sin(phi)*Math.sin(phiseed) +
+ Math.cos(phi)*Math.cos(phiseed)) +
+ Math.cos(theta)*Math.cos(thetaseed);
+ cond = (dotp > Math.cos(_radius));
+ break;
+ }
+
+ if(cond)
+ {
+ // 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)
+ {
+ BasicCluster clus = new BasicCluster();
+ clus.setPropertyCalculator(_clusterPropertyCalculator);
+ 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 ClusterESort());
+ // 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(out.get(i), out.get(j));
+ if (dTheta<2*_radius)
+ {
+ resolve(out.get(i), out.get(j), decoder);
+ }
+ }
+ }
+ }
+
+// System.out.println("found "+out.size()+ " clusters");
+ return out;
+ }
/**
- * Processes an Event to find CalorimeterClusters
+ * Calculate the angle between two Clusters
*
- * @param event The Event to process
+ * @param c1 First Cluster
+ * @param c2 Second Cluster
+ * @return The angle between the two clusters
*/
- public void process(EventHeader event)
+
+ public double dTheta(BasicCluster c1, BasicCluster c2)
+ {
+ _clusterPropertyCalculator.calculateProperties(c1.getCalorimeterHits());
+ Lorentz4Vector v1 = _clusterPropertyCalculator.vector();
+ _clusterPropertyCalculator.calculateProperties(c2.getCalorimeterHits());
+ Lorentz4Vector v2 = _clusterPropertyCalculator.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 Cluster
+ * @param c2 Second Cluster
+ */
+
+ public void resolve(BasicCluster c1, BasicCluster c2, CalorimeterIDDecoder decoder)
{
- 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);
+ // do not recalculate cluster axis until all reshuffling is done
+ // do not want the cones to shift
+ // this behavior may change in the future
- String name = event.getMetaData(collection).getName();
- if (clusters.size() > 0) event.put(name+"FCClusters",clusters);
+ _clusterPropertyCalculator.calculateProperties(c1.getCalorimeterHits());
+ Lorentz4Vector v1 = _clusterPropertyCalculator.vector();
+ double phi1 = v1.phi();
+ double theta1 = v1.theta();
+ _clusterPropertyCalculator.calculateProperties(c2.getCalorimeterHits());
+ Lorentz4Vector v2 = _clusterPropertyCalculator.vector();
+ double phi2 = v2.phi();
+ double theta2 = v2.theta();
+ List<CalorimeterHit> cells1 = c1.getCalorimeterHits();
+ List<CalorimeterHit> cells2 = c2.getCalorimeterHits();
+ 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.calculateProperties();
+ c2.calculateProperties();
}
-
+
public String toString()
{
- return "FixedConeClusterer with radius "+_radius+" seed Energy "+_Eseed+" minimum energy "+_Emin;
+ return "FixedConeClusterer with radius "+_radius+" seed Energy "+_seedEnergy+" minimum energy "+_minEnergy+"distance metric "+_dm;
}
}
lcsim/src/org/lcsim/recon/cluster/fixedcone
diff -N CreateFCClusters.java
--- CreateFCClusters.java 8 Jul 2005 17:15:06 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,295 +0,0 @@
-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 CreateFCClusters
-{
- private double _radius;
- private double _Eseed;
- private double _Emin;
- private int _numLayers;
- private double _samplingFraction;
- private double[] _layerEnergy;
-
- CalorimeterIDDecoder _decoder;
-
- private static final double PI=Math.PI;
- private static final double TWOPI = 2.*PI;
- public FCClusterPropertyCalculator cpc;
-
- /**
- * 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 CreateFCClusters(double radius, double seed, double minE)
- {
- _radius = radius;
- // overwrite later with sampling fraction correction
- _Eseed = seed;
- _Emin = minE;
- }
- public void setRadius(double radius)
- {
- _radius = radius;
- }
- public void setSeed(double seed)
- {
- _Eseed = seed;
- }
- public void setMinE(double minE)
- {
- _Emin = minE;
- }
-
-
- /**
- * Make clusters from the input list
- */
- public List<BasicCluster> makeClusters(List<CalorimeterHit> in, CalorimeterIDDecoder decoder)
- {
- List<BasicCluster> out = new ArrayList<BasicCluster>();
- cpc = new FCClusterPropertyCalculator(decoder);
-
- 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()>_Eseed)
- {
- 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()>_Emin)
- {
- BasicCluster clus = new BasicCluster();
- clus.setPropertyCalculator(cpc);
- 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 ClusterESort());
- // 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((BasicCluster)out.get(i), (BasicCluster)out.get(j));
- if (dTheta<2*_radius)
- {
-// resolve((BasicCluster)out.get(i), (BasicCluster)out.get(j));
- }
- }
- }
- }
-
-// System.out.println("found "+out.size()+ " clusters");
- return out;
- }
-
-
- /**
- * Calculate the angle between two Clusters
- *
- * @param c1 First Cluster
- * @param c2 Second Cluster
- * @return The angle between the two clusters
- */
-
- public double dTheta(BasicCluster c1, BasicCluster c2)
- {
- cpc.calculateProperties(c1.getCalorimeterHits());
- Lorentz4Vector v1 = cpc.vector();
- cpc.calculateProperties(c2.getCalorimeterHits());
- Lorentz4Vector v2 = cpc.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 Cluster
- * @param c2 Second Cluster
- */
-
- public void resolve(BasicCluster c1, BasicCluster c2)
- {
- // do not recalculate cluster axis until all reshuffling is done
- // do not want the cones to shift
- // this behavior may change in the future
-
- cpc.calculateProperties(c1.getCalorimeterHits());
- Lorentz4Vector v1 = cpc.vector();
- double phi1 = v1.phi();
- double theta1 = v1.theta();
- cpc.calculateProperties(c2.getCalorimeterHits());
- Lorentz4Vector v2 = cpc.vector();
- double phi2 = v2.phi();
- double theta2 = v2.theta();
- List<CalorimeterHit> cells1 = c1.getCalorimeterHits();
- List<CalorimeterHit> cells2 = c2.getCalorimeterHits();
- 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)
- {
- c1.removeHit(p2);
- c2.addHit(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)
- {
- c2.removeHit(p2);
- c1.addHit(p2);
- }
- }
- }
-
-
- public String toString()
- {
- return "FixedConeClusterer with radius "+_radius+" seed Energy "+_Eseed+" minimum energy "+_Emin;
- }
-}
lcsim/src/org/lcsim/recon/cluster/nn
diff -u -r1.5 -r1.6
--- NearestNeighborCluster.java 8 Jul 2005 17:15:07 -0000 1.5
+++ NearestNeighborCluster.java 18 Jul 2005 03:06:19 -0000 1.6
@@ -13,72 +13,72 @@
*/
public class NearestNeighborCluster extends BasicCluster
{
- /**
- * Construct a NearestNeighborCluster. Note that the constructor actually performs the
- * clustering, given a seed hit.
- * @param decoder The CellID decoder appropriate for this collection of hits.
- * @param hitmap The list of hit calorimeter cells available for clustering.
- * @param hit The seed calorimeter hit for this cluster.
- * @param key The key for the seed calorimeter hit.
- * @param threshold The energy threshold that must be met or exceeded in order
- * to be considered in the cluster.
- */
-
- public NearestNeighborCluster(CalorimeterIDDecoder decoder, Map<Long, CalorimeterHit> hitmap, CalorimeterHit hit, Long key, int dU, int dV, int dLayer, double threshold)
- {
- // start by adding this hit to the cluster
- addHit(hit);
-
- // set the decoder to this hit
- decoder.setID(key);
- // remove this hit from the map so it can't be used again
- hitmap.remove(key);
-
- // loop over the hits in the cluster and add all its neighbors
- // note that hits.size() grows as we add cells
- for(int i=0; i<hits.size(); ++i)
- {
- CalorimeterHit c = hits.get(i);
- decoder.setID(c.getCellID());
-
- long[] neighbors = decoder.getNeighbourIDs(dLayer, dU, dV);
- // loop over all neighboring cell Ids
- for (int j=0; j<neighbors.length; ++j)
- {
- CalorimeterHit h = hitmap.get(neighbors[j]);
- // is the neighboring cell id hit?
- // if so, does it meet or exceed threshold?
- if (h != null)
+ /**
+ * Construct a NearestNeighborCluster. Note that the constructor actually performs the
+ * clustering, given a seed hit.
+ * @param decoder The CellID decoder appropriate for this collection of hits.
+ * @param hitmap The list of hit calorimeter cells available for clustering.
+ * @param hit The seed calorimeter hit for this cluster.
+ * @param key The key for the seed calorimeter hit.
+ * @param threshold The energy threshold that must be met or exceeded in order
+ * to be considered in the cluster.
+ */
+
+ public NearestNeighborCluster(CalorimeterIDDecoder decoder, Map<Long, CalorimeterHit> hitmap, CalorimeterHit hit, Long key, int dU, int dV, int dLayer, double threshold)
+ {
+ // start by adding this hit to the cluster
+ addHit(hit);
+
+ // set the decoder to this hit
+ decoder.setID(key);
+ // remove this hit from the map so it can't be used again
+ hitmap.remove(key);
+
+ // loop over the hits in the cluster and add all its neighbors
+ // note that hits.size() grows as we add cells
+ for(int i=0; i<hits.size(); ++i)
+ {
+ CalorimeterHit c = hits.get(i);
+ decoder.setID(c.getCellID());
+
+ long[] neighbors = decoder.getNeighbourIDs(dLayer, dU, dV);
+ // loop over all neighboring cell Ids
+ for (int j=0; j<neighbors.length; ++j)
{
- if (h.getEnergy() >= threshold)
- addHit(h);
- hitmap.remove(neighbors[j]);
+ CalorimeterHit h = hitmap.get(neighbors[j]);
+ // is the neighboring cell id hit?
+ // if so, does it meet or exceed threshold?
+ if (h != null)
+ {
+ if (h.getEnergy() >= threshold)
+ addHit(h);
+ hitmap.remove(neighbors[j]);
+ }
}
- }
- }
- // calculateDerivedQuantities();
- }
-
- /**
- * Construct a NearestNeighborCluster. Note that the constructor actually performs the
- * clustering, given a seed hit, without considering an energy threshold.
- * @param decoder The CellID decoder appropriate for this collection of hits.
- * @param hitmap The list of hit calorimeter cells available for clustering.
- * @param hit The seed calorimeter hit for this cluster.
- * @param key The key for the seed calorimeter hit.
- */
- public NearestNeighborCluster(CalorimeterIDDecoder decoder, Map<Long, CalorimeterHit> hitmap, CalorimeterHit hit, Long key, int dU, int dV, int dLayer)
- {
- this(decoder, hitmap, hit, key, dU, dV, dLayer, 0.0);
- }
-
- /**
- * String representation of this object.
- * @return A String representation of this object
- */
-
- public String toString()
- {
- return "NearestNeighborCluster with "+hits.size()+ " calorimeter hits";
- }
+ }
+ // calculateDerivedQuantities();
+ }
+
+ /**
+ * Construct a NearestNeighborCluster. Note that the constructor actually performs the
+ * clustering, given a seed hit, without considering an energy threshold.
+ * @param decoder The CellID decoder appropriate for this collection of hits.
+ * @param hitmap The list of hit calorimeter cells available for clustering.
+ * @param hit The seed calorimeter hit for this cluster.
+ * @param key The key for the seed calorimeter hit.
+ */
+ public NearestNeighborCluster(CalorimeterIDDecoder decoder, Map<Long, CalorimeterHit> hitmap, CalorimeterHit hit, Long key, int dU, int dV, int dLayer)
+ {
+ this(decoder, hitmap, hit, key, dU, dV, dLayer, 0.0);
+ }
+
+ /**
+ * String representation of this object.
+ * @return A String representation of this object
+ */
+
+ public String toString()
+ {
+ return "NearestNeighborCluster with "+hits.size()+ " calorimeter hits";
+ }
}
lcsim/src/org/lcsim/recon/cluster/nn
diff -u -r1.6 -r1.7
--- NearestNeighborClusterer.java 8 Jul 2005 17:15:07 -0000 1.6
+++ NearestNeighborClusterer.java 18 Jul 2005 03:06:19 -0000 1.7
@@ -14,96 +14,81 @@
import org.lcsim.util.Driver;
/**
- * A simple driver which applies a nearest neighbor clustering algorithm
- * to calorimeter cell collections contained in the event.
- *
+ * Class to create a list of of Nearest Neighbor clusters
* @author Norman A. Graf with additions of threshold code by Eric J. Benavidez ([log in to unmask])
*/
-public class NearestNeighborClusterer extends Driver
+public class NearestNeighborClusterer
{
- // 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 CreateNNClusters cnnc;
-
- /**
- * Creates a new instance of NearestNeighborClusterer with a domain of (1,1,1)
- * and no threshold.
- * @param ncells The minimum number of cells required to constitute a cluster.
- */
- public NearestNeighborClusterer(int ncells)
- {
- this(1, 1, 1, ncells, 0);
- }
-
- /**
- * Creates a new instance of NearestNeighborClusterer
- * @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 NearestNeighborClusterer(int dU, int dV, int dLayer, int ncells,
- double threshold)
- {
- _dU = dU;
- _dV = dV;
- _dLayer = dLayer;
- _minNcells = ncells;
- _thresh = threshold;
- cnnc = new CreateNNClusters(_dU,_dV,_dLayer,_minNcells,_thresh);
- }
-
- /**
- * Creates a new NearestNeighborClusterer 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 NearestNeighborClusterer(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 = cnnc.findClusters(collection,decoder);
-
- String name = event.getMetaData(collection).getName();
- if (clusters.size() > 0)
- event.put(name+"NNClusters",clusters);
- }
- }
-
+ // 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;
+
+ /**
+ */
+ public NearestNeighborClusterer()
+ {
+ this(1, 1, 1, 1, 0.);
+ }
+ /**
+ */
+ public NearestNeighborClusterer(int dU, int dV, int dLayer, int ncells, double threshold)
+ {
+ _dU = dU;
+ _dV = dV;
+ _dLayer = dLayer;
+ _minNcells = ncells;
+ _thresh = threshold;
+ }
+ public void setGrid(int dU, int dV, int dLayer)
+ {
+ _dU = dU;
+ _dV = dV;
+ _dLayer = dLayer;
+ }
+ public void setMinNcells(int ncells)
+ {
+ _minNcells = ncells;
+ }
+ public void setThreshold(double threshold)
+ {
+ _thresh = threshold;
+ }
+ public List<BasicCluster> findClusters(List<CalorimeterHit> hitlist, CalorimeterIDDecoder decoder)
+ {
+ // create a map of cells keyed on their index
+ List<BasicCluster> clusters = new ArrayList();
+ Map<Long, CalorimeterHit> hitmap = new HashMap<Long, CalorimeterHit>();
+ for (CalorimeterHit hit : hitlist)
+ {
+ long id = hit.getCellID();
+ hitmap.put(id, hit);
+ }
+ while (!hitmap.isEmpty())
+ {
+ Long k = hitmap.keySet().iterator().next();
+ CalorimeterHit hit = hitmap.get(k);
+ // start a cluster, constructor will aggregate remaining hits unto itself
+ // NB : Must pass the threshold into the constructor, otherwise threshold
+ // will not be used.
+ NearestNeighborCluster nnclus = new NearestNeighborCluster(decoder, hitmap, hit, k, _dU, _dV, _dLayer, _thresh);
+ if(nnclus.getCalorimeterHits().size()>_minNcells)
+ {
+ clusters.add(nnclus);
+ }
+ }
+ return clusters;
+ }
+
public String toString()
{
return "NearestNeighborClusterer "+_dU+" "+_dV+" "+_dLayer+" minimum of "+_minNcells+" cells";
- }
-
-}
+ }
+}
\ No newline at end of file
lcsim/src/org/lcsim/recon/cluster/nn
diff -N CreateNNClusters.java
--- CreateNNClusters.java 8 Jul 2005 17:15:07 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,90 +0,0 @@
-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;
-
-/**
- * Class to create a list of of Nearest Neighbor clusters
- * @author Norman A. Graf with additions of threshold code by Eric J. Benavidez ([log in to unmask])
- */
-public class CreateNNClusters
-{
- // 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;
-
- /**
- */
- public CreateNNClusters()
- {
- this(1, 1, 1, 1, 0.);
- }
- /**
- */
- public CreateNNClusters(int dU, int dV, int dLayer, int ncells, double threshold)
- {
- _dU = dU;
- _dV = dV;
- _dLayer = dLayer;
- _minNcells = ncells;
- _thresh = threshold;
- }
- public void setGrid(int dU, int dV, int dLayer)
- {
- _dU = dU;
- _dV = dV;
- _dLayer = dLayer;
- }
- public void setMinNcells(int ncells)
- {
- _minNcells = ncells;
- }
- public void setThreshold(double threshold)
- {
- _thresh = threshold;
- }
- public List<BasicCluster> findClusters(List<CalorimeterHit> hitlist, CalorimeterIDDecoder decoder)
- {
- // create a map of cells keyed on their index
- List<BasicCluster> clusters = new ArrayList();
- Map<Long, CalorimeterHit> hitmap = new HashMap<Long, CalorimeterHit>();
- for (CalorimeterHit hit : hitlist)
- {
- long id = hit.getCellID();
- hitmap.put(id, hit);
- }
- while (!hitmap.isEmpty())
- {
- Long k = hitmap.keySet().iterator().next();
- CalorimeterHit hit = hitmap.get(k);
- // start a cluster, constructor will aggregate remaining hits unto itself
- // NB : Must pass the threshold into the constructor, otherwise threshold
- // will not be used.
- NearestNeighborCluster nnclus = new NearestNeighborCluster(decoder, hitmap, hit, k, _dU, _dV, _dLayer, _thresh);
- if(nnclus.getCalorimeterHits().size()>_minNcells)
- {
- clusters.add(nnclus);
- }
- }
- return clusters;
- }
-
-}
lcsim/src/org/lcsim/recon/cluster/util
diff -u -r1.2 -r1.3
--- BasicCluster.java 8 Jul 2005 17:15:08 -0000 1.2
+++ BasicCluster.java 18 Jul 2005 03:06:19 -0000 1.3
@@ -10,287 +10,287 @@
/**
* Default implementation of Cluster Interface for Simulation
- *
+ *
* @author cassell
*/
public class BasicCluster implements Cluster
{
- protected List<CalorimeterHit> hits = new ArrayList<CalorimeterHit>();
- protected List<Cluster> clusters = new ArrayList<Cluster>();
- protected List<Subdetector> detectors = new ArrayList<Subdetector>();
- public double raw_energy;
- public double corrected_energy;
- protected List<Double> subdetector_raw_energies = new ArrayList<Double>();
- protected List<Double> subdetector_corrected_energies = new ArrayList<Double>();
- protected List<Double> hit_energies = new ArrayList<Double>();
- protected ClusterPropertyCalculator cluster_property_calculator = new DefaultClusterPropertyCalculator();
- protected boolean needsPropertyCalculation = true;
- protected double[] position;
- protected double[] positionError;
- protected double iphi;
- protected double itheta;
- protected double[] directionError;
- protected double[] shapeParameters;
-
- /**
- * Add a CalorimeterHit to the cluster
- */
- public void addHit(CalorimeterHit hit)
- {
- hits.add(hit);
- raw_energy += hit.getEnergy();
- /*
- Subdetector d = hit.getSubdetector();
- int detector_index = detectors.indexOf(d);
- if(detector_index <0)
- {
- detectors.add(d);
- detector_index = detectors.indexOf(d);
- subdetector_raw_energies.add(0.);
- subdetector_corrected_energies.add(0.);
- }
- double sampling_fraction = d.getSamplingFraction();
- */
- int detector_index = 0;
- double sampling_fraction = 1.;
- if(subdetector_raw_energies.size() < 1)
- {
- subdetector_raw_energies.add(0.);
- subdetector_corrected_energies.add(0.);
- }
- double hce = hit.getEnergy()/sampling_fraction;
- corrected_energy += hce;
- hit_energies.add(hce);
- subdetector_raw_energies.set(detector_index,subdetector_raw_energies.get(detector_index) + hit.getEnergy());
- subdetector_corrected_energies.set(detector_index,subdetector_corrected_energies.get(detector_index) + hce);
- needsPropertyCalculation = true;
- }
-
- /**
- * Remove a CalorimeterHit from the cluster
- */
- public void removeHit(CalorimeterHit hit)
- {
- int indx = hits.indexOf(hit);
- hits.remove(hit);
- raw_energy -= hit.getEnergy();
- /*
- Subdetector d = hit.getSubdetector();
- int detector_index = detectors.indexOf(d);
- if(detector_index <0)
- {
- detectors.add(d);
- detector_index = detectors.indexOf(d);
- subdetector_raw_energies.add(0.);
- subdetector_corrected_energies.add(0.);
- }
- double sampling_fraction = d.getSamplingFraction();
- */
- int detector_index = 0;
- double sampling_fraction = 1.;
- double hce = hit.getEnergy()/sampling_fraction;
- corrected_energy -= hce;
- hit_energies.remove(indx);
- subdetector_raw_energies.set(detector_index,subdetector_raw_energies.get(detector_index) - hit.getEnergy());
- subdetector_corrected_energies.set(detector_index,subdetector_corrected_energies.get(detector_index) - hce);
- needsPropertyCalculation = true;
- }
- /**
- * Add a Cluster to the cluster
- */
- public void addCluster(Cluster cluster)
- {
- clusters.add(cluster);
- List<CalorimeterHit> chits = cluster.getCalorimeterHits();
- if(chits.size() > 0)
- {
- for(int i=0;i<chits.size();i++)
- {
- addHit(chits.get(i));
- }
- needsPropertyCalculation = true;
- }
- else
- {
- corrected_energy += cluster.getEnergy();
- double[] sde = cluster.getSubdetectorEnergies();
- for(int i=0;i<sde.length;i++)
- {
- subdetector_corrected_energies.set(i,subdetector_corrected_energies.get(i) + sde[i]);
- }
- }
- }
-
- /**
- * Return the hits comprising the cluster as per interface
- */
- public List<CalorimeterHit> getCalorimeterHits()
- {
- return hits;
- }
-
- /**
- * Return the clusters comprising the cluster as per interface
- */
- public List<Cluster> getClusters()
- {
- return clusters;
- }
-
- /**
- * Return the corrected Energy of the cluster as per interface
- */
- public double getEnergy()
- {
- return corrected_energy;
- }
-
- /**
- * Return the corrected energy contribution of each hit
- * comprising the cluster as per interface
- */
- public double[] getHitContributions()
- {
- double[] earray = new double[hit_energies.size()];
- for(int i=0;i<hit_energies.size();i++)
- {
- earray[i] = hit_energies.get(i);
- }
- return earray;
- }
-
- /**
- * Return the phi direction of the cluster as per interface
- */
- public double getIPhi()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return iphi;
- }
-
- /**
- * Return the theta direction of the cluster as per interface
- */
- public double getITheta()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return itheta;
- }
-
- /**
- * Return the direction error of the cluster as per interface
- */
- public double[] getDirectionError()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return directionError;
- }
-
- /**
- * Return the position of the cluster as per interface
- */
- public double[] getPosition()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return position;
- }
-
- /**
- * Return the position error of the cluster as per interface
- */
- public double[] getPositionError()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return positionError;
- }
-
- /**
- * Return the shape parameters of the cluster as per interface
- */
- public double[] getShape()
- {
- if(needsPropertyCalculation)
- {
- calculateProperties();
- }
- return shapeParameters;
- }
-
- /**
- * Return the subdetector energy contributions
- * of the cluster as per interface
- */
- public double[] getSubdetectorEnergies()
- {
- double[] sdearray = new double[subdetector_corrected_energies.size()];
- for(int i=0;i<subdetector_corrected_energies.size();i++)
- {
- sdearray[i] = subdetector_corrected_energies.get(i);
- }
- return sdearray;
- }
-
- /**
- * Return a bit mask defining the type of cluster as per interface
- */
- public int getType()
- {
- return 0;
- }
-
- /**
- * Return the sum of the raw energies from the hits in the cluster
- */
- public double getRawEnergy()
- {
- return raw_energy;
- }
-
- /**
- * Allow for recalculation of cluster properties, or
- * prevent recalculation after adding to cluster
- */
- public void setNeedsPropertyCalculation(boolean tf)
- {
- needsPropertyCalculation = tf;
- }
-
- /**
- * Override the default ClusterPropertyCalculator
- */
- public void setPropertyCalculator(ClusterPropertyCalculator cpc)
- {
- cluster_property_calculator = cpc;
- }
-
- /**
- * Calculate the cluster properties from the hits
- */
- public void calculateProperties()
- {
- cluster_property_calculator.calculateProperties(hits);
- position = cluster_property_calculator.getPosition();
- positionError = cluster_property_calculator.getPositionError();
- iphi = cluster_property_calculator.getIPhi();
- itheta = cluster_property_calculator.getITheta();
- directionError = cluster_property_calculator.getDirectionError();
- shapeParameters = cluster_property_calculator.getShapeParameters();
- needsPropertyCalculation = false;
- }
-
+ protected List<CalorimeterHit> hits = new ArrayList<CalorimeterHit>();
+ protected List<Cluster> clusters = new ArrayList<Cluster>();
+ protected List<Subdetector> detectors = new ArrayList<Subdetector>();
+ public double raw_energy;
+ public double corrected_energy;
+ protected List<Double> subdetector_raw_energies = new ArrayList<Double>();
+ protected List<Double> subdetector_corrected_energies = new ArrayList<Double>();
+ protected List<Double> hit_energies = new ArrayList<Double>();
+ protected ClusterPropertyCalculator cluster_property_calculator = new DefaultClusterPropertyCalculator();
+ protected boolean needsPropertyCalculation = true;
+ protected double[] position;
+ protected double[] positionError;
+ protected double iphi;
+ protected double itheta;
+ protected double[] directionError;
+ protected double[] shapeParameters;
+
+ /**
+ * Add a CalorimeterHit to the cluster
+ */
+ public void addHit(CalorimeterHit hit)
+ {
+ hits.add(hit);
+ raw_energy += hit.getEnergy();
+ /*
+ Subdetector d = hit.getSubdetector();
+ int detector_index = detectors.indexOf(d);
+ if(detector_index <0)
+ {
+ detectors.add(d);
+ detector_index = detectors.indexOf(d);
+ subdetector_raw_energies.add(0.);
+ subdetector_corrected_energies.add(0.);
+ }
+ double sampling_fraction = d.getSamplingFraction();
+ */
+ int detector_index = 0;
+ double sampling_fraction = 1.;
+ if(subdetector_raw_energies.size() < 1)
+ {
+ subdetector_raw_energies.add(0.);
+ subdetector_corrected_energies.add(0.);
+ }
+ double hce = hit.getEnergy()/sampling_fraction;
+ corrected_energy += hce;
+ hit_energies.add(hce);
+ subdetector_raw_energies.set(detector_index,subdetector_raw_energies.get(detector_index) + hit.getEnergy());
+ subdetector_corrected_energies.set(detector_index,subdetector_corrected_energies.get(detector_index) + hce);
+ needsPropertyCalculation = true;
+ }
+
+ /**
+ * Remove a CalorimeterHit from the cluster
+ */
+ public void removeHit(CalorimeterHit hit)
+ {
+ int indx = hits.indexOf(hit);
+ hits.remove(hit);
+ raw_energy -= hit.getEnergy();
+ /*
+ Subdetector d = hit.getSubdetector();
+ int detector_index = detectors.indexOf(d);
+ if(detector_index <0)
+ {
+ detectors.add(d);
+ detector_index = detectors.indexOf(d);
+ subdetector_raw_energies.add(0.);
+ subdetector_corrected_energies.add(0.);
+ }
+ double sampling_fraction = d.getSamplingFraction();
+ */
+ int detector_index = 0;
+ double sampling_fraction = 1.;
+ double hce = hit.getEnergy()/sampling_fraction;
+ corrected_energy -= hce;
+ hit_energies.remove(indx);
+ subdetector_raw_energies.set(detector_index,subdetector_raw_energies.get(detector_index) - hit.getEnergy());
+ subdetector_corrected_energies.set(detector_index,subdetector_corrected_energies.get(detector_index) - hce);
+ needsPropertyCalculation = true;
+ }
+ /**
+ * Add a Cluster to the cluster
+ */
+ public void addCluster(Cluster cluster)
+ {
+ clusters.add(cluster);
+ List<CalorimeterHit> chits = cluster.getCalorimeterHits();
+ if(chits.size() > 0)
+ {
+ for(int i=0;i<chits.size();i++)
+ {
+ addHit(chits.get(i));
+ }
+ needsPropertyCalculation = true;
+ }
+ else
+ {
+ corrected_energy += cluster.getEnergy();
+ double[] sde = cluster.getSubdetectorEnergies();
+ for(int i=0;i<sde.length;i++)
+ {
+ subdetector_corrected_energies.set(i,subdetector_corrected_energies.get(i) + sde[i]);
+ }
+ }
+ }
+
+ /**
+ * Return the hits comprising the cluster as per interface
+ */
+ public List<CalorimeterHit> getCalorimeterHits()
+ {
+ return hits;
+ }
+
+ /**
+ * Return the clusters comprising the cluster as per interface
+ */
+ public List<Cluster> getClusters()
+ {
+ return clusters;
+ }
+
+ /**
+ * Return the corrected Energy of the cluster as per interface
+ */
+ public double getEnergy()
+ {
+ return corrected_energy;
+ }
+
+ /**
+ * Return the corrected energy contribution of each hit
+ * comprising the cluster as per interface
+ */
+ public double[] getHitContributions()
+ {
+ double[] earray = new double[hit_energies.size()];
+ for(int i=0;i<hit_energies.size();i++)
+ {
+ earray[i] = hit_energies.get(i);
+ }
+ return earray;
+ }
+
+ /**
+ * Return the phi direction of the cluster as per interface
+ */
+ public double getIPhi()
+ {
+ if(needsPropertyCalculation)
+ {
+ calculateProperties();
+ }
+ return iphi;
+ }
+
+ /**
+ * Return the theta direction of the cluster as per interface
+ */
+ public double getITheta()
+ {
+ if(needsPropertyCalculation)
+ {
+ calculateProperties();
+ }
+ return itheta;
+ }
+
+ /**
+ * Return the direction error of the cluster as per interface
+ */
+ public double[] getDirectionError()
+ {
+ if(needsPropertyCalculation)
+ {
+ calculateProperties();
+ }
+ return directionError;
+ }
+
+ /**
+ * Return the position of the cluster as per interface
+ */
+ public double[] getPosition()
+ {
+ if(needsPropertyCalculation)
+ {
+ calculateProperties();
+ }
+ return position;
+ }
+
+ /**
+ * Return the position error of the cluster as per interface
+ */
+ public double[] getPositionError()
+ {
+ if(needsPropertyCalculation)
+ {
+ calculateProperties();
+ }
+ return positionError;
+ }
+
+ /**
+ * Return the shape parameters of the cluster as per interface
+ */
+ public double[] getShape()
+ {
+ if(needsPropertyCalculation)
+ {
+ calculateProperties();
+ }
+ return shapeParameters;
+ }
+
+ /**
+ * Return the subdetector energy contributions
+ * of the cluster as per interface
+ */
+ public double[] getSubdetectorEnergies()
+ {
+ double[] sdearray = new double[subdetector_corrected_energies.size()];
+ for(int i=0;i<subdetector_corrected_energies.size();i++)
+ {
+ sdearray[i] = subdetector_corrected_energies.get(i);
+ }
+ return sdearray;
+ }
+
+ /**
+ * Return a bit mask defining the type of cluster as per interface
+ */
+ public int getType()
+ {
+ return 0;
+ }
+
+ /**
+ * Return the sum of the raw energies from the hits in the cluster
+ */
+ public double getRawEnergy()
+ {
+ return raw_energy;
+ }
+
+ /**
+ * Allow for recalculation of cluster properties, or
+ * prevent recalculation after adding to cluster
+ */
+ public void setNeedsPropertyCalculation(boolean tf)
+ {
+ needsPropertyCalculation = tf;
+ }
+
+ /**
+ * Override the default ClusterPropertyCalculator
+ */
+ public void setPropertyCalculator(ClusterPropertyCalculator cpc)
+ {
+ cluster_property_calculator = cpc;
+ }
+
+ /**
+ * Calculate the cluster properties from the hits
+ */
+ public void calculateProperties()
+ {
+ cluster_property_calculator.calculateProperties(hits);
+ position = cluster_property_calculator.getPosition();
+ positionError = cluster_property_calculator.getPositionError();
+ iphi = cluster_property_calculator.getIPhi();
+ itheta = cluster_property_calculator.getITheta();
+ directionError = cluster_property_calculator.getDirectionError();
+ shapeParameters = cluster_property_calculator.getShapeParameters();
+ needsPropertyCalculation = false;
+ }
+
}
lcsim/src/org/lcsim/recon/cluster/util
diff -u -r1.2 -r1.3
--- TensorClusterPropertyCalculator.java 8 Jul 2005 17:15:08 -0000 1.2
+++ TensorClusterPropertyCalculator.java 18 Jul 2005 03:06:19 -0000 1.3
@@ -9,279 +9,280 @@
/**
* Implementation of ClusterPropertyCalculator Interface using
- * inertia Tensor calculation to derive parameters
- *
+ * inertia Tensor calculation to derive parameters
+ *
* @author cassell
*/
public class TensorClusterPropertyCalculator implements ClusterPropertyCalculator
{
- protected double[] position;
- protected double[] positionError;
- protected double iphi;
- protected double itheta;
- protected double[] directionError;
- protected double[] shapeParameters;
- protected double[] mm_NE;
- protected double[] mm_CE;
- protected double[][] mm_PA;
-
- /**
- * Calculate the cluster properties from the hits
- */
- public TensorClusterPropertyCalculator()
- {
- position = new double[3];
- positionError = new double[6];
- iphi = 0.;
- itheta = 0.;
- directionError = new double[6];
- shapeParameters = new double[6];
- }
- public void calculateProperties(List<CalorimeterHit> hits)
- {
- mm_NE = new double[3];
- mm_CE = new double[3];
- mm_PA = new double[3][3];
- for(int i=0;i<3;++i)
- {
- mm_NE[i] = 0.;
- mm_CE[i] = 0.;
- for(int j=0;j<3;++j){mm_PA[i][j]= 0.;}
- }
- double Etot = 0.0;
- double Exx = 0.0;
- double Eyy = 0.0;
- double Ezz = 0.0;
- double Exy = 0.0;
- double Eyz = 0.0;
- double Exz = 0.0;
- double CEx = 0.0;
- double CEy = 0.0;
- double CEz = 0.0;
- double CEr = 0.0;
- double E1 = 0.0;
- double E2 = 0.0;
- double E3 = 0.0;
- double NE1 = 0.0;
- double NE2 = 0.0;
- double NE3 = 0.0;
- double Tr = 0.0;
- double M = 0.0;
- double Det = 0.0;
- int nhits = hits.size();
- for(int i=0;i<hits.size();i++)
- {
- CalorimeterHit hit = hits.get(i);
- // CalorimeterIDDecoder decoder = hit.getDecoder();
- // decoder.setID(hit.getCellID());
- // double[] pos = new double[3];
- // pos[0] = decoder.getX();
- // pos[1] = decoder.getY();
- // pos[2] = decoder.getZ();
- double[] pos = hit.getPosition();
- // Subdetector d = hit.getSubdetector();
- // double sampling_fraction = d.getSamplingFraction();
- double sampling_fraction = 1.;
- double E = hit.getEnergy()/sampling_fraction;
- Etot += E;
- CEx += E*pos[0];
- CEy += E*pos[1];
- CEz += E*pos[2];
- Exx += E*(pos[1]*pos[1] + pos[2]*pos[2]);
- Eyy += E*(pos[0]*pos[0] + pos[2]*pos[2]);
- Ezz += E*(pos[1]*pos[1] + pos[0]*pos[0]);
- Exy -= E*pos[0]*pos[1];
- Eyz -= E*pos[1]*pos[2];
- Exz -= E*pos[0]*pos[2];
- }
- CEx = CEx/Etot;
- CEy = CEy/Etot;
- CEz = CEz/Etot;
- double CErSq = CEx*CEx + CEy*CEy + CEz*CEz;
- CEr = Math.sqrt(CErSq);
- // now go to center of energy coords.
- if (nhits > 3 )
- {
- Exx = Exx - Etot*(CErSq - CEx*CEx);
- Eyy = Eyy - Etot*(CErSq - CEy*CEy);
- Ezz = Ezz - Etot*(CErSq - CEz*CEz);
- Exy = Exy + Etot*CEx*CEy;
- Eyz = Eyz + Etot*CEy*CEz;
- Exz = Exz + Etot*CEz*CEx;
-
- //
- Tr = Exx + Eyy + Ezz;
- double Dxx = Eyy*Ezz - Eyz*Eyz;
- double Dyy = Ezz*Exx - Exz*Exz;
- double Dzz = Exx*Eyy - Exy*Exy;
- M = Dxx + Dyy + Dzz;
- double Dxy = Exy*Ezz - Exz*Eyz;
- double Dxz = Exy*Eyz - Exz*Eyy;
- Det = Exx*Dxx - Exy*Dxy + Exz*Dxz;
- double xt = Tr*Tr - 3*M;
- double sqrtxt = Math.sqrt(xt);
- // eqn to solve for eigenvalues is x**3 - x**2*Tr + x*M - Det = 0
- // crosses y axis at -Det and inflection points are
-
- double mE1 = 0.;
- double mE2 = 0.;
- double mE3 = 0.;
- double a = (3*M - Tr*Tr)/3.;
- double b = (-2*Tr*Tr*Tr + 9*Tr*M -27*Det)/27.;
- double test = b*b/4. + a*a*a/27.;
- if(test >= 0.01)
- {
- System.out.println("AbstractCluster: Only 1 real root!!!");
- System.out.println(" nhits = " + nhits + "\n");
- System.out.println(" a,b,test = " + a + " " + b + " " + test + "\n");
- }
- else
- {
- double temp;
- if(test >= 0.)temp = 1.;
- else temp = Math.sqrt(b*b*27./(-a*a*a*4.));
- if(b > 0.)temp = -temp;
- double phi = Math.acos(temp);
- double temp1 = 2.*Math.sqrt(-a/3.);
- mE1 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3.);
- mE2 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3. + 2.*Math.PI/3.);
- mE3 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3. + 4.*Math.PI/3.);
- }
- if(mE1 < mE2)
- {
- if(mE1 < mE3)
- {
- E1 = mE1;
- if(mE2 < mE3)
- {
- E2 = mE2;
- E3 = mE3;
- }
- else
- {
- E2 = mE3;
- E3 = mE2;
- }
- }
- else
- {
- E1 = mE3;
- E2 = mE1;
- E3 = mE2;
- }
- }
- else
- {
- if(mE2 < mE3)
- {
- E1 = mE2;
- if(mE1 < mE3)
- {
- E2 = mE1;
- E3 = mE3;
- }
- else
- {
- E2 = mE3;
- E3 = mE1;
- }
- }
- else
- {
- E1 = mE3;
- E2 = mE2;
- E3 = mE1;
- }
- }
-
- NE1 = E1/Etot;
- NE2 = E2/Etot;
- NE3 = E3/Etot;
- double[] EV = new double[3];
- EV[0] = E1;
- EV[1] = E2;
- EV[2] = E3;
- // Now calculate principal axes
- for( int i = 0 ; i < 3 ; i++ )
- {
- double[] C = new double[3];
- C[0] = 1.0;
- C[2] = (Exy*Exy + (Eyy - EV[i])*(EV[i] - Exx))/
- ((Eyy - EV[i])*Exz - Eyz*Exy);
- C[1] = (EV[i] - Exx - Exz*C[2])/Exy;
- double norm = Math.sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
- mm_PA[i][0] = C[0]/norm;
- mm_PA[i][1] = C[1]/norm;
- mm_PA[i][2] = C[2]/norm;
- }
- }
- mm_NE[0] = NE1;
- mm_NE[1] = NE2;
- mm_NE[2] = NE3;
- mm_CE[0] = CEx;
- mm_CE[1] = CEy;
- mm_CE[2] = CEz;
- for(int i=0;i<6;i++)
- {
- positionError[i] = 0.;
- directionError[i] = 0.;
- shapeParameters[i] = 0.;
- }
- for(int i=0;i<3;i++)
- {
- position[i] = mm_CE[i];
- shapeParameters[i] = mm_NE[i];
- }
- if(nhits > 3)
- {
- double dr = Math.sqrt( (position[0]+mm_PA[0][0])*(position[0]+mm_PA[0][0]) +
- (position[1]+mm_PA[0][1])*(position[1]+mm_PA[0][1]) +
- (position[2]+mm_PA[0][2])*(position[2]+mm_PA[0][2]) ) -
- Math.sqrt( (position[0])*(position[0]) +
- (position[1])*(position[1]) +
- (position[2])*(position[2]) ) ;
- double sign = 1.;
- if(dr < 0.)sign = -1.;
- itheta = Math.acos(sign*mm_PA[0][2]);
- iphi = Math.atan2(sign*mm_PA[0][1],sign*mm_PA[0][0]);
- }
- else
- {
- itheta = 999.;
- iphi = 999.;
- }
- }
- public double[] getPosition()
- {
- return position;
- }
- public double[] getPositionError()
- {
- return positionError;
- }
- public double getIPhi()
- {
- return iphi;
- }
- public double getITheta()
- {
- return itheta;
- }
- public double[] getDirectionError()
- {
- return directionError;
- }
- public double[] getShapeParameters()
- {
- return shapeParameters;
- }
- public double[] getNormalizedEigenvalues()
- {
- return mm_NE;
- }
- public double[][] getPrincipleAxis()
- {
- return mm_PA;
- }
-
+ protected double[] position;
+ protected double[] positionError;
+ protected double iphi;
+ protected double itheta;
+ protected double[] directionError;
+ protected double[] shapeParameters;
+ protected double[] mm_NE;
+ protected double[] mm_CE;
+ protected double[][] mm_PA;
+
+ /**
+ * Calculate the cluster properties from the hits
+ */
+ public TensorClusterPropertyCalculator()
+ {
+ position = new double[3];
+ positionError = new double[6];
+ iphi = 0.;
+ itheta = 0.;
+ directionError = new double[6];
+ shapeParameters = new double[6];
+ }
+ public void calculateProperties(List<CalorimeterHit> hits)
+ {
+ mm_NE = new double[3];
+ mm_CE = new double[3];
+ mm_PA = new double[3][3];
+ for(int i=0;i<3;++i)
+ {
+ mm_NE[i] = 0.;
+ mm_CE[i] = 0.;
+ for(int j=0;j<3;++j)
+ {mm_PA[i][j]= 0.;}
+ }
+ double Etot = 0.0;
+ double Exx = 0.0;
+ double Eyy = 0.0;
+ double Ezz = 0.0;
+ double Exy = 0.0;
+ double Eyz = 0.0;
+ double Exz = 0.0;
+ double CEx = 0.0;
+ double CEy = 0.0;
+ double CEz = 0.0;
+ double CEr = 0.0;
+ double E1 = 0.0;
+ double E2 = 0.0;
+ double E3 = 0.0;
+ double NE1 = 0.0;
+ double NE2 = 0.0;
+ double NE3 = 0.0;
+ double Tr = 0.0;
+ double M = 0.0;
+ double Det = 0.0;
+ int nhits = hits.size();
+ for(int i=0;i<hits.size();i++)
+ {
+ CalorimeterHit hit = hits.get(i);
+ // CalorimeterIDDecoder decoder = hit.getDecoder();
+ // decoder.setID(hit.getCellID());
+ // double[] pos = new double[3];
+ // pos[0] = decoder.getX();
+ // pos[1] = decoder.getY();
+ // pos[2] = decoder.getZ();
+ double[] pos = hit.getPosition();
+ // Subdetector d = hit.getSubdetector();
+ // double sampling_fraction = d.getSamplingFraction();
+ double sampling_fraction = 1.;
+ double E = hit.getEnergy()/sampling_fraction;
+ Etot += E;
+ CEx += E*pos[0];
+ CEy += E*pos[1];
+ CEz += E*pos[2];
+ Exx += E*(pos[1]*pos[1] + pos[2]*pos[2]);
+ Eyy += E*(pos[0]*pos[0] + pos[2]*pos[2]);
+ Ezz += E*(pos[1]*pos[1] + pos[0]*pos[0]);
+ Exy -= E*pos[0]*pos[1];
+ Eyz -= E*pos[1]*pos[2];
+ Exz -= E*pos[0]*pos[2];
+ }
+ CEx = CEx/Etot;
+ CEy = CEy/Etot;
+ CEz = CEz/Etot;
+ double CErSq = CEx*CEx + CEy*CEy + CEz*CEz;
+ CEr = Math.sqrt(CErSq);
+ // now go to center of energy coords.
+ if (nhits > 3 )
+ {
+ Exx = Exx - Etot*(CErSq - CEx*CEx);
+ Eyy = Eyy - Etot*(CErSq - CEy*CEy);
+ Ezz = Ezz - Etot*(CErSq - CEz*CEz);
+ Exy = Exy + Etot*CEx*CEy;
+ Eyz = Eyz + Etot*CEy*CEz;
+ Exz = Exz + Etot*CEz*CEx;
+
+ //
+ Tr = Exx + Eyy + Ezz;
+ double Dxx = Eyy*Ezz - Eyz*Eyz;
+ double Dyy = Ezz*Exx - Exz*Exz;
+ double Dzz = Exx*Eyy - Exy*Exy;
+ M = Dxx + Dyy + Dzz;
+ double Dxy = Exy*Ezz - Exz*Eyz;
+ double Dxz = Exy*Eyz - Exz*Eyy;
+ Det = Exx*Dxx - Exy*Dxy + Exz*Dxz;
+ double xt = Tr*Tr - 3*M;
+ double sqrtxt = Math.sqrt(xt);
+ // eqn to solve for eigenvalues is x**3 - x**2*Tr + x*M - Det = 0
+ // crosses y axis at -Det and inflection points are
+
+ double mE1 = 0.;
+ double mE2 = 0.;
+ double mE3 = 0.;
+ double a = (3*M - Tr*Tr)/3.;
+ double b = (-2*Tr*Tr*Tr + 9*Tr*M -27*Det)/27.;
+ double test = b*b/4. + a*a*a/27.;
+ if(test >= 0.01)
+ {
+ System.out.println("AbstractCluster: Only 1 real root!!!");
+ System.out.println(" nhits = " + nhits + "\n");
+ System.out.println(" a,b,test = " + a + " " + b + " " + test + "\n");
+ }
+ else
+ {
+ double temp;
+ if(test >= 0.)temp = 1.;
+ else temp = Math.sqrt(b*b*27./(-a*a*a*4.));
+ if(b > 0.)temp = -temp;
+ double phi = Math.acos(temp);
+ double temp1 = 2.*Math.sqrt(-a/3.);
+ mE1 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3.);
+ mE2 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3. + 2.*Math.PI/3.);
+ mE3 = Tr/3. + 2.*Math.sqrt(-a/3.)*Math.cos(phi/3. + 4.*Math.PI/3.);
+ }
+ if(mE1 < mE2)
+ {
+ if(mE1 < mE3)
+ {
+ E1 = mE1;
+ if(mE2 < mE3)
+ {
+ E2 = mE2;
+ E3 = mE3;
+ }
+ else
+ {
+ E2 = mE3;
+ E3 = mE2;
+ }
+ }
+ else
+ {
+ E1 = mE3;
+ E2 = mE1;
+ E3 = mE2;
+ }
+ }
+ else
+ {
+ if(mE2 < mE3)
+ {
+ E1 = mE2;
+ if(mE1 < mE3)
+ {
+ E2 = mE1;
+ E3 = mE3;
+ }
+ else
+ {
+ E2 = mE3;
+ E3 = mE1;
+ }
+ }
+ else
+ {
+ E1 = mE3;
+ E2 = mE2;
+ E3 = mE1;
+ }
+ }
+
+ NE1 = E1/Etot;
+ NE2 = E2/Etot;
+ NE3 = E3/Etot;
+ double[] EV = new double[3];
+ EV[0] = E1;
+ EV[1] = E2;
+ EV[2] = E3;
+ // Now calculate principal axes
+ for( int i = 0 ; i < 3 ; i++ )
+ {
+ double[] C = new double[3];
+ C[0] = 1.0;
+ C[2] = (Exy*Exy + (Eyy - EV[i])*(EV[i] - Exx))/
+ ((Eyy - EV[i])*Exz - Eyz*Exy);
+ C[1] = (EV[i] - Exx - Exz*C[2])/Exy;
+ double norm = Math.sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
+ mm_PA[i][0] = C[0]/norm;
+ mm_PA[i][1] = C[1]/norm;
+ mm_PA[i][2] = C[2]/norm;
+ }
+ }
+ mm_NE[0] = NE1;
+ mm_NE[1] = NE2;
+ mm_NE[2] = NE3;
+ mm_CE[0] = CEx;
+ mm_CE[1] = CEy;
+ mm_CE[2] = CEz;
+ for(int i=0;i<6;i++)
+ {
+ positionError[i] = 0.;
+ directionError[i] = 0.;
+ shapeParameters[i] = 0.;
+ }
+ for(int i=0;i<3;i++)
+ {
+ position[i] = mm_CE[i];
+ shapeParameters[i] = mm_NE[i];
+ }
+ if(nhits > 3)
+ {
+ double dr = Math.sqrt( (position[0]+mm_PA[0][0])*(position[0]+mm_PA[0][0]) +
+ (position[1]+mm_PA[0][1])*(position[1]+mm_PA[0][1]) +
+ (position[2]+mm_PA[0][2])*(position[2]+mm_PA[0][2]) ) -
+ Math.sqrt( (position[0])*(position[0]) +
+ (position[1])*(position[1]) +
+ (position[2])*(position[2]) ) ;
+ double sign = 1.;
+ if(dr < 0.)sign = -1.;
+ itheta = Math.acos(sign*mm_PA[0][2]);
+ iphi = Math.atan2(sign*mm_PA[0][1],sign*mm_PA[0][0]);
+ }
+ else
+ {
+ itheta = 999.;
+ iphi = 999.;
+ }
+ }
+ public double[] getPosition()
+ {
+ return position;
+ }
+ public double[] getPositionError()
+ {
+ return positionError;
+ }
+ public double getIPhi()
+ {
+ return iphi;
+ }
+ public double getITheta()
+ {
+ return itheta;
+ }
+ public double[] getDirectionError()
+ {
+ return directionError;
+ }
+ public double[] getShapeParameters()
+ {
+ return shapeParameters;
+ }
+ public double[] getNormalizedEigenvalues()
+ {
+ return mm_NE;
+ }
+ public double[][] getPrincipleAxis()
+ {
+ return mm_PA;
+ }
+
}
lcsim/src/org/lcsim/recon/cluster/util
diff -N FCClusterPropertyCalculator.java
--- FCClusterPropertyCalculator.java 8 Jul 2005 17:15:08 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,330 +0,0 @@
-package org.lcsim.recon.cluster.util;
-
-import java.util.List;
-import org.lcsim.spacegeom.PrincipalAxesLineFitter;
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.geometry.CalorimeterIDDecoder;
-import org.lcsim.util.fourvec.Lorentz4Vector;
-import org.lcsim.util.fourvec.Momentum4Vector;
-
-/**
- * A class encapsulating the behavior of a calorimeter cluster.
- *
- * @author Norman A. Graf
- * @version 1.0
- */
-public class FCClusterPropertyCalculator implements ClusterPropertyCalculator
-{
- private CalorimeterIDDecoder _decoder;
- private Lorentz4Vector _vec;
- // second moment of the cluster
- private double _width;
- private double[] _layerEnergy;
- private double _clusterEnergy;
- private double[] _layerWidth;
- private double[] _centroid;
- private double[] _directionCosines;
- private double _samplingFraction;
- private CalorimeterHit _hottestCell;
- private double _highestCellEnergy;
- private boolean _isEndCap;
- private boolean _isNorth;
- private double _chisq = 99999.;
- private int layers;
- private double _itheta;
- private double _iphi;
-
- /**
- * Fully qualified constructor
- *
- * @param decoder The CalorimeterIDDecoder used for indexing
- * @param vec The Lorentz four-vector
- * @param members A Vector of CalorimeterHits
- * @param layers The number of layers in the Calorimeter
- * @param samplingFraction The conversion from measured to deposited energy
- */
- public FCClusterPropertyCalculator(CalorimeterIDDecoder decoder)
- {
- _decoder = decoder;
- layers = 64;
- }
- public void calculateProperties(List<CalorimeterHit> hits)
- {
- _vec = calculateVec(hits);
- _layerEnergy = new double[layers];
- _layerWidth = new double[layers];
- // the array of cell (x,y,z) coordinates
- double[][] points = new double[3][hits.size()];
- int npoints=0;
- _highestCellEnergy = 0.;
-
- for(CalorimeterHit h : hits )
- {
- // Subdetector d = hit.getSubdetector();
- // double _sampling_fraction = d.getSamplingFraction();
- double _sampling_fraction = 1.;
- _decoder.setID(h.getCellID());
- // this is a change...
- // incorporate sampling fractions to get correct cluster energy
- double e = h.getEnergy()/_samplingFraction;
- if (e>_highestCellEnergy)
- {
- _highestCellEnergy = e;
- _hottestCell = h;
- // for now let highest energy cells determine location
-// _isEndCap = cell.isEndcap();
-// _isNorth = cell.isNorth();
- }
- // calculate the energy-weighted cluster width (second moment)
- double dtheta=_vec.theta() - _decoder.getTheta();
- double dphi=_vec.phi() - _decoder.getPhi();
- // phi-wrap at 0:2pi?
- if (dphi>Math.PI)
- {
- dphi-=2.*Math.PI;
- }
- if (dphi<-Math.PI)
- {
- dphi+=2.*Math.PI;
- }
- double dRw=(dtheta*dtheta + dphi*dphi)*e;
- _width+=dRw;
- // increment the energy deposited in this layer
- _layerEnergy[_decoder.getLayer()]+=e;
- _clusterEnergy+=e;
- // increment the width (second moment of energy) in this layer
- _layerWidth[_decoder.getLayer()]+=dRw;
- //store the hit (x,y,z) coordinates
- points[0][npoints] = _decoder.getX();
- points[1][npoints] = _decoder.getY();
- points[2][npoints] = _decoder.getZ();
- npoints++;
- }
- // normalize the second moments
- _width /= _clusterEnergy;
- for(int i=0; i<layers; ++i)
- {
- _layerWidth[i]/=_clusterEnergy;
- }
-
- // fit a straight line through the cells and store the results
- PrincipalAxesLineFitter lf = new PrincipalAxesLineFitter();
- lf.fit(points);
- _centroid = lf.centroid();
- _directionCosines = lf.dircos();
- double dr = Math.sqrt( (_centroid[0]+_directionCosines[0])*(_centroid[0]+_directionCosines[0]) +
- (_centroid[1]+_directionCosines[1])*(_centroid[1]+_directionCosines[1]) +
- (_centroid[2]+_directionCosines[2])*(_centroid[2]+_directionCosines[2]) ) -
- Math.sqrt( (_centroid[0])*(_centroid[0]) +
- (_centroid[1])*(_centroid[1]) +
- (_centroid[2])*(_centroid[2]) ) ;
- double sign = 1.;
- if(dr < 0.)sign = -1.;
- _itheta = Math.acos(_directionCosines[2]);
- _iphi = Math.atan2(_directionCosines[1],_directionCosines[0]);
-
- // finish up the cluster (base class method)
-// calculateDerivedQuantities();
- }
-
- /**
- * Calculate the cluster four-momentum.
- * The Lorentz four-vector is derived from the cluster cells.
- *
- */
- public Lorentz4Vector calculateVec(List<CalorimeterHit> hits)
- {
- Lorentz4Vector sum = new Momentum4Vector();
- double[] sums = {0.,0.,0.};
- double wtsum = 0.;
- for(int i=0;i<hits.size();i++)
- {
- CalorimeterHit hit = hits.get(i);
- // Subdetector d = hit.getSubdetector();
- // double sampling_fraction = d.getSamplingFraction();
- double sampling_fraction = 1.;
- double[] pos = new double[3];
- _decoder.setID(hit.getCellID());
- pos[0] = _decoder.getX();
- pos[1] = _decoder.getY();
- pos[2] = _decoder.getZ();
- double wt = hit.getEnergy()/sampling_fraction;
- wtsum += wt;
- for(int j=0;j<3;j++)
- {
- sums[j] += wt*pos[j];
- }
- }
- sum.plusEquals(new Momentum4Vector(sums[0], sums[1], sums[2], wtsum));
- return sum;
- }
-
- /**
- * The cluster width (energy second moment).
- *
- * @return The cluster width
- */
- public double width()
- {
- return _width;
- }
-
- /**
- * The cluster four-momentum
- *
- * @return The Lorentz four-vector
- */
- public Lorentz4Vector vector()
- {
- return _vec;
- }
-
- /**
- * The constituent cells
- *
- * @return Vector of the CalorimeterHits constituting the cluster.
- */
- /**
- * The cluster energy deposited in a specific layer
- *
- * @return The cluster energy in layer <b>layer</b>
- */
- public double layerEnergy(int layer)
- {
- return _layerEnergy[layer];
- }
-
- /**
- * The cluster layer energies
- *
- * @return The array of cluster energies in each layer.
- */
- public double[] layerEnergies()
- {
- return _layerEnergy;
- }
-
-
- /**
- * The cluster energy corrected for sampling fractions
- *
- * @return The cluster energy
- */
- public double clusterEnergy()
- {
- return _clusterEnergy;
- }
-
- /**
- * The energy of the highest energy cell in this cluster
- *
- * @return The energy of the highest energy cell in this cluster corrected by the sampling fraction.
- */
- public double highestCellEnergy()
- {
- return _highestCellEnergy;
- }
-
- /**
- * The CalorimeterHit in this cluster with the highest energy
- *
- * @return The CalorimeterHit in this cluster with the highest energy
- */
- public CalorimeterHit hottestCell()
- {
- return _hottestCell;
- }
-
- /**
- * The cluster width (energy second moment) deposited in a specific layer
- *
- * @return The cluster width in layer <b>layer</b>
- */
- public double layerWidth(int layer)
- {
- return _layerWidth[layer];
- }
-
- /**
- * The unweighted spatial centroid (x,y,z) of the cluster line fit
- *
- * @return The unweighted spatial centroid (x,y,z) of the cluster
- */
- public double[] centroid()
- {
- return _centroid;
- }
-
- /**
- * The direction cosines of the cluster line fit
- *
- * @return The direction cosines of the cluster
- */
- public double[] directionCosines()
- {
- return _directionCosines;
- }
-
-
- /**
- * Returns topological position of cluster.
- *
- * @return true if in EndCap
- */
- public boolean isEndCap()
- {
- return _isEndCap;
- }
-
-
- /**
- * Returns topological position of cluster.
- *
- * @return true if in "North" EndCap
- */
- public boolean isNorth()
- {
- return _isNorth;
- }
-
- public void setChisq(double chisq)
- {
- _chisq = chisq;
- }
-
- public double chisq()
- {
- return _chisq;
- }
- public double[] getPosition()
- {
- return _centroid;
- }
- public double[] getPositionError()
- {
- double[] positionError = {0.,0.,0.,0.,0.,0.};
- return positionError;
- }
- public double getIPhi()
- {
- return _iphi;
- }
- public double getITheta()
- {
- return _itheta;
- }
- public double[] getDirectionError()
- {
- double[] directionError = {0.,0.,0.,0.,0.,0.};
- return directionError;
- }
- public double[] getShapeParameters()
- {
- double[] shapeParameters = {0.,0.,0.,0.,0.,0.};
- shapeParameters[0] = _width;
- return shapeParameters;
- }
-
-
-}
-
CVSspam 0.2.8
