6 modified files
hps-java/src/main/java/org/lcsim/hps/analysis/ecal
diff -u -r1.10 -r1.11
--- HPSEcalPlotsDriver.java 2 Sep 2011 21:16:20 -0000 1.10
+++ HPSEcalPlotsDriver.java 2 Sep 2011 21:17:49 -0000 1.11
@@ -1,475 +1,475 @@
-package org.lcsim.hps.analysis.ecal;
-
-import hep.aida.ICloud1D;
-import hep.aida.ICloud2D;
-import hep.aida.IHistogram1D;
-import hep.aida.IHistogram2D;
-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.Map.Entry;
-import java.util.Set;
-
-import org.lcsim.detector.identifier.IIdentifier;
-import org.lcsim.detector.identifier.IIdentifierHelper;
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.Cluster;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.util.ParticleTypeClassifier;
-import org.lcsim.hps.recon.ecal.HPSEcalCluster;
-import org.lcsim.hps.recon.ecal.HPSRawCalorimeterHit;
-import org.lcsim.hps.util.ClockSingleton;
-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.10 2011/09/02 21:16:20 meeg Exp $
- */
-public class HPSEcalPlotsDriver extends Driver {
-
- String ecalCollectionName = "EcalHits";
- String ecalClusterCollectionName = "EcalClusters";
- AIDA aida = AIDA.defaultInstance();
- // MCParticle plots.
- ICloud1D primaryEPlot;
- ICloud1D fsCountPlot;
- ICloud1D fsEPlot;
- ICloud1D fsGammaEPlot;
- ICloud1D fsElectronEPlot;
- ICloud1D fsPositronEPlot;
- // 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;
- ICloud1D crystalXPlot;
- ICloud1D crystalYPlot;
- //ICloud2D cellXYPlot;
- IHistogram2D crystalXYPlot;
- // Cluster plots.
- IHistogram1D nclusPlot;
- ICloud1D clusEPlot;
- ICloud1D clusTotEPlot;
- ICloud1D leadClusEPlot;
- ICloud1D leadClus2EPlot;
- //ICloud1D clusResTop3Plot;
- IHistogram1D clusHitPlot;
- ICloud1D clusSeedEPlot;
- ICloud1D clusSeedDistPlot;
- ICloud2D leadClusAndPrimaryPlot;
- IHistogram1D clusNHits;
- IHistogram2D hitXYPlot;
- int numTriggers = 0;
-
- 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 setEcalCollectionName(String ecalCollectionName) {
- this.ecalCollectionName = ecalCollectionName;
- }
-
- public void setEcalClusterCollectionName(String ecalClusterCollectionName) {
- this.ecalClusterCollectionName = ecalClusterCollectionName;
- }
-
- public void startOfData() {
- fsCountPlot = aida.cloud1D("MCParticle: Number of Final State Particles");
- fsCountPlot.annotation().addItem("xAxisLabel", "Number of FS Particles");
-
- fsEPlot = aida.cloud1D("MCParticle: FS Particle E");
- fsEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
-
- fsGammaEPlot = aida.cloud1D("MCParticle: FS Gamma E");
- fsGammaEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
-
- fsElectronEPlot = aida.cloud1D("MCParticle: FS Electron E");
- fsElectronEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
-
- fsPositronEPlot = aida.cloud1D("MCParticle: FS Positron E");
- fsPositronEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
-
- 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]");
-
- 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]");
-
- crystalXPlot = aida.cloud1D(ecalCollectionName + " : X Field Value");
- crystalXPlot.annotation().addItem("xAxisLabel", "Number of Entries");
-
- crystalYPlot = aida.cloud1D(ecalCollectionName + " : Y Field Value");
- crystalYPlot.annotation().addItem("xAxisLabel", "Number of Entries");
-
- crystalXYPlot = aida.histogram2D(
- ecalCollectionName + " : X & Y ID Values",
- 46, -23., 23., 5, 1., 6.);
-
- 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");
-
- 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", "Distance [mm]");
- clusSeedDistPlot.annotation().addItem("yAxisScale", "log");
-
- clusSeedEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Energy");
- clusSeedEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
-
- clusNHits = aida.histogram1D(ecalClusterCollectionName + " : Number of Hits per Cluster", 13, 0, 13);
- clusNHits.annotation().addItem("xAxisLabel", "Number of Hits");
-
- leadClusAndPrimaryPlot = aida.cloud2D(ecalClusterCollectionName + " : Lead Cluster E vs Highest Primary Particle E");
- }
-
- public void process(EventHeader event) {
- // MCParticles
- List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
-
- // Final State particles.
- List<MCParticle> fsParticles = makeGenFSParticleList(mcparticles);
-
- //System.out.println("fsParticles="+fsParticles.size());
- fsCountPlot.fill(fsParticles.size());
-
- for (MCParticle fs : fsParticles) {
- double fsE = fs.getEnergy();
- int fsPdg = fs.getPDGID();
- fsEPlot.fill(fsE);
- if (ParticleTypeClassifier.isElectron(fsPdg)) {
- fsElectronEPlot.fill(fsE);
- } else if (ParticleTypeClassifier.isPositron(fsPdg)) {
- fsPositronEPlot.fill(fsE);
- } else if (ParticleTypeClassifier.isPhoton(fsPdg)) {
- fsGammaEPlot.fill(fsE);
- }
- }
-
- // 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();
-
- // primary particle with most E
- MCParticle primary = getPrimary(mcparticles);
- double primaryE = primary.getEnergy();
- primaryEPlot.fill(primaryE);
-
- // event energy
- double eventE = getPrimaryE(mcparticles);
- eventEPlot.fill(eventE);
-
- List<HPSEcalCluster> clusters = event.get(HPSEcalCluster.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;
- boolean quadrants[] = new boolean[4];
- boolean trigger = false;
- for (HPSEcalCluster clus : clusters) {
- clusNHits.fill(clus.getCalorimeterHits().size());
-
- double e = clus.getEnergy();
- clusEPlot.fill(e);
- clusE += e;
- HPSRawCalorimeterHit seedHit = (HPSRawCalorimeterHit) clus.getSeedHit();
- //double maxe = 0;
- for (CalorimeterHit hit : clus.getCalorimeterHits()) {
- if (hitClusMap.containsKey(hit)) {
- int nshared = hitClusMap.get(hit);
- ++nshared;
-
- 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();
- clusSeedDistPlot.fill(Math.sqrt(y * y + z * z));
-
- // Seed E.
- clusSeedEPlot.fill(seedHit.getRawEnergy());
-
- quadrants[seedHit.getQuadrant() - 1] = true;
- }
-
- //TODO: change the quadrants back for correct X coordinate
- if ((quadrants[0] && quadrants[2]) || (quadrants[1] && quadrants[3])) {
- trigger = true;
- numTriggers++;
- }
-
- // 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());
- }
-
- // 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();
-
- // Get ID helper.
- IIdentifierHelper helper =
- event.getMetaData(hits).getIDDecoder().getSubdetector().getDetectorElement().getIdentifierHelper();
-
-
-
- // Check unique IDs.
- Set<Long> ids = new HashSet<Long>();
-
- // n hits
- hitCountPlot.fill(nhits);
-
- int nhits100MeV = 0;
- int nhitsOver100MeV = 0;
-
- if (trigger && numTriggers<=100) {
- hitXYPlot = aida.histogram2D(
- ecalCollectionName + " : hit E, event " + String.format("%07d", ClockSingleton.getClock()),
- 47, -23.5, 23.5, 11, -5.5, 5.5);
- }
-
- // 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++;
- }
-
- // 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());
-
- // Add to ECal energy sum.
- esum += eraw;
-
- // X and Y identifier values.
- IIdentifier id = hit.getIdentifier();
- int ix = helper.unpack(id).getValue(helper.getFieldIndex("ix"));
- int iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
- crystalXPlot.fill(ix);
- crystalYPlot.fill(iy);
- crystalXYPlot.fill(ix, iy);
- if (trigger && numTriggers<=100)
- hitXYPlot.fill(ix, iy, eraw);
- }
-
- hitUnder100MeVPlot.fill(nhits100MeV);
- 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());
- }
-
- 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> makeGenFSParticleList(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 endOfData() {
- System.out.printf("Trigger count: %d\n", numTriggers);
- }
-}
+package org.lcsim.hps.analysis.ecal;
+
+import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+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.Map.Entry;
+import java.util.Set;
+
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.hps.recon.ecal.HPSRawCalorimeterHit;
+import org.lcsim.hps.util.ClockSingleton;
+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.11 2011/09/02 21:17:49 meeg Exp $
+ */
+public class HPSEcalPlotsDriver extends Driver {
+
+ String ecalCollectionName = "EcalHits";
+ String ecalClusterCollectionName = "EcalClusters";
+ AIDA aida = AIDA.defaultInstance();
+ // MCParticle plots.
+ ICloud1D primaryEPlot;
+ ICloud1D fsCountPlot;
+ ICloud1D fsEPlot;
+ ICloud1D fsGammaEPlot;
+ ICloud1D fsElectronEPlot;
+ ICloud1D fsPositronEPlot;
+ // 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;
+ ICloud1D crystalXPlot;
+ ICloud1D crystalYPlot;
+ //ICloud2D cellXYPlot;
+ IHistogram2D crystalXYPlot;
+ // Cluster plots.
+ IHistogram1D nclusPlot;
+ ICloud1D clusEPlot;
+ ICloud1D clusTotEPlot;
+ ICloud1D leadClusEPlot;
+ ICloud1D leadClus2EPlot;
+ //ICloud1D clusResTop3Plot;
+ IHistogram1D clusHitPlot;
+ ICloud1D clusSeedEPlot;
+ ICloud1D clusSeedDistPlot;
+ ICloud2D leadClusAndPrimaryPlot;
+ IHistogram1D clusNHits;
+ IHistogram2D hitXYPlot;
+ int numTriggers = 0;
+
+ 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 setEcalCollectionName(String ecalCollectionName) {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setEcalClusterCollectionName(String ecalClusterCollectionName) {
+ this.ecalClusterCollectionName = ecalClusterCollectionName;
+ }
+
+ public void startOfData() {
+ fsCountPlot = aida.cloud1D("MCParticle: Number of Final State Particles");
+ fsCountPlot.annotation().addItem("xAxisLabel", "Number of FS Particles");
+
+ fsEPlot = aida.cloud1D("MCParticle: FS Particle E");
+ fsEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
+
+ fsGammaEPlot = aida.cloud1D("MCParticle: FS Gamma E");
+ fsGammaEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
+
+ fsElectronEPlot = aida.cloud1D("MCParticle: FS Electron E");
+ fsElectronEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
+
+ fsPositronEPlot = aida.cloud1D("MCParticle: FS Positron E");
+ fsPositronEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
+
+ 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]");
+
+ 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]");
+
+ crystalXPlot = aida.cloud1D(ecalCollectionName + " : X Field Value");
+ crystalXPlot.annotation().addItem("xAxisLabel", "Number of Entries");
+
+ crystalYPlot = aida.cloud1D(ecalCollectionName + " : Y Field Value");
+ crystalYPlot.annotation().addItem("xAxisLabel", "Number of Entries");
+
+ crystalXYPlot = aida.histogram2D(
+ ecalCollectionName + " : X & Y ID Values",
+ 46, -23., 23., 5, 1., 6.);
+
+ 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");
+
+ 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", "Distance [mm]");
+ clusSeedDistPlot.annotation().addItem("yAxisScale", "log");
+
+ clusSeedEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Energy");
+ clusSeedEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusNHits = aida.histogram1D(ecalClusterCollectionName + " : Number of Hits per Cluster", 13, 0, 13);
+ clusNHits.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ leadClusAndPrimaryPlot = aida.cloud2D(ecalClusterCollectionName + " : Lead Cluster E vs Highest Primary Particle E");
+ }
+
+ public void process(EventHeader event) {
+ // MCParticles
+ List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
+
+ // Final State particles.
+ List<MCParticle> fsParticles = makeGenFSParticleList(mcparticles);
+
+ //System.out.println("fsParticles="+fsParticles.size());
+ fsCountPlot.fill(fsParticles.size());
+
+ for (MCParticle fs : fsParticles) {
+ double fsE = fs.getEnergy();
+ int fsPdg = fs.getPDGID();
+ fsEPlot.fill(fsE);
+ if (ParticleTypeClassifier.isElectron(fsPdg)) {
+ fsElectronEPlot.fill(fsE);
+ } else if (ParticleTypeClassifier.isPositron(fsPdg)) {
+ fsPositronEPlot.fill(fsE);
+ } else if (ParticleTypeClassifier.isPhoton(fsPdg)) {
+ fsGammaEPlot.fill(fsE);
+ }
+ }
+
+ // 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();
+
+ // primary particle with most E
+ MCParticle primary = getPrimary(mcparticles);
+ double primaryE = primary.getEnergy();
+ primaryEPlot.fill(primaryE);
+
+ // event energy
+ double eventE = getPrimaryE(mcparticles);
+ eventEPlot.fill(eventE);
+
+ List<HPSEcalCluster> clusters = event.get(HPSEcalCluster.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;
+ boolean quadrants[] = new boolean[4];
+ boolean trigger = false;
+ for (HPSEcalCluster clus : clusters) {
+ clusNHits.fill(clus.getCalorimeterHits().size());
+
+ double e = clus.getEnergy();
+ clusEPlot.fill(e);
+ clusE += e;
+ HPSRawCalorimeterHit seedHit = (HPSRawCalorimeterHit) clus.getSeedHit();
+ //double maxe = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ if (hitClusMap.containsKey(hit)) {
+ int nshared = hitClusMap.get(hit);
+ ++nshared;
+
+ 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();
+ clusSeedDistPlot.fill(Math.sqrt(y * y + z * z));
+
+ // Seed E.
+ clusSeedEPlot.fill(seedHit.getRawEnergy());
+
+ quadrants[seedHit.getQuadrant() - 1] = true;
+ }
+
+ //TODO: change the quadrants back for correct X coordinate
+ if ((quadrants[0] && quadrants[2]) || (quadrants[1] && quadrants[3])) {
+ trigger = true;
+ numTriggers++;
+ }
+
+ // 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());
+ }
+
+ // 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();
+
+ // Get ID helper.
+ IIdentifierHelper helper =
+ event.getMetaData(hits).getIDDecoder().getSubdetector().getDetectorElement().getIdentifierHelper();
+
+
+
+ // Check unique IDs.
+ Set<Long> ids = new HashSet<Long>();
+
+ // n hits
+ hitCountPlot.fill(nhits);
+
+ int nhits100MeV = 0;
+ int nhitsOver100MeV = 0;
+
+ if (trigger && numTriggers<=100) {
+ hitXYPlot = aida.histogram2D(
+ ecalCollectionName + " : hit E, event " + String.format("%07d", ClockSingleton.getClock()),
+ 47, -23.5, 23.5, 11, -5.5, 5.5);
+ }
+
+ // 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++;
+ }
+
+ // 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());
+
+ // Add to ECal energy sum.
+ esum += eraw;
+
+ // X and Y identifier values.
+ IIdentifier id = hit.getIdentifier();
+ int ix = helper.unpack(id).getValue(helper.getFieldIndex("ix"));
+ int iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
+ crystalXPlot.fill(ix);
+ crystalYPlot.fill(iy);
+ crystalXYPlot.fill(ix, iy);
+ if (trigger && numTriggers<=100)
+ hitXYPlot.fill(ix, iy, eraw);
+ }
+
+ hitUnder100MeVPlot.fill(nhits100MeV);
+ 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());
+ }
+
+ 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> makeGenFSParticleList(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 endOfData() {
+ System.out.printf("Trigger count: %d\n", numTriggers);
+ }
+}
hps-java/src/main/java/org/lcsim/hps/recon/ecal
diff -u -r1.14 -r1.15
--- HPSEcal1BitClusterer.java 2 Sep 2011 21:16:20 -0000 1.14
+++ HPSEcal1BitClusterer.java 2 Sep 2011 21:17:49 -0000 1.15
@@ -1,320 +1,320 @@
-package org.lcsim.hps.recon.ecal;
-
-import java.util.ArrayList;
-import java.util.Collection;
-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.HPSEcal3.NeighborMap;
-import org.lcsim.geometry.subdetector.HPSEcal3;
-import org.lcsim.util.Driver;
-import org.lcsim.util.lcio.LCIOConstants;
-import org.lcsim.recon.cluster.util.Clusterer;
-
-/**
- * Creates clusters from CalorimeterHits in the HPSEcal detector.
- *
- * The clustering algorithm is from JLab Hall B 6 GeV DVCS Trigger Design doc.
- *
- * @author Jeremy McCormick <[log in to unmask]>
- * @author Tim Nelson <[log in to unmask]>
- * @author Sho Uemura <[log in to unmask]>
- * @version $Id: HPSEcal1BitClusterer.java,v 1.14 2011/09/02 21:16:20 meeg Exp $
- */
-public class HPSEcal1BitClusterer extends Driver implements Clusterer {
-
- HPSEcal3 ecal;
- String ecalName;
- String ecalCollectionName;
- String clusterCollectionName = "EcalClusters";
- // Threshold E for discriminator.
- double hitEMin = .05;
- // Threshold hit count for cluster. Cluster must have more than this many hits.
- int clusterThreshold = 5;
- // Map of crystals to their neighbors.
- NeighborMap neighborMap = null;
-
- public HPSEcal1BitClusterer() {
- }
-
- public void setClusterCollectionName(String clusterCollectionName) {
- this.clusterCollectionName = clusterCollectionName;
- }
-
- public void setEcalCollectionName(String ecalCollectionName) {
- this.ecalCollectionName = ecalCollectionName;
- }
-
- public void setHitEMin(double hitEMin) {
- this.hitEMin = hitEMin;
- }
-
- public void setClusterThreshold(int clusterThreshold) {
- this.clusterThreshold = clusterThreshold;
- }
-
- public void setEcalName(String ecalName) {
- this.ecalName = ecalName;
- }
-
- public void startOfData() {
- if (ecalCollectionName == null) {
- throw new RuntimeException("The parameter ecalCollectionName was not set!");
- }
-
- if (ecalName == null) {
- throw new RuntimeException("The parameter ecalName was not set!");
- }
- }
-
- public void detectorChanged(Detector detector) {
- // Get the Subdetector.
- ecal = (HPSEcal3) detector.getSubdetector(ecalName);
-
- // Cache ref to neighbor map.
- neighborMap = ecal.getNeighborMap();
-
- //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());
-
- //System.out.println(neighborMap.toString());
- }
-
- public void process(EventHeader event) {
- //System.out.println(this.getClass().getCanonicalName() + " - process");
-
- // Get the list of raw ECal hits.
- List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
- if (hits == null)
- throw new RuntimeException("Event is missing ECal raw hits collection!");
-
- // Put Cluster collection into event.
- int flag = 1 << LCIOConstants.CLBIT_HITS;
- event.put(clusterCollectionName, createClusters(hits), Cluster.class, flag);
- }
-
- public List<Cluster> createClusters(List<CalorimeterHit> hits) {
- // Hit map.
- Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
-
- // Make a hit map for quick lookup by ID.
- for (CalorimeterHit hit : hits) {
- hitMap.put(hit.getCellID(), hit);
- }
-
- return createClusters(hitMap);
- }
-
- public List<Cluster> createClusters(Map<Long, CalorimeterHit> map) {
- // New Cluster list to be added to event.
- List<Cluster> clusters = new ArrayList<Cluster>();
-
- // Get the decoder for the ECal IDs.
- IDDecoder dec = ecal.getIDDecoder();
-
- Map<Long, Integer> hitCounts = new HashMap<Long, Integer>();
-
- Collection<CalorimeterHit> hitCollection = map.values();
- // Loop over ECal hits to count hit towers
- for (CalorimeterHit hit : hitCollection) {
- // Cut on min seed E.
- if (hit.getRawEnergy() < hitEMin) {
- continue;
- }
-
- // Get neighbor crystal IDs.
- Set<Long> neighbors = neighborMap.get(hit.getCellID());
-
- if (neighbors == null) {
- throw new RuntimeException("Oops! Set of neighbors is null!");
- }
-
- Integer hitCount = hitCounts.get(hit.getCellID());
- if (hitCount == null)
- hitCounts.put(hit.getCellID(), 1);
- else
- hitCounts.put(hit.getCellID(), hitCount + 1);
-
- // Loop over neighbors to make hit list for cluster.
- for (Long neighborId : neighbors) {
- hitCount = hitCounts.get(neighborId);
- if (hitCount == null)
- hitCounts.put(neighborId, 1);
- else
- hitCounts.put(neighborId, hitCount + 1);
- }
- }
-
- //for each crystal with a nonzero hit count, test for cluster
- for (Long possibleCluster : hitCounts.keySet()) {
- //System.out.printf("%d, %d hits\n",possibleCluster,hitCounts.get(possibleCluster));
- if (hitCounts.get(possibleCluster) <= clusterThreshold) {
- continue;
- }
-
- // Get neighbor crystal IDs.
- Set<Long> neighbors = neighborMap.get(possibleCluster);
-
- if (neighbors == null) {
- throw new RuntimeException("Oops! Set of neighbors is null!");
- }
-
- //Apply peak detector scheme.
- // Set the ID.
- dec.setID(possibleCluster);
- // Get ID field values.
- int x1 = dec.getValue("ix");
- int y1 = dec.getValue("iy");
- Integer hitCount = hitCounts.get(possibleCluster);
-
- //System.out.printf("Possible cluster: x: %d, y: %d, hits: %d\n",x1,y1,hitCount);
- boolean isCluster = true;
- for (Long neighborId : neighbors) {
- // Set the ID.
- dec.setID(neighborId);
- Integer neighborHitCount = hitCounts.get(neighborId);
- if (neighborHitCount == null)
- continue;
- // Get ID field values.
- int x2 = dec.getValue("ix");
- int y2 = dec.getValue("iy");
- if (y1>0) {
- if (x1 > 0) {
- //quadrant 1
- if (x1 == 1) {
- //special case: left edge of quadrant
- if (x1 > x2 && y1 < y2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (x1 > x2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else if (x1 < x2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (y1 < y2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else {
- if (hitCount > neighborHitCount)
- continue;
- }
- } else {
- if (x1 > x2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else if (x1 < x2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (y1 < y2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else {
- if (hitCount > neighborHitCount)
- continue;
- }
- }
- } else {
- //quadrant 2
- if (y1 > y2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else if (y1 < y2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (x1 > x2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else {
- if (hitCount > neighborHitCount)
- continue;
- }
- }
- } else {
- if (x1 < 0) {
- //quadrant 3
- if (x1 == 1) {
- //special case: left edge of quadrant
- if (x1 < x2 && y1 > y2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (x1 < x2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else if (x1 > x2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (y1 > y2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else {
- if (hitCount > neighborHitCount)
- continue;
- }
- } else {
- if (x1 < x2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else if (x1 > x2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (y1 > y2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else {
- if (hitCount > neighborHitCount)
- continue;
- }
- }
- } else {
- //quadrant 4
- if (y1 < y2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else if (y1 > y2) {
- if (hitCount > neighborHitCount)
- continue;
- } else if (x1 < x2) {
- if (hitCount >= neighborHitCount)
- continue;
- } else {
- if (hitCount > neighborHitCount)
- continue;
- }
- }
- }
- isCluster = false;
- break;
- }
-
- if (isCluster) {
- //System.out.printf("Cluster: x: %d, y: %d, side:%d, hits: %d\n",x1,y1,side1,hitCount);
- HPSEcalCluster cluster = new HPSEcalCluster();
- CalorimeterHit hit = map.get(possibleCluster);
- if (hit != null && hit.getRawEnergy() > hitEMin) {
- cluster.addSeedHit(hit);
- }
-
- for (Long neighborId : neighbors) {
- // Find the neighbor hit in the event if it exists.
- hit = map.get(neighborId);
- if (hit != null && hit.getRawEnergy() > hitEMin) {
- cluster.addHit(hit);
- }
- }
- clusters.add(cluster);
- }
- }
- return clusters;
- }
+package org.lcsim.hps.recon.ecal;
+
+import java.util.ArrayList;
+import java.util.Collection;
+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.HPSEcal3.NeighborMap;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.util.Driver;
+import org.lcsim.util.lcio.LCIOConstants;
+import org.lcsim.recon.cluster.util.Clusterer;
+
+/**
+ * Creates clusters from CalorimeterHits in the HPSEcal detector.
+ *
+ * The clustering algorithm is from JLab Hall B 6 GeV DVCS Trigger Design doc.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Tim Nelson <[log in to unmask]>
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: HPSEcal1BitClusterer.java,v 1.15 2011/09/02 21:17:49 meeg Exp $
+ */
+public class HPSEcal1BitClusterer extends Driver implements Clusterer {
+
+ HPSEcal3 ecal;
+ String ecalName;
+ String ecalCollectionName;
+ String clusterCollectionName = "EcalClusters";
+ // Threshold E for discriminator.
+ double hitEMin = .05;
+ // Threshold hit count for cluster. Cluster must have more than this many hits.
+ int clusterThreshold = 5;
+ // Map of crystals to their neighbors.
+ NeighborMap neighborMap = null;
+
+ public HPSEcal1BitClusterer() {
+ }
+
+ public void setClusterCollectionName(String clusterCollectionName) {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ public void setEcalCollectionName(String ecalCollectionName) {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setHitEMin(double hitEMin) {
+ this.hitEMin = hitEMin;
+ }
+
+ public void setClusterThreshold(int clusterThreshold) {
+ this.clusterThreshold = clusterThreshold;
+ }
+
+ public void setEcalName(String ecalName) {
+ this.ecalName = ecalName;
+ }
+
+ public void startOfData() {
+ if (ecalCollectionName == null) {
+ throw new RuntimeException("The parameter ecalCollectionName was not set!");
+ }
+
+ if (ecalName == null) {
+ throw new RuntimeException("The parameter ecalName was not set!");
+ }
+ }
+
+ public void detectorChanged(Detector detector) {
+ // Get the Subdetector.
+ ecal = (HPSEcal3) detector.getSubdetector(ecalName);
+
+ // Cache ref to neighbor map.
+ neighborMap = ecal.getNeighborMap();
+
+ //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());
+
+ //System.out.println(neighborMap.toString());
+ }
+
+ public void process(EventHeader event) {
+ //System.out.println(this.getClass().getCanonicalName() + " - process");
+
+ // Get the list of raw ECal hits.
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ if (hits == null)
+ throw new RuntimeException("Event is missing ECal raw hits collection!");
+
+ // Put Cluster collection into event.
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, createClusters(hits), Cluster.class, flag);
+ }
+
+ public List<Cluster> createClusters(List<CalorimeterHit> hits) {
+ // Hit map.
+ Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
+
+ // Make a hit map for quick lookup by ID.
+ for (CalorimeterHit hit : hits) {
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ return createClusters(hitMap);
+ }
+
+ public List<Cluster> createClusters(Map<Long, CalorimeterHit> map) {
+ // New Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ // Get the decoder for the ECal IDs.
+ IDDecoder dec = ecal.getIDDecoder();
+
+ Map<Long, Integer> hitCounts = new HashMap<Long, Integer>();
+
+ Collection<CalorimeterHit> hitCollection = map.values();
+ // Loop over ECal hits to count hit towers
+ for (CalorimeterHit hit : hitCollection) {
+ // Cut on min seed E.
+ if (hit.getRawEnergy() < hitEMin) {
+ continue;
+ }
+
+ // Get neighbor crystal IDs.
+ Set<Long> neighbors = neighborMap.get(hit.getCellID());
+
+ if (neighbors == null) {
+ throw new RuntimeException("Oops! Set of neighbors is null!");
+ }
+
+ Integer hitCount = hitCounts.get(hit.getCellID());
+ if (hitCount == null)
+ hitCounts.put(hit.getCellID(), 1);
+ else
+ hitCounts.put(hit.getCellID(), hitCount + 1);
+
+ // Loop over neighbors to make hit list for cluster.
+ for (Long neighborId : neighbors) {
+ hitCount = hitCounts.get(neighborId);
+ if (hitCount == null)
+ hitCounts.put(neighborId, 1);
+ else
+ hitCounts.put(neighborId, hitCount + 1);
+ }
+ }
+
+ //for each crystal with a nonzero hit count, test for cluster
+ for (Long possibleCluster : hitCounts.keySet()) {
+ //System.out.printf("%d, %d hits\n",possibleCluster,hitCounts.get(possibleCluster));
+ if (hitCounts.get(possibleCluster) <= clusterThreshold) {
+ continue;
+ }
+
+ // Get neighbor crystal IDs.
+ Set<Long> neighbors = neighborMap.get(possibleCluster);
+
+ if (neighbors == null) {
+ throw new RuntimeException("Oops! Set of neighbors is null!");
+ }
+
+ //Apply peak detector scheme.
+ // Set the ID.
+ dec.setID(possibleCluster);
+ // Get ID field values.
+ int x1 = dec.getValue("ix");
+ int y1 = dec.getValue("iy");
+ Integer hitCount = hitCounts.get(possibleCluster);
+
+ //System.out.printf("Possible cluster: x: %d, y: %d, hits: %d\n",x1,y1,hitCount);
+ boolean isCluster = true;
+ for (Long neighborId : neighbors) {
+ // Set the ID.
+ dec.setID(neighborId);
+ Integer neighborHitCount = hitCounts.get(neighborId);
+ if (neighborHitCount == null)
+ continue;
+ // Get ID field values.
+ int x2 = dec.getValue("ix");
+ int y2 = dec.getValue("iy");
+ if (y1>0) {
+ if (x1 > 0) {
+ //quadrant 1
+ if (x1 == 1) {
+ //special case: left edge of quadrant
+ if (x1 > x2 && y1 < y2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (x1 > x2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else if (x1 < x2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (y1 < y2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else {
+ if (hitCount > neighborHitCount)
+ continue;
+ }
+ } else {
+ if (x1 > x2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else if (x1 < x2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (y1 < y2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else {
+ if (hitCount > neighborHitCount)
+ continue;
+ }
+ }
+ } else {
+ //quadrant 2
+ if (y1 > y2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else if (y1 < y2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (x1 > x2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else {
+ if (hitCount > neighborHitCount)
+ continue;
+ }
+ }
+ } else {
+ if (x1 < 0) {
+ //quadrant 3
+ if (x1 == 1) {
+ //special case: left edge of quadrant
+ if (x1 < x2 && y1 > y2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (x1 < x2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else if (x1 > x2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (y1 > y2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else {
+ if (hitCount > neighborHitCount)
+ continue;
+ }
+ } else {
+ if (x1 < x2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else if (x1 > x2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (y1 > y2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else {
+ if (hitCount > neighborHitCount)
+ continue;
+ }
+ }
+ } else {
+ //quadrant 4
+ if (y1 < y2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else if (y1 > y2) {
+ if (hitCount > neighborHitCount)
+ continue;
+ } else if (x1 < x2) {
+ if (hitCount >= neighborHitCount)
+ continue;
+ } else {
+ if (hitCount > neighborHitCount)
+ continue;
+ }
+ }
+ }
+ isCluster = false;
+ break;
+ }
+
+ if (isCluster) {
+ //System.out.printf("Cluster: x: %d, y: %d, side:%d, hits: %d\n",x1,y1,side1,hitCount);
+ HPSEcalCluster cluster = new HPSEcalCluster();
+ CalorimeterHit hit = map.get(possibleCluster);
+ if (hit != null && hit.getRawEnergy() > hitEMin) {
+ cluster.addSeedHit(hit);
+ }
+
+ for (Long neighborId : neighbors) {
+ // Find the neighbor hit in the event if it exists.
+ hit = map.get(neighborId);
+ if (hit != null && hit.getRawEnergy() > hitEMin) {
+ cluster.addHit(hit);
+ }
+ }
+ clusters.add(cluster);
+ }
+ }
+ return clusters;
+ }
}
\ No newline at end of file
hps-java/src/main/java/org/lcsim/hps/recon/ecal
diff -u -r1.12 -r1.13
--- HPSEcalClusterer.java 2 Sep 2011 21:16:20 -0000 1.12
+++ HPSEcalClusterer.java 2 Sep 2011 21:17:49 -0000 1.13
@@ -1,175 +1,175 @@
-package org.lcsim.hps.recon.ecal;
-
-import java.util.ArrayList;
-import java.util.Collection;
-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.subdetector.HPSEcal3.NeighborMap;
-import org.lcsim.geometry.subdetector.HPSEcal3;
-import org.lcsim.recon.cluster.util.Clusterer;
-import org.lcsim.util.Driver;
-import org.lcsim.util.lcio.LCIOConstants;
-
-/**
- * This Driver creates clusters from the CalorimeterHits of an
- * {@link org.lcsim.geometry.subdetectur.HPSEcal3} detector.
- *
- * The clustering algorithm is from pages 83 and 84 of the HPS Proposal.
- *
- * @author Jeremy McCormick <[log in to unmask]>
- * @author Tim Nelson <[log in to unmask]>
- *
- * @version $Id: HPSEcalClusterer.java,v 1.12 2011/09/02 21:16:20 meeg Exp $
- */
-public class HPSEcalClusterer extends Driver implements Clusterer {
-
- HPSEcal3 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;
- // Map of crystals to their neighbors.
- NeighborMap neighborMap = null;
-
- 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 (ecalName == null)
- throw new RuntimeException("The parameter ecalName was not set!");
- }
-
- public void detectorChanged(Detector detector) {
- // Get the Subdetector.
- ecal = (HPSEcal3) detector.getSubdetector(ecalName);
-
- // Cache ref to neighbor map.
- neighborMap = ecal.getNeighborMap();
-
- //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());
-
- //System.out.println(neighborMap.toString());
- }
-
- public void process(EventHeader event) {
- //System.out.println(this.getClass().getCanonicalName() + " - process");
-
- // Get the list of raw ECal hits.
- List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
- if (hits == null)
- throw new RuntimeException("Event is missing ECal raw hits collection!");
-
- // Put Cluster collection into event.
- int flag = 1 << LCIOConstants.CLBIT_HITS;
- event.put(clusterCollectionName, createClusters(hits), Cluster.class, flag);
- }
-
- public List<Cluster> createClusters(List<CalorimeterHit> hits) {
- // Hit map.
- Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
-
- // Make a hit map for quick lookup by ID.
- for (CalorimeterHit hit : hits) {
- hitMap.put(hit.getCellID(), hit);
- }
-
- return createClusters(hitMap);
- }
-
- public List<Cluster> createClusters(Map<Long, CalorimeterHit> map) {
-
- // New Cluster list to be added to event.
- List<Cluster> clusters = new ArrayList<Cluster>();
-
- Collection<CalorimeterHit> hits = map.values();
-
- // Loop over ECal hits to find cluster seeds.
- for (CalorimeterHit hit : hits) {
- // Cut on min seed E.
- if (hit.getRawEnergy() < seedEMin)
- continue;
-
- // Get neighbor crystal IDs.
- Set<Long> neighbors = neighborMap.get(hit.getCellID());
-
- if (neighbors == null)
- throw new RuntimeException("Oops! Set of neighbors is null!");
-
- // List for neighboring hits.
- List<CalorimeterHit> neighborHits = new ArrayList<CalorimeterHit>();
-
- // Loop over neighbors to make hit list for cluster.
- boolean isSeed = true;
- for (Long neighborId : neighbors) {
- // Find the neighbor hit in the event if it exists.
- CalorimeterHit neighborHit = map.get(neighborId);
-
- // Was this neighbor cell hit?
- if (neighborHit != null) {
- // Check if neighbor cell has more energy.
- if (neighborHit.getRawEnergy() > hit.getRawEnergy()) {
- // Neighbor has more energy, so cell is not a seed.
- isSeed = false;
- break;
- }
-
- // Add to cluster if above min E.
- if (neighborHit.getRawEnergy() >= addEMin) {
- neighborHits.add(neighborHit);
- }
- }
- }
-
- // Did we find a seed?
- if (isSeed) {
- // Make a cluster from the hit list.
- HPSEcalCluster cluster = new HPSEcalCluster();
- cluster.addSeedHit(hit);
- for (CalorimeterHit clusHit : neighborHits) {
- cluster.addHit(clusHit);
- }
- clusters.add(cluster);
- }
- }
- return clusters;
- }
+package org.lcsim.hps.recon.ecal;
+
+import java.util.ArrayList;
+import java.util.Collection;
+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.subdetector.HPSEcal3.NeighborMap;
+import org.lcsim.geometry.subdetector.HPSEcal3;
+import org.lcsim.recon.cluster.util.Clusterer;
+import org.lcsim.util.Driver;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ * This Driver creates clusters from the CalorimeterHits of an
+ * {@link org.lcsim.geometry.subdetectur.HPSEcal3} detector.
+ *
+ * The clustering algorithm is from pages 83 and 84 of the HPS Proposal.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @author Tim Nelson <[log in to unmask]>
+ *
+ * @version $Id: HPSEcalClusterer.java,v 1.13 2011/09/02 21:17:49 meeg Exp $
+ */
+public class HPSEcalClusterer extends Driver implements Clusterer {
+
+ HPSEcal3 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;
+ // Map of crystals to their neighbors.
+ NeighborMap neighborMap = null;
+
+ 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 (ecalName == null)
+ throw new RuntimeException("The parameter ecalName was not set!");
+ }
+
+ public void detectorChanged(Detector detector) {
+ // Get the Subdetector.
+ ecal = (HPSEcal3) detector.getSubdetector(ecalName);
+
+ // Cache ref to neighbor map.
+ neighborMap = ecal.getNeighborMap();
+
+ //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());
+
+ //System.out.println(neighborMap.toString());
+ }
+
+ public void process(EventHeader event) {
+ //System.out.println(this.getClass().getCanonicalName() + " - process");
+
+ // Get the list of raw ECal hits.
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ if (hits == null)
+ throw new RuntimeException("Event is missing ECal raw hits collection!");
+
+ // Put Cluster collection into event.
+ int flag = 1 << LCIOConstants.CLBIT_HITS;
+ event.put(clusterCollectionName, createClusters(hits), Cluster.class, flag);
+ }
+
+ public List<Cluster> createClusters(List<CalorimeterHit> hits) {
+ // Hit map.
+ Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
+
+ // Make a hit map for quick lookup by ID.
+ for (CalorimeterHit hit : hits) {
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ return createClusters(hitMap);
+ }
+
+ public List<Cluster> createClusters(Map<Long, CalorimeterHit> map) {
+
+ // New Cluster list to be added to event.
+ List<Cluster> clusters = new ArrayList<Cluster>();
+
+ Collection<CalorimeterHit> hits = map.values();
+
+ // Loop over ECal hits to find cluster seeds.
+ for (CalorimeterHit hit : hits) {
+ // Cut on min seed E.
+ if (hit.getRawEnergy() < seedEMin)
+ continue;
+
+ // Get neighbor crystal IDs.
+ Set<Long> neighbors = neighborMap.get(hit.getCellID());
+
+ if (neighbors == null)
+ throw new RuntimeException("Oops! Set of neighbors is null!");
+
+ // List for neighboring hits.
+ List<CalorimeterHit> neighborHits = new ArrayList<CalorimeterHit>();
+
+ // Loop over neighbors to make hit list for cluster.
+ boolean isSeed = true;
+ for (Long neighborId : neighbors) {
+ // Find the neighbor hit in the event if it exists.
+ CalorimeterHit neighborHit = map.get(neighborId);
+
+ // Was this neighbor cell hit?
+ if (neighborHit != null) {
+ // Check if neighbor cell has more energy.
+ if (neighborHit.getRawEnergy() > hit.getRawEnergy()) {
+ // Neighbor has more energy, so cell is not a seed.
+ isSeed = false;
+ break;
+ }
+
+ // Add to cluster if above min E.
+ if (neighborHit.getRawEnergy() >= addEMin) {
+ neighborHits.add(neighborHit);
+ }
+ }
+ }
+
+ // Did we find a seed?
+ if (isSeed) {
+ // Make a cluster from the hit list.
+ HPSEcalCluster cluster = new HPSEcalCluster();
+ cluster.addSeedHit(hit);
+ for (CalorimeterHit clusHit : neighborHits) {
+ cluster.addHit(clusHit);
+ }
+ clusters.add(cluster);
+ }
+ }
+ return clusters;
+ }
}
\ No newline at end of file
hps-java/src/main/java/org/lcsim/hps/recon/ecal
diff -u -r1.5 -r1.6
--- HPSEcalSimpleReadoutDriver.java 2 Sep 2011 21:16:20 -0000 1.5
+++ HPSEcalSimpleReadoutDriver.java 2 Sep 2011 21:17:49 -0000 1.6
@@ -1,58 +1,58 @@
-package org.lcsim.hps.recon.ecal;
-
-import java.util.ArrayList;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Map;
-
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.geometry.IDDecoder;
-
-/**
- * Performs readout of ECal hits.
- * No time evolution - this just integrates all hits in a cycle.
- *
- * @author Sho Uemura <[log in to unmask]>
- * @version $Id: HPSEcalSimpleReadoutDriver.java,v 1.5 2011/09/02 21:16:20 meeg Exp $
- */
-public class HPSEcalSimpleReadoutDriver extends HPSEcalReadoutDriver {
-
- //buffer for deposited energy
- Map<Long, Double> eDepMap = null;
-
- public List<HPSRawCalorimeterHit> readHits() {
- IDDecoder dec = ecal.getIDDecoder();
- List<HPSRawCalorimeterHit> hitList = new ArrayList<HPSRawCalorimeterHit>();
- for (Long cellID : eDepMap.keySet()) {
- dec.setID(cellID);
-// int ix = dec.getValue("ix");
-// int iy = dec.getValue("iy");
-// //temporary hack to disable crystals and flip X coordinate
-// int side = dec.getValue("side");
-// if (iy == 1 && ix*side >= -10 && ix*side <= -2)
-// continue;
- if (eDepMap.get(cellID) > threshold)
- hitList.add(new HPSRawCalorimeterHit(eDepMap.get(cellID), dec.getPosition(), 0.0, cellID, hitType));
- }
- //reset hit integration
- eDepMap = new HashMap<Long, Double>();
- return hitList;
- }
-
- public void putHits(List<CalorimeterHit> hits) {
- //fill the readout buffers
- for (CalorimeterHit hit : hits) {
- Double eDep = eDepMap.get(hit.getCellID());
- if (eDep == null) {
- eDepMap.put(hit.getCellID(), hit.getRawEnergy());
- } else {
- eDepMap.put(hit.getCellID(), eDep + hit.getRawEnergy());
- }
- }
- }
-
- public void initReadout() {
- //initialize buffers
- eDepMap = new HashMap<Long, Double>();
- }
+package org.lcsim.hps.recon.ecal;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+
+/**
+ * Performs readout of ECal hits.
+ * No time evolution - this just integrates all hits in a cycle.
+ *
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: HPSEcalSimpleReadoutDriver.java,v 1.6 2011/09/02 21:17:49 meeg Exp $
+ */
+public class HPSEcalSimpleReadoutDriver extends HPSEcalReadoutDriver {
+
+ //buffer for deposited energy
+ Map<Long, Double> eDepMap = null;
+
+ public List<HPSRawCalorimeterHit> readHits() {
+ IDDecoder dec = ecal.getIDDecoder();
+ List<HPSRawCalorimeterHit> hitList = new ArrayList<HPSRawCalorimeterHit>();
+ for (Long cellID : eDepMap.keySet()) {
+ dec.setID(cellID);
+// int ix = dec.getValue("ix");
+// int iy = dec.getValue("iy");
+// //temporary hack to disable crystals and flip X coordinate
+// int side = dec.getValue("side");
+// if (iy == 1 && ix*side >= -10 && ix*side <= -2)
+// continue;
+ if (eDepMap.get(cellID) > threshold)
+ hitList.add(new HPSRawCalorimeterHit(eDepMap.get(cellID), dec.getPosition(), 0.0, cellID, hitType));
+ }
+ //reset hit integration
+ eDepMap = new HashMap<Long, Double>();
+ return hitList;
+ }
+
+ public void putHits(List<CalorimeterHit> hits) {
+ //fill the readout buffers
+ for (CalorimeterHit hit : hits) {
+ Double eDep = eDepMap.get(hit.getCellID());
+ if (eDep == null) {
+ eDepMap.put(hit.getCellID(), hit.getRawEnergy());
+ } else {
+ eDepMap.put(hit.getCellID(), eDep + hit.getRawEnergy());
+ }
+ }
+ }
+
+ public void initReadout() {
+ //initialize buffers
+ eDepMap = new HashMap<Long, Double>();
+ }
}
\ No newline at end of file
hps-java/src/main/java/org/lcsim/hps/recon/ecal
diff -u -r1.6 -r1.7
--- HPSEcalTimeEvolutionReadoutDriver.java 2 Sep 2011 21:16:20 -0000 1.6
+++ HPSEcalTimeEvolutionReadoutDriver.java 2 Sep 2011 21:17:49 -0000 1.7
@@ -1,84 +1,84 @@
-package org.lcsim.hps.recon.ecal;
-
-import java.util.ArrayList;
-import java.util.HashMap;
-import java.util.List;
-import java.util.Map;
-
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.geometry.IDDecoder;
-import org.lcsim.hps.util.ClockSingleton;
-import org.lcsim.hps.util.RingBuffer;
-
-/**
- * Performs readout of ECal hits.
- * Simulates time evolution of preamp output pulse.
- *
- * @author Sho Uemura <[log in to unmask]>
- * @version $Id: HPSEcalTimeEvolutionReadoutDriver.java,v 1.6 2011/09/02 21:16:20 meeg Exp $
- */
-public class HPSEcalTimeEvolutionReadoutDriver extends HPSEcalReadoutDriver {
-
- //buffer for deposited energy
- Map<Long, RingBuffer> eDepMap = null;
- //length of ring buffer (in readout cycles)
- int bufferLength = 20;
- //shaper time constant in ns
- double t0 = 18.0;
-
- public HPSEcalTimeEvolutionReadoutDriver() {
- }
-
- public void setBufferLength(int bufferLength) {
- this.bufferLength = bufferLength;
- eDepMap = new HashMap<Long, RingBuffer>();
- }
-
- public List<HPSRawCalorimeterHit> readHits() {
- IDDecoder dec = ecal.getIDDecoder();
- List<HPSRawCalorimeterHit> hitList = new ArrayList<HPSRawCalorimeterHit>();
- for (Long cellID : eDepMap.keySet()) {
- RingBuffer eDepBuffer = eDepMap.get(cellID);
- if (eDepBuffer.currentValue() > threshold) {
- dec.setID(cellID);
- hitList.add(new HPSRawCalorimeterHit(eDepBuffer.currentValue(), dec.getPosition(), 0.0, cellID, hitType));
- }
- eDepBuffer.step();
- }
- return hitList;
- }
-
- public void putHits(List<CalorimeterHit> hits) {
- //fill the readout buffers
- for (CalorimeterHit hit : hits) {
- IDDecoder dec = ecal.getIDDecoder();
- dec.setID(hit.getCellID());
-// int ix = dec.getValue("ix");
-// int iy = dec.getValue("iy");
-// //temporary hack to disable crystals and flip X coordinate
-// int side = dec.getValue("side");
-// if (iy == 1 && ix * side >= -10 && ix * side <= -2)
-// continue;
-
- RingBuffer eDepBuffer = eDepMap.get(hit.getCellID());
- if (eDepBuffer == null) {
- eDepBuffer = new RingBuffer(bufferLength);
- eDepMap.put(hit.getCellID(), eDepBuffer);
- }
- for (int i = 0; i < bufferLength; i++) {
- eDepBuffer.addToCell(i, hit.getRawEnergy() * pulseAmplitude(i * ClockSingleton.getDt() * readoutCycle - hit.getTime()));
- }
- }
- }
-
- public void initReadout() {
- //initialize buffers
- eDepMap = new HashMap<Long, RingBuffer>();
- }
-
- private double pulseAmplitude(double time) {
- if (time < 0.0)
- return 0.0;
- return (time / t0) * Math.exp(1.0 - time / t0);
- }
+package org.lcsim.hps.recon.ecal;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.hps.util.ClockSingleton;
+import org.lcsim.hps.util.RingBuffer;
+
+/**
+ * Performs readout of ECal hits.
+ * Simulates time evolution of preamp output pulse.
+ *
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: HPSEcalTimeEvolutionReadoutDriver.java,v 1.7 2011/09/02 21:17:49 meeg Exp $
+ */
+public class HPSEcalTimeEvolutionReadoutDriver extends HPSEcalReadoutDriver {
+
+ //buffer for deposited energy
+ Map<Long, RingBuffer> eDepMap = null;
+ //length of ring buffer (in readout cycles)
+ int bufferLength = 20;
+ //shaper time constant in ns
+ double t0 = 18.0;
+
+ public HPSEcalTimeEvolutionReadoutDriver() {
+ }
+
+ public void setBufferLength(int bufferLength) {
+ this.bufferLength = bufferLength;
+ eDepMap = new HashMap<Long, RingBuffer>();
+ }
+
+ public List<HPSRawCalorimeterHit> readHits() {
+ IDDecoder dec = ecal.getIDDecoder();
+ List<HPSRawCalorimeterHit> hitList = new ArrayList<HPSRawCalorimeterHit>();
+ for (Long cellID : eDepMap.keySet()) {
+ RingBuffer eDepBuffer = eDepMap.get(cellID);
+ if (eDepBuffer.currentValue() > threshold) {
+ dec.setID(cellID);
+ hitList.add(new HPSRawCalorimeterHit(eDepBuffer.currentValue(), dec.getPosition(), 0.0, cellID, hitType));
+ }
+ eDepBuffer.step();
+ }
+ return hitList;
+ }
+
+ public void putHits(List<CalorimeterHit> hits) {
+ //fill the readout buffers
+ for (CalorimeterHit hit : hits) {
+ IDDecoder dec = ecal.getIDDecoder();
+ dec.setID(hit.getCellID());
+// int ix = dec.getValue("ix");
+// int iy = dec.getValue("iy");
+// //temporary hack to disable crystals and flip X coordinate
+// int side = dec.getValue("side");
+// if (iy == 1 && ix * side >= -10 && ix * side <= -2)
+// continue;
+
+ RingBuffer eDepBuffer = eDepMap.get(hit.getCellID());
+ if (eDepBuffer == null) {
+ eDepBuffer = new RingBuffer(bufferLength);
+ eDepMap.put(hit.getCellID(), eDepBuffer);
+ }
+ for (int i = 0; i < bufferLength; i++) {
+ eDepBuffer.addToCell(i, hit.getRawEnergy() * pulseAmplitude(i * ClockSingleton.getDt() * readoutCycle - hit.getTime()));
+ }
+ }
+ }
+
+ public void initReadout() {
+ //initialize buffers
+ eDepMap = new HashMap<Long, RingBuffer>();
+ }
+
+ private double pulseAmplitude(double time) {
+ if (time < 0.0)
+ return 0.0;
+ return (time / t0) * Math.exp(1.0 - time / t0);
+ }
}
\ No newline at end of file
hps-java/src/main/java/org/lcsim/hps/recon/ecal
diff -u -r1.2 -r1.3
--- HPSRawCalorimeterHit.java 2 Sep 2011 21:16:19 -0000 1.2
+++ HPSRawCalorimeterHit.java 2 Sep 2011 21:17:49 -0000 1.3
@@ -1,49 +1,49 @@
-package org.lcsim.hps.recon.ecal;
-
-import org.lcsim.event.base.BaseCalorimeterHit;
-import org.lcsim.geometry.IDDecoder;
-
-/**
- * An implementation of CalorimeterHit, with a constructor that sets rawEnergy
- * for use in ECalReadout
- *
- * @author Sho Uemura
- * @version $Id: HPSRawCalorimeterHit.java,v 1.2 2011/09/02 21:16:19 meeg Exp $
- */
-public class HPSRawCalorimeterHit extends BaseCalorimeterHit {
-
- /**
- * Fully qualified constructor that sets rawEnergy
- *
- * @param energy Raw energy for this cell
- * @param position Global Cartesian coordinate for this cell
- * @param time Time of energy deposition
- * @param id Cell ID
- * @param type Type
- */
- public HPSRawCalorimeterHit(double energy, double[] position, double time, long id, int type) {
- this.rawEnergy = energy;
- this.position = position;
- this.time = time;
- this.id = id;
- this.type = type;
- }
-
- public int getQuadrant() {
- IDDecoder dec = getIDDecoder();
- dec.setID(id);
- if (dec.getValue("ix")>0) {
- if (dec.getValue("iy")>0) {
- return 1;
- } else {
- return 4;
- }
- } else {
- if (dec.getValue("iy")>0) {
- return 2;
- } else {
- return 3;
- }
- }
- }
-}
+package org.lcsim.hps.recon.ecal;
+
+import org.lcsim.event.base.BaseCalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+
+/**
+ * An implementation of CalorimeterHit, with a constructor that sets rawEnergy
+ * for use in ECalReadout
+ *
+ * @author Sho Uemura
+ * @version $Id: HPSRawCalorimeterHit.java,v 1.3 2011/09/02 21:17:49 meeg Exp $
+ */
+public class HPSRawCalorimeterHit extends BaseCalorimeterHit {
+
+ /**
+ * Fully qualified constructor that sets rawEnergy
+ *
+ * @param energy Raw energy for this cell
+ * @param position Global Cartesian coordinate for this cell
+ * @param time Time of energy deposition
+ * @param id Cell ID
+ * @param type Type
+ */
+ public HPSRawCalorimeterHit(double energy, double[] position, double time, long id, int type) {
+ this.rawEnergy = energy;
+ this.position = position;
+ this.time = time;
+ this.id = id;
+ this.type = type;
+ }
+
+ public int getQuadrant() {
+ IDDecoder dec = getIDDecoder();
+ dec.setID(id);
+ if (dec.getValue("ix")>0) {
+ if (dec.getValue("iy")>0) {
+ return 1;
+ } else {
+ return 4;
+ }
+ } else {
+ if (dec.getValue("iy")>0) {
+ return 2;
+ } else {
+ return 3;
+ }
+ }
+ }
+}
CVSspam 0.2.8