hps-java/src/main/java/org/lcsim/hps/analysis/ecal
diff -N HPSEcalPlotsDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HPSEcalPlotsDriver.java 16 May 2011 19:14:29 -0000 1.1
@@ -0,0 +1,518 @@
+package org.lcsim.contrib.jeremym.hps;
+
+import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
+import hep.aida.IHistogram1D;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import java.util.Map.Entry;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.units.SystemOfUnits;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * Diagnostic plots for HPS ECal.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @version $Id: HPSEcalPlotsDriver.java,v 1.1 2011/05/16 19:14:29 jeremy Exp $
+ */
+// TODO ZY plot can use h2d with small bins in Y to avoid weird binning effect
+// TODO ZY of cluster positions (algo from Tim to find pos)
+// TODO ZY of cluster seed
+
+public class HPSEcalPlotsDriver extends Driver
+{
+ String ecalCollectionName = "EcalHits";
+ String ecalClusterCollectionName = "EcalClusters";
+
+ AIDA aida = AIDA.defaultInstance();
+
+ // MCParticle plots.
+ ICloud1D primaryEPlot;
+ ICloud1D fsPlot;
+
+ // CalHit plots.
+ IHistogram1D hitEPlot;
+ ICloud1D ecalEPlot;
+ ICloud1D eventEPlot;
+ ICloud2D hitXZPlot;
+ ICloud2D hitYZPlot;
+ ICloud1D hitUnder100MeVPlot;
+ IHistogram1D hitOver100MeVPlot;
+ ICloud1D maxHitEPlot;
+ IHistogram1D maxTimePlot;
+ ICloud1D timePlot;
+ ICloud1D hitCountPlot;
+ ICloud1D idCountPlot;
+
+ // Cluster plots.
+ IHistogram1D nclusPlot;
+ ICloud1D clusEPlot;
+ ICloud1D clusTotEPlot;
+ ICloud1D leadClusEPlot;
+ ICloud1D leadClus2EPlot;
+ ICloud1D clusResTop3Plot;
+ IHistogram1D clusHitPlot;
+ ICloud1D clusSeedEPlot;
+ ICloud1D clusSeedDistPlot;
+ ICloud2D leadClusAndPrimaryPlot;
+
+ //ICloud1D hit50MeVPlot;
+ //ICloud1D hit30MeVPlot;
+
+ //ICloud1D residualAllPlot;
+ //ICloud1D residualPlot;
+ //ICloud1D residualMaxHitPlot;
+ //ICloud1D residualTop2Plot;
+ //ICloud1D residualTop3Plot;
+
+ class MCParticleEComparator implements Comparator<MCParticle>
+ {
+ public int compare(MCParticle p1, MCParticle p2)
+ {
+ double e1 = p1.getEnergy();
+ double e2 = p2.getEnergy();
+ if (e1 < e2)
+ {
+ return -1;
+ }
+ else if (e1 == e2)
+ {
+ return 0;
+ }
+ else
+ {
+ return 1;
+ }
+ }
+ }
+
+ class ClusterEComparator implements Comparator<Cluster>
+ {
+ public int compare(Cluster o1, Cluster o2)
+ {
+ if (o1.getEnergy() < o2.getEnergy())
+ {
+ return -1;
+ }
+ else if (o1.getEnergy() > o2.getEnergy())
+ {
+ return 1;
+ }
+ else
+ {
+ return 0;
+ }
+ }
+ }
+
+ public void startOfData()
+ {
+ timePlot = aida.cloud1D(ecalCollectionName + " : Hit Time");
+ timePlot.annotation().addItem("yAxisScale", "log");
+ timePlot.annotation().addItem("xAxisLabel", "Time [ns]");
+
+ maxTimePlot = aida.histogram1D(ecalCollectionName + " : Max Time", 200, 0, 1000);
+ maxTimePlot.annotation().addItem("yAxisScale", "log");
+ maxTimePlot.annotation().addItem("xAxisLabel", "Time [ns]");
+
+ maxHitEPlot = aida.cloud1D(ecalCollectionName + " : Max Hit E in Event");
+ maxHitEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ hitEPlot = aida.histogram1D(ecalCollectionName + " : Hit Energy", 200, 0, 3500);
+ hitEPlot.annotation().addItem("yAxisScale", "log");
+ hitEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ hitCountPlot = aida.cloud1D(ecalCollectionName + " : Hit Count");
+ hitCountPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ idCountPlot = aida.cloud1D(ecalCollectionName + " : Uniq Hit IDs");
+ idCountPlot.annotation().addItem("xAxisLabel", "Number of Unique IDs in Event");
+
+ ecalEPlot = aida.cloud1D(ecalCollectionName + " : Total E in Event");
+ ecalEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ primaryEPlot = aida.cloud1D("MCParticle: Highest Primary E in Event");
+ primaryEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ eventEPlot = aida.cloud1D("MCParticle: Total Gen FS Electron E in Event");
+ eventEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ //residualPlot = aida.cloud1D(ecalCollectionName + " : Residual with Highest E Particle [GeV]");
+ //residualPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ //residualAllPlot = aida.cloud1D(ecalCollectionName + " : Resdual with All Final State Electrons [GeV]");
+ //residualAllPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ //residualMaxHitPlot = aida.cloud1D(ecalCollectionName + " : Residual of Max Hit and Highest E Electron");
+ //residualMaxHitPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ hitYZPlot = aida.cloud2D(ecalCollectionName + " : Y vs Z");
+ hitYZPlot.annotation().addItem("xAxisLabel", "Y [mm]");
+ hitYZPlot.annotation().addItem("yAxisLabel", "Z [mm]");
+
+ hitXZPlot = aida.cloud2D(ecalCollectionName + " : X vs Z");
+ hitXZPlot.annotation().addItem("xAxisLabel", "X [mm]");
+ hitXZPlot.annotation().addItem("yAxisLabel", "Z [mm]");
+
+ //residualTop2Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 2 FS Particle E");
+ //residualTop2Plot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ //residualTop3Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 3 FS Particle E");
+ //residualTop3Plot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ fsPlot = aida.cloud1D("MCParticle: Number of Final State Particles");
+ fsPlot.annotation().addItem("xAxisLabel", "Number of FS Particles");
+
+ clusEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster E");
+ clusEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ nclusPlot = aida.histogram1D(ecalClusterCollectionName + " : Number of Clusters", 20, 0, 20);
+ nclusPlot.annotation().addItem("xAxisLabel", "Number of Clusters");
+
+ //hit50MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 50 MeV");
+ //hit50MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ //hit30MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 30 MeV");
+ //hit30MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ hitUnder100MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 100 MeV");
+ hitUnder100MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ hitOver100MeVPlot = aida.histogram1D(ecalCollectionName + " : Hits with E >= 100 MeV", 20, 0, 20);
+ hitUnder100MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ leadClusEPlot = aida.cloud1D(ecalClusterCollectionName + " : Leading Cluster E");
+ leadClusEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ leadClus2EPlot = aida.cloud1D(ecalClusterCollectionName + " : Second Leading Cluster E");
+ leadClus2EPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusTotEPlot = aida.cloud1D(ecalClusterCollectionName + " : Total Clus E in Event");
+ clusTotEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusResTop3Plot = aida.cloud1D(ecalClusterCollectionName + " : Total Clus E Residual with Top 3 Particles");
+ clusResTop3Plot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusHitPlot = aida.histogram1D(ecalClusterCollectionName + " : Number of Clusters per Hit", 5, 1, 6);
+ clusHitPlot.annotation().addItem("xAxisLabel", "Number of Clusters");
+
+ clusSeedDistPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Hit Distance from Beam Axis");
+ clusSeedDistPlot.annotation().addItem("xAxisLabel", "Number of Clusters");
+ clusSeedDistPlot.annotation().addItem("yAxisScale", "log");
+
+ clusSeedEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Hit Distance from Beam Axis");
+ clusSeedEPlot.annotation().addItem("xAxisLabel", "Distance [mm]");
+
+ leadClusAndPrimaryPlot = aida.cloud2D(ecalClusterCollectionName + " : Lead Cluster E vs Highest Primary Particle E");
+ }
+
+ public void process(EventHeader event)
+ {
+ // sum and max vars
+ double esum = 0;
+ double emax = 0;
+ double tmax = 0;
+
+ // Loop over hits from ECal collection.
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ int nhits = hits.size();
+
+ // MCParticles
+ List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
+
+ // Final State particles.
+ List<MCParticle> fsParticles = getFinalStateParticles(mcparticles);
+
+ //System.out.println("fsParticles="+fsParticles.size());
+ fsPlot.fill(fsParticles.size());
+
+ // Sort MCParticles on energy.
+ Collections.sort(fsParticles, new MCParticleEComparator());
+
+ // Energy of top two FS particles.
+ double e2 = fsParticles.get(0).getEnergy() + fsParticles.get(1).getEnergy();
+
+ // Energy of top three FS particles.
+ double e3 = e2 + fsParticles.get(2).getEnergy();
+
+ // Check unique IDs.
+ Set<Long> ids = new HashSet<Long>();
+
+ // n hits
+ hitCountPlot.fill(nhits);
+
+ //int naccept = 0;
+
+ //int nhits50MeV = 0;
+ //int nhits30MeV = 0;
+ int nhits100MeV = 0;
+ int nhitsOver100MeV = 0;
+
+ // Loop over ECal hits.
+ for (CalorimeterHit hit : hits)
+ {
+ // get raw E
+ double eraw = hit.getRawEnergy();
+
+ if (eraw >= 0.1)
+ {
+ nhitsOver100MeV++;
+ }
+
+ if (eraw < 0.1)
+ {
+ nhits100MeV++;
+ }
+
+ //if (eraw < 0.05)
+ //{
+ // nhits50MeV++;
+ //}
+
+ //if (eraw < 0.03)
+ //{
+ // nhits30MeV++;
+ //}
+
+ // time
+ timePlot.fill(hit.getTime());
+
+ // YZ
+ hitYZPlot.fill(hit.getPosition()[1], hit.getPosition()[2]);
+
+ // XZ
+ hitXZPlot.fill(hit.getPosition()[0], hit.getPosition()[2]);
+
+ // hit E
+ hitEPlot.fill(eraw / SystemOfUnits.MeV);
+
+ // check E max
+ if (eraw > emax)
+ emax = eraw;
+
+ // check T max
+ if (hit.getTime() > tmax)
+ tmax = hit.getTime();
+
+ if (ids.contains(hit.getCellID()))
+ throw new RuntimeException("Duplicate cell ID: " + hit.getCellID());
+
+ ids.add(hit.getCellID());
+
+ // incr E sum
+ esum += eraw;
+ }
+
+ hitUnder100MeVPlot.fill(nhits100MeV);
+ //hit50MeVPlot.fill(nhits50MeV);
+ //hit30MeVPlot.fill(nhits30MeV);
+ hitOver100MeVPlot.fill(nhitsOver100MeV);
+
+ // total E in Cal
+ ecalEPlot.fill(esum);
+
+ // max hit E in event
+ maxHitEPlot.fill(emax);
+
+ // max hit time
+ maxTimePlot.fill(tmax);
+
+ // number of unique hit ids
+ idCountPlot.fill(ids.size());
+
+ // primary particle with most E
+ MCParticle primary = getPrimary(mcparticles);
+ double primaryE = primary.getEnergy();
+ primaryEPlot.fill(primaryE);
+
+ // residual using highest E primary
+ //double singleResidual = esum - primaryE;
+ //residualPlot.fill(singleResidual);
+
+ // event energy
+ double eventE = getPrimaryE(mcparticles);
+ eventEPlot.fill(eventE);
+
+ // event residual
+ //double residual = esum - eventE;
+ //residualAllPlot.fill(residual);
+
+ // hottest hit residual with highest E particle
+ //residualMaxHitPlot.fill(emax - primaryE);
+
+ // Residual with two
+ //double res2 = esum - e2;
+ //residualTop2Plot.fill(res2);
+
+ //double res3 = esum - e3;
+ //residualTop3Plot.fill(res3);
+
+ List<Cluster> clusters = event.get(Cluster.class, ecalClusterCollectionName);
+ if (clusters == null)
+ throw new RuntimeException("Missing cluster collection!");
+ Collections.sort(clusters, new ClusterEComparator());
+
+ nclusPlot.fill(clusters.size());
+
+ // Leading cluster E.
+ if (clusters.size() > 0)
+ {
+ Cluster leadClus = clusters.get(clusters.size()-1);
+ leadClusEPlot.fill(leadClus.getEnergy());
+
+ leadClusAndPrimaryPlot.fill(primaryE, leadClus.getEnergy());
+ }
+
+ // Second leading cluster E.
+ if (clusters.size() > 1)
+ {
+ leadClus2EPlot.fill(clusters.get(clusters.size()-2).getEnergy());
+ }
+
+ Map<CalorimeterHit,Integer> hitClusMap = new HashMap<CalorimeterHit,Integer>();
+
+ double clusE = 0;
+ Set<CalorimeterHit> naughty = new HashSet<CalorimeterHit>();
+ for (Cluster clus : clusters)
+ {
+ double e = clus.getEnergy();
+ clusEPlot.fill(e);
+ clusE += e;
+ CalorimeterHit seedHit = null;
+ double maxe = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits())
+ {
+ if (hitClusMap.containsKey(hit))
+ {
+ int nshared = hitClusMap.get(hit);
+ ++nshared;
+
+ // debug
+ if (nshared >= 4)
+ naughty.add(hit);
+ //
+
+ hitClusMap.put(hit, nshared);
+ }
+ else
+ {
+ hitClusMap.put(hit, 1);
+ }
+
+ if (hit.getRawEnergy() > maxe)
+ {
+ seedHit = hit;
+ }
+ }
+
+ // Seed distance from X axis.
+ Hep3Vector pos = seedHit.getPositionVec();
+ double y = pos.y();
+ double z = pos.z();
+ clusSeedEPlot.fill(Math.sqrt(y*y+z*z));
+
+ // Seed E.
+ clusSeedEPlot.fill(seedHit.getRawEnergy());
+ }
+
+ // debug of clusters with too many hits
+ for (CalorimeterHit duggan : naughty)
+ {
+ for (Cluster clus : clusters)
+ {
+ if (clus.getCalorimeterHits().contains(duggan))
+ {
+ System.out.println("naughty clus w/ " + clus.getCalorimeterHits().size() + " hits and E = " + clus.getEnergy());
+ for (CalorimeterHit clusHit : clus.getCalorimeterHits())
+ {
+ IDDecoder dec = clusHit.getIDDecoder();
+ dec.setID(clusHit.getCellID());
+ int ix = dec.getValue("ix");
+ int iy = dec.getValue("iy");
+ int side = dec.getValue("side");
+ System.out.println(" x, y, side, E: " + ix + ", " + iy + ", " + side + ", " + clusHit.getRawEnergy());
+ }
+ }
+ }
+ }
+
+ // Total E in all clusters.
+ clusTotEPlot.fill(clusE);
+
+ // Residual of cluster total E and E from top 3 primary particles.
+ clusResTop3Plot.fill(clusE - e3);
+
+ for (Entry<CalorimeterHit, Integer> clusHit : hitClusMap.entrySet())
+ {
+ clusHitPlot.fill(clusHit.getValue());
+ }
+ }
+
+ private double getPrimaryE(List<MCParticle> particles)
+ {
+ double electronE = 0;
+ for (MCParticle particle : particles)
+ {
+ if (particle.getGeneratorStatus() == MCParticle.FINAL_STATE
+ && Math.abs(particle.getPDGID()) == 11)
+ {
+ electronE += particle.getEnergy();
+ }
+ }
+ return electronE;
+ }
+
+ private MCParticle getPrimary(List<MCParticle> particles)
+ {
+ double maxE = 0;
+ MCParticle primary = null;
+ for (MCParticle particle : particles)
+ {
+ if (particle.getGeneratorStatus() == MCParticle.FINAL_STATE
+ && particle.getEnergy() > maxE)
+ {
+ maxE = particle.getEnergy();
+ primary = particle;
+ }
+ }
+ return primary;
+ }
+
+ private static List<MCParticle> getFinalStateParticles(List<MCParticle> mcparticles)
+ {
+ List<MCParticle> fsParticles = new ArrayList<MCParticle>();
+ for (MCParticle mcparticle : mcparticles)
+ {
+ if (mcparticle.getGeneratorStatus() == MCParticle.FINAL_STATE)
+ {
+ fsParticles.add(mcparticle);
+ }
+ }
+ return fsParticles;
+ }
+
+ public void setEcalCollectionName(String ecalCollectionName)
+ {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setEcalClusterCollectionName(String ecalClusterCollectionName)
+ {
+ this.ecalClusterCollectionName = ecalClusterCollectionName;
+ }
+}
hps-java/src/main/java/org/lcsim/hps/recon/ecal
diff -N HPSEcalClusterer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HPSEcalClusterer.java 16 May 2011 19:14:29 -0000 1.1
@@ -0,0 +1,194 @@
+package org.lcsim.contrib.jeremym.hps;
+
+import java.util.ArrayList;
+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.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.geometry.subdetector.HPSEcal;
+import org.lcsim.geometry.subdetector.HPSEcal.XYSide;
+import org.lcsim.geometry.util.IDEncoder;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.util.Driver;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ * Creates clusters from CalorimeterHits in the HPSEcal detector.
+ *
+ * Clustering algorithm is from pages 83 and 84 of the HPS Proposal.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @version $Id: HPSEcalClusterer.java,v 1.1 2011/05/16 19:14:29 jeremy Exp $
+ */
+public class HPSEcalClusterer extends Driver
+{
+ HPSEcal ecal;
+
+ String ecalCollectionName;
+ String ecalName;
+ String clusterCollectionName = "EcalClusters";
+
+ // Minimum E for cluster seed.
+ double seedEMin = .05;
+
+ // Minimum E to add hit to cluster.
+ double addEMin = .03;
+
+ // Odd or even number of crystals in X.
+ boolean oddX;
+
+ public HPSEcalClusterer()
+ {}
+
+ public void setClusterCollectionName(String clusterCollectionName)
+ {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ public void setSeedEMin(double seedEMin)
+ {
+ this.seedEMin = seedEMin;
+ }
+
+ public void setAddEMin(double addEMin)
+ {
+ this.addEMin = addEMin;
+ }
+
+ public void setEcalCollectionName(String ecalCollectionName)
+ {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setEcalName(String ecalName)
+ {
+ this.ecalName = ecalName;
+ }
+
+ public void startOfData()
+ {
+ if (ecalCollectionName == null)
+ throw new RuntimeException("The parameter ecalCollectionName was not set!");
+
+ if (ecalCollectionName == null)
+ throw new RuntimeException("The parameter ecalName was not set!");
+
+ System.out.println(this.getClass().getCanonicalName());
+ System.out.println(" seedEMin="+seedEMin);
+ System.out.println(" addEMin="+addEMin);
+ System.out.println();
+ }
+
+ public void detectorChanged(Detector detector)
+ {
+ ecal = (HPSEcal)detector.getSubdetector(ecalName);
+
+ System.out.println(ecal.getName());
+ System.out.println(" nx="+ecal.nx());
+ System.out.println(" ny="+ecal.ny());
+ System.out.println(" beamgap="+ecal.beamGap());
+ System.out.println(" dface="+ecal.distanceToFace());
+ }
+
+ public void process(EventHeader event)
+ {
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ if (hits == null)
+ throw new RuntimeException("Event is missing ECal hits collection!");
+
+ IDDecoder dec = ecal.getIDDecoder();
+
+ // Hit map.
+ Map<Long,CalorimeterHit> hitMap = new HashMap<Long,CalorimeterHit>();
+
+ // Make map of x, y, and side to hit.
+ for (CalorimeterHit hit : hits)
+ {
+ dec.setID(hit.getCellID());
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ // Setup an encoder from the decoder.
+ IDEncoder enc = new IDEncoder(dec.getIDDescription());
+
+ // Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ // Loop over ECal hits.
+ for (CalorimeterHit hit : hits)
+ {
+ // Cut on min seed E.
+ if (hit.getRawEnergy() < seedEMin)
+ continue;
+
+ // Set decoder's ID.
+ dec.setID(hit.getCellID());
+
+ // Get ID values.
+ int ix = dec.getValue("ix");
+ int iy = dec.getValue("iy");
+ int side = dec.getValue("side");
+
+ // Get neighboring crystals to hit.
+ Set<XYSide> neighbors = ecal.getNeighbors(ix, iy, side);
+
+ //System.out.println("neighbors of " + ix + ", " + iy + ", " + side);
+
+ int[] buffer = new int[dec.getFieldCount()];
+ dec.getValues(buffer);
+
+ // Loop over neighbors.
+ List<CalorimeterHit> neighborHits = new ArrayList<CalorimeterHit>();
+ boolean isSeed = true;
+ for (XYSide neighbor : neighbors)
+ {
+ buffer[dec.getFieldIndex("ix")] = neighbor.x();
+ buffer[dec.getFieldIndex("iy")] = neighbor.y();
+ buffer[dec.getFieldIndex("side")] = neighbor.side();
+
+ // Create ID for this crystal.
+ // FIXME This is redundant and inefficient.
+ long neighborId = enc.setValues(buffer);
+
+ // Find neighbor hit if it exists.
+ CalorimeterHit neighborHit = hitMap.get(neighborId);
+
+ // Got a hit?
+ if (neighborHit != null)
+ {
+ // Neighbor is hotter so skip this hit.
+ if (neighborHit.getRawEnergy() > hit.getRawEnergy())
+ {
+ isSeed = false;
+ break;
+ }
+
+ // Add to cluster if above min E.
+ if (neighborHit.getRawEnergy() >= addEMin)
+ neighborHits.add(neighborHit);
+ }
+ }
+
+ if (isSeed)
+ {
+ // Make cluster.
+ BasicCluster cluster = new BasicCluster();
+ cluster.addHit(hit);
+ for (CalorimeterHit clusHit : neighborHits)
+ {
+ cluster.addHit(clusHit);
+ }
+ clusters.add(cluster);
+ }
+ }
+
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, clusters, Cluster.class, flag);
+ }
+}
\ No newline at end of file
