8 added files
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N AnalysisUtils.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AnalysisUtils.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,306 @@
+/*
+ * 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.recon.tracking.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.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.event.Track;
+import org.lcsim.recon.tracking.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 {
+ //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() {
+ }
+
+ /**
+ * 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;
+ }
+ /**
+ * 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)).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(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
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N EfficiencyHistograms.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EfficiencyHistograms.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,95 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.analysis;
+
+import hep.aida.IHistogram1D;
+import java.util.Map;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author Pelham Keahey
+ */
+public class EfficiencyHistograms {
+ private AIDA aida = AIDA.defaultInstance();
+ private AnalysisUtils util = new AnalysisUtils();
+ //map retrived from a util method that gives you a map of mcparticles with an integer, the integer is either
+ //0: no track found but should have been
+ //1: track found and meets cuts
+ private Map<MCParticle,Integer> effimap;
+ //used to cut down on time spent on recalcuting variables when the event hasn't even changed (variables affected by event)
+ private EventHeader currentEvent=null;
+ //histograms
+ IHistogram1D thetaeff,pteff,dcaeff,z0eff,phi0eff,tanLeff,omegaeff;
+
+ public EfficiencyHistograms(){
+ //set dir and create histogram shells
+ //histograms to be displayed
+ aida.tree().mkdir("Efficiency");
+ thetaeff = aida.histogramFactory().createHistogram1D("Efficiency/Theta","", 90, 0., 180., "type=efficiency");
+ pteff = aida.histogramFactory().createHistogram1D("Efficiency/Pt","", 100, 0., 50., "type=efficiency");
+ aida.tree().mkdir("Efficiency/HelixParam");
+ dcaeff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/DCA","", 120, -4., 4., "type=efficiency");
+ z0eff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/z0","", 120, -4., 4., "type=efficiency");
+ phi0eff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/phi0","", 120, 0., 360., "type=efficiency");
+ tanLeff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/tanL","", 120, -10., 10., "type=efficiency");
+ omegaeff = aida.histogramFactory().createHistogram1D("Efficiency/HelixParam/Omega","", 120, -.02, .02, "type=efficiency");
+
+ }
+
+ public void drawEfficiencyHistograms(MCParticle mcp,EventHeader event){
+ //calculator used to get helix parameters given an MCParticle and a BField
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
+ //method to prevent recalculating variables based on the event, when the event is still the same
+ if(currentEvent==null || !currentEvent.equals(event)) {
+ effimap = util.FindEfficiency(event);
+ currentEvent = event;
+ }
+ //helix parameters from the calculator above (used for cuts later)
+ double theta = calc.getTheta();
+ double degreetheta = 180*theta/Math.PI;
+ double pt = calc.getMCTransverseMomentum();
+ double dca = calc.getDCA();
+ double z0 = calc.getZ0();
+ //check to see if the efficiency map has the mcparicle and get the associated integer (1 or 0)
+ if(effimap.containsKey(mcp)){
+ int wgt = effimap.get(mcp);
+ //cuts for tanL: pt,dca and z0
+ if(pt>1.1 && Math.abs(dca) < 9 && Math.abs(z0) < 9){
+ thetaeff.fill(degreetheta, wgt);
+ tanLeff.fill(calc.getSlopeSZPlane(), wgt);
+ aida.histogram1D("Efficiency/MCType/tanL/"+mcp.getType().getName(),120, -10., 10., "type=efficiency").fill(calc.getSlopeSZPlane(), wgt);
+ }
+ //Pt Cuts: Cos(theta) DCA z0
+ if (Math.cos(theta) < 0.985 && Math.abs(dca) < 9 && Math.abs(z0) < 9) {
+ pteff.fill(pt,wgt);
+ aida.histogram1D("Efficiency/MCType/Pt/"+mcp.getType().getName(), 100, 0., 50., "type=efficiency").fill(pt, wgt);
+ }
+ //Omega Cut: Pt
+ if(pt>1.1){
+ omegaeff.fill(calc.getMCOmega(), wgt);
+ aida.histogram1D("Efficiency/MCType/Omega/"+mcp.getType().getName(), 120, -.02, .02, "type=efficiency").fill(calc.getMCOmega(), wgt);
+ }
+ //Z0 Cuts: Pt DCA Cos(theta)
+ if(pt>1.1 && Math.abs(dca) < 9 && Math.cos(theta) < 0.985){
+ z0eff.fill(calc.getZ0(), wgt);
+ aida.histogram1D("Efficiency/MCType/z0/"+mcp.getType().getName(), 120, -4., 4., "type=efficiency").fill(calc.getZ0(), wgt);
+ }
+ //DCA Cuts: Pt Cos(theta) z0
+ if(pt>1.1 && Math.cos(theta) < 0.985 && Math.abs(z0) < 9){
+ dcaeff.fill(calc.getDCA(), wgt);
+ aida.histogram1D("Efficiency/MCType/DCA/"+mcp.getType().getName(), 120, -4., 4., "type=efficiency").fill(calc.getDCA(), wgt);
+ }
+ //Phi0 Cuts: Pt Cos(theta) z0 DCA
+ if(pt>1.1 && Math.cos(theta) < 0.985 && Math.abs(z0) < 9 && Math.abs(dca) < 9){
+ phi0eff.fill(180*calc.getPhi0()/Math.PI, wgt);
+ aida.histogram1D("Efficiency/MCType/phi0/"+mcp.getType().getName(), 120, 0., 360., "type=efficiency").fill(calc.getDCA(), wgt);
+ }
+ }
+ }
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N HelixParamCalculator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixParamCalculator.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,205 @@
+/*
+ * HelixParameterCalculator.java
+ *
+ * Created on July 7th, 2008, 11:09 AM
+ *
+ *
+ */
+package org.lcsim.recon.tracking.seedtracker.analysis;
+
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.EventHeader;
+
+/**
+ * Class used for calculating MC particle track paramters
+ * @author Pelham Keahey
+ *
+ */
+public class HelixParamCalculator {
+
+ MCParticle mcp;
+ private double R,BField,theta,arclength;
+ //Some varibles are usesd in other calculations, thus they are global
+ /*
+ xc, yc --coordinates
+ mcdca -- MC distance of closest approach
+ mcphi0 --azimuthal angle
+ tanL -- Slope SZ plane ds/dz
+ x0,y0 are the position of the particle at the dca
+ */
+ private double xc,yc,mcdca,mcphi0,tanL;
+ private double x0,y0;
+ /**
+ * Constructor that is fed a magnetic field and MCPARTICLE
+ * @param mcpc
+ * @param cBField
+ */
+ public HelixParamCalculator(MCParticle mcpc,double cBField)
+ {
+ //mc and event global varibles used for calculation
+ mcp=mcpc;
+
+ //Returns the MagneticField at point (0,0,0), assumes constant magfield
+ BField = cBField;
+
+ //Calculate theta, the of the helix projected into an SZ plane, from the z axis
+ 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;
+ theta = Math.acos(cth);
+
+ //Calculate Radius of the Helix
+ R = ((mcp.getCharge())*(mcp.getMomentum().magnitude()*Math.sin(theta))/(.0003*BField));
+
+ //Slope in the Dz/Ds sense, tanL Calculation
+ tanL = mcp.getPZ()/(Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()));
+
+ //Distance of closest approach Calculation
+ xc = mcp.getOriginX() + R * Math.sin(Math.atan2(mcp.getPY(),mcp.getPX()));
+ yc = mcp.getOriginY() - R * Math.cos(Math.atan2(mcp.getPY(),mcp.getPX()));
+ double xcyc = Math.sqrt(xc*xc + yc*yc);
+ if(mcp.getCharge()>0)
+ {
+ mcdca = R - xcyc;
+ }
+ else
+ {
+ mcdca = R + xcyc;
+ }
+
+
+ //azimuthal calculation of the momentum at the DCA, phi0, Calculation
+ mcphi0 = Math.atan2(xc/(R-mcdca), -yc/(R-mcdca));
+ if(mcphi0<0)
+ {
+ mcphi0 += 2*Math.PI;
+ }
+ //z0 Calculation, z position of the particle at dca
+ x0 = -mcdca*Math.sin(mcphi0);
+ y0 = mcdca*Math.sin(mcphi0);
+ arclength = (((mcp.getOriginX()-x0)*Math.cos(mcphi0))+((mcp.getOriginY()-y0)*Math.sin(mcphi0)));
+
+ }
+ /**
+ * Calculates the B-Field from event
+ * @param mcpc
+ * @param eventc
+ */
+ public HelixParamCalculator(MCParticle mcpc,EventHeader eventc)
+ {
+ this(mcpc,eventc.getDetector().getFieldMap().getField(new BasicHep3Vector(0.,0.,0.)).z());
+ }
+ /**
+ * Return the magneticfield at point 0,0,0
+ * @return double BField
+ */
+ public double getMagField()
+ {
+ return BField;
+ }
+ /**
+ * Return the radius of the Helix track
+ * @return double R
+ */
+ public double getRadius()
+ {
+ return R;
+ }
+ /**
+ * Return the theta angle for the projection of the helix in the SZ plane
+ * from the z axis
+ * @return double theta
+ */
+ public double getTheta()
+ {
+ return theta;
+ }
+ /**
+ * Return the particle's momentum
+ * @return double mcp momentum
+ */
+ public double getMCMomentum()
+ {
+ return mcp.getMomentum().magnitude();
+ }
+ /**
+ * Return the curvature (omega)
+ * @return double omega
+ */
+ public double getMCOmega()
+ {
+ return mcp.getCharge()/((mcp.getMomentum().magnitude()*Math.sin(theta))/(.0003*BField));
+ }
+ /**
+ * Return the transvers momentum of the MC particle, Pt
+ * @return double Pt
+ */
+ public double getMCTransverseMomentum()
+ {
+ return (mcp.getMomentum().magnitude())*Math.sin(theta);
+ }
+ /**
+ * Return the slope of the helix in the SZ plane, tan(lambda)
+ * @return double tanL
+ */
+ public double getSlopeSZPlane()
+ {
+ return tanL;
+ }
+ /**
+ * Return the distance of closest approach
+ * @return double mcdca
+ */
+ public double getDCA()
+ {
+ return mcdca;
+ }
+ /**
+ * Return the azimuthal angle of the momentum when at the position of closest approach
+ * @return double mcphi0
+ */
+ public double getPhi0()
+ {
+ return mcphi0;
+ }
+ /**
+ * Return the z position at the distance of closest approach
+ * @return double z0 position
+ */
+ public double getZ0()
+ {
+ double x0 = -mcdca*Math.sin(mcphi0);
+ double y0 = mcdca*Math.sin(mcphi0);
+ double s = (((mcp.getOriginX()-x0)*Math.cos(mcphi0))+((mcp.getOriginY()-y0)*Math.sin(mcphi0)));
+ return mcp.getOriginZ()-(s*tanL);
+ }
+ /**
+ * Return the arclength of the helix from the ORIGIN TO THE DCA
+ * @return double arclength
+ */
+ public double getArcLength()
+ {
+ return arclength;
+ }
+ /**
+ * Return the x position of the particle when at the dca
+ * @return double arclength
+ */
+ public double getX0()
+ {
+ return x0;
+ }
+ /**
+ * Return the y position of the particle at the dca
+ * @return double arclength
+ */
+ public double getY0()
+ {
+ return y0;
+ }
+
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N HelixParamHistograms.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixParamHistograms.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,177 @@
+package org.lcsim.recon.tracking.seedtracker.analysis;
+
+
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.recon.tracking.seedtracker.analysis.HelixParamCalculator;
+
+/**
+ *
+ * @author Pelham Keahey
+ * Thursday, July 3rd 2008
+ */
+public class HelixParamHistograms {
+ //helix from mc partcle for layer hit 1st
+ //simtracker hit where is nrg
+ //pos from tracker hit 2nd
+ //compare fitted helix to mcparticle hits r phi z
+ private AIDA aida = AIDA.defaultInstance();
+ //Histrogram Variables
+ private double datom,dtl,ddca,dphi0,dz,dm,dtm;
+ private double mcom,mctl,mcdca,mcphi0,mcz,mcm,mctm;
+ //Error Variables
+ private double curverror,sloperror,dcaerror,phi0error,z0error;
+
+ //Constructor for Histrogram Drawer, Calculations done inside as well
+ public HelixParamHistograms(HelicalTrackFit fit,MCParticle mcp,EventHeader event)
+ {
+
+
+ HelixParamCalculator calc = new HelixParamCalculator(mcp,event);
+ //************ Data Calculations for MOMENTUM, OMEGA, TRANS. MOMENTUM *********
+ dm = fit.p(calc.getMagField());
+ datom = fit.curvature();
+ dtm = fit.pT(calc.getMagField());
+
+
+ //************* MCParticles Calculations for MOMENTUM, OMEGA, TRANS. MOMENTUM ******
+ mcm = calc.getMCMomentum();
+ mcom = calc.getMCOmega();
+ mctm = calc.getMCTransverseMomentum();
+
+
+ //************ Tan Lamda Calculations tanL tsnL tan L *************
+ dtl = fit.slope();
+ mctl = calc.getSlopeSZPlane();
+
+
+ //************ DCA Calculations DCA DCA DCA DCA DCA DCA *************
+ ddca = fit.dca();
+ mcdca = calc.getDCA();
+
+
+ //******* phi0 Calculation phi0 phi0 phi0 phi0 ******************
+ dphi0 = fit.phi0();
+ mcphi0 = calc.getPhi0();
+
+
+ //******* z0 Calculations z0 z0 z0 z0 z0 z0 z0 ************************
+ dz = fit.z0();
+ mcz= calc.getZ0();
+
+
+ //Get error for calculating the pull
+ curverror = fit.getCurveError();
+ sloperror = fit.getSlopeError();
+ dcaerror = fit.getDcaError();
+ phi0error = fit.getPhi0Error();
+ z0error = fit.getZ0Error();
+
+
+ }
+ public void DrawDataPositive()
+ {
+
+ aida.histogram1D("HelixParam/1Dat-Positive/Omega", 300, -1*Math.pow(10, -3), Math.pow(10,-3)).fill(datom);
+ aida.histogram1D("HelixParam/1Dat-Positive/tanL", 200, -8., 8.).fill(dtl);
+ aida.histogram1D("HelixParam/1Dat-Positive/DCA", 200, -.1, .1).fill(ddca);
+ aida.histogram1D("HelixParam/1Dat-Positive/phi0", 200, -1., 3*Math.PI).fill(dphi0);
+ aida.histogram1D("HelixParam/1Dat-Positive/z0" , 200, -.1, .1).fill(dz);
+ aida.histogram1D("HelixParam/1Dat-Positive/Momentum", 200, 0., 200.).fill(dm);
+ aida.histogram1D("HelixParam/1Dat-Positive/Trans.Momentum", 200, 0., 150.).fill(dtm);
+ }
+ public void DrawDataNegative()
+ {
+ aida.histogram1D("HelixParam/2Dat-Negative/Omega", 300, -1*Math.pow(10, -3), Math.pow(10,-3)).fill(datom);
+ aida.histogram1D("HelixParam/2Dat-Negative/tanL", 200, -8., 8.).fill(dtl);
+ aida.histogram1D("HelixParam/2Dat-Negative/DCA", 200, -.1, .1).fill(ddca);
+ aida.histogram1D("HelixParam/2Dat-Negative/phi0", 200, -1., 3*Math.PI).fill(dphi0);
+ aida.histogram1D("HelixParam/2Dat-Negative/z0" , 200, -.1, .1).fill(dz);
+ aida.histogram1D("HelixParam/2Dat-Negative/Momentum", 200, 0., 200.).fill(dm);
+ aida.histogram1D("HelixParam/2Dat-Negative/Trans.Momentum", 200, 0., 150.).fill(dtm);
+ }
+ public void DrawMCPositive()
+ {
+
+ aida.histogram1D("HelixParam/3MC-Positive/Omega", 300, -1*Math.pow(10, -3), Math.pow(10,-3)).fill(mcom);
+ aida.histogram1D("HelixParam/3MC-Positive/tanL", 200, -8., 8.).fill(mctl);
+ aida.histogram1D("HelixParam/3MC-Positive/DCA", 200, -.01, .01).fill(mcdca);
+ aida.histogram1D("HelixParam/3MC-Positive/phi0", 200, -1., 3*Math.PI).fill(mcphi0);
+ aida.histogram1D("HelixParam/3MC-Positive/z0" , 200, -.001, .001).fill(mcz);
+ aida.histogram1D("HelixParam/3MC-Positive/Momentum", 200, 0., 200.).fill(mcm);
+ aida.histogram1D("HelixParam/3MC-Positive/Trans.Momentum", 200, 0., 150.).fill(mctm);
+ }
+ public void DrawMCNegative()
+ {
+ aida.histogram1D("HelixParam/4MC-Negative/Omega", 300, -1*Math.pow(10, -3), Math.pow(10,-3)).fill(mcom);
+ aida.histogram1D("HelixParam/4MC-Negative/tanL", 200, -8., 8.).fill(mctl);
+ aida.histogram1D("HelixParam/4MC-Negative/DCA", 200, -.01, .01).fill(mcdca);
+ aida.histogram1D("HelixParam/4MC-Negative/phi0", 200, -1., 3*Math.PI).fill(mcphi0);
+ aida.histogram1D("HelixParam/4MC-Negative/z0" , 200, -.001, .001).fill(mcz);
+ aida.histogram1D("HelixParam/4MC-Negative/Momentum", 200, 0., 200.).fill(mcm);
+ aida.histogram1D("HelixParam/4MC-Negative/Trans.Momentum", 200, 0., 150.).fill(mctm);
+
+ }
+ public void DrawResidualPositive()
+ {
+ aida.histogram1D("HelixParam/Residual-Positive Charges/Momentum" , 200, -10., 10.).fill(dm-mcm);
+ aida.histogram1D("HelixParam/Residual-Positive Charges/Trans.Momentum", 200, -5., 5.).fill(dtm-mctm);
+ aida.histogram1D("HelixParam/Residual-Positive Charges/Omega" , 300, -5*Math.pow(10, -6), 5*Math.pow(10,-6)).fill(datom-mcom);
+ aida.histogram1D("HelixParam/Residual-Positive Charges/tanL" , 300, -.01, .01).fill(dtl-mctl);
+ aida.histogram1D("HelixParam/Residual-Positive Charges/DCA" , 200, -.1, .1).fill(ddca-mcdca);
+ aida.histogram1D("HelixParam/Residual-Positive Charges/phi0" , 200, -.002, .002).fill(dphi0-mcphi0);
+ aida.histogram1D("HelixParam/Residual-Positive Charges/z0" , 200, -.1, .1).fill(dz-mcz);
+
+ }
+ public void DrawResidualNegative()
+ {
+ aida.histogram1D("HelixParam/Residual-Negative Charges/Momentum" , 200, -10., 10.).fill(dm-mcm);
+ aida.histogram1D("HelixParam/Residual-Negative Charges/Trans.Momentum", 200, -5., 5.).fill(dtm-mctm);
+ aida.histogram1D("HelixParam/Residual-Negative Charges/Omega" , 300, -5*Math.pow(10, -6), 5*Math.pow(10,-6)).fill(datom-mcom);
+ aida.histogram1D("HelixParam/Residual-Negative Charges/tanL" , 300, -.01, .01).fill(dtl-mctl);
+ aida.histogram1D("HelixParam/Residual-Negative Charges/DCA" , 200, -.1, .1).fill(ddca-mcdca);
+ aida.histogram1D("HelixParam/Residual-Negative Charges/phi0" , 200, -.002, .002).fill(dphi0-mcphi0);
+ aida.histogram1D("HelixParam/Residual-Negative Charges/z0" , 200, -.1, .1).fill(dz-mcz);
+
+ }
+ public void DrawResidual()
+ {
+ aida.histogram1D("HelixParam/Residual-both/Momentum" , 200, -10., 10.).fill(dm-mcm);
+ aida.histogram1D("HelixParam/Residual-both/Trans.Momentum", 200, -5., 5.).fill(dtm-mctm);
+ aida.histogram1D("HelixParam/Residual-both/Omega" , 300, -5*Math.pow(10, -6), 5*Math.pow(10,-6)).fill(datom-mcom);
+ aida.histogram1D("HelixParam/Residual-both/tanL" , 300, -.01, .01).fill(dtl-mctl);
+ aida.histogram1D("HelixParam/Residual-both/DCA" , 200, -.1, .1).fill(ddca-mcdca);
+ aida.histogram1D("HelixParam/Residual-both/phi0" , 200, -.002, .002).fill(dphi0-mcphi0);
+ aida.histogram1D("HelixParam/Residual-both/z0" , 200, -.1, .1).fill(dz-mcz);
+
+ }
+ public void DrawPullPositive()
+ {
+ aida.histogram1D("HelixParam/Pull-Positive Charges/Omega", 300, -10, 10).fill((datom-mcom)/curverror);
+ aida.histogram1D("HelixParam/Pull-Positive Charges/tanL" , 300, -11, 11).fill((dtl-mctl)/sloperror);
+ aida.histogram1D("HelixParam/Pull-Positive Charges/dca" , 300, -10, 10).fill((ddca-mcdca)/dcaerror);
+ aida.histogram1D("HelixParam/Pull-Positive Charges/phi0" , 300, -11, 11).fill((dphi0-mcphi0)/phi0error);
+ aida.histogram1D("HelixParam/Pull-Positive Charges/z0" , 300, -50, 50).fill((dz-mcz)/z0error);
+
+ }
+ public void DrawPullNegative()
+ {
+ aida.histogram1D("HelixParam/Pull-Negative Charges/Omega", 300, -10, 10).fill((datom-mcom)/curverror);
+ aida.histogram1D("HelixParam/Pull-Negative Charges/tanL" , 300, -11, 11).fill((dtl-mctl)/sloperror);
+ aida.histogram1D("HelixParam/Pull-Negative Charges/dca" , 300, -10, 10).fill((ddca-mcdca)/dcaerror);
+ aida.histogram1D("HelixParam/Pull-Negative Charges/phi0" , 300, -11, 11).fill((dphi0-mcphi0)/phi0error);
+ aida.histogram1D("HelixParam/Pull-Negative Charges/z0" , 300, -50, 50).fill((dz-mcz)/z0error);
+ }
+ public void DrawPull()
+ {
+ aida.histogram1D("HelixParam/Pull-Both/Omega", 300, -10, 10).fill((datom-mcom)/curverror);
+ aida.histogram1D("HelixParam/Pull-Both/tanL" , 300, -11, 11).fill((dtl-mctl)/sloperror);
+ aida.histogram1D("HelixParam/Pull-Both/dca " , 300, -10, 10).fill((ddca-mcdca)/dcaerror);
+ aida.histogram1D("HelixParam/Pull-Both/phi0" , 300, -11, 11).fill((dphi0-mcphi0)/phi0error);
+ aida.histogram1D("HelixParam/Pull-Both/z0 " , 300, -50, 50).fill((dz-mcz)/z0error);
+ }
+
+
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N HistogramAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HistogramAnalysisDriver.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,143 @@
+/*
+ * AnalysisDriver.java
+ *
+ * Created on February 11, 2008, 11:47 AM
+ *
+ */
+
+package org.lcsim.recon.tracking.seedtracker.analysis;
+
+
+import org.lcsim.recon.tracking.seedtracker.analysis.TrackHitHistograms;
+import org.lcsim.recon.tracking.seedtracker.analysis.HelixParamHistograms;
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.analysis.AnalysisUtils;
+import org.lcsim.recon.tracking.seedtracker.analysis.EfficiencyHistograms;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.recon.tracking.seedtracker.analysis.HelixParamCalculator;
+import org.lcsim.recon.tracking.seedtracker.analysis.PurityHistograms;
+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 HistogramAnalysisDriver extends Driver {
+ private AIDA aida = AIDA.defaultInstance();
+ private double b=0,c=0,d=0;
+ /** Creates a new instance of AnalysisDriver */
+ public HistogramAnalysisDriver() {
+
+ }
+
+ /**
+ * Process the current event
+ * @param event EventHeader for this event
+ */
+ PurityHistograms draw3 = new PurityHistograms();
+ EfficiencyHistograms draw4 = new EfficiencyHistograms();
+ public void process(EventHeader event){
+
+
+ List<Track> tracklist = event.getTracks();
+ Map<MCParticle, Track> trkmap = new HashMap<MCParticle, Track>();
+
+ for (Track track : tracklist) {
+ Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ if (track instanceof SeedTrack) {
+ SeedCandidate seed = ((SeedTrack) track).getSeedCandidate();
+ for (HelicalTrackHit hit : seed.getHits()) {
+
+
+
+
+ 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;
+
+ if (trkmap.containsKey(mcp)) {
+ Track trk = trkmap.get(mcp);
+ if (trk instanceof SeedTrack) {
+ SeedTrack strk = (SeedTrack) trk;
+ HelicalTrackFit fit = strk.getSeedCandidate().getHelix();
+
+ //Class to draw histrogram
+ HelixParamHistograms draw = new HelixParamHistograms(fit,mcp,event);
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
+ AnalysisUtils util = new AnalysisUtils();
+
+ //---Helix Parameter Histrograms------------Helix Parameter Histrograms---\\
+ if(mcp.getCharge()>0)
+ {
+ draw.DrawResidualPositive();
+ draw.DrawPullPositive();
+ }
+ else
+ {
+ draw.DrawResidualNegative();
+ draw.DrawPullNegative();
+ }
+ draw.DrawResidual();
+ draw.DrawPull();
+
+
+ TrackHitHistograms draw2 = new TrackHitHistograms(fit, mcp,event);
+
+ //---MCParticle vs HelicalTrackFit----------------------------MCParticle vs HelicalTrackFit---
+ draw2.drawMCParticlevHelicalTrackHit();
+
+ //---SimTracker vs HelicalTrackFit---------------------------SimTracker vs HelicalTrackFit---
+ draw2.drawSimTrackervHelicalTrack();
+
+ //---Purity---------------------------Purity---
+ draw3.drawPurityHistograms(mcp, trk);
+ double a = util.FindTrackPurity(trk);
+ c++;
+ d += Math.sqrt((a*(1-a))/c);
+ if(a==1)b++;
+
+
+ }
+ }
+ //--Efficiency---------------------------Efficiency---
+ draw4.drawEfficiencyHistograms(mcp,event);
+ }
+ System.out.println((b/c)*100+"+-"+(d/c)*100+"%");
+ }
+}
+
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N PurityHistograms.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PurityHistograms.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,54 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.analysis;
+
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author Pelham Keahey
+ */
+public class PurityHistograms {
+ private AIDA aida = AIDA.defaultInstance();
+ private AnalysisUtils util = new AnalysisUtils();
+ /**
+ *
+ * @param mcp MCParticle currently being "looked" at
+ * @param trk the mcp's track
+ */
+ public void drawPurityHistograms(MCParticle mcp,Track trk){
+ //used to calculate helix parameters given a MCParticle and BField
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, 5);
+ if (trk instanceof SeedTrack) {
+ //create a seed track
+ SeedTrack strk = (SeedTrack) trk;
+ //get total size of the hits used to create the track
+ int hitsize = strk.getSeedCandidate().getHits().size();
+ //get strat. name (used for histograms to analyze purity by strat.)
+ String stratname = strk.getStrategy().getName();
+ //find the purity of the track
+ double purity = util.FindTrackPurity(trk);
+ //shows the fake tracks (fake track is a track with a purity less than .5)
+ if(purity<.5)aida.histogram1D("Purity/Purity<.5 only", 200, 0, 1.1).fill(purity);
+ //1D histograms
+ aida.histogram1D("Purity/ByStrategy/Strat: "+stratname, 200,0,1.1).fill(purity);
+ aida.histogram1D("Purity/Purity", 200, 0, 1.1).fill(purity);
+ //only shows the purity when it is not perfect
+ if(purity!=1)aida.histogram1D("Purity/Purity<1 only", 200, 0, 1.1).fill(purity);
+ aida.histogram1D("Purity/ByHits/Purity to hits: "+hitsize, 300, 0, 1.1).fill(purity);
+ //2D histograms
+ aida.histogram1D("Purity/Hits v Purity", 200, 0, 20,"type=efficiency").fill(hitsize,purity);
+ aida.cloud2D("Purity/Hits v. Purity (2d)").setTitle("");
+ aida.cloud2D("Purity/Hits v. Purity (2d)").fill(hitsize,purity);
+ aida.histogram1D("Purity/Hits v 1-Purity", 200, 0, 20,"type=efficiency").fill(hitsize,1-purity);
+ aida.histogram1D("Purity/Theta v Purity", 100, 10, 100, "type=efficiency").fill(180*calc.getTheta()/Math.PI,purity);
+ aida.histogram1D("Purity/Pt v Purity", 200, 0, 70,"type=efficiency").fill(calc.getMCTransverseMomentum(),purity);
+ }
+ }
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N SeedTrackDiag.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SeedTrackDiag.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,135 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.analysis;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Iterator;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedLayer;
+import org.lcsim.recon.tracking.seedtracker.SeedLayer.SeedType;
+import org.lcsim.recon.tracking.seedtracker.analysis.AnalysisUtils;
+import org.lcsim.recon.tracking.seedtracker.diagnostic.EmptySeedTrackerDiagnostics;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.line.SlopeInterceptLineFit;
+import org.lcsim.fit.zsegment.ZSegmentFit;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author Pelham Keahey
+ */
+public class SeedTrackDiag extends EmptySeedTrackerDiagnostics{
+ AIDA aida = AIDA.defaultInstance();
+ //List of MCParticles
+ List<MCParticle> mcheck = new ArrayList<MCParticle>();
+ AnalysisUtils util = new AnalysisUtils();
+ @Override
+ public void setEvent(EventHeader event){
+ super.setEvent(event);
+
+
+ }
+
+ @Override
+ /*
+ * Will print out all the information about tracks in the event that should have been found ,but weren't that met the cuts.
+ */
+ public void fireFitterFitMade(SeedCandidate seed, HelicalTrackFit hfit, CircleFit cfit, SlopeInterceptLineFit lfit, ZSegmentFit zfit, FitStatus fitStatus, boolean success) {
+
+ List<MCParticle> temp = new ArrayList<MCParticle>();
+ if(fitStatus != FitStatus.Success )return;
+ if (success) return;
+ Set<MCParticle> mcps = new HashSet<MCParticle>();
+ Map<MCParticle,List<HelicalTrackHit>> hitmcmap = new HashMap<MCParticle,List<HelicalTrackHit>>();
+ List<String> idlist = new ArrayList<String>();
+ for (SeedLayer l : seed.getSeedStrategy().getLayers(SeedType.Seed)){
+ idlist.add(l.getDetName()+l.getLayer()+l.getBarrelEndcapFlag().toString());
+ }
+ ArrayList<HelicalTrackHit> seedlayerhits = new ArrayList<HelicalTrackHit>();
+ if(seedlayerhits.size()>0)System.out.println("Oh Snap");
+ for(HelicalTrackHit hito : seed.getHits()){
+ mcps.addAll(hito.getMCParticles());
+ if(idlist.contains(hito.getLayerIdentifier())){
+ seedlayerhits.add(hito);
+ for(MCParticle mcp : hito.getMCParticles()){
+ hitmcmap.put(mcp, seedlayerhits);
+ temp.add(mcp);
+ }
+ }
+ }
+ for (HelicalTrackHit h : seed.getHits()){
+ Iterator<MCParticle> iter = mcps.iterator();
+ while (iter.hasNext()){
+ MCParticle next = iter.next();
+ if (!h.getMCParticles().contains(next)){
+ //System.out.println(mcps.toString());
+ iter.remove();
+ }
+ }
+ }
+ if(mcps.isEmpty())return;
+ /*//Check method to see if you are getting the right mc particles in the list
+ List<MCParticle> testlist = util.FindMCParticleInAllHits(seed.getHits());
+ if(testlist.containsAll(temp)){System.out.println("true");}*/
+ /* for (HelicalTrackHit h : seed.getHits()) {
+ //hitmcmap.get(h).addAll(h.getMCParticles());
+ mcps.addAll(h.getMCParticles());
+ }
+
+
+
+ }*/
+ double costheta = 1/(Math.sqrt(1+((Math.pow(1/hfit.slope(), 2)))));
+ if(hfit.pT(5)>1.1 && Math.abs(hfit.dca()) < 9 && Math.abs(hfit.z0()) < 9 && Math.abs(costheta) < 0.985 ){
+ for(MCParticle mcp : mcps){
+ boolean flag = true;
+ if(hitmcmap.containsKey(mcp)){
+ if(mcheck.contains(mcp)){
+ flag = false;
+ break;
+ }
+ else{
+ mcheck.add(mcp);
+ }
+ if(flag){
+ System.out.println("From SeedTrackerDig Tool");
+ System.out.println("P:"+mcp.getMomentum().toString());
+ System.out.println("Endpoint:"+mcp.getEndPoint().toString());
+ System.out.print("SeedlayerHits: ");
+ for(HelicalTrackHit seedhit : hitmcmap.get(mcp)){
+ System.out.print(seedhit.getLayerIdentifier()+"(time: "+seedhit.getTime()+")");
+ System.out.print("("+seedhit.getPosition()[0]+",");
+ System.out.print(seedhit.getPosition()[1]+",");
+ System.out.println(+seedhit.getPosition()[2]+")");
+ }
+ System.out.println("chisq0: "+hfit.chisq()[0]);
+ System.out.println("chisq1: "+hfit.chisq()[1]);
+ System.out.println("chisqtotal: "+hfit.chisqtot());
+ System.out.println();
+ aida.histogram1D("Efficiency/NoNHit/Has Fit Missed chsq:",200,0,2).fill(1);
+ }
+ }
+ }
+
+
+ /*for(HelicalTrackHit h : seed.getHits()){
+ aida.histogram1D("Efficiency/NoNHit/scattering/"+h.getLayerIdentifier(),200,-10,10).fill(+h.dr());
+ aida.histogram1D("Efficiency/NoNHit/scattering/"+h.getLayerIdentifier(),200,-10,10).fill(+h.drphi());
+ }*/
+ }
+ }
+ }
+
lcsim/src/org/lcsim/recon/tracking/seedtracker/analysis
diff -N TrackHitHistograms.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackHitHistograms.java 27 Aug 2008 17:58:15 -0000 1.1
@@ -0,0 +1,302 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.analysis;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.EventHeader;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.recon.tracking.seedtracker.analysis.HelixParamCalculator;
+import org.lcsim.util.aida.AIDA;
+import java.util.List;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+
+/**
+ *
+ * @author Pelham Keahey
+ */
+public class TrackHitHistograms {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private double phi,r,z;
+ private double phih,rh,zh;
+ private double simr,simz,simphi;
+ HelicalTrackFit fit;
+ MCParticle mcp;
+ EventHeader event;
+ List<HelicalTrackHit> hitlist;
+
+ /**
+ *
+ * @param pfit HelicalTrackfit
+ * @param pmcp MCParticle
+ * @param pevent EventHeader
+ * @param phit List of Hits
+ */
+ public TrackHitHistograms(HelicalTrackFit pfit, MCParticle pmcp,EventHeader pevent)
+ {
+ fit=pfit;
+ mcp=pmcp;
+ event = pevent;
+ hitlist = event.get(HelicalTrackHit.class,"HelicalTrackHits");
+
+
+ }/**
+ *
+ * Calculates the hit locations in polar coordinates and displays residual histograms between the hit locations from the mcparticle
+ * and the helicaltrackhits
+ */
+ public void drawMCParticlevHelicalTrackHit()
+ {
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, event);
+
+ for(HelicalTrackHit hit : hitlist)
+ {
+
+ //Calculations for MC particle
+
+ //find x0,y0 coordinates
+ double x0 = calc.getX0();
+ double y0 = calc.getY0();
+
+ //find the pathlength of the helix directly from the hit
+ double pathlength=0;
+ double closer=999999999,difference;
+
+ List<Double> listd;
+ HelixUtils util = new HelixUtils();
+ if(hit.BarrelEndcapFlag().isBarrel()){
+ listd = util.PathToCylinder(fit, hit.r(), 5000, 10);
+ closer = 999999999;
+ if(!listd.isEmpty()){
+ pathlength = listd.get(0);
+ }
+ for(double temp : listd){
+ BasicHep3Vector tempvec = (BasicHep3Vector) util.PointOnHelix(fit, temp);
+ difference = Math.abs(hit.z()-tempvec.z());
+ if(difference<closer){
+ closer = difference;
+ pathlength = temp;
+ }
+ }
+
+ }
+ else{
+ pathlength = util.PathToZPlane(fit, hit.z());
+ }
+
+ //varible for easier calculation temp = C*s/2
+ double temp = ((calc.getMCOmega()*pathlength)/2);
+
+ //find x,y coordinates
+ double x = x0 + (((pathlength*Math.sin(temp))/temp)*Math.cos(calc.getPhi0()-temp));
+ double y = y0 + (((pathlength*Math.sin(temp))/temp)*Math.sin(calc.getPhi0()-temp));
+
+ //Find MC phi, MC radius (r), MC z coordinate EDIT
+ phi = Math.atan2(y, x);
+ if(phi<0)
+ {
+ phi+=2*Math.PI;
+ }
+ r = Math.sqrt(x*x + y*y);
+ z = calc.getZ0() + pathlength*calc.getSlopeSZPlane();
+
+ //Calculations directly from hit
+
+ //Calculations directly from HelicalTrackHit
+ phih = hit.phi();
+ rh = hit.r();
+ zh = hit.z();
+
+ if(hit.getLayerIdentifier().contains("Tracker"))
+ {
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Tracker/Barrel/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -25, 25).fill((r*phi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Tracker/Barrel/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Tracker/Barrel/Z/Z-"+hit.getLayerIdentifier() , 200, -55,55).fill(z-zh);
+ }
+ }
+ if(hit.getLayerIdentifier().contains("Tracker"))
+ {
+ if(hit.BarrelEndcapFlag().isEndcap())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Tracker/EndCap/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -25, 25).fill((r*phi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Tracker/EndCap/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Tracker/EndCap/Z/Z-"+hit.getLayerIdentifier() , 200, -5,5).fill(z-zh);
+ }
+ }
+ if(hit.getLayerIdentifier().contains("Vertex"))
+ {
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Vertex/Barrel/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -.5, .5).fill((r*phi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Vertex/Barrel/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Vertex/Barrel/Z/Z-"+hit.getLayerIdentifier() , 200, -55,55).fill(z-zh);
+ }
+ }
+ if(hit.getLayerIdentifier().contains("Vertex"))
+ {
+ if(hit.BarrelEndcapFlag().isEndcap())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Vertex/EndCap/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -30, 30).fill((r*phi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Vertex/EndCap/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("HitTypeAnalysis/MCParticleVSHelicalTrackFit/Residual-Vertex/EndCap/Z/Z-"+hit.getLayerIdentifier() , 200, -55,55).fill(z-zh);
+ }
+ }
+ }
+ }
+ /*
+ * Calculates the hits in polar and displays the residual histograms of simtrackerhits from helicaltrackhits
+ */
+ public void drawSimTrackervHelicalTrack()
+ {
+
+ List<SimTrackerHit> simhitlist = new ArrayList<SimTrackerHit>();
+
+ for(List<SimTrackerHit> l : event.get(SimTrackerHit.class))
+ {
+ simhitlist.addAll(l);
+ }
+
+ double [] alld;
+ double distancex,distancey,distancez;
+ double closest=999999999,testclose;
+ double zstripdiff;
+ double closestdiff;
+ double smallres;
+ double us,resunmeas;
+
+
+ for(HelicalTrackHit hit : hitlist)
+ {
+ if(hit instanceof HelicalTrackCross)
+ {
+ smallres = 999999999;
+ HelicalTrackCross cross = (HelicalTrackCross) hit;
+
+ for(HelicalTrackStrip strip : cross.getStrips())
+ {
+ SimTrackerHit closesthit =null;
+ closestdiff = 999999999;
+ for(SimTrackerHit simhit : simhitlist)
+ {
+ BasicHep3Vector simclose = new BasicHep3Vector(simhit.getPoint());
+ zstripdiff = Math.abs(simclose.z()- strip.origin().z());
+ if(zstripdiff<closestdiff)
+ {
+ closestdiff = zstripdiff;
+ closesthit = simhit;
+ }
+
+ }
+ BasicHep3Vector closesthitvec = new BasicHep3Vector(closesthit.getPoint());
+ BasicHep3Vector diffvector = new BasicHep3Vector();
+ diffvector = (BasicHep3Vector) VecOp.sub(closesthitvec, strip.origin());
+ us = VecOp.dot(diffvector, strip.u());
+ resunmeas = strip.umeas() - us;
+ if(resunmeas<smallres)
+ {
+ smallres = resunmeas;
+ }
+ aida.histogram1D("HitTypeAnalysis/Unmeasured/EndCaps/"+Math.round(strip.origin().z())+":ZPOS Unmeasured Strip ", 200, -1,1).fill(smallres);
+
+ }
+ }
+ }
+ for(HelicalTrackHit hit : hitlist)
+ {
+ SimTrackerHit hitholder = null;
+ closest = 999999999;
+ for(SimTrackerHit simhit : simhitlist)
+ {
+
+ //Search for the closest points
+ alld = simhit.getPoint();
+ distancex = hit.x()-alld[0];
+ distancey = hit.y()-alld[1];
+ distancez = hit.z()-alld[2];
+ testclose = Math.sqrt((distancex*distancex) + (distancey*distancey)+(distancez*distancez));
+ if(testclose<closest)
+ {
+ closest=testclose;
+ hitholder = simhit;
+ }
+
+
+ }
+ if(hitholder==null)break;
+ alld = hitholder.getPoint();
+ simr = Math.sqrt((alld[0]*alld[0])+(alld[1]*alld[1]));
+ simz = alld[2];
+ simphi = Math.atan2(alld[1], alld[0]);
+ if(simphi<0)
+ {
+ simphi+=2*Math.PI;
+ }
+
+ rh = hit.r();
+ zh = hit.z();
+ phih = hit.phi();
+
+ if(hit.getLayerIdentifier().contains("Tracker"))
+ {
+
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Tracker/Barrel/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -1,1).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Tracker/Barrel/Residual/R/R"+hit.getLayerIdentifier(), 200, -.5,.5).fill(simr-rh);
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Tracker/Barrel/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -55,55).fill(simz-zh);
+ }
+
+ else
+ {
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Tracker/Endcap/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -10, 10).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Tracker/Endcap/Residual/R/R"+hit.getLayerIdentifier(), 200, -3,3).fill(simr-rh);
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Tracker/Endcap/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -3,3).fill(simz-zh);
+ }
+ }
+ else
+ {
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Vertex/Barrel/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -1,1).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Vertex/Barrel/Residual/R/R"+hit.getLayerIdentifier(), 200, -1,1).fill(simr-rh);
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Vertex/Barrel/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -1,1).fill(simz-zh);
+ }
+ else
+ {
+
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Vertex/Endcap/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -1, 1).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Vertex/Endcap/Residual/R/R"+hit.getLayerIdentifier(), 200, -1,1).fill(simr-rh);
+ aida.histogram1D("HitTypeAnalysis/SimTrackerVSHelicalTrack/Vertex/Endcap/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -1,1).fill(simz-zh);
+ }
+ }
+ }
+ }
+
+
+
+
+}
CVSspam 0.2.8