java/trunk/analysis/src/main/java/org/hps/analysis/ecal
--- java/trunk/analysis/src/main/java/org/hps/analysis/ecal/HPSEcalTriggerPlotsDriver.java 2014-04-08 01:05:06 UTC (rev 466)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/ecal/HPSEcalTriggerPlotsDriver.java 2014-04-08 06:41:10 UTC (rev 467)
@@ -2,17 +2,14 @@
import hep.aida.IHistogram2D;
+import java.util.ArrayList;
import java.util.List;
import org.hps.readout.ecal.TriggerDriver;
import org.hps.recon.ecal.ECalUtils;
-import org.hps.recon.ecal.HPSCalorimeterHit;
import org.hps.recon.ecal.HPSEcalCluster;
-import org.lcsim.detector.identifier.IIdentifier;
-import org.lcsim.detector.identifier.IIdentifierHelper;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.EventHeader;
-import org.lcsim.geometry.IDDecoder;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -24,22 +21,26 @@
* Exp $
*/
public class HPSEcalTriggerPlotsDriver extends Driver {
-
+ // LCSim collection names.
String ecalCollectionName = "EcalHits";
String clusterCollectionName = "EcalClusters";
+
+ // Histogram factory.
AIDA aida = AIDA.defaultInstance();
- IHistogram2D hitXYPlot;
- IHistogram2D hitXYPlot50;
- IHistogram2D hitXYPlot100;
- IHistogram2D hitXYPlot200;
- IHistogram2D hitXYPlot500;
- IHistogram2D hitXYPlot1000;
+
+ // Energy cuts for hit histograms (in MeV).
+ double energyCut[] = { 0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
+
+ // Energy cut histograms.
+ IHistogram2D hitXYPlot[] = new IHistogram2D[energyCut.length];
+
+ // Additional histograms.
IHistogram2D crystalDeadTime;
IHistogram2D clusterHitXYPlot;
IHistogram2D seedHitXYPlot;
IHistogram2D triggerClusterHitXYPlot;
IHistogram2D triggerSeedHitXYPlot;
- IDDecoder dec = null;
+
private int coincidenceWindow = 2;
private double tp = 14.0;
private double threshold = 50 * 10 * 0.15 * ECalUtils.MeV;
@@ -51,133 +52,104 @@
public void setClusterCollectionName(String clusterCollectionName) {
this.clusterCollectionName = clusterCollectionName;
}
-
- @Override
+
public void startOfData() {
- hitXYPlot = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Hits",
- 46, -23, 23, 11, -5.5, 5.5);
- hitXYPlot50 = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Hits above 50 MeV",
- 46, -23, 23, 11, -5.5, 5.5);
- hitXYPlot100 = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Hits above 100 MeV",
- 46, -23, 23, 11, -5.5, 5.5);
- hitXYPlot200 = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Hits above 200 MeV",
- 46, -23, 23, 11, -5.5, 5.5);
- hitXYPlot500 = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Hits above 500 MeV",
- 46, -23, 23, 11, -5.5, 5.5);
- hitXYPlot1000 = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Hits above 1000 MeV",
- 46, -23, 23, 11, -5.5, 5.5);
- crystalDeadTime = aida.histogram2D(
- "Trigger plots: " + ecalCollectionName + " : Crystal dead time",
- 46, -23, 23, 11, -5.5, 5.5);
- clusterHitXYPlot = aida.histogram2D(
- "Trigger plots: " + clusterCollectionName + " : Crystals in clusters",
- 47, -23.5, 23.5, 11, -5.5, 5.5);
- seedHitXYPlot = aida.histogram2D(
- "Trigger plots: " + clusterCollectionName + " : Seed hits",
- 47, -23.5, 23.5, 11, -5.5, 5.5);
- triggerClusterHitXYPlot = aida.histogram2D(
- "Trigger plots: " + clusterCollectionName + " : Crystals in clusters, with trigger",
- 47, -23.5, 23.5, 11, -5.5, 5.5);
- triggerSeedHitXYPlot = aida.histogram2D(
- "Trigger plots: " + clusterCollectionName + " : Seed hits, with trigger",
- 47, -23.5, 23.5, 11, -5.5, 5.5);
+ // Initialize a hit histogram for each declared energy.
+ for(int e = 0; e < energyCut.length; e++) {
+ hitXYPlot[e] = aida.histogram2D("Trigger Plots: " + ecalCollectionName +
+ " : Hits above " + energyCut[e] + " MeV", 46, -23, 23, 11, -5.5, 5.5);
+ }
+ // Initialize the remaining plots.
+ crystalDeadTime = aida.histogram2D("Trigger Plots: " + ecalCollectionName +
+ " : Crystal dead time", 46, -23, 23, 11, -5.5, 5.5);
+ clusterHitXYPlot = aida.histogram2D("Trigger Plots: " + clusterCollectionName +
+ " : Crystals in clusters", 47, -23.5, 23.5, 11, -5.5, 5.5);
+ seedHitXYPlot = aida.histogram2D("Trigger Plots: " + clusterCollectionName +
+ " : Seed hits", 47, -23.5, 23.5, 11, -5.5, 5.5);
+ triggerClusterHitXYPlot = aida.histogram2D("Trigger Plots: " + clusterCollectionName +
+ " : Crystals in clusters, with trigger", 47, -23.5, 23.5, 11, -5.5, 5.5);
+ triggerSeedHitXYPlot = aida.histogram2D("Trigger Plots: " + clusterCollectionName +
+ " : Seed hits, with trigger", 47, -23.5, 23.5, 11, -5.5, 5.5);
}
-
- @Override
+
public void process(EventHeader event) {
- List<HPSEcalCluster> clusters = event.get(HPSEcalCluster.class, clusterCollectionName);
- if (clusters == null) {
- throw new RuntimeException("Missing cluster collection!");
- }
-
- List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
- if (hits == null) {
- throw new RuntimeException("Missing hit collection!");
- }
-
- // Get ID helper.
- IIdentifierHelper helper =
- event.getMetaData(hits).getIDDecoder().getSubdetector().getDetectorElement().getIdentifierHelper();
-
+ // If the current event has the indicated hit collection,
+ // use it as the hit list.
+ List<CalorimeterHit> hits;
+ if(event.hasCollection(CalorimeterHit.class, ecalCollectionName)) {
+ hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ }
+ // If it does not, then use an empty list to avoid crashing.
+ else { hits = new ArrayList<CalorimeterHit>(0); }
+
+ // If the current event has the indicated cluster collection,
+ // use it as the cluster list.
+ List<HPSEcalCluster> clusters;
+ if(event.hasCollection(HPSEcalCluster.class, clusterCollectionName)) {
+ clusters = event.get(HPSEcalCluster.class, clusterCollectionName);
+ }
+ // If it does not, then use an empty list to avoid crashing.
+ else { clusters = new ArrayList<HPSEcalCluster>(0); }
+
+ // Populate hit plots.
for (CalorimeterHit hit : hits) {
+ // Get the hit crystal position.
int ix = hit.getIdentifierFieldValue("ix");
int iy = hit.getIdentifierFieldValue("iy");
- hitXYPlot.fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
- if (hit.getRawEnergy() > 50.0 * ECalUtils.MeV) {
- hitXYPlot50.fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
- if (hit.getRawEnergy() > 100.0 * ECalUtils.MeV) {
- hitXYPlot100.fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
- if (hit.getRawEnergy() > 200.0 * ECalUtils.MeV) {
- hitXYPlot200.fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
- if (hit.getRawEnergy() > 500.0 * ECalUtils.MeV) {
- hitXYPlot500.fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
- if (hit.getRawEnergy() > 1000.0 * ECalUtils.MeV) {
- hitXYPlot1000.fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
- }
- }
- }
- }
+ double energy = hit.getRawEnergy();
+
+ // Loop through the energy plots and fill them if the hit
+ // is over the current energy threshold/
+ for(int e = 0; e < energyCut.length; e++) {
+ if(energy > energyCut[e] * ECalUtils.MeV) {
+ hitXYPlot[e].fill(ix - 0.5 * Math.signum(ix), iy, 1.0 / coincidenceWindow);
+ }
}
+
+ // Generate the dead time plot.
double deadTime = 0;
for (int time = 0; time < 500; time++) {
- if (hit.getRawEnergy() * pulseAmplitude(time) > threshold) {
- deadTime += 1e-6; //units of milliseconds
- } else if (time > 2 * tp || deadTime != 0) {
- break;
- }
+ if (hit.getRawEnergy() * pulseAmplitude(time) > threshold) { deadTime += 1e-6; } // units of milliseconds
+ else if (time > 2 * tp || deadTime != 0) { break; }
}
crystalDeadTime.fill(ix - 0.5 * Math.signum(ix), iy, deadTime / coincidenceWindow);
}
-
- for (HPSEcalCluster clus : clusters) {
- HPSCalorimeterHit seedHit = (HPSCalorimeterHit) clus.getSeedHit();
- IIdentifier id = seedHit.getIdentifier();
- int ix = helper.unpack(id).getValue(helper.getFieldIndex("ix"));
- int iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
-// dec = seedHit.getIDDecoder();
-// dec.setID(seedHit.getCellID());
+
+ // Check the trigger bit.
+ boolean trigger = TriggerDriver.triggerBit();
+
+ // Populate cluster based plots.
+ for (HPSEcalCluster cluster : clusters) {
+ // Get the cluster's seed hit position.
+ CalorimeterHit seed = cluster.getSeedHit();
+ int ix = seed.getIdentifierFieldValue("ix");
+ int iy = seed.getIdentifierFieldValue("iy");
+
+ // Fill the seed hit plot.
seedHitXYPlot.fill(ix, iy);
- for (CalorimeterHit hit : clus.getCalorimeterHits()) {
- id = hit.getIdentifier();
- ix = helper.unpack(id).getValue(helper.getFieldIndex("ix"));
- iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
-// dec = hit.getIDDecoder();
-// dec.setID(hit.getCellID());
+
+ // If the trigger bit is set, add the seed to the trigger
+ // plot histogram as well.
+ if(trigger) { triggerSeedHitXYPlot.fill(ix, iy); }
+
+ // Populate the component hit histogram.
+ for (CalorimeterHit hit : cluster.getCalorimeterHits()) {
+ // Get the component hit location.
+ ix = hit.getIdentifierFieldValue("ix");
+ iy = hit.getIdentifierFieldValue("iy");
+
+ // Add it to the plot.
clusterHitXYPlot.fill(ix, iy);
+
+ // If the trigger bit is set, add the component to
+ // the trigger cluster histogram too.
+ if(trigger) { triggerClusterHitXYPlot.fill(ix, iy); }
}
}
-
- if (TriggerDriver.triggerBit()) {
- for (HPSEcalCluster clus : clusters) {
- HPSCalorimeterHit seedHit = (HPSCalorimeterHit) clus.getSeedHit();
- IIdentifier id = seedHit.getIdentifier();
- int ix = helper.unpack(id).getValue(helper.getFieldIndex("ix"));
- int iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
-// dec = seedHit.getIDDecoder();
-// dec.setID(seedHit.getCellID());
- triggerSeedHitXYPlot.fill(ix, iy);
- for (CalorimeterHit hit : clus.getCalorimeterHits()) {
- id = hit.getIdentifier();
- ix = helper.unpack(id).getValue(helper.getFieldIndex("ix"));
- iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
-// dec = hit.getIDDecoder();
-// dec.setID(hit.getCellID());
- triggerClusterHitXYPlot.fill(ix, iy);
- }
- }
- }
}
-
+
private double pulseAmplitude(double time) {
- if (time <= 0.0) {
- return 0.0;
- }
+ if (time <= 0.0) { return 0.0; }
return (time / tp) * Math.exp(1.0 - time / tp);
}
}
