hps-java/src/main/java/org/lcsim/hps/analysis/ecal
diff -u -r1.6 -r1.7
--- HPSEcalPlotsDriver.java 27 Jul 2011 23:33:44 -0000 1.6
+++ HPSEcalPlotsDriver.java 29 Aug 2011 23:07:02 -0000 1.7
@@ -23,6 +23,9 @@
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;
@@ -31,32 +34,29 @@
* Diagnostic plots for HPS ECal.
*
* @author Jeremy McCormick <[log in to unmask]>
- * @version $Id: HPSEcalPlotsDriver.java,v 1.6 2011/07/27 23:33:44 jeremy Exp $
+ * @version $Id: HPSEcalPlotsDriver.java,v 1.7 2011/08/29 23:07:02 meeg Exp $
*/
-public class HPSEcalPlotsDriver extends Driver
-{
+public class HPSEcalPlotsDriver extends Driver {
+
String ecalCollectionName = "EcalHits";
String ecalClusterCollectionName = "EcalClusters";
-
AIDA aida = AIDA.defaultInstance();
-
// MCParticle plots.
ICloud1D primaryEPlot;
- ICloud1D fsCountPlot;
- ICloud1D fsEPlot;
+ ICloud1D fsCountPlot;
+ ICloud1D fsEPlot;
ICloud1D fsGammaEPlot;
ICloud1D fsElectronEPlot;
ICloud1D fsPositronEPlot;
-
// CalHit plots.
IHistogram1D hitEPlot;
- ICloud1D ecalEPlot;
- ICloud1D eventEPlot;
+ ICloud1D ecalEPlot;
+ ICloud1D eventEPlot;
ICloud2D hitXZPlot;
ICloud2D hitYZPlot;
ICloud1D hitUnder100MeVPlot;
IHistogram1D hitOver100MeVPlot;
- ICloud1D maxHitEPlot;
+ ICloud1D maxHitEPlot;
IHistogram1D maxTimePlot;
ICloud1D timePlot;
ICloud1D hitCountPlot;
@@ -65,11 +65,10 @@
ICloud1D crystalYPlot;
//ICloud2D cellXYPlot;
IHistogram2D crystalXYPlot;
-
// Cluster plots.
IHistogram1D nclusPlot;
ICloud1D clusEPlot;
- ICloud1D clusTotEPlot;
+ ICloud1D clusTotEPlot;
ICloud1D leadClusEPlot;
ICloud1D leadClus2EPlot;
//ICloud1D clusResTop3Plot;
@@ -78,416 +77,399 @@
ICloud1D clusSeedDistPlot;
ICloud2D leadClusAndPrimaryPlot;
IHistogram1D clusNHits;
-
- class MCParticleEComparator implements Comparator<MCParticle>
- {
- public int compare(MCParticle p1, MCParticle p2)
- {
+ 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)
- {
+ if (e1 < e2) {
return -1;
- }
- else if (e1 == e2)
- {
+ } else if (e1 == e2) {
return 0;
- }
- else
- {
+ } else {
return 1;
}
- }
+ }
}
-
- class ClusterEComparator implements Comparator<Cluster>
- {
- public int compare(Cluster o1, Cluster o2)
- {
- if (o1.getEnergy() < o2.getEnergy())
- {
+
+ 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())
- {
+ } else if (o1.getEnergy() > o2.getEnergy()) {
return 1;
- }
- else
- {
+ } else {
return 0;
}
- }
+ }
}
-
- public void setEcalCollectionName(String ecalCollectionName)
- {
+
+ public void setEcalCollectionName(String ecalCollectionName) {
this.ecalCollectionName = ecalCollectionName;
}
-
- public void setEcalClusterCollectionName(String ecalClusterCollectionName)
- {
+
+ public void setEcalClusterCollectionName(String ecalClusterCollectionName) {
this.ecalClusterCollectionName = ecalClusterCollectionName;
- }
-
- public void startOfData()
- {
+ }
+
+ 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 = 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",
+ 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]");
-
+ 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("xAxisLabel", "Distance [mm]");
clusSeedDistPlot.annotation().addItem("yAxisScale", "log");
-
- clusSeedEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Hit Distance from Beam Axis");
- clusSeedEPlot.annotation().addItem("xAxisLabel", "Distance [mm]");
-
+
+ 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)
- {
- // 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();
-
+
+ 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)
- {
+
+ 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))
- {
+ if (ParticleTypeClassifier.isElectron(fsPdg)) {
+ fsElectronEPlot.fill(fsE);
+ } else if (ParticleTypeClassifier.isPositron(fsPdg)) {
fsPositronEPlot.fill(fsE);
- }
- else if (ParticleTypeClassifier.isPhoton(fsPdg))
- {
+ } 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;
+ }
+
+ 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) {
+ hitXYPlot = aida.histogram2D(
+ ecalCollectionName + " : hit E, event " + String.format("%05d", ClockSingleton.getClock()),
+ 47, -23.5, 23.5, 11, -5.5, 5.5);
+ }
+
// Loop over ECal hits.
- for (CalorimeterHit hit : hits)
- {
+ for (CalorimeterHit hit : hits) {
// get raw E
double eraw = hit.getRawEnergy();
-
- if (eraw >= 0.1)
- {
+
+ if (eraw >= 0.1) {
nhitsOver100MeV++;
- }
-
- if (eraw < 0.1)
- {
+ }
+
+ if (eraw < 0.1) {
nhits100MeV++;
}
-
+
// time
timePlot.fill(hit.getTime());
-
+
// YZ
- hitYZPlot.fill(hit.getPosition()[1], hit.getPosition()[2]);
-
+ 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)
+ 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"));
+ int iy = helper.unpack(id).getValue(helper.getFieldIndex("iy"));
+ int side = helper.unpack(id).getValue(helper.getFieldIndex("side"));
crystalXPlot.fill(ix);
crystalYPlot.fill(iy);
crystalXYPlot.fill(ix, iy);
+ if (trigger)
+ hitXYPlot.fill(ix, iy * side, 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());
-
- // 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<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;
- for (Cluster clus : clusters)
- {
- clusNHits.fill(clus.getCalorimeterHits().size());
-
- 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;
-
- 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());
- }
-
- // 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)
- {
+ private double getPrimaryE(List<MCParticle> particles) {
double electronE = 0;
- for (MCParticle particle : particles)
- {
- if (particle.getGeneratorStatus() == MCParticle.FINAL_STATE
- && Math.abs(particle.getPDGID()) == 11)
- {
+ 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)
- {
+
+ 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)
- {
+ for (MCParticle particle : particles) {
+ if (particle.getGeneratorStatus() == MCParticle.FINAL_STATE
+ && particle.getEnergy() > maxE) {
maxE = particle.getEnergy();
primary = particle;
}
}
- return primary;
+ return primary;
}
-
- private static List<MCParticle> makeGenFSParticleList(List<MCParticle> mcparticles)
- {
+
+ private static List<MCParticle> makeGenFSParticleList(List<MCParticle> mcparticles) {
List<MCParticle> fsParticles = new ArrayList<MCParticle>();
- for (MCParticle mcparticle : mcparticles)
- {
- if (mcparticle.getGeneratorStatus() == MCParticle.FINAL_STATE)
- {
+ 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);
+ }
}
