hps-java/src/main/java/org/lcsim/hps/users/omoreno
diff -u -r1.7 -r1.8
--- SvtTrackRecoEfficiency.java 21 Nov 2012 08:33:18 -0000 1.7
+++ SvtTrackRecoEfficiency.java 18 Dec 2012 08:53:00 -0000 1.8
@@ -22,6 +22,7 @@
import hep.physics.vec.VecOp;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
//--- lcsim ---//
@@ -51,7 +52,7 @@
/**
*
* @author Omar Moreno <[log in to unmask]>
- * @version $Id: SvtTrackRecoEfficiency.java,v 1.7 2012/11/21 08:33:18 omoreno Exp $
+ * @version $Id: SvtTrackRecoEfficiency.java,v 1.8 2012/12/18 08:53:00 omoreno Exp $
*/
public class SvtTrackRecoEfficiency extends Driver {
@@ -76,6 +77,7 @@
String fittedRawTrackerHitCollectionName = "SVTFittedRawTrackerHits";
String trackCollectionName = "MatchedTracks";
String stereoHitCollectionName = "RotatedHelicalTrackHits";
+ String siTrackerHitCollectionName = "StripClusterer_SiTrackerHitStrip1D";
int eventNumber = 0;
int plotterIndex, histo1DIndex, histo2DIndex;
@@ -99,58 +101,47 @@
/**
* Enable/Disable debug
*/
- public void setDebug(boolean debug)
- {
+ public void setDebug(boolean debug){
this.debug = debug;
}
/**
* Set the name of the file to output efficiency data to
*/
- public void setEfficiencyOutputFile(String efficiencyOutputFile)
- {
+ public void setEfficiencyOutputFile(String efficiencyOutputFile){
this.efficiencyOutputFile = efficiencyOutputFile;
}
/**
* Set the name of the file to output momentum data to
*/
- public void setMomentumOutputFile(String momentumOutputFile)
- {
+ public void setMomentumOutputFile(String momentumOutputFile){
this.momentumOutputFile = momentumOutputFile;
}
/**
* Set the required number of layers an MC particle must hit
*/
- public void setTotalLayersHit(int totalLayersHit)
- {
+ public void setTotalLayersHit(int totalLayersHit){
if(totalLayersHit%2 == 1) throw new RuntimeException(this.getClass().getSimpleName() + ": Total number of layers hit must be even");
this.totalLayersHit = totalLayersHit;
}
/**
+ * Print a debug statement
*
- * @param strategyResource
- */
- public void setStrategyResourcePath(String strategyResourcePath)
- {
- this.strategyResourcePath = strategyResourcePath;
- }
-
- /**
- * Print debug
+ * @param message : debug message
*/
- private void printDebug(String message)
- {
+ private void printDebug(String message){
+ if(debug){
System.out.println(this.getClass().getSimpleName() + ": " + message);
+ }
}
/**
*
*/
- protected void detectorChanged(Detector detector)
- {
+ protected void detectorChanged(Detector detector){
super.detectorChanged(detector);
// setup AIDA
@@ -179,66 +170,27 @@
PlotUtils.setup1DRegion(plotters.get(plotterIndex), "Tracking Efficiency", 0, "Momentum [GeV]", histo1D.get(histo1DIndex));
histo1DIndex++;
plotterIndex++;
+
plotters.add(PlotUtils.setupPlotter("Momentum", 0, 0));
histo1D.add(aida.histogram1D("Momentum", 60, 0, 6));
PlotUtils.setup1DRegion(plotters.get(plotterIndex), "Momentum", 0, "Momentum [GeV]", histo1D.get(histo1DIndex));
plotterIndex++;
histo1DIndex++;
- plotters.add(PlotUtils.setupPlotter("Stereo Hits", 0, 0));
- histo1D.add(aida.histogram1D("Stereo Hits", 5, 0, 5));
- PlotUtils.setup1DRegion(plotters.get(plotterIndex), "Stereo Hits", 0, "Number of Stereo Hits", histo1D.get(histo1DIndex));
- plotterIndex++;
- histo1DIndex++;
-
- plotters.add(PlotUtils.setupPlotter("Hit Positions - Missed SimTrackerHits", 5, 4));
- for(SiSensor sensor : SvtUtils.getInstance().getSensors()){
- String title = SvtUtils.getInstance().getDescription(sensor) + " - Hit Positions - Missed SimTrackerHits";
- if(SvtUtils.getInstance().isTopLayer(sensor)){
- histo2D.add(aida.histogram2D(title, 50, -50, 50, 50, 0, 50));
- } else {
- histo2D.add(aida.histogram2D(title, 50, -50, 50, 50, -50, 0));
- }
- PlotUtils.setup2DRegion(plotters.get(plotterIndex), title, PlotUtils.getPlotterRegion(sensor), "x [mm]", "y [mm]", histo2D.get(histo2DIndex));
- histo2DIndex++;
- }
- plotterIndex++;
-
- plotters.add(PlotUtils.setupPlotter("t0 - Missed SimTrackerHits", 5, 4));
- for(SiSensor sensor : SvtUtils.getInstance().getSensors()){
- String title = SvtUtils.getInstance().getDescription(sensor) + " - t0 - Missed SimTrackerHits";
- histo1D.add(aida.histogram1D(title, 50, 0, 100));
- PlotUtils.setup1DRegion(plotters.get(plotterIndex), title, PlotUtils.getPlotterRegion(sensor), "t0 [ns]", histo1D.get(histo1DIndex));
- histo1DIndex++;
- }
- plotterIndex++;
-
- plotters.add(PlotUtils.setupPlotter("Hit Positions - Found SimTrackerHits", 5, 4));
- for(SiSensor sensor : SvtUtils.getInstance().getSensors()){
- String title = SvtUtils.getInstance().getDescription(sensor) + " - Hit Positions - Found SimTrackerHits";
- if(SvtUtils.getInstance().isTopLayer(sensor)){
- histo2D.add(aida.histogram2D(title, 50, -50, 50, 50, 0, 50));
- } else {
- histo2D.add(aida.histogram2D(title, 50, -50, 50, 50, -50, 0));
- }
- PlotUtils.setup2DRegion(plotters.get(plotterIndex), title, PlotUtils.getPlotterRegion(sensor), "x [mm]", "y [mm]", histo2D.get(histo2DIndex));
- histo2DIndex++;
- }
- plotterIndex++;
-
- plotters.add(PlotUtils.setupPlotter("t0 - Found SimTrackerHits", 5, 4));
- for(SiSensor sensor : SvtUtils.getInstance().getSensors()){
- String title = SvtUtils.getInstance().getDescription(sensor) + " - t0 - Found SimTrackerHits";
- histo1D.add(aida.histogram1D(title, 50, 0, 100));
- PlotUtils.setup1DRegion(plotters.get(plotterIndex), title, PlotUtils.getPlotterRegion(sensor), "t0 [ns]", histo1D.get(histo1DIndex));
- histo1DIndex++;
- }
- plotterIndex++;
}
for(IPlotter plotter : plotters){
plotter.show();
}
}
+
+ private String samplesToString(short[] samples){
+ String sampleList = "[ ";
+ for(short sample : samples){
+ sampleList += Short.toString(sample) + ", ";
+ }
+ sampleList += "]";
+ return sampleList;
+ }
/**
* Dflt Ctor
@@ -246,34 +198,45 @@
public SvtTrackRecoEfficiency(){}
@Override
- protected void process(EventHeader event)
- {
+ protected void process(EventHeader event){
+ // For now, only look at events with a single track
if(event.get(Track.class, trackCollectionName).size() > 1) return;
-
eventNumber++;
// If the event doesn't contain SimTrackerHits, skip the event
if(!event.hasCollection(SimTrackerHit.class, simTrackerHitCollectionName)) return;
List<SimTrackerHit> simTrackerHits = event.get(SimTrackerHit.class, simTrackerHitCollectionName);
- if(debug)
- this.printDebug("Event " + eventNumber + " contains " + simTrackerHits.size() + " SimTrackerHits");
-
- // Add the SimTrackerHits to its respective sensor readout
- /*for(SimTrackerHit simHitTrackerHit : simTrackerHits){
- ((SiSensor) simHitTrackerHit.getDetectorElement()).getReadout().addHit(simHitTrackerHit);
- }*/
+ this.printDebug("\nEvent " + eventNumber + " contains " + simTrackerHits.size() + " SimTrackerHits");
+ // Loop through all SimTrackerHits and confirm that a corresponding RawTrackerHit was created
+ for(SimTrackerHit simTrackHit : simTrackerHits){
+
+ this.printDebug("SimTrackerHit Layer Number: " + simTrackHit.getLayerNumber());
+ }
// Get the list of RawTrackerHits and add them to the sensor readout
List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, rawTrackerHitCollectionName);
+ String volume;
for(RawTrackerHit rawHit : rawHits){
+ if(SvtUtils.getInstance().isTopLayer((SiSensor)rawHit.getDetectorElement())){
+ volume = "Top Volume ";
+ } else {
+ volume = "Bottom Volume ";
+ }
+ this.printDebug(volume + "RawTrackerHit Channel #: " + rawHit.getIdentifierFieldValue("strip") + " Layer Number: " + rawHit.getLayerNumber()
+ + " Samples: " + samplesToString(rawHit.getADCValues()));
((SiSensor) rawHit.getDetectorElement()).getReadout().addHit(rawHit);
}
- // Get the list of HPSFittedRawTrackerHits and add them to the sensor readout
- List<HPSFittedRawTrackerHit> fittedHits = event.get(HPSFittedRawTrackerHit.class, fittedRawTrackerHitCollectionName);
- for(HPSFittedRawTrackerHit fittedHit : fittedHits){
- ((SiSensor) fittedHit.getRawTrackerHit().getDetectorElement()).getReadout().addHit(fittedHit);
+ if(event.hasCollection(SiTrackerHit.class, siTrackerHitCollectionName)){
+ List<SiTrackerHit> hitlist = event.get(SiTrackerHit.class, siTrackerHitCollectionName);
+ for(SiTrackerHit siTrackerHit : hitlist){
+ this.printDebug("Cluster layer: " + SvtUtils.getInstance().getLayerNumber(siTrackerHit.getSensor()));
+ this.printDebug("Cluster is comprised by the following raw hits:");
+ for(RawTrackerHit rawHit : siTrackerHit.getRawHits()){
+ this.printDebug("RawTrackerHit Channel #: " + rawHit.getIdentifierFieldValue("strip") + " Layer Number: " + rawHit.getLayerNumber());
+ }
+ }
}
// Get the MC Particles associated with the SimTrackerHits
@@ -289,11 +252,11 @@
// Get the magnetic field
Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
- this.printDebug("BField: " + event.getDetector().getFieldMap().getField(IP).y());
+ //this.printDebug("BField: " + event.getDetector().getFieldMap().getField(IP).y());
// Check if the MC particle track should be found by the tracking algorithm
- // Note: Only require 4 of the 5 SVT layers to be hit; Layers must be in the same SVT volume
findable = new FindableTrack(event, simTrackerHits);
+
// Use an iterator to avoid ConcurrentModificationException
Iterator<MCParticle> mcParticleIterator = mcParticles.iterator();
trackIsFindable = false;
@@ -317,20 +280,14 @@
}
}
- // If a track is findable, check if a track was actually found otherwise return
- if(!event.hasCollection(Track.class, trackCollectionName) || !trackIsFindable){
- if(trackIsFindable)
- throw new RuntimeException("Track is not suppose to be 'findable");
- return;
- }
+ // Check if the event contains a Track collection, otherwise return
+ if(!event.hasCollection(Track.class, trackCollectionName)) return;
List<Track> tracks = event.get(Track.class, trackCollectionName);
- if(debug)
- this.printDebug("Event " + eventNumber + " contains " + tracks.size() + " Tracks");
+ this.printDebug("Event " + eventNumber + " contains " + tracks.size() + " Tracks");
// Relate a stereo hits to a SimTrackerHit; This is a required argument by TrackAnalysis
List<HelicalTrackHit> stereoHits = event.get(HelicalTrackHit.class, stereoHitCollectionName);
- if(debug)
- this.printDebug("Event " + eventNumber + " contains " + stereoHits.size() + " HelicalTrackHits");
+ this.printDebug("Event " + eventNumber + " contains " + stereoHits.size() + " HelicalTrackHits");
RelationalTable<HelicalTrackHit, MCParticle> hitToMC = stereoHitToMC(stereoHits, simTrackerHits);
// Check if an MC particle is related to a found track
@@ -350,9 +307,8 @@
if(!mcParticles.isEmpty() && trackingEfficiencyPlots){
// If the list still contains MC Particles, a matching track wasn't found
- if(debug)
- this.printDebug("No matching track found");
-
+ this.printDebug("No matching track found");
+
// Check that all stereoHits were correctly assigned to an MCParticle
for(MCParticle mcParticle : mcParticles){
@@ -410,151 +366,11 @@
bottomStereoLayerHit[(stereoHit.Layer() - 1)/2] = true;
}
}
-
- // Loop through all of the layers and check if it contains both a SimTrackerHit pair
- // and a stereo hit
- SiSensor axialSensor = null;
- SiSensor stereoSensor = null;
- String title = null;
- for(int index = 0; index < layerHit.length; index++){
- if(isTopTrack){
- axialSensor = SvtUtils.getInstance().getTopSensor(2*index + 1, 0);
- stereoSensor = SvtUtils.getInstance().getTopSensor(2*index + 2, 0);
-
- if(axialSensor.getReadout().getHits(SimTrackerHit.class).isEmpty()
- || stereoSensor.getReadout().getHits(SimTrackerHit.class).isEmpty()){
- this.printDebug("One of the sensors does not have a SimTrackerHit associated with it; moving on to next layer");
- continue;
- }
-
- this.printDebug(SvtUtils.getInstance().getDescription(axialSensor) + ": SimTrackerHit: " + axialSensor.getReadout().getHits(SimTrackerHit.class).size());
- this.printDebug(SvtUtils.getInstance().getDescription(stereoSensor) + ": SimTrackerHit: " + stereoSensor.getReadout().getHits(SimTrackerHit.class).size());
- if(layerHit[index]){
- this.printDebug("Top Layer " + (index+1) + " was hit");
- }
-
- for(SimTrackerHit simHit : axialSensor.getReadout().getHits(SimTrackerHit.class)){
- if(layerHit[index] && !topStereoLayerHit[index]){
- this.printDebug("Top Layer " + (index+1) + " is missing a stereo hit");
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - Hit Positions - Missed SimTrackerHits";
- } else if(layerHit[index] && topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - Hit Positions - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
- aida.histogram2D(title).fill(simHit.getPositionVec().x(), simHit.getPositionVec().y());
- }
-
- for(HPSFittedRawTrackerHit fittedRawHit : axialSensor.getReadout().getHits(HPSFittedRawTrackerHit.class)){
- if(layerHit[index] && !topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - t0 - Missed SimTrackerHits";
- } else if(layerHit[index] && topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - t0 - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
-
- aida.histogram1D(title).fill(fittedRawHit.getT0());
- }
-
- for(SimTrackerHit simHit : stereoSensor.getReadout().getHits(SimTrackerHit.class)){
- if(layerHit[index] && !topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - Hit Positions - Missed SimTrackerHits";
- } else if(layerHit[index] && topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - Hit Positions - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
- aida.histogram2D(title).fill(simHit.getPositionVec().x(), simHit.getPositionVec().y());
- }
-
- for(HPSFittedRawTrackerHit fittedRawHit : stereoSensor.getReadout().getHits(HPSFittedRawTrackerHit.class)){
- if(layerHit[index] && !topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - t0 - Missed SimTrackerHits";
- } else if(layerHit[index] && topStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - t0 - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
-
- aida.histogram1D(title).fill(fittedRawHit.getT0());
- }
-
- } else {
- axialSensor = SvtUtils.getInstance().getBottomSensor(2*index + 2, 0);
- stereoSensor = SvtUtils.getInstance().getBottomSensor(2*index + 1, 0);
-
- if(axialSensor.getReadout().getHits(SimTrackerHit.class).isEmpty()
- || stereoSensor.getReadout().getHits(SimTrackerHit.class).isEmpty()){
- this.printDebug("One of the sensors does not have a SimTrackerHit associated with it; moving on to next layer");
- continue;
- }
-
- this.printDebug(SvtUtils.getInstance().getDescription(axialSensor) + ": SimTrackerHit: " + axialSensor.getReadout().getHits(SimTrackerHit.class).size());
- this.printDebug(SvtUtils.getInstance().getDescription(stereoSensor) + ": SimTrackerHit: " + stereoSensor.getReadout().getHits(SimTrackerHit.class).size());
- if(layerHit[index]){
- this.printDebug("Bottom Layer " + (index+1) + " was hit");
- }
-
- for(SimTrackerHit simHit : axialSensor.getReadout().getHits(SimTrackerHit.class)){
- if(debug)
- this.printDebug(SvtUtils.getInstance().getDescription(axialSensor) + ": Position: " + simHit.getPositionVec().toString());
-
- if(layerHit[index] && !bottomStereoLayerHit[index]){
- this.printDebug("Bottom Layer " + (index+1) + " is missing a stereo hit");
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - Hit Positions - Missed SimTrackerHits";
- } else if(layerHit[index] && bottomStereoLayerHit[index]) {
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - Hit Positions - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
- aida.histogram2D(title).fill(simHit.getPositionVec().x(), simHit.getPositionVec().y());
- }
-
- for(HPSFittedRawTrackerHit fittedRawHit : axialSensor.getReadout().getHits(HPSFittedRawTrackerHit.class)){
- if(layerHit[index] && !bottomStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - t0 - Missed SimTrackerHits";
- } else if(layerHit[index] && bottomStereoLayerHit[index]) {
- title = SvtUtils.getInstance().getDescription(axialSensor) + " - t0 - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
- aida.histogram1D(title).fill(fittedRawHit.getT0());
- }
-
- for(SimTrackerHit simHit : stereoSensor.getReadout().getHits(SimTrackerHit.class)){
- if(layerHit[index] && !bottomStereoLayerHit[index]){
- this.printDebug("Bottom Layer " + (index+1) + " is missing a stereo hit");
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - Hit Positions - Missed SimTrackerHits";
- } else if(layerHit[index] && bottomStereoLayerHit[index]) {
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - Hit Positions - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
- aida.histogram2D(title).fill(simHit.getPositionVec().x(), simHit.getPositionVec().y());
- }
-
- for(HPSFittedRawTrackerHit fittedRawHit : stereoSensor.getReadout().getHits(HPSFittedRawTrackerHit.class)){
- if(layerHit[index] && !bottomStereoLayerHit[index]){
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - t0 - Missed SimTrackerHits";
- } else if(layerHit[index] && bottomStereoLayerHit[index]) {
- title = SvtUtils.getInstance().getDescription(stereoSensor) + " - t0 - Found SimTrackerHits";
- } else throw new RuntimeException("Histogram title was not found");
- aida.histogram1D(title).fill(fittedRawHit.getT0());
- }
- }
- }
-
+
aida.histogram1D("Tracking Efficiency").fill(mcParticle.getMomentum().magnitude(), 0);
aida.histogram1D("Momentum").fill(mcParticle.getMomentum().magnitude());
- aida.histogram1D("Stereo Hits").fill(hitToMC.allTo(mcParticle).size());
-
- // Check if all HelicalTrackHits were properly assigned to an MC particle
- Set<HelicalTrackHit> mcpStereoHits = hitToMC.allTo(mcParticle);
- Iterator<HelicalTrackHit> stereoHitIterator = stereoHits.iterator();
- while(stereoHitIterator.hasNext()){
- HelicalTrackHit stereoHit = stereoHitIterator.next();
- if(mcpStereoHits.contains(stereoHit)){
- this.printDebug("MC particle: " + mcParticle.getPDGID() + ": Matched stereo hit position: " + stereoHit.getCorrectedPosition().toString());
- stereoHitIterator.remove();
- }
- }
-
}
-
- // If the list of HelicalTrackHits is not empty, there exist unmatched hits
- /*if(!stereoHits.isEmpty()){
- for(HelicalTrackHit stereoHit : stereoHits){
- this.printDebug("Unmatched stereo hit position: " + stereoHit.getCorrectedPosition().toString());
- }
- }*/
}
}