 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/users/omoreno on MAIN | |||
| SvtTrackRecoEfficiency.java | +244 | -90 | 1.1 -> 1.2 |
Update to use FindableTrack
diff -u -r1.1 -r1.2 --- SvtTrackRecoEfficiency.java 18 Sep 2012 23:14:28 -0000 1.1 +++ SvtTrackRecoEfficiency.java 25 Sep 2012 01:57:48 -0000 1.2 @@ -1,50 +1,81 @@
package org.lcsim.hps.users.omoreno; //--- java ---//
+import hep.aida.ICloud2D; +import hep.aida.IHistogram1D; +import hep.aida.IHistogram2D; +import hep.aida.IPlotter; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; + +import java.util.HashMap; +import java.util.Iterator;
import java.util.List; import java.util.Map;
+import java.util.Set;
import java.util.TreeMap; import java.util.ArrayList; //--- lcsim ---// import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
import org.lcsim.detector.tracker.silicon.SiSensor; import org.lcsim.event.EventHeader; import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
import org.lcsim.event.SimTrackerHit; import org.lcsim.event.RawTrackerHit; import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseRelationalTable; +import org.lcsim.fit.helicaltrack.HelicalTrackCross; +import org.lcsim.fit.helicaltrack.HelicalTrackHit; +import org.lcsim.fit.helicaltrack.HelicalTrackStrip; +import org.lcsim.geometry.Detector;
-//--- hps-java ---//
+import org.lcsim.hps.monitoring.AIDAFrame;
import org.lcsim.hps.recon.tracking.SvtUtils;
+//--- hps-java ---// +import org.lcsim.hps.recon.tracking.TrackAnalysis;
import org.lcsim.hps.recon.tracking.TrackUtils;
+import org.lcsim.hps.recon.tracking.FindableTrack;
/** * * @author Omar Moreno <[log in to unmask]>
- * @version $Id: SvtTrackRecoEfficiency.java,v 1.1 2012/09/18 23:14:28 omoreno Exp $
+ * @version $Id: SvtTrackRecoEfficiency.java,v 1.2 2012/09/25 01:57:48 omoreno Exp $
*/
public class SvtTrackRecoEfficiency extends Driver {
- // Map to store SimTrackerHits for a given event - TreeMap<Integer, List<SimTrackerHit>> eventToSimTrackerHits = new TreeMap<Integer, List<SimTrackerHit>>();
+ private AIDA aida; + private List<IPlotter> plotters = new ArrayList<IPlotter>(); + private List<IHistogram1D> histo1D = new ArrayList<IHistogram1D>(); + private List<SimTrackerHit> simTrackerHits = new ArrayList<SimTrackerHit>(); + private Map<Track, TrackAnalysis> trkToTrkAnalysis = new HashMap<Track, TrackAnalysis>(); +
TrackUtils trkUtils = new TrackUtils();
+ FindableTrack findable = null; + TrackAnalysis trkAnalysis = null; + RelationalTable<SimTrackerHit, MCParticle> simHitToMcParticle;
// Collection Names
String simTrackerHitCollectionName = "TrackerHits";
String rawTrackerHitCollectionName = "SVTRawTrackerHits";
String trackCollectionName = "MatchedTracks";
+ String stereoHitCollectionName = "RotatedHelicalTrackHits";
int eventNumber = 0;
- // List<Integer> topSimTrackerHits = new ArrayList<Integer>(10); - // List<Integer> bottomSimTrackerHits = new ArrayList<Integer>(10);
+ int plotterIndex, histo1DIndex;
int[] topSimTrackerHits;
int[] bottomSimTrackerHits;
+ double findableTracks, foundTracks;
double topPossibleTracks, bottomPossibleTracks, possibleTracks;
double totalTopTracks, totalBottomTracks, totalTracks;
-
+ int totalLayersHit = 10;
boolean debug = false;
+ boolean trackingEfficiencyPlots = true; + boolean trackMatch = false; + boolean trackIsFindable = false;
/**
* Enable/Disable debug
@@ -52,6 +83,35 @@
public void setDebug(boolean debug){
this.debug = debug;
}
+
+ /**
+ *
+ */
+ protected void detectorChanged(Detector detector)
+ {
+ super.detectorChanged(detector);
+
+ // setup AIDA
+ aida = AIDA.defaultInstance();
+ aida.tree().cd("/");
+
+ if(trackingEfficiencyPlots){
+ plotters.add(PlotUtils.setupPlotter("Tracking Efficiency", 0, 0));
+ histo1D.add(aida.histogram1D("Tracking Efficiency", 60, 0, 6));
+ 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++;
+ }
+
+ for(IPlotter plotter : plotters){
+ plotter.show();
+ }
+ }
/**
* Dflt Ctor
@@ -59,117 +119,211 @@
public SvtTrackRecoEfficiency(){
}
- public void process(EventHeader event){
+ @Override
+ protected void process(EventHeader event)
+ {
eventNumber++;
- // If the event contains SimTrackerHits store them
+ // If the event contains SimTrackerHits store them for later use
if(event.hasCollection(SimTrackerHit.class, simTrackerHitCollectionName)){
- List<SimTrackerHit> simTrackerHits = event.get(SimTrackerHit.class, simTrackerHitCollectionName); - eventToSimTrackerHits.put(eventNumber, simTrackerHits); - System.out.println(this.getClass().getSimpleName() + ": Event: " + eventNumber + " Number of SimTrackerHits: " + simTrackerHits.size()); -
+ simTrackerHits.addAll(event.get(SimTrackerHit.class, simTrackerHitCollectionName));
+ if(debug){
+ System.out.println(this.getClass().getSimpleName() + ": Event: " + eventNumber + " Number of SimTrackerHits: " + simTrackerHits.size());
+ System.out.print(this.getClass().getSimpleName() + ": MC Particles: ");
+ for(MCParticle mcParticle : event.getMCParticles()){
+ System.out.print(mcParticle.getPDGID() + " ");
+ }
+ System.out.print("\n");
+ }
}
// The SimTrackerHits should not be stored for more than 200 events
- if(!eventToSimTrackerHits.isEmpty()){
- Integer firstEvent = eventToSimTrackerHits.firstKey();
- if(firstEvent.intValue() + 200 < eventNumber){
- eventToSimTrackerHits.remove(firstEvent);
- }
+ if(eventNumber%500 == 0 && !simTrackerHits.isEmpty()){
+ simTrackerHits.clear();
}
- // If the event doesn't contain RawTrackerHits skip it
+ // Skip the event if it doesn't contain RawTrackerHits; Only interested in triggered events
if(!event.hasCollection(RawTrackerHit.class, rawTrackerHitCollectionName)) return;
-
- if(eventToSimTrackerHits.isEmpty()){
- throw new RuntimeException(this.getClass().getSimpleName() + ": There are no SimTrackerHits associated with the RawTrackerHits");
- }
-
- List<SimTrackerHit> simTrackerHits = eventToSimTrackerHits.get(eventToSimTrackerHits.firstKey());
+
+
+ // If the collection of SimTrackerHits is empty, something is wrong
+ //if(simTrackerHits.isEmpty()){
+ // throw new RuntimeException(this.getClass().getSimpleName() + ": There are no SimTrackerHits associated with the RawTrackerHits");
+ //}
+
+ System.out.println(this.getClass().getSimpleName() + ": Number of SimTrackerHits: " + simTrackerHits.size());
+
+ // Get the MC Particles associated with the SimTrackerHits
+ List<MCParticle> mcParticles = new ArrayList<MCParticle>();
+ System.out.print(this.getClass().getSimpleName() + ": MC Particles: ");
for(SimTrackerHit simTrackerHit : simTrackerHits){
- MCParticle particle = simTrackerHit.getMCParticle();
- System.out.println(this.getClass().getSimpleName() + ": MC Particle: " + particle.getPDGID());
- System.out.println(this.getClass().getSimpleName() + ": SimTrackerHit: " + SvtUtils.getInstance().getDescription((SiSensor) simTrackerHit.getDetectorElement()));
- }
-
- this.countHitsPerLayer(simTrackerHits);
-
- // Check if a track should be found
- if(this.hasConsecutiveHits(topSimTrackerHits)){
- topPossibleTracks++;
- possibleTracks++;
- System.out.println(this.getClass().getSimpleName() + ": Found Top Possible Track");
- } else if(this.hasConsecutiveHits(bottomSimTrackerHits)){
- bottomPossibleTracks++;
- System.out.println(this.getClass().getSimpleName() + ": Found Bottom Possible Track");
- possibleTracks++;
- } else {
- return;
- }
-
- // Get the RawTrackerHits from the event
-// List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, rawTrackerHitCollectionName);
-// for(RawTrackerHit rawHit : rawHits){
-// System.out.println(this.getClass().getSimpleName() + ": RawTrackerHit: " + SvtUtils.getInstance().getDescription((SiSensor) rawHit.getDetectorElement()));
-// }
-
- // If consecutive hits were found, check if a track was also found
- if(!event.hasCollection(Track.class, trackCollectionName)) return;
-
+ System.out.print(simTrackerHit.getMCParticle().getPDGID() + " ");
+ if(mcParticles.contains(simTrackerHit.getMCParticle())) continue;
+ mcParticles.add(simTrackerHit.getMCParticle());
+ }
+ System.out.print("\n");
+
+ // 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
+ findable = new FindableTrack(event, simTrackerHits);
+ Iterator<MCParticle> mcParticleIterator = mcParticles.iterator();
+ trackIsFindable = false;
+ while(mcParticleIterator.hasNext()){
+ MCParticle mcParticle = mcParticleIterator.next();
+ if(findable.isTrackFindable(mcParticle, 8)){
+
+ // Check that all SimTrackerHits are within the same detector volume
+ Set<SimTrackerHit> trackerHits = findable.getSimTrackerHits(mcParticle);
+ if(this.isSameSvtVolume(trackerHits)){
+ if(debug) System.out.println(this.getClass().getSimpleName() + ": Track is findable ...");
+ findableTracks++;
+ trackIsFindable = true;
+ }
+ } else {
+ mcParticleIterator.remove();
+ }
+ }
+
+ // If a track is findable, check if a track was actually found otherwise return
+ if(!event.hasCollection(Track.class, trackCollectionName) || !trackIsFindable) return;
List<Track> tracks = event.get(Track.class, trackCollectionName);
-
+
+ if(mcParticles.isEmpty() && !tracks.isEmpty()){
+ throw new RuntimeException(this.getClass().getSimpleName() + ": Tracks have no associated MC particle");
+ }
+
+ // Relate a stereo hits to a SimTrackerHit; This is a required argument by TrackAnalysis
+ RelationalTable<HelicalTrackHit, SimTrackerHit> hitToMC =
+ stereoHitToMC(event.get(HelicalTrackHit.class, stereoHitCollectionName), simTrackerHits);
+
+ // Check if an MC particle is related to a found track
for(Track track : tracks){
- trkUtils.setTrack(track);
- if(trkUtils.getZ0() > 0){
- totalTopTracks++;
- System.out.println(this.getClass().getSimpleName() + ": Found Top Track");
- }
- else{
- totalBottomTracks++;
- System.out.println(this.getClass().getSimpleName() + ": Found Bottom Track");
+ trkAnalysis = new TrackAnalysis(track, hitToMC);
+ if(mcParticles.contains(trkAnalysis.getMCParticle())){
+ System.out.println(this.getClass().getSimpleName() + ": Track match found");
+ foundTracks++;
+ if(trackingEfficiencyPlots){
+ aida.histogram1D("Tracking Efficiency").fill(trkAnalysis.getMCParticle().getMomentum().magnitude(), 1);
+ aida.histogram1D("Momentum").fill(trkAnalysis.getMCParticle().getMomentum().magnitude());
+ }
+ mcParticles.remove(trkAnalysis.getMCParticle());
+ }
+ }
+
+
+ if(trackingEfficiencyPlots){
+ // If the list still contains MC Particles, a matching track wasn't found
+ for(MCParticle mcParticle : mcParticles){
+ aida.histogram1D("Tracking Efficiency").fill(mcParticle.getMomentum().magnitude(), 0);
+ aida.histogram1D("Momentum").fill(mcParticle.getMomentum().magnitude());
}
- totalTracks++;
}
}
-
- /**
- *
- */
- private boolean hasConsecutiveHits(int[] topSimTrackerHits2){
- for(int index = 0; index < topSimTrackerHits2.length; index++){
- if(topSimTrackerHits2[index] == 0) return false;
- }
- return true;
+
+ private boolean isSameSvtVolume(Set<SimTrackerHit> simTrackerHits){
+ int volumeIndex = 0;
+ for(SimTrackerHit simTrackerHit : simTrackerHits){
+ if(SvtUtils.getInstance().isTopLayer((SiSensor) simTrackerHit.getDetectorElement())) volumeIndex++;
+ else volumeIndex--;
+ }
+ return Math.abs(volumeIndex) == simTrackerHits.size();
}
- /**
- *
- */
- private void countHitsPerLayer(List<SimTrackerHit> simTrackerHits){
- topSimTrackerHits = new int[10];
- bottomSimTrackerHits = new int[10];
-
- SiSensor sensor = null;
+ private RelationalTable<HelicalTrackHit, SimTrackerHit> stereoHitToMC(List<HelicalTrackHit> stereoHits, List<SimTrackerHit> simTrackerHits){
+ RelationalTable hitToMC = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ Map<Integer, List<SimTrackerHit>> layerToSimTrackerHit = new HashMap<Integer, List<SimTrackerHit>>();
+ Map<Integer, List<HelicalTrackStrip>> layerToHelicalTrackStrip = new HashMap<Integer, List<HelicalTrackStrip>>();
+
+ // Sort the SimTrackerHits by Layer
for(SimTrackerHit simTrackerHit : simTrackerHits){
- // The only hits of interest are those created by an electron (or positron)
- if(Math.abs(simTrackerHit.getMCParticle().getPDGID()) != 11) continue;
- sensor = (SiSensor) simTrackerHit.getDetectorElement();
- if(SvtUtils.getInstance().isTopLayer(sensor)){
- topSimTrackerHits[SvtUtils.getInstance().getLayerNumber(sensor) - 1]++;
- } else {
- bottomSimTrackerHits[SvtUtils.getInstance().getLayerNumber(sensor) - 1]++;
+ if(!layerToSimTrackerHit.containsKey(simTrackerHit.getLayer()))
+ layerToSimTrackerHit.put(simTrackerHit.getLayer(), new ArrayList<SimTrackerHit>());
+ layerToSimTrackerHit.get(simTrackerHit.getLayer()).add(simTrackerHit);
+ }
+
+ for(HelicalTrackHit stereoHit : stereoHits){
+ HelicalTrackCross cross = (HelicalTrackCross) stereoHit;
+ for(HelicalTrackStrip strip : cross.getStrips()){
+ if(layerToSimTrackerHit.get(strip.layer()).size() == 1){
+ hitToMC.add(stereoHit, layerToSimTrackerHit.get(strip.layer()).get(0).getMCParticle());
+ layerToSimTrackerHit.remove(strip.layer());
+ }
+ else if(layerToSimTrackerHit.get(strip.layer()).size() > 1){
+ System.out.println(this.getClass().getSimpleName() + ": Layer with multiple hits found.");
+ }
+ }
}
+
+ return hitToMC;
+ }
+
+ public boolean isTrackFindable(MCParticle mcParticle, int nLayers){
+
+ if(nLayers%2 == 1) throw new RuntimeException(this.getClass().getSimpleName() + ": The required number of layers hit must be even");
+
+ System.out.println(this.getClass().getSimpleName() + ": The number of relations in the table are: " + simHitToMcParticle.size());
+
+ System.out.println(this.getClass().getSimpleName() + ": Number of SimTrackerHits associated with the MC Particle: " + simHitToMcParticle.allTo(mcParticle).size());
+
+ // Get the list of SimTrackerHits associated with the MC particle
+ Set<SimTrackerHit> simHits = simHitToMcParticle.allTo(mcParticle);
+
+ System.out.println("Number of simtracker hits " + simHits.size());
+
+ // Find the layers hit
+ boolean[] layerHit = new boolean[totalLayersHit];
+ for(SimTrackerHit simHit : simHits){
+ layerHit[simHit.getLayer()-1] = true;
}
+
+ int nLayersHit = 0;
+ // Check how many pairs of layers were hit
+ for(int index = 0; index < totalLayersHit; index += 2){
+ System.out.println(this.getClass().getSimpleName() + ": " + layerHit[index] + " " + layerHit[index+1]);
+ if(layerHit[index] && layerHit[index+1]){
+ System.out.println(this.getClass().getSimpleName() + ": Layer was hit!");
+ nLayersHit += 2;
+ }
+ }
+
+ System.out.println(this.getClass().getSimpleName() + ": Number of layers hit: " + nLayersHit);
+
+ return nLayers == nLayersHit;
+ }
+
+ private Hep3Vector getClusterPosition(HelicalTrackStrip strip){
+ Hep3Vector origin = strip.origin();
+ Hep3Vector u = strip.u();
+ double umeas = strip.umeas();
+ Hep3Vector uvec = VecOp.mult(umeas, u);
+ return VecOp.add(origin, uvec);
}
@Override
public void endOfData(){
+
+
+ if(trackingEfficiencyPlots){
+ double[] efficiency = new double[60];
+ for(int index = 0; index < 60; index++){
+ if(aida.histogram1D("Tracking Efficiency").binEntries(index) == 0) efficiency[index] = 0;
+ else {
+ efficiency[index] = aida.histogram1D("Tracking Efficiency").binHeight(index)/aida.histogram1D("Tracking Efficiency").binEntries(index);
+ }
+ System.out.println("Tracking Efficiency - Bin " + index + " " + efficiency[index]);
+ }
+ aida.histogram1D("Tracking Efficiency").reset();
+ for(int index = 0; index < 60; index++){
+ aida.histogram1D("Tracking Efficiency").fill(index, efficiency[index]);
+ }
+ }
+
System.out.println("%===============================================================%");
System.out.println("%============== Track Reconstruction Efficiencies ==============%");
System.out.println("%===============================================================%\n%");
- if(possibleTracks > 0){
- System.out.println("% Total Track Reconstruction Efficiency: " + (totalTracks/possibleTracks)*100 + "%");
+ if(findableTracks > 0){
+ System.out.println("% Total Track Reconstruction Efficiency: " + foundTracks + "/" + findableTracks + "=" + (foundTracks/findableTracks)*100 + "%");
}
if(topPossibleTracks > 0){
System.out.println("% Top Track Reconstruction Efficiency: " + (totalTopTracks/topPossibleTracks)*100 + "%");
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1