lcsim-contrib/src/main/java/org/lcsim/contrib/jeremym/hps
diff -N HPSEcalPlotsDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HPSEcalPlotsDriver.java 28 Apr 2011 00:34:39 -0000 1.1
@@ -0,0 +1,344 @@
+package org.lcsim.contrib.jeremym.hps;
+
+import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
+import hep.aida.IHistogram1D;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Set;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.units.SystemOfUnits;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * Diagnostic plots for HPS ECal.
+ *
+ * @author Jeremy McCormick <[log in to unmask]>
+ * @version $Id: HPSEcalPlotsDriver.java,v 1.1 2011/04/28 00:34:39 jeremy Exp $
+ */
+// TODO Sort Ecal hits on energy and plot top 2,3,4,5
+// TODO Add plot of percentage of total FS particle E in highest hit and top 2 (?) hits
+// TODO ZY plot can use h2d with small bins in Y to avoid weird binning effect
+
+public class HPSEcalPlotsDriver extends Driver
+{
+ String ecalCollectionName = "EcalHits";
+ String ecalClusterCollectionName = "EcalClusters";
+
+ AIDA aida = AIDA.defaultInstance();
+
+ IHistogram1D hitEPlot;
+ ICloud1D ecalEPlot;
+ ICloud1D primaryEPlot;
+ ICloud1D eventEPlot;
+ ICloud1D residualPlot;
+ ICloud1D residualAllPlot;
+ ICloud1D residualMaxHitPlot;
+ ICloud1D maxHitEPlot;
+ ICloud2D hitYZPlot;
+ IHistogram1D maxTimePlot;
+ ICloud1D timePlot;
+ ICloud1D hitCountPlot;
+ ICloud1D idCountPlot;
+ ICloud1D residualTop2Plot;
+ ICloud1D residualTop3Plot;
+ ICloud1D lowEHitPlot;
+ ICloud2D hitXZPlot;
+ IHistogram1D acceptHitPlot;
+ ICloud1D fsPlot;
+ ICloud1D clusEPlot;
+ ICloud1D nclusPlot;
+
+ double hitEnergyMin = 0;
+
+ static 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;
+ }
+ }
+ }
+
+ public void startOfData()
+ {
+ timePlot = aida.cloud1D(ecalCollectionName + " : Hit Time [ns]");
+ timePlot.annotation().addItem("yAxisScale", "log");
+
+ maxTimePlot = aida.histogram1D(ecalCollectionName + " : Max Time [ns]", 200, 0, 1000);
+ maxTimePlot.annotation().addItem("yAxisScale", "log");
+
+ maxHitEPlot = aida.cloud1D(ecalCollectionName + " : Max Hit E in Event [GeV]");
+
+ hitEPlot = aida.histogram1D(ecalCollectionName + " : Hit Energy [MeV]", 200, 0, 3500);
+ hitEPlot.annotation().addItem("yAxisScale", "log");
+
+ hitCountPlot = aida.cloud1D(ecalCollectionName + " : Hit Count");
+
+ idCountPlot = aida.cloud1D(ecalCollectionName + " : Uniq Hit IDs");
+
+ ecalEPlot = aida.cloud1D(ecalCollectionName + " : E Event [GeV]");
+
+ primaryEPlot = aida.cloud1D("MCParticle: Highest Primary Energy in Event [GeV]");
+
+ eventEPlot = aida.cloud1D("MCParticle: Sum of FS Electron Energies in Event [GeV]");
+
+ residualPlot = aida.cloud1D(ecalCollectionName + " : E Residual with Highest E Particle [GeV]");
+
+ residualAllPlot = aida.cloud1D(ecalCollectionName + " : E Residual with All Final State Electrons [GeV]");
+
+ residualMaxHitPlot = aida.cloud1D(ecalCollectionName + " : Residual of Max Hit and Highest E Electron");
+
+ hitYZPlot = aida.cloud2D(ecalCollectionName + " : Y vs Z");
+
+ hitXZPlot = aida.cloud2D(ecalCollectionName + " : X vs Z");
+
+ residualTop2Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 2 FS Particle E");
+
+ residualTop3Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 3 FS Particle E");
+
+ lowEHitPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 100 MeV");
+
+ acceptHitPlot = aida.histogram1D(ecalCollectionName + " : Hits with E >= 100 MeV", 15, 0, 15);
+
+ fsPlot = aida.cloud1D("MCParticle: Number of Final State Particlesin Event");
+
+ clusEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster E");
+
+ nclusPlot = aida.cloud1D(ecalClusterCollectionName + " : Number of Clusters");
+ }
+
+ private static List<MCParticle> getFinalStateParticles(List<MCParticle> mcparticles)
+ {
+ List<MCParticle> fsParticles = new ArrayList<MCParticle>();
+ for (MCParticle mcparticle : mcparticles)
+ {
+ if (mcparticle.getGeneratorStatus() == MCParticle.FINAL_STATE)
+ {
+ fsParticles.add(mcparticle);
+ }
+ }
+ return fsParticles;
+ }
+
+ public void process(EventHeader event)
+ {
+ // sum and max vars
+ double esum = 0;
+ double emax = 0;
+ double tmax = 0;
+
+ // Loop over hits from ECal collection.
+ List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName);
+ int nhits = hits.size();
+
+ // MCParticles
+ List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
+
+ // Final State particles.
+ List<MCParticle> fsParticles = getFinalStateParticles(mcparticles);
+
+ //System.out.println("fsParticles="+fsParticles.size());
+ fsPlot.fill(fsParticles.size());
+
+ // Sort MCParticles on energy.
+ Collections.sort(fsParticles, new MCParticleEComparator());
+
+ // Energy of top two FS particles.
+ double e2 = fsParticles.get(0).getEnergy() + fsParticles.get(1).getEnergy();
+
+ // Energy of top three FS particles.
+ double e3 = e2 + fsParticles.get(2).getEnergy();
+
+ // Check unique IDs.
+ Set<Long> ids = new HashSet<Long>();
+
+ // n hits
+ hitCountPlot.fill(nhits);
+
+ int nlowe = 0;
+ int naccept = 0;
+
+ // Loop over ECal hits.
+
+ //System.out.println("nhits="+hits.size());
+
+ for (CalorimeterHit hit : hits)
+ {
+ // get raw E
+ double eraw = hit.getRawEnergy();
+
+ // Keep track of # hits under energy threshold.
+ if (eraw < hitEnergyMin)
+ {
+ ++nlowe;
+ //continue;
+ }
+ else
+ {
+ ++naccept;
+ }
+
+ // time
+ timePlot.fill(hit.getTime());
+
+ // YZ
+ hitYZPlot.fill(hit.getPosition()[1], hit.getPosition()[2]);
+
+ // XZ
+ hitXZPlot.fill(hit.getPosition()[0], hit.getPosition()[2]);
+
+ // hit E
+ hitEPlot.fill(eraw / SystemOfUnits.MeV);
+
+ // check E max
+ if (eraw > emax)
+ emax = eraw;
+
+ // check T max
+ if (hit.getTime() > tmax)
+ tmax = hit.getTime();
+
+ if (ids.contains(hit.getCellID()))
+ throw new RuntimeException("Duplicate cell ID: " + hit.getCellID());
+
+ ids.add(hit.getCellID());
+
+ // incr E sum
+ esum += eraw;
+
+ }
+
+ //System.out.println("nskipped="+nlowe);
+ //System.out.println("naccept="+naccept);
+
+ //System.out.println("tmax: " + tmax);
+
+ // Low E hits.
+ lowEHitPlot.fill(nlowe);
+
+ // Hits passed E cut.
+ acceptHitPlot.fill(naccept);
+
+ // total E in Cal
+ ecalEPlot.fill(esum);
+
+ // max hit E in event
+ maxHitEPlot.fill(emax);
+
+ // max hit time
+ maxTimePlot.fill(tmax);
+
+ // number of unique hit ids
+ idCountPlot.fill(ids.size());
+
+ // primary particle with most E
+ MCParticle primary = getPrimary(mcparticles);
+ double primaryE = primary.getEnergy();
+ primaryEPlot.fill(primaryE);
+
+ // residual using highest E primary
+ double singleResidual = esum - primaryE;
+ residualPlot.fill(singleResidual);
+
+ // event energy
+ double eventE = getPrimaryE(mcparticles);
+ eventEPlot.fill(eventE);
+
+ // event residual
+ double residual = esum - eventE;
+ residualAllPlot.fill(residual);
+
+ // hottest hit residual with highest E particle
+ residualMaxHitPlot.fill(emax - primaryE);
+
+ // Residual with two
+ double res2 = emax - e2;
+ residualTop2Plot.fill(res2);
+
+ double res3 = emax - e3;
+ residualTop3Plot.fill(res3);
+
+ List<Cluster> clusters = event.get(Cluster.class, ecalClusterCollectionName);
+ if (clusters == null)
+ throw new RuntimeException("Missing cluster collection!");
+
+ nclusPlot.fill(clusters.size());
+
+ for (Cluster clus : clusters)
+ {
+ clusEPlot.fill(clus.getEnergy());
+ }
+ }
+
+ 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;
+ }
+
+ public void setEcalCollectionName(String ecalCollectionName)
+ {
+ this.ecalCollectionName = ecalCollectionName;
+ }
+
+ public void setEcalClusterCollectionName(String ecalClusterCollectionName)
+ {
+ this.ecalClusterCollectionName = ecalClusterCollectionName;
+ }
+
+ /**
+ * Set minimum hit energy in GeV.
+ * @param hitEnergyMin
+ */
+ public void setHitEnergyMin(double hitEnergyMin)
+ {
+ //System.out.println("hitEnergyMin="+hitEnergyMin);
+ this.hitEnergyMin = hitEnergyMin;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/jeremym/hps
diff -N HPSECalPlotsDriver.java
--- HPSECalPlotsDriver.java 22 Apr 2011 01:10:07 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,189 +0,0 @@
-package org.lcsim.contrib.jeremym.hps;
-
-
-import hep.aida.ICloud1D;
-import hep.aida.ICloud2D;
-import hep.aida.IHistogram1D;
-
-import java.util.HashSet;
-import java.util.List;
-import java.util.Set;
-
-import org.lcsim.event.CalorimeterHit;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.units.SystemOfUnits;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-
-/**
- * @author Jeremy McCormick <[log in to unmask]>
- * @version $Id: HPSECalPlotsDriver.java,v 1.1 2011/04/22 01:10:07 jeremy Exp $
- */
-public class HPSECalPlotsDriver extends Driver
-{
- String ecalCollectionName = "EcalHits";
- AIDA aida = AIDA.defaultInstance();
-
- IHistogram1D hitEPlot;
- ICloud1D ecalEPlot;
- ICloud1D primaryEPlot;
- ICloud1D eventEPlot;
- ICloud1D residualPlot;
- ICloud1D residualAllPlot;
- ICloud1D residualMaxHitPlot;
- ICloud1D maxHitEPlot;
- ICloud2D hitYZPlot;
- IHistogram1D maxTimePlot;
- ICloud1D timePlot;
- ICloud1D hitCountPlot;
- ICloud1D idCountPlot;
-
- public void startOfData()
- {
- timePlot = aida.cloud1D(ecalCollectionName + " : Hit Time [ns]");
- timePlot.annotation().addItem("yAxisScale", "log");
-
- maxTimePlot = aida.histogram1D(ecalCollectionName + " : Max Time [ns]", 200, 0, 1000);
- maxTimePlot.annotation().addItem("yAxisScale", "log");
-
- maxHitEPlot = aida.cloud1D(ecalCollectionName + " : Max Hit E in Event [GeV]");
-
- hitEPlot = aida.histogram1D(ecalCollectionName + " : Hit Energy [MeV]", 200, 0, 3500);
- hitEPlot.annotation().addItem("yAxisScale", "log");
-
- hitCountPlot = aida.cloud1D(ecalCollectionName + " : Hit Count");
-
- idCountPlot = aida.cloud1D(ecalCollectionName + " : Uniq Hit IDS");
-
- ecalEPlot = aida.cloud1D(ecalCollectionName + " : E Event [GeV]");
-
- primaryEPlot = aida.cloud1D("Highest Primary Particle Energy in Event [GeV]");
-
- eventEPlot = aida.cloud1D("Sum of Final State Electron Energies in Event [GeV]");
-
- residualPlot = aida.cloud1D(ecalCollectionName + " : E Residual with Highest E Particle [GeV]");
-
- residualAllPlot = aida.cloud1D(ecalCollectionName + " : E Residual with All Final State Electrons [GeV]");
-
- residualMaxHitPlot = aida.cloud1D(ecalCollectionName + " : Residual of Max Hit and Electron E");
-
- hitYZPlot = aida.cloud2D(ecalCollectionName + " : Y vs Z");
- }
-
- 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();
-
- Set<Long> ids = new HashSet<Long>();
-
- // n hits
- hitCountPlot.fill(nhits);
-
- for (CalorimeterHit hit : hits)
- {
- // time
- timePlot.fill(hit.getTime());
-
- // YZ
- hitYZPlot.fill(hit.getPosition()[1], hit.getPosition()[2]);
-
- // get raw E
- double eraw = hit.getRawEnergy();
-
- // 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();
-
- ids.add(hit.getCellID());
-
- // incr E sum
- esum += eraw;
-
- }
-
- //System.out.println("tmax: " + tmax);
-
- // 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(event.get(MCParticle.class).get(0));
- double primaryE = primary.getEnergy();
- primaryEPlot.fill(primaryE);
-
- // residual using highest E primary
- double singleResidual = esum - primaryE;
- residualPlot.fill(singleResidual);
-
- // event energy
- double eventE = getPrimaryE(event.get(MCParticle.class).get(0));
- eventEPlot.fill(eventE);
-
- // event residual
- double residual = esum - eventE;
- residualAllPlot.fill(residual);
-
- // hottest hit residual with highest E particle
- residualMaxHitPlot.fill(emax - primaryE);
-
- }
-
- 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;
- }
-
- public void setEcalCollectionName(String ecalCollectionName)
- {
- this.ecalCollectionName = ecalCollectionName;
- }
-}
