8 removed files
lcsim/src/org/lcsim/contrib/seedtracker/analysis
diff -N AnalysisUtils.java
--- AnalysisUtils.java 20 Aug 2008 01:34:21 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,306 +0,0 @@
-/*
- * 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;
-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
- * @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/contrib/seedtracker/analysis
diff -N EfficiencyHistograms.java
--- EfficiencyHistograms.java 15 Aug 2008 22:22:03 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,95 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.lcsim.contrib.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/contrib/seedtracker/analysis
diff -N HelixParamCalculator.java
--- HelixParamCalculator.java 6 Aug 2008 22:21:56 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,205 +0,0 @@
-/*
- * HelixParameterCalculator.java
- *
- * Created on July 7th, 2008, 11:09 AM
- *
- *
- */
-package org.lcsim.contrib.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/contrib/seedtracker/analysis
diff -N HelixParamHistograms.java
--- HelixParamHistograms.java 15 Aug 2008 22:59:00 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,177 +0,0 @@
-package org.lcsim.contrib.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.contrib.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/contrib/seedtracker/analysis
diff -N HistogramAnalysisDriver.java
--- HistogramAnalysisDriver.java 15 Aug 2008 22:59:20 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,143 +0,0 @@
-/*
- * AnalysisDriver.java
- *
- * Created on February 11, 2008, 11:47 AM
- *
- */
-
-package org.lcsim.contrib.seedtracker.analysis;
-
-
-import org.lcsim.contrib.seedtracker.analysis.TrackHitHistograms;
-import org.lcsim.contrib.seedtracker.analysis.HelixParamHistograms;
-import java.util.List;
-import java.util.Map;
-import java.util.HashMap;
-
-import org.lcsim.contrib.seedtracker.SeedTrack;
-import org.lcsim.contrib.seedtracker.SeedCandidate;
-import org.lcsim.contrib.seedtracker.analysis.AnalysisUtils;
-import org.lcsim.contrib.seedtracker.analysis.EfficiencyHistograms;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.contrib.seedtracker.analysis.HelixParamCalculator;
-import org.lcsim.contrib.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/contrib/seedtracker/analysis
diff -N PurityHistograms.java
--- PurityHistograms.java 15 Aug 2008 22:48:31 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,54 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.lcsim.contrib.seedtracker.analysis;
-
-import org.lcsim.contrib.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/contrib/seedtracker/analysis
diff -N SeedTrackDiag.java
--- SeedTrackDiag.java 13 Aug 2008 18:35:26 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,135 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.lcsim.contrib.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.contrib.seedtracker.SeedCandidate;
-import org.lcsim.contrib.seedtracker.SeedLayer;
-import org.lcsim.contrib.seedtracker.SeedLayer.SeedType;
-import org.lcsim.contrib.seedtracker.analysis.AnalysisUtils;
-import org.lcsim.contrib.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/contrib/seedtracker/analysis
diff -N TrackHitHistograms.java
--- TrackHitHistograms.java 15 Aug 2008 22:59:08 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,302 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.lcsim.contrib.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.contrib.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