lcsim/src/org/lcsim/contrib/seedtracker/analysis
diff -u -r1.1 -r1.2
--- AnalysisUtils.java 4 Aug 2008 18:43:20 -0000 1.1
+++ AnalysisUtils.java 12 Aug 2008 23:41:59 -0000 1.2
@@ -9,20 +9,34 @@
package org.lcsim.contrib.seedtracker.analysis;
+import hep.physics.vec.BasicHep3Vector;
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 org.lcsim.contrib.seedtracker.SeedCandidate;
+import org.lcsim.event.Track;
+import org.lcsim.contrib.seedtracker.SeedTrack;
+import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
/**
* This class contains utilities useful for seedtracker performance studies
- * @author Richard Partridge
+ * @authors Richard Partridge,Pelham Keahey
* @version 1.0
*/
public class AnalysisUtils {
-
+ //used to cache a current event in order to prevent needless looping
+ private static EventHeader currentEvent=null;
+ //trkmap returned in the FindTrackMap method
+ private static Map<MCParticle,Track> trkmap;
/** Creates a new instance of AnalysisUtils */
public AnalysisUtils() {
}
@@ -77,4 +91,216 @@
// Done looping over all MCParticle candidates - return the list with any matches that were found
return mcplist;
}
-}
\ No newline at end of file
+ /**
+ * Finds the purity of the given track
+ * @param trk find purity of this track
+ * @return a double purity
+ */
+ public static double FindTrackPurity(Track trk){
+ double purity=0;
+ //MCParticle and Int map, int will be the number of times this mcparticle shows up in the track
+ Map<MCParticle,Integer> occurrences = new HashMap<MCParticle,Integer>();
+ double mostoccured;
+ int index;
+ //creates a seedcandiate for the track
+ if (trk instanceof SeedTrack) {
+ SeedTrack strk = (SeedTrack) trk;
+ for (HelicalTrackHit h : strk.getSeedCandidate().getHits()){
+ //Create List of MCParticles in the SeedCandidate
+ List<MCParticle> mcpl = h.getMCParticles();
+ //loop through these mcparticles and count how many times that particle occurs
+ for(MCParticle part : mcpl){
+ if(occurrences.containsKey(part)){
+ index = occurrences.get(part);
+ occurrences.put(part,index+1);
+ }
+ else {
+ occurrences.put(part, 1);
+ }
+ }
+
+ }
+ //If only one MCParticle exsists in the map, then purity if perfect
+ if(occurrences.size()==1){
+ purity =1;
+ }
+ //otherwise save the value of the most occuring mcparticle against the total size of the track's hitlist
+ else{
+ mostoccured = -1;
+ for(MCParticle p : occurrences.keySet()){
+ int occured = occurrences.get(p);
+ if(occured>mostoccured){
+ mostoccured = occured;
+ }
+ }
+ purity = mostoccured / strk.getSeedCandidate().getHits().size();
+ }
+ }
+ return purity;
+ }
+ /**
+ * @param event
+ * @return a map of MCParticles and Integers. The int is either 1 or 0, 1 being track found 0, track not found
+ * For the MCParticles in the map are only particles that can be expected to have a found track based upon cuts in: pt, cos(theta),
+ * z0 and dca
+ * THIS METHOD USES DEFAULT CUTS FOR MCParticle:
+ * pt>1.1 && Math.abs(dca) < 9 && Math.abs(z0) < 9 && Math.cos(theta) < 0.985
+ */
+ public static Map<MCParticle, Integer> FindEfficiency(EventHeader event){
+ return FindEfficiencyParam(event, 1.1 , 9.0 , 9.0 , 0.985);
+ }
+ /**
+ * @param event
+ * @param ptcut min transverse mom.
+ * @param dcacut max dca
+ * @param z0cut max z0
+ * @param costhetacut max cos(theta)
+ * @return Returns a map of MCParticles and Integers. The int is either 1 or 0, 1 being track found 0, track not found
+ * For the MCParticles in the map are only particles that can be expected to have a found track based upon cuts in: pt, cos(theta),
+ * z0 and dca
+ * THIS METHOD ALLOWS CUTS TO BE SET
+ */
+ public static Map<MCParticle,Integer> FindEfficiencyParam(EventHeader event,double ptcut,double dcacut,double z0cut,double costhetacut){
+ //prevent from creating a trkmap over and over again if its the same event (saves a lot of time)
+ if(currentEvent==null || !currentEvent.equals(event)) {
+ trkmap = FindTrackMap(event);
+ currentEvent = event;
+ }
+ //cut parameters
+ double ptmin = ptcut;
+ double dcamax = dcacut;
+ double z0max = z0cut;
+ double costhetamax = costhetacut;
+
+ //Map that will be returned giving the particle and if a track was found for it
+ Map<MCParticle,Integer> effimap = new HashMap<MCParticle,Integer>();
+
+ //Loop through all the MCParticles in the event
+ for(MCParticle mcp : event.getMCParticles()){
+ //0: if no track 1: if track found
+ int wgt=0;
+ //used to calcute the helix params of mcp to be applied to cuts
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, event.getDetector().getFieldMap().getField(new BasicHep3Vector(0.,0.,0.)).z());
+ double theta = calc.getTheta();
+ double pt = calc.getMCTransverseMomentum();
+ double dca = calc.getDCA();
+ double z0 = calc.getZ0();
+
+ //Map of mcp with a string that will be a list of layers the particle hit
+ Map<MCParticle,List<String>> hitlayermap = new HashMap<MCParticle,List<String>>();
+ List<HelicalTrackHit> listjmk = event.get(HelicalTrackHit.class,"HelicalTrackHits");
+
+ //Sort the list of hits in order of there time of impact
+ Collections.sort(listjmk, new Comparator() {
+ public int compare(Object o1, Object o2) {
+ HelicalTrackHit h1 = (HelicalTrackHit) o1;
+ HelicalTrackHit h2 = (HelicalTrackHit) o2;
+
+ double t1, t2;
+ try{
+ t1 = ((RawTrackerHit)h1.getRawHits().get(0)).getSimTrackerHit().get(0).getTime();
+ t2 = ((RawTrackerHit)h2.getRawHits().get(0)).getSimTrackerHit().get(0).getTime();
+ } catch(NullPointerException npe) {
+ t1 = h1.getTime();
+ t2 = h2.getTime();
+ }
+
+ return Double.compare(t1, t2);
+ }
+ });
+
+ //Loop through the hits
+ for(HelicalTrackHit hito : listjmk){
+ //get the mcp from the hitlist
+ List<MCParticle> mplist = hito.getMCParticles();
+ for(MCParticle mp : mplist){
+ //if the mcp is already in the map add to its hit order list
+ if(hitlayermap.containsKey(mp)){
+ hitlayermap.get(mp).add(hito.getLayerIdentifier());
+ }
+ //if the mcp isn't already in the list, create a new list of strings to add to the map
+ else{
+ ArrayList<String> layers = new ArrayList<String>();
+ layers.add(hito.getLayerIdentifier());
+ hitlayermap.put(mp, layers);
+ ArrayList<HelicalTrackHit> trackhits = new ArrayList<HelicalTrackHit>();
+ trackhits.add(hito);
+ }
+ }
+ }
+
+ //check for duplicate MCParticles, also used to determine the number of hits
+ HashSet<String> dupemcp = new HashSet<String>();
+ if(hitlayermap.containsKey(mcp)){
+ for(String e : hitlayermap.get(mcp)){
+ dupemcp.add(e);
+ }
+ }
+
+ //To be considered a mcp must have at least seven hit and actually...and be there in the map
+ if(hitlayermap.containsKey(mcp) && dupemcp.size()>=7){
+ //if a track was found, the efficiency is 1, put it into the effimap
+ if(trkmap.containsKey(mcp)){
+ wgt=1;
+ effimap.put(mcp, wgt);
+ }
+ //if no track was found, but meets all the parameters put it into the map with wgt=0
+ else if(pt>ptmin && Math.abs(dca) < dcamax && Math.abs(z0) < z0max && Math.cos(theta) < costhetamax ){
+ effimap.put(mcp, wgt);
+ }
+ }
+ }
+ //return the efficiency map
+ return effimap;
+ }
+ /**
+ *
+ * @param event
+ * @return a track map of MCParticles with its associated track
+ */
+ public static Map<MCParticle, Track> FindTrackMap(EventHeader event){
+ //get list of tracks
+ List<Track> tracklist = event.getTracks();
+ //the map to be returned
+ Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
+ //loop through the tracks
+ for (Track track : tracklist) {
+ Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ if (track instanceof SeedTrack) {
+ //create a seed canidate and get the hits
+ SeedCandidate seed = ((SeedTrack) track).getSeedCandidate();
+ //loop through the seedcandidate hits
+ for (HelicalTrackHit hit : seed.getHits()) {
+ //get a list of mc particles from the seed hits
+ List<MCParticle> mclist = hit.getMCParticles();
+ for (MCParticle mcp : mclist) {
+ //add one to the entries if the mcp is already in the map
+ if (mcmap.containsKey(mcp)) {
+ int entries = mcmap.get(mcp) + 1;
+ mcmap.put(mcp, entries);
+ }
+ //else just put 1 in the map
+ else
+ mcmap.put(mcp, 1);
+ }
+ }
+ MCParticle mcmatch = null;
+ int nmatch = 0;
+ //loop through the mcmap map and get the greatest int value within
+ for (Map.Entry<MCParticle, Integer> me : mcmap.entrySet()) {
+ if (me.getValue() > nmatch) {
+ nmatch = me.getValue();
+ mcmatch = me.getKey();
+ }
+ }
+ //if the track already is in the map...thats bad, otherwise add the mcmatch (MCParticle) key and its track
+ if (trkmap.containsKey(mcmatch)) System.out.println("more than one track associated with an MCParticle!!");
+ else trkmap.put(mcmatch, track);
+ }
+ }
+ //return the map
+ return trkmap;
+ }
+
+}
+
\ No newline at end of file
