lcsim/src/org/lcsim/contrib/Pelham/Example1
diff -N TrackHitHistograms.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackHitHistograms.java 25 Jul 2008 18:29:28 -0000 1.1
@@ -0,0 +1,322 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Pelham.Example1;
+import java.util.Map;
+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.fit.helicaltrack.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;
+import java.util.HashMap;
+
+/**
+ *
+ * @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
+ * @param pmcp
+ * @param pevent
+ * @param phit
+ */
+ public TrackHitHistograms(HelicalTrackFit pfit, MCParticle pmcp, EventHeader pevent,List<HelicalTrackHit> phit)
+ {
+ fit=pfit;
+ event = pevent;
+ mcp=pmcp;
+ hitlist = phit;
+
+
+ }
+ public void drawMCParticlevHelicalTrackHit()
+ {
+ HelixParamCalculator calc = new HelixParamCalculator(mcp, event);
+ double [] pars = new double[5];
+ pars[0]=calc.getDCA();
+ pars[1]=calc.getPhi0();
+ pars[2]=calc.getMCOmega();
+ pars[3]=calc.getZ0();
+ pars[4]=calc.getSlopeSZPlane();
+
+ for(HelicalTrackHit hit : hitlist)
+ {
+ Map<HelicalTrackHit, Double> lol = new HashMap<HelicalTrackHit, Double>();
+ lol.put(hit, calc.getArcLength());
+
+ HelicalTrackFit test = new HelicalTrackFit(pars,null, null, null, lol, null);
+ //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;
+ double closer=999999999,difference;
+ double arclength = test.PathMap().get(hit);
+
+ List<Double> listd;
+ HelixUtils util = new HelixUtils();
+ if(hit.BarrelEndcapFlag().isBarrel()){
+ listd = util.PathToCylinder(fit, hit.r(), 5000, 10);
+ closer = 999999999;
+ 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;
+ //System.out.println(arclength+":"+pathlength+"fitmap vs pathcylinderbarrel");
+ }
+ }
+
+ }
+ else{
+ pathlength = util.PathToZPlane(fit, hit.z());
+ //System.out.println(arclength+":"+pathlength+"fitmap vs pathcylinderendcap");
+ }
+
+
+ /*
+ //find point closest on barrel
+ util.PathToCylinder(fit, r, simr, j);
+ //for endcap for pathlength to given z value
+ util.PathToZPlane(fit, z);
+ util.PointOnHelix(fit, r);
+
+ HelicalTrackCross cross = new HelicalTrackCross();
+ //HelicalTrackCross From cross you can get two strips
+ //from a strip
+ //u is measuemn
+ //unmeasured coordinite hitutils*/
+
+
+
+ //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().startsWith("Tracker"))
+ {
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Tracker/Barrel/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -25, 25).fill((r*phi)-(rh*phih));
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Tracker/Barrel/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Tracker/Barrel/Z/Z-"+hit.getLayerIdentifier() , 200, -55,55).fill(z-zh);
+ }
+ }
+ if(hit.getLayerIdentifier().startsWith("Tracker"))
+ {
+ if(hit.BarrelEndcapFlag().isEndcap())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Tracker/EndCap/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -25, 25).fill((r*phi)-(rh*phih));
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Tracker/EndCap/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Tracker/EndCap/Z/Z-"+hit.getLayerIdentifier() , 200, -5,5).fill(z-zh);
+ }
+ }
+ if(hit.getLayerIdentifier().startsWith("Vertex"))
+ {
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Vertex/Barrel/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -.5, .5).fill((r*phi)-(rh*phih));
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Vertex/Barrel/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Vertex/Barrel/Z/Z-"+hit.getLayerIdentifier() , 200, -55,55).fill(z-zh);
+ }
+ }
+ if(hit.getLayerIdentifier().startsWith("Vertex"))
+ {
+ if(hit.BarrelEndcapFlag().isEndcap())
+ {
+
+ //Residual Histograms
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Vertex/EndCap/r*phi/r*phi-"+hit.getLayerIdentifier(), 200, -30, 30).fill((r*phi)-(rh*phih));
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Vertex/EndCap/R/R-"+hit.getLayerIdentifier() , 200, -.1,.1).fill(r-rh);
+ aida.histogram1D("MCParticleVSHelicalTrackFit/Residual-Vertex/EndCap/Z/Z-"+hit.getLayerIdentifier() , 200, -55,55).fill(z-zh);
+ }
+ }
+ }
+ }
+ 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("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;
+ }
+
+
+ }
+
+ 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().startsWith("Tracker"))
+ {
+
+ if(hit.BarrelEndcapFlag().isBarrel())
+ {
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("SimTrackerVSHelicalTrack/Tracker/Barrel/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -1,1).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("SimTrackerVSHelicalTrack/Tracker/Barrel/Residual/R/R"+hit.getLayerIdentifier(), 200, -.5,.5).fill(simr-rh);
+ aida.histogram1D("SimTrackerVSHelicalTrack/Tracker/Barrel/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -55,55).fill(simz-zh);
+ }
+
+ else
+ {
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("SimTrackerVSHelicalTrack/Tracker/Endcap/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -10, 10).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("SimTrackerVSHelicalTrack/Tracker/Endcap/Residual/R/R"+hit.getLayerIdentifier(), 200, -3,3).fill(simr-rh);
+ aida.histogram1D("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("SimTrackerVSHelicalTrack/Vertex/Barrel/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -1,1).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("SimTrackerVSHelicalTrack/Vertex/Barrel/Residual/R/R"+hit.getLayerIdentifier(), 200, -1,1).fill(simr-rh);
+ aida.histogram1D("SimTrackerVSHelicalTrack/Vertex/Barrel/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -1,1).fill(simz-zh);
+ }
+ else
+ {
+
+ //Residual Histograms: SimTrackerHit - HelicalTrackHit
+ aida.histogram1D("SimTrackerVSHelicalTrack/Vertex/Endcap/Residual/r*phi/R*Phi"+hit.getLayerIdentifier(), 200, -1, 1).fill((simr*simphi)-(rh*phih));
+ aida.histogram1D("SimTrackerVSHelicalTrack/Vertex/Endcap/Residual/R/R"+hit.getLayerIdentifier(), 200, -1,1).fill(simr-rh);
+ aida.histogram1D("SimTrackerVSHelicalTrack/Vertex/Endcap/Residual/Z/Z"+hit.getLayerIdentifier(), 200, -1,1).fill(simz-zh);
+ }
+ }
+ }
+ }
+
+
+
+
+}