diff -N TestRunTrackReconEfficiency.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ TestRunTrackReconEfficiency.java	25 Apr 2013 21:03:40 -0000	1.1
@@ -0,0 +1,317 @@

+package org.lcsim.hps.users.omoreno;
+
+//--- java ---//
+import java.util.ArrayList;
+import java.util.List;
+
+//--- hep ---//
+import hep.aida.IHistogram1D;
+import hep.aida.IPlotter;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+//--- lcsim ---//
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.event.Track;
+import org.lcsim.event.GenericObject;
+
+//--- hps-java ---//
+import org.lcsim.hps.recon.tracking.SvtTrackExtrapolator;
+import org.lcsim.hps.recon.tracking.TrackUtils;
+import org.lcsim.hps.recon.ecal.HPSEcalCluster;
+
+/**
+ * Class to calculate track reconstruction efficiency using a tag and probe
+ * method. 
+ * 
+ * @author Omar Moreno <[log in to unmask]>
+ * @version $Id: TestRunTrackReconEfficiency.java,v 1.1 2013/04/25 21:03:40 omoreno Exp $
+ */
+public class TestRunTrackReconEfficiency  extends Driver {
+    
+    private AIDA aida;
+    private List<IPlotter>     plotters = new ArrayList<IPlotter>();
+    private List<IHistogram1D> histo1D = new ArrayList<IHistogram1D>();
+    TrackUtils trkUtil = new TrackUtils();
+    List<Track> topTracks; 
+    List<Track> botTracks;
+    SvtTrackExtrapolator extrapolator = new SvtTrackExtrapolator(); 
+    
+    double eventNumber     = 0;
+    double nOppositeVolume = 0;
+    double nWithinWindow   = 0;
+    double nAboveThreshold = 0;
+    double nTrigClusterTrackMatch = 0;
+    double findableTracks, findableTopTracks,findableBottomTracks;
+    double totalTracks, totalTopTracks, totalBottomTracks; 
+    double thresholdEnergy = 0;
+    
+    boolean debug = false;  
+    boolean topTrackIsFindable, bottomTrackIsFindable;
+    boolean topTrigger = false;
+     
+    // Collection Names
+    String stereoHitCollectionName = "HelicalTrackHits";
+    String trackCollectionName = "MatchedTracks";
+    String ecalClustersCollectionName = "EcalClusters";
+    String triggerDataCollectionName = "TriggerBank";
+    
+    /**
+     * Dflt Ctor
+     */
+    public TestRunTrackReconEfficiency(){};
+
+    /**
+     * Enable/disble debug mode
+     */
+    public void setDebug(boolean debug) {
+        this.debug = debug;
+    }
+    
+    /**
+     * 
+     */
+    public void setThresholdEnergy(double thresholdEnergy){
+    	this.thresholdEnergy = thresholdEnergy;
+    }
+     
+    /**
+     * 
+     */
+    protected void detectorChanged(Detector detector){
+        super.detectorChanged(detector);
+    }
+    
+    /**
+     * 
+     */
+    protected void process(EventHeader event)
+    {
+        eventNumber++;
+        
+        // First check if the event contains a collection of Ecal clusters.  If
+        // it doesn't skip the event.
+        if(!event.hasCollection(HPSEcalCluster.class, ecalClustersCollectionName)) return;
+       
+        // Get the list of Ecal clusters in the event
+        List<HPSEcalCluster> ecalClusters = event.get(HPSEcalCluster.class, ecalClustersCollectionName);
+        
+        // Only look at events which have two Ecal cluster 
+        if(ecalClusters.size() != 2) return;
+
+        // Check that the Ecal clusters are in opposite Ecal volumes. If 
+       	// they don't, skip the event.
+       	if(!hasClustersInOppositeVolumes(ecalClusters)){
+       		this.printDebug("Ecal clusters are not in opposite volumes");
+       		return;
+       	}
+       	nOppositeVolume++;
+      
+       	// Check that the Ecal clusters lie within some pre-defined window. If
+       	// they don't, skip the event. 
+        if(!isClusterWithinWindow(ecalClusters.get(0)) || !isClusterWithinWindow(ecalClusters.get(1))){
+        		this.printDebug("Ecal cluster falls outside of window.");
+        		return;        	
+        }
+        nWithinWindow++;
+       
+        // Check that the Ecal clusters are above the threshold energy.  If 
+        // they don't, skip the event.
+        if(!isClusterAboveEnergyThreshold(ecalClusters.get(0)) || !isClusterAboveEnergyThreshold(ecalClusters.get(1))){
+        		this.printDebug("Ecal cluster energies are below threshold.");
+        		return;        	        	
+        }
+       	nAboveThreshold++;
+ 
+       	// Check that the difference between the Ecal cluster energies is 
+        // reasonable
+       
+        // Check if the event contains a collection of tracks.  If it doesn't,
+        // move on to the next event.
+        if(!event.hasCollection(Track.class, trackCollectionName)){
+        	this.printDebug("Event doesn't contain a collection of tracks!");
+        	return;
+        }
+
+        // Get the list of tracks from the event
+        List<Track> tracks = event.get(Track.class, trackCollectionName);
+
+        // If there are no tracks in the collection, move on to the next event. 
+        if(tracks.isEmpty()){
+        	this.printDebug("Event doesn't contain any tracks!");
+        	return;  
+        }
+       
+        // Sort the tracks by SVT volume
+    	topTracks = new ArrayList<Track>();
+    	botTracks = new ArrayList<Track>();
+        for(Track track : tracks){
+    		if(track.getTrackStates().get(0).getZ0() > 0)  topTracks.add(track);
+    		else if(track.getTrackStates().get(0).getZ0() < 0) botTracks.add(track);
+    	}
+        
+       	// Get the trigger information from the event
+       	List<GenericObject> triggerData = event.get(GenericObject.class, triggerDataCollectionName);
+       	GenericObject triggerDatum = triggerData.get(0);
+       	if(triggerDatum.getIntVal(4) > 0){
+       		this.printDebug("Ecal triggered by top cluster");
+       		topTrigger = true;
+       	} else if(triggerDatum.getIntVal(5) > 0){
+       		this.printDebug("Ecal triggered by bottom cluster");       		
+       		topTrigger = false;
+       	}
+       	
+       	// Match a track to the trigger cluster
+       	HPSEcalCluster matchedCluster = null; 
+       	for(HPSEcalCluster ecalCluster : ecalClusters){
+       		if(ecalCluster.getPosition()[1] > 0 && topTrigger){
+       			if(!isClusterMatchedToTrack(ecalCluster, topTracks)){
+       				this.printDebug("Trigger cluster-track match was not found.");
+       				return;
+       			}
+       			matchedCluster = ecalCluster; 
+       			findableBottomTracks++;
+       			break;
+       		} else if( ecalCluster.getPosition()[1] < 0 && !topTrigger){
+       			if(!isClusterMatchedToTrack(ecalCluster, botTracks)){
+       				this.printDebug("Trigger cluster-track match was not found.");
+       				return;
+       			}
+       			matchedCluster = ecalCluster;
+       			findableTopTracks++;
+       			break;
+       		}
+       	}
+       	if(matchedCluster != null) ecalClusters.remove(matchedCluster);
+       	nTrigClusterTrackMatch++;
+       	
+        // If the cluster passes all requirements, then there is likely a track
+        // associated with it
+        findableTracks++;
+        
+        // Now check if a track is associated with the non-trigger cluster
+        if(topTrigger){
+       			if(!isClusterMatchedToTrack(ecalClusters.get(0), botTracks)){
+       				this.printDebug("Non trigger cluster-track match was not found.");
+       				return;
+       			}
+       			totalBottomTracks++;
+        } else if(!topTrigger){
+       			if(!isClusterMatchedToTrack(ecalClusters.get(0), topTracks)){
+       				this.printDebug("Non trigger cluster-track match was not found.");
+       				return;
+       			}        	       
+       			totalTopTracks++;
+        } 
+        ++totalTracks; 
+    }    
+
+    /**
+     * Print a debug message if they are enabled.
+     * 
+     * @param debugMessage : message to be printed
+     */ 
+    private void printDebug(String debugMessage){
+        if(debug){
+            System.out.println(this.getClass().getSimpleName() + ": " + debugMessage);
+        }
+    }
+    
+    /**
+     * 
+     */
+    private boolean isClusterWithinWindow(HPSEcalCluster clusterPosition){
+    	return true;
+    }
+    
+    /**
+     * 
+     */
+    private boolean isClusterAboveEnergyThreshold(HPSEcalCluster ecalCluster){
+    	if(ecalCluster.getEnergy() > thresholdEnergy) return true;
+    	return false;
+    }
+     
+    /**
+     * 
+     */
+    private boolean hasClustersInOppositeVolumes(List<HPSEcalCluster> ecalClusters){
+    	this.printPosition(ecalClusters.get(0).getPosition());
+    	this.printPosition(ecalClusters.get(1).getPosition());
+    	if((ecalClusters.get(0).getPosition()[1] > 0 && ecalClusters.get(1).getPosition()[1] < 0)
+    			|| (ecalClusters.get(0).getPosition()[1] < 0 && ecalClusters.get(1).getPosition()[1] > 0)){
+    		return true; 	
+    	}
+    	return false;
+    }
+    
+    /**
+     * 
+     */
+    private boolean isClusterMatchedToTrack(HPSEcalCluster cluster, List<Track> tracks){
+    	Hep3Vector clusterPos = new BasicHep3Vector(cluster.getPosition());
+    	double rMax = Double.MAX_VALUE;
+    	Track matchedTrack = null; 
+    	for(Track track : tracks){
+    		extrapolator.setTrack(track);
+    		
+    		Hep3Vector trkPosAtShowerMax = extrapolator.extrapolateTrack(clusterPos.z());
+    		if(Double.isNaN(trkPosAtShowerMax.x()) || Double.isNaN(trkPosAtShowerMax.y())){
+    			this.printDebug("Invalid track position");
+    			return false; 
+    		}
+    		this.printDebug("Track position at shower max: " + trkPosAtShowerMax.toString());
+ 
+    		// Find the distance between the track position at shower
+    		// max and the cluster position
+    		double r = VecOp.sub(trkPosAtShowerMax, clusterPos).magnitude();
+    		this.printDebug("Distance between Ecal cluster and track position at shower max: " + r + " mm");
+        	
+    		// Check if the track is the closest to the cluster.  If it is, then
+    		// save the track and contineu looping over all other tracks
+    		if (r < rMax /*&& r <= maxTrackClusterDistance*/) {
+    			rMax = r;
+    			matchedTrack = track;
+    		}
+    	}
+    	if(matchedTrack != null) return true;
+    	return false;
+    }
+    
+    /**
+     * 
+     */
+    private void printPosition(double[] position){
+    	this.printDebug("[ " + position[0] + ", " + position[1] + ", " + position[2] + " ]");
+    }
+    
+    
+    @Override
+    public void endOfData(){ 
+        System.out.println("%===================================================================% \n");
+        if(nOppositeVolume > 0){
+        	System.out.println("Total events passing opposite volume requirement: " + nOppositeVolume + " / " + eventNumber + " = " + (nOppositeVolume/eventNumber)*100 + "%");
+        }
+        if(nAboveThreshold > 0){
+        	System.out.println("Total events with both clusters above energy threshold: " + nAboveThreshold + " / " + eventNumber + " = " + (nAboveThreshold/eventNumber)*100 + "%");
+        }
+        if(nTrigClusterTrackMatch > 0){
+        	System.out.println("Total events with a trigger cluster-track match: " + nTrigClusterTrackMatch + " / " + eventNumber + " = " + (nTrigClusterTrackMatch/eventNumber)*100 + "%");
+        }
+        if(findableTracks > 0){
+            System.out.println("% Total Track Reconstruction Efficiency: " + totalTracks + " / " + findableTracks + " = " + (totalTracks / findableTracks) * 100 + "%");
+        }
+        if(findableTopTracks > 0){
+            System.out.println("% Total Top Track Reconstruction Efficiency: " + totalTopTracks + " / " + findableTopTracks + " = " + (totalTopTracks / findableTopTracks)* 100 + "%");
+        }
+        if(findableBottomTracks > 0){
+            System.out.println("% Total Bottom Track Reconstruction Efficiency: " + totalBottomTracks + " / " + findableBottomTracks + " = " + (totalBottomTracks / findableBottomTracks) * 100 + "%");
+        }
+        System.out.println("\n%===================================================================% \n");
+    }
+    
+}