lcsim/sandbox/Partridge
diff -N RPAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RPAnalysisDriver.java 4 Sep 2008 22:50:44 -0000 1.1
@@ -0,0 +1,150 @@
+/*
+ * AnalysisDriver.java
+ *
+ * Created on February 11, 2008, 11:47 AM
+ *
+ */
+
+package org.lcsim.contrib.seedtracker.example;
+
+import hep.aida.IHistogram1D;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+
+import org.lcsim.constants.Constants;
+import org.lcsim.contrib.seedtracker.SeedTrack;
+import org.lcsim.contrib.seedtracker.SeedCandidate;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+
+
+/**
+ *
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class AnalysisDriver extends Driver {
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogram1D h1;
+ private IHistogram1D h2;
+
+ /** Creates a new instance of AnalysisDriver */
+ public AnalysisDriver() {
+// h1 = aida.histogramFactory().createHistogram1D("pT Efficiency", "", 100, 0., 50., "type=efficiency");
+ h2 = aida.histogramFactory().createHistogram1D("theta Efficiency", "", 90, 0., 180., "type=efficiency");
+
+ }
+
+ /**
+ * Process the current event
+ * @param event EventHeader for this event
+ */
+ public void process(EventHeader event) {
+ List<Track> tracklist = event.getTracks();
+ Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
+ for (Track track : tracklist) {
+ List<TrackerHit> hitlist = track.getTrackerHits();
+ Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ if (track instanceof SeedTrack) {
+ SeedCandidate seed = ((SeedTrack) track).getSeedCandidate();
+ int nhits = seed.getHits().size();
+ double theta = 180. * Math.acos(seed.getHelix().cth()) / Math.PI;
+ aida.cloud1D("Polar Angle for tracks with "+nhits+" hits").fill(180.*Math.acos(seed.getHelix().cth())/Math.PI);
+ if (nhits > 8) aida.cloud1D("Polar Angle for Tracks with 9 or more hits").fill(theta);
+ if (nhits < 9) aida.cloud1D("Polar Angle for Tracks with 8 or fewer hits").fill(theta);
+
+ for (HelicalTrackHit hit : seed.getHits()) {
+ double x = hit.getCorrectedPosition().x();
+ double y = hit.getCorrectedPosition().y();
+ HelicalTrackFit helix = seed.getHelix();
+ if (hit.getMCParticles().size() > 0) {
+ MCParticle particle = hit.getMCParticles().get(0);
+ Hep3Vector p = particle.getMomentum();
+ double phi = Math.atan2(p.y(), p.x());
+ double R = Math.sqrt(p.x()*p.x() + p.y()*p.y()) / (5. * Constants.fieldConversion);
+ double RS = R;
+ if (particle.getCharge() < 0.) RS = -1. * RS;
+ Hep3Vector start = particle.getOrigin();
+ double xc = start.x() + RS * Math.sin(phi);
+ double yc = start.y() - RS * Math.cos(phi);
+// double xc = helix.xc();
+// double yc = helix.yc();
+// double R = Math.abs(helix.R());
+ double RHit = Math.sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
+ double drphi_ms = 0.;
+ if (helix.ScatterMap().containsKey(hit)) drphi_ms = helix.ScatterMap().get(hit).drphi();
+ double drphi_res = hit.drphi();
+ double drphi = Math.sqrt(drphi_ms*drphi_ms + drphi_res*drphi_res);
+// aida.cloud1D("r-phi residual for layer "+hit.getLayerIdentifier()).fill(RHit-R);
+// aida.cloud1D("Pull for layer "+hit.getLayerIdentifier()).fill((RHit - R)/drphi);
+// aida.cloud1D("Hit resolution for layer "+hit.getLayerIdentifier()).fill(drphi_res);
+// aida.cloud2D("Hit MS resolution for layer "+hit.getLayerIdentifier()).fill(drphi_ms, helix.pT(5.));
+ }
+ List<MCParticle> mclist = hit.getMCParticles();
+ for (MCParticle mcp : mclist) {
+ if (mcmap.containsKey(mcp)) {
+ int entries = mcmap.get(mcp) + 1;
+ mcmap.put(mcp, entries);
+ } else
+ mcmap.put(mcp, 1);
+ }
+ }
+ MCParticle mcmatch = null;
+ int nmatch = 0;
+ for (Map.Entry<MCParticle, Integer> me : mcmap.entrySet()) {
+ if (me.getValue() > nmatch) {
+ nmatch = me.getValue();
+ mcmatch = me.getKey();
+ }
+ }
+ if (trkmap.containsKey(mcmatch)) System.out.println("more than one track associated with an MCParticle!!");
+ else trkmap.put(mcmatch, track);
+ }
+ }
+
+ List<MCParticle> mclist = event.getMCParticles();
+ for (MCParticle mcp : mclist) {
+ if (mcp.getCharge() == 0) continue;
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) continue;
+ double px = mcp.getPX();
+ double py = mcp.getPY();
+ double pz = mcp.getPZ();
+ double pt = Math.sqrt(px*px + py*py);
+ double p = Math.sqrt(pt*pt + pz*pz);
+ double cth = pz / p;
+ double theta = 180. * Math.acos(cth) / Math.PI;
+ double wgt = 0.;
+ if (trkmap.containsKey(mcp)) wgt = 1.;
+ if (pt > 1.1) {
+ aida.profile1D("Efficiency vs theta", 90, 0., 180.).fill(theta, wgt);
+ h2.fill(theta, wgt);
+ aida.histogram1D("MC angle", 90, 0., 180.).fill(theta);
+ }
+ if (Math.abs(cth) < 0.985) {
+ aida.profile1D("Efficiency vs pT", 100, 0., 50.).fill(pt, wgt);
+ aida.histogram1D("pT Efficiency", 100, 0., 50., "type=efficiency").fill(pt,wgt);
+ aida.histogram1D("MC pT", 100, 0., 50.).fill(pt);
+ }
+ if (!trkmap.containsKey(mcp)) {
+ List<HelicalTrackHit> hits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+ if (hits.size() > 6) {
+ System.out.println("Failed to find track. Found "+hits.size()+" hits @ theta = "
+ +180.*Math.acos(mcp.getMomentum().z() / mcp.getMomentum().magnitude())/Math.PI);
+ for (HelicalTrackHit hit : hits) {
+ System.out.println("Hit in "+hit.getLayerIdentifier()+" at x= "+hit.x()+" y= "+hit.y()+" z= "+hit.z());
+ }
+ }
+ }
+ }
+ return;
+ }
+}
\ No newline at end of file
lcsim/sandbox/Partridge
diff -N RPAnalysisUtils.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RPAnalysisUtils.java 4 Sep 2008 22:50:44 -0000 1.1
@@ -0,0 +1,388 @@
+/*
+ * AnalysisUtils.java
+ *
+ * Created on August 4, 2008, 11:00 AM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.seedtracker.analysis;
+
+import hep.physics.vec.BasicHep3Vector;
+import java.util.ArrayList;
+<<<<<<< AnalysisUtils.java
+import java.util.HashMap;
+=======
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.HashSet;
+>>>>>>> 1.2
+import java.util.List;
+import java.util.Map;
+import org.lcsim.contrib.seedtracker.SeedLayer;
+import org.lcsim.contrib.seedtracker.SeedLayer.SeedType;
+import org.lcsim.contrib.seedtracker.SeedStrategy;
+import org.lcsim.event.EventHeader;
+
+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
+ * @authors Richard Partridge,Pelham Keahey
+ * @version 1.0
+ */
+public class AnalysisUtils {
+<<<<<<< AnalysisUtils.java
+ private static String hitname = "HelicalTrackHits";
+
+=======
+ //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;
+>>>>>>> 1.2
+ /** Creates a new instance of AnalysisUtils */
+ public AnalysisUtils() {
+ }
+
+ /**
+ * Finds those MCParticles that are associated with every hit in a list of
+ * HelicalTrackHits. A list of size 0 is returned if there are no MCParticles
+ * that are associated with all hits.
+ * @param hitlist List of hits to be checked
+ * @return List of MCParticles that are associated with all hits
+ */
+ public static List<MCParticle> FindMCParticleInAllHits(List<HelicalTrackHit> hitlist) {
+
+ // Create a list for MCParticles that are in all hits
+ List<MCParticle> mcplist = new ArrayList<MCParticle>();
+
+ // 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()) {
+
+ // The boolean good will remain true as long as the MCParticle is in each hit checked
+ boolean good = true;
+ // Loop over the rest of the hits in the input list
+ for (int i = 1; i < hitlist.size(); i++) {
+
+ // The boolean match will be set true if we find the matching MCParticle in the current hit
+ boolean match = false;
+ // Loop over MCParticles in the current hit
+ for (MCParticle mcp2 : hitlist.get(i).getMCParticles()) {
+ // Check for a match
+ if (mcp.equals(mcp2)) {
+ // Found a match - signal it and we can quit checking
+ match = true;
+ break;
+ }
+ }
+
+ // If this hit did not have a match, the MCParticle being checked is not in all hits
+ if (match == false) {
+ good = false;
+ break;
+ }
+ }
+
+ // If we finish looping through all hits and find the MCParticle being checked is in all
+ // hits, add it to the list of MCParticles that are in all hits
+ if (good) mcplist.add(mcp);
+ }
+
+ // Done looping over all MCParticle candidates - return the list with any matches that were found
+ return mcplist;
+ }
+<<<<<<< AnalysisUtils.java
+
+ public static List<MCParticle> FindableMCParticles(EventHeader event, List<SeedStrategy> strategylist) {
+
+ // Create a list of findable MCParticles
+ List<MCParticle> mclist = new ArrayList<MCParticle>();
+
+ // Get the helical track hits from the event
+ List<HelicalTrackHit> hitcol = (List<HelicalTrackHit>) event.get(hitname);
+
+ // Create a map between the MCParticles and a list of hits containing that MCParticle
+ Map<MCParticle, List<HelicalTrackHit>> mcmap = new HashMap<MCParticle, List<HelicalTrackHit>>();
+
+ // Loop over the hits and the MCParticles in those hits
+ for (HelicalTrackHit hit : hitcol) {
+ for (MCParticle mcp : hit.getMCParticles()) {
+
+ // If this is the first time we come across this MCParticle, create a new list of hits
+ if (!mcmap.containsKey(mcp)) mcmap.put(mcp, new ArrayList<HelicalTrackHit>());
+ // Add this hit to the list
+ mcmap.get(mcp).add(hit);
+ }
+ }
+
+ // Loop over the MCParticles that have one or more HelicalTrackHit
+ for (MCParticle mcp : mcmap.keySet()) {
+
+ // Loop over strategies until we come across a findable strategy
+ boolean findable = false;
+ for (SeedStrategy strategy : strategylist) {
+ if (CheckHits(strategy, mcmap.get(mcp))) {
+
+ }
+ }
+ }
+
+ private static boolean CheckHits(SeedStrategy strategy, List<HelicalTrackHit> hitcol) {
+ int nseed = 0;
+ int nconfirm = 0;
+ int nextend = 0;
+ boolean goodstrategy = true;
+ List<SeedLayer> lyrlist = strategy.getLayerList();
+ for (SeedLayer lyr : lyrlist) {
+ boolean match = false;
+ for (HelicalTrackHit hit : hitcol) {
+ if (MatchLayers(hit, lyr)) {
+ match = true;
+ break;
+ }
+ }
+ if (match) {
+ SeedType type = lyr.getType();
+ if (type.equals(SeedType.Seed)) nseed++;
+ if (type.equals(SeedType.Confirm)) nconfirm++;
+ if (type.equals(SeedType.Extend)) nextend++;
+ }
+ }
+ int ntot = nseed + nconfirm + nextend;
+ return nseed == 3 && nconfirm >= strategy.getMinConfirm() && ntot >= strategy.getMinHits();
+ }
+
+ }
+ private static boolean MatchLayers(HelicalTrackHit hit, SeedLayer lyr) {
+ if (!hit.Detector().equals(lyr.getDetName())) return false;
+ if (hit.Layer() != lyr.getLayer()) return false;
+ return hit.BarrelEndcapFlag().equals(lyr.getBarrelEndcapFlag());
+ }
+}=======
+ /**
+ * 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;
+ }
+
+}
+ >>>>>>> 1.2
