hps-java/src/main/java/org/lcsim/hps/users/omoreno
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");
+ }
+
+}