lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -u -r1.1 -r1.2
--- AnalysisUtils.java 27 Aug 2008 17:58:15 -0000 1.1
+++ AnalysisUtils.java 11 Sep 2008 23:41:41 -0000 1.2
@@ -24,6 +24,7 @@
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
@@ -33,10 +34,10 @@
* @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;
+ //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() {
}
@@ -55,7 +56,7 @@
// If there are 0 hits, return an empty list
if (hitlist.size() == 0) return mcplist;
-
+
// Loop over the MCParticles in the first hit
for (MCParticle mcp : hitlist.get(0).getMCParticles()) {
@@ -99,7 +100,7 @@
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>();
+ Map<MCParticle,Integer> occurrences = new HashMap<MCParticle,Integer>();
double mostoccured;
int index;
//creates a seedcandiate for the track
@@ -110,13 +111,12 @@
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(occurrences.containsKey(part)){
+ index = occurrences.get(part);
+ occurrences.put(part,index+1);
+ } else {
+ occurrences.put(part, 1);
+ }
}
}
@@ -133,18 +133,18 @@
mostoccured = occured;
}
}
- purity = mostoccured / strk.getSeedCandidate().getHits().size();
- }
- }
- return purity;
- }
+ 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
+ * 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);
@@ -163,20 +163,20 @@
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;
+ 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()){
+ 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
@@ -185,76 +185,66 @@
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)).getSimTrackerHits().get(0).getTime();
- t2 = ((RawTrackerHit)h2.getRawHits().get(0)).getSimTrackerHits().get(0).getTime();
- } catch(NullPointerException npe) {
- t1 = h1.getTime();
- t2 = h2.getTime();
+ return Double.compare(SimTime(h1), SimTime(h2));
+ }
+ });
+
+ //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());
}
-
- return Double.compare(t1, t2);
+ //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);
}
- });
-
- //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;
+ }
+
+ //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
*/
@@ -265,7 +255,7 @@
Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
//loop through the tracks
for (Track track : tracklist) {
- Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ 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();
@@ -278,7 +268,7 @@
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);
@@ -302,5 +292,35 @@
return trkmap;
}
-}
-
\ No newline at end of file
+ private static double SimTime(HelicalTrackHit hit) {
+
+ // If there are no raw hits, we probably have virtual segmentation hits -> get time from hit
+ if (hit.getRawHits().size() == 0) return hit.getTime();
+
+ // We have raw tracker hits -> get time from SimTrackerHits
+ List<RawTrackerHit> rawlist = (List<RawTrackerHit>) hit.getRawHits();
+
+ // Default is to return a time of 0
+ double time = 0.;
+
+ // Loop over the RawTrackerHits
+ for (RawTrackerHit raw : rawlist) {
+ List<SimTrackerHit> simlist = raw.getSimTrackerHits();
+
+ // Find the average time of this raw hit
+ double rawtime = 0.;
+ if (simlist.size() > 0) {
+ for (SimTrackerHit simhit : simlist) {
+ rawtime += simhit.getTime();
+ }
+ rawtime /= simlist.size();
+ }
+
+ // Average the times of the raw hits
+ time += rawtime;
+ }
+ time /= rawlist.size();
+
+ return time;
+ }
+}
\ No newline at end of file