lcsim-contrib/src/main/java/org/lcsim/contrib/jeremym/hps
diff -u -r1.2 -r1.3
--- HPSEcalPlotsDriver.java 28 Apr 2011 01:56:37 -0000 1.2
+++ HPSEcalPlotsDriver.java 29 Apr 2011 00:07:45 -0000 1.3
@@ -3,18 +3,23 @@
import hep.aida.ICloud1D;
import hep.aida.ICloud2D;
import hep.aida.IHistogram1D;
+import hep.physics.vec.Hep3Vector;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
+import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
+import java.util.Map;
import java.util.Set;
+import java.util.Map.Entry;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.event.Cluster;
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
+import org.lcsim.geometry.IDDecoder;
import org.lcsim.units.SystemOfUnits;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -23,11 +28,11 @@
* Diagnostic plots for HPS ECal.
*
* @author Jeremy McCormick <[log in to unmask]>
- * @version $Id: HPSEcalPlotsDriver.java,v 1.2 2011/04/28 01:56:37 jeremy Exp $
+ * @version $Id: HPSEcalPlotsDriver.java,v 1.3 2011/04/29 00:07:45 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
+// TODO ZY of cluster positions (algo from Tim to find pos)
+// TODO ZY of cluster seed
public class HPSEcalPlotsDriver extends Driver
{
@@ -35,33 +40,47 @@
String ecalClusterCollectionName = "EcalClusters";
AIDA aida = AIDA.defaultInstance();
+
+ // MCParticle plots.
+ ICloud1D primaryEPlot;
+ ICloud1D fsPlot;
+ // CalHit plots.
IHistogram1D hitEPlot;
ICloud1D ecalEPlot;
- ICloud1D primaryEPlot;
- ICloud1D eventEPlot;
- ICloud1D residualPlot;
- ICloud1D residualAllPlot;
- ICloud1D residualMaxHitPlot;
- ICloud1D maxHitEPlot;
+ ICloud1D eventEPlot;
+ ICloud2D hitXZPlot;
ICloud2D hitYZPlot;
+ ICloud1D hitUnder100MeVPlot;
+ IHistogram1D hitOver100MeVPlot;
+ ICloud1D maxHitEPlot;
IHistogram1D maxTimePlot;
ICloud1D timePlot;
ICloud1D hitCountPlot;
ICloud1D idCountPlot;
- ICloud1D residualTop2Plot;
- ICloud1D residualTop3Plot;
- ICloud2D hitXZPlot;
- IHistogram1D acceptHitPlot;
- ICloud1D fsPlot;
+
+ // Cluster plots.
+ IHistogram1D nclusPlot;
ICloud1D clusEPlot;
- ICloud1D nclusPlot;
- ICloud1D hit50MeVPlot;
- ICloud1D hit30MeVPlot;
- ICloud1D hit100MeVPlot;
- ICloud1D hitOver100MeVPlot;
+ ICloud1D clusTotEPlot;
+ ICloud1D leadClusEPlot;
+ ICloud1D leadClus2EPlot;
+ ICloud1D clusResTop3Plot;
+ IHistogram1D clusHitPlot;
+ ICloud1D clusSeedEPlot;
+ ICloud1D clusSeedDistPlot;
+ ICloud2D leadClusAndPrimaryPlot;
+
+ //ICloud1D hit50MeVPlot;
+ //ICloud1D hit30MeVPlot;
- static class MCParticleEComparator implements Comparator<MCParticle>
+ //ICloud1D residualAllPlot;
+ //ICloud1D residualPlot;
+ //ICloud1D residualMaxHitPlot;
+ //ICloud1D residualTop2Plot;
+ //ICloud1D residualTop3Plot;
+
+ class MCParticleEComparator implements Comparator<MCParticle>
{
public int compare(MCParticle p1, MCParticle p2)
{
@@ -82,6 +101,25 @@
}
}
+ class ClusterEComparator implements Comparator<Cluster>
+ {
+ public int compare(Cluster o1, Cluster o2)
+ {
+ if (o1.getEnergy() < o2.getEnergy())
+ {
+ return -1;
+ }
+ else if (o1.getEnergy() > o2.getEnergy())
+ {
+ return 1;
+ }
+ else
+ {
+ return 0;
+ }
+ }
+ }
+
public void startOfData()
{
timePlot = aida.cloud1D(ecalCollectionName + " : Hit Time");
@@ -114,14 +152,14 @@
eventEPlot = aida.cloud1D("MCParticle: Total Gen FS Electron E in Event");
eventEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
- residualPlot = aida.cloud1D(ecalCollectionName + " : Residual with Highest E Particle [GeV]");
- residualPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+ //residualPlot = aida.cloud1D(ecalCollectionName + " : Residual with Highest E Particle [GeV]");
+ //residualPlot.annotation().addItem("xAxisLabel", "E [GeV]");
- residualAllPlot = aida.cloud1D(ecalCollectionName + " : Resdual with All Final State Electrons [GeV]");
- residualAllPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+ //residualAllPlot = aida.cloud1D(ecalCollectionName + " : Resdual with All Final State Electrons [GeV]");
+ //residualAllPlot.annotation().addItem("xAxisLabel", "E [GeV]");
- residualMaxHitPlot = aida.cloud1D(ecalCollectionName + " : Residual of Max Hit and Highest E Electron");
- residualMaxHitPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+ //residualMaxHitPlot = aida.cloud1D(ecalCollectionName + " : Residual of Max Hit and Highest E Electron");
+ //residualMaxHitPlot.annotation().addItem("xAxisLabel", "E [GeV]");
hitYZPlot = aida.cloud2D(ecalCollectionName + " : Y vs Z");
hitYZPlot.annotation().addItem("xAxisLabel", "Y [mm]");
@@ -131,14 +169,11 @@
hitXZPlot.annotation().addItem("xAxisLabel", "X [mm]");
hitXZPlot.annotation().addItem("yAxisLabel", "Z [mm]");
- residualTop2Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 2 FS Particle E");
- residualTop2Plot.annotation().addItem("xAxisLabel", "E [GeV]");
+ //residualTop2Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 2 FS Particle E");
+ //residualTop2Plot.annotation().addItem("xAxisLabel", "E [GeV]");
- residualTop3Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 3 FS Particle E");
- residualTop3Plot.annotation().addItem("xAxisLabel", "E [GeV]");
-
- acceptHitPlot = aida.histogram1D(ecalCollectionName + " : Hits with E >= 100 MeV", 15, 0, 15);
- acceptHitPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+ //residualTop3Plot = aida.cloud1D(ecalCollectionName + " : Residual with Top 3 FS Particle E");
+ //residualTop3Plot.annotation().addItem("xAxisLabel", "E [GeV]");
fsPlot = aida.cloud1D("MCParticle: Number of Final State Particles");
fsPlot.annotation().addItem("xAxisLabel", "Number of FS Particles");
@@ -146,35 +181,46 @@
clusEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster E");
clusEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
- nclusPlot = aida.cloud1D(ecalClusterCollectionName + " : Number of Clusters");
+ nclusPlot = aida.histogram1D(ecalClusterCollectionName + " : Number of Clusters", 20, 0, 20);
nclusPlot.annotation().addItem("xAxisLabel", "Number of Clusters");
- hit50MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 50 MeV");
- hit50MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+ //hit50MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 50 MeV");
+ //hit50MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
- hit30MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 30 MeV");
- hit30MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+ //hit30MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 30 MeV");
+ //hit30MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
- hit100MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 100 MeV");
- hit100MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+ hitUnder100MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E < 100 MeV");
+ hitUnder100MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
- hitOver100MeVPlot = aida.cloud1D(ecalCollectionName + " : Hits with E >= 100 MeV");
- hit100MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
- }
-
- 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;
+ hitOver100MeVPlot = aida.histogram1D(ecalCollectionName + " : Hits with E >= 100 MeV", 20, 0, 20);
+ hitUnder100MeVPlot.annotation().addItem("xAxisLabel", "Number of Hits");
+
+ leadClusEPlot = aida.cloud1D(ecalClusterCollectionName + " : Leading Cluster E");
+ leadClusEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ leadClus2EPlot = aida.cloud1D(ecalClusterCollectionName + " : Second Leading Cluster E");
+ leadClus2EPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusTotEPlot = aida.cloud1D(ecalClusterCollectionName + " : Total Clus E in Event");
+ clusTotEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusResTop3Plot = aida.cloud1D(ecalClusterCollectionName + " : Total Clus E Residual with Top 3 Particles");
+ clusResTop3Plot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ clusHitPlot = aida.histogram1D(ecalClusterCollectionName + " : Number of Clusters per Hit", 5, 1, 6);
+ clusHitPlot.annotation().addItem("xAxisLabel", "Number of Clusters");
+
+ clusSeedDistPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Hit Distance from Beam Axis");
+ clusSeedDistPlot.annotation().addItem("xAxisLabel", "Number of Clusters");
+ clusSeedDistPlot.annotation().addItem("yAxisScale", "log");
+
+ clusSeedEPlot = aida.cloud1D(ecalClusterCollectionName + " : Cluster Seed Hit Distance from Beam Axis");
+ clusSeedEPlot.annotation().addItem("xAxisLabel", "Distance [mm]");
+
+ leadClusAndPrimaryPlot = aida.cloud2D(ecalClusterCollectionName + " : Lead Cluster E vs Highest Primary Particle E");
}
-
+
public void process(EventHeader event)
{
// sum and max vars
@@ -210,10 +256,10 @@
// n hits
hitCountPlot.fill(nhits);
- int naccept = 0;
+ //int naccept = 0;
- int nhits50MeV = 0;
- int nhits30MeV = 0;
+ //int nhits50MeV = 0;
+ //int nhits30MeV = 0;
int nhits100MeV = 0;
int nhitsOver100MeV = 0;
@@ -227,18 +273,21 @@
{
nhitsOver100MeV++;
}
- else if (eraw < 0.1)
+
+ if (eraw < 0.1)
{
nhits100MeV++;
}
- else if (eraw < 0.05)
- {
- nhits50MeV++;
- }
- else if (eraw < 0.03)
- {
- nhits30MeV++;
- }
+
+ //if (eraw < 0.05)
+ //{
+ // nhits50MeV++;
+ //}
+
+ //if (eraw < 0.03)
+ //{
+ // nhits30MeV++;
+ //}
// time
timePlot.fill(hit.getTime());
@@ -269,14 +318,11 @@
esum += eraw;
}
- hit100MeVPlot.fill(nhits100MeV);
- hit50MeVPlot.fill(nhits50MeV);
- hit30MeVPlot.fill(nhits30MeV);
+ hitUnder100MeVPlot.fill(nhits100MeV);
+ //hit50MeVPlot.fill(nhits50MeV);
+ //hit30MeVPlot.fill(nhits30MeV);
hitOver100MeVPlot.fill(nhitsOver100MeV);
-
- // Hits passed E cut.
- acceptHitPlot.fill(naccept);
-
+
// total E in Cal
ecalEPlot.fill(esum);
@@ -295,38 +341,127 @@
primaryEPlot.fill(primaryE);
// residual using highest E primary
- double singleResidual = esum - primaryE;
- residualPlot.fill(singleResidual);
+ //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);
+ //double residual = esum - eventE;
+ //residualAllPlot.fill(residual);
// hottest hit residual with highest E particle
- residualMaxHitPlot.fill(emax - primaryE);
+ //residualMaxHitPlot.fill(emax - primaryE);
// Residual with two
- double res2 = emax - e2;
- residualTop2Plot.fill(res2);
+ //double res2 = esum - e2;
+ //residualTop2Plot.fill(res2);
- double res3 = emax - e3;
- residualTop3Plot.fill(res3);
+ //double res3 = esum - e3;
+ //residualTop3Plot.fill(res3);
List<Cluster> clusters = event.get(Cluster.class, ecalClusterCollectionName);
if (clusters == null)
throw new RuntimeException("Missing cluster collection!");
+ Collections.sort(clusters, new ClusterEComparator());
nclusPlot.fill(clusters.size());
+ // Leading cluster E.
+ if (clusters.size() > 0)
+ {
+ Cluster leadClus = clusters.get(clusters.size()-1);
+ leadClusEPlot.fill(leadClus.getEnergy());
+
+ leadClusAndPrimaryPlot.fill(primaryE, leadClus.getEnergy());
+ }
+
+ // Second leading cluster E.
+ if (clusters.size() > 1)
+ {
+ leadClus2EPlot.fill(clusters.get(clusters.size()-2).getEnergy());
+ }
+
+ Map<CalorimeterHit,Integer> hitClusMap = new HashMap<CalorimeterHit,Integer>();
+
+ double clusE = 0;
+ Set<CalorimeterHit> naughty = new HashSet<CalorimeterHit>();
for (Cluster clus : clusters)
{
- clusEPlot.fill(clus.getEnergy());
+ double e = clus.getEnergy();
+ clusEPlot.fill(e);
+ clusE += e;
+ CalorimeterHit seedHit = null;
+ double maxe = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits())
+ {
+ if (hitClusMap.containsKey(hit))
+ {
+ int nshared = hitClusMap.get(hit);
+ ++nshared;
+
+ // debug
+ if (nshared >= 4)
+ naughty.add(hit);
+ //
+
+ hitClusMap.put(hit, nshared);
+ }
+ else
+ {
+ hitClusMap.put(hit, 1);
+ }
+
+ if (hit.getRawEnergy() > maxe)
+ {
+ seedHit = hit;
+ }
+ }
+
+ // Seed distance from X axis.
+ Hep3Vector pos = seedHit.getPositionVec();
+ double y = pos.y();
+ double z = pos.z();
+ clusSeedEPlot.fill(Math.sqrt(y*y+z*z));
+
+ // Seed E.
+ clusSeedEPlot.fill(seedHit.getRawEnergy());
}
- }
+
+ // debug of clusters with too many hits
+ for (CalorimeterHit duggan : naughty)
+ {
+ for (Cluster clus : clusters)
+ {
+ if (clus.getCalorimeterHits().contains(duggan))
+ {
+ System.out.println("naughty clus w/ " + clus.getCalorimeterHits().size() + " hits and E = " + clus.getEnergy());
+ for (CalorimeterHit clusHit : clus.getCalorimeterHits())
+ {
+ IDDecoder dec = clusHit.getIDDecoder();
+ dec.setID(clusHit.getCellID());
+ int ix = dec.getValue("ix");
+ int iy = dec.getValue("iy");
+ int side = dec.getValue("side");
+ System.out.println(" x, y, side, E: " + ix + ", " + iy + ", " + side + ", " + clusHit.getRawEnergy());
+ }
+ }
+ }
+ }
+
+ // Total E in all clusters.
+ clusTotEPlot.fill(clusE);
+
+ // Residual of cluster total E and E from top 3 primary particles.
+ clusResTop3Plot.fill(clusE - e3);
+
+ for (Entry<CalorimeterHit, Integer> clusHit : hitClusMap.entrySet())
+ {
+ clusHitPlot.fill(clusHit.getValue());
+ }
+ }
private double getPrimaryE(List<MCParticle> particles)
{
@@ -358,6 +493,19 @@
return primary;
}
+ private static List<MCParticle> getFinalStateParticles(List<MCParticle> mcparticles)
+ {
+ List<MCParticle> fsParticles = new ArrayList<MCParticle>();
+ for (MCParticle mcparticle : mcparticles)
+ {
+ if (mcparticle.getGeneratorStatus() == MCParticle.FINAL_STATE)
+ {
+ fsParticles.add(mcparticle);
+ }
+ }
+ return fsParticles;
+ }
+
public void setEcalCollectionName(String ecalCollectionName)
{
this.ecalCollectionName = ecalCollectionName;
