hps-java/src/main/java/org/lcsim/hps/recon/tracking
diff -u -r1.5 -r1.6
--- TrackUtils.java 28 Sep 2012 22:07:28 -0000 1.5
+++ TrackUtils.java 9 Oct 2012 01:16:11 -0000 1.6
@@ -1,46 +1,63 @@
package org.lcsim.hps.recon.tracking;
//--- hep ---//
+import hep.physics.matrix.SymmetricMatrix;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.HashMap;
import java.util.List;
+import java.util.Map;
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.solids.Box;
+import org.lcsim.detector.solids.Point3D;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.SiSensor;
//--- org.lcsim ---//
import org.lcsim.event.Track;
import org.lcsim.event.TrackState;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
-import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.*;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
/**
*
* @author Omar Moreno <[log in to unmask]>
- * @version $Id: TrackUtils.java,v 1.5 2012/09/28 22:07:28 phansson Exp $
+ * @version $Id: TrackUtils.java,v 1.6 2012/10/09 01:16:11 phansson Exp $
* TODO: Switch to JLab coordinates
*/
public class TrackUtils {
- Track track;
+ boolean _debug = false;
+ boolean isTrackSet = false;
double[] trackParameters;
/**
*
*/
public TrackUtils(){
- track = null;
+ isTrackSet = false;
}
/**
*
*/
- public void setTrack(Track track){
- this.track = track;
+ public void setTrack(Track track){
trackParameters = new double[9];
this.setTrackParameters(track);
+ isTrackSet = true;
+ }
+
+ /**
+ *
+ */
+ public void setTrack(double[] pars){
+ trackParameters = new double[9];
+ this.setTrackParameters(pars);
+ isTrackSet = true;
}
/**
@@ -111,7 +128,7 @@
public double getXOnHelixProjection(double y){
// Check if a track has been set
- if(track == null) throw new RuntimeException("Track has not been set!");
+ if(!isTrackSet) throw new RuntimeException("Track has not been set!");
// Find the x position
return this.getXC()
@@ -124,11 +141,15 @@
public double getYOnHelixProjection(double x){
// Check if a track has been set
- if(track == null) throw new RuntimeException("Track has not been set!");
-
+ if(!isTrackSet) throw new RuntimeException("Track has not been set!");
+ if(_debug) System.out.println(this.getClass().getSimpleName() + ": getYOnHelixProjection for x=" + x + " R=" + this.getR() + " YC=" + this.getYC() + " XC=" + this.getXC());
// Find the y position
- return this.getYC()
- + Math.signum(this.getR())*Math.sqrt(this.getR()*this.getR() - Math.pow(x-this.getXC(), 2));
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": => y = " + Math.sqrt(this.getR()*this.getR() - Math.pow(x-this.getXC(), 2)));
+ System.out.println(this.getClass().getSimpleName() + ": => y = " + (this.getYC() + Math.signum(this.getR())*Math.sqrt(this.getR()*this.getR() - Math.pow(x-this.getXC(), 2))));
+ }
+
+ return this.getYC() + Math.signum(this.getR())*Math.sqrt(this.getR()*this.getR() - Math.pow(x-this.getXC(), 2));
}
/**
@@ -150,7 +171,7 @@
*
*/
public void printTrackParameters(){
- System.out.println(" DOCA: " + this.getDoca() +
+ System.out.println(this.getClass().getSimpleName() + ": DOCA: " + this.getDoca() +
" phi0: " + this.getPhi0() +
" R: " + this.getR() +
" x0: " + this.getX0() +
@@ -159,8 +180,7 @@
" tan(Lambda): " + this.getTanLambda() +
" xc: " + this.getXC() +
" yc: " + this.getYC() );
-
- }
+ }
/**
@@ -178,16 +198,37 @@
trackParameters[6] = track.getTrackParameter(4); // tan(Lambda)
trackParameters[7] = (this.getR() - this.getDoca())*Math.sin(this.getPhi0()); // xc
trackParameters[8] = -(this.getR() - this.getDoca())*Math.cos(this.getPhi0()); // yc
+ if(_debug) this.printTrackParameters();
+ }
+
+ /**
+ *
+ */
+ private void setTrackParameters(double[] pars){
+
+ // All track parameters are in LCSim coordinates system
+ trackParameters[0] = pars[0]; // DOCA
+ trackParameters[1] = pars[1]; // phi0
+ trackParameters[2] = 1/pars[2]; // R
+ trackParameters[3] = -this.getDoca()*Math.sin(this.getPhi0()); // x0
+ trackParameters[4] = this.getDoca()*Math.cos(this.getPhi0()); // y0
+ trackParameters[5] = pars[3]; // z0
+ trackParameters[6] = pars[4]; // tan(Lambda)
+ trackParameters[7] = (this.getR() - this.getDoca())*Math.sin(this.getPhi0()); // xc
+ trackParameters[8] = -(this.getR() - this.getDoca())*Math.cos(this.getPhi0()); // yc
+ if(_debug) this.printTrackParameters();
+
}
/**
*
*/
public Hep3Vector extrapolateHelixToXPlane(double x){
- double y = this.getYOnHelixProjection(x);
+ if(_debug) System.out.println(this.getClass().getSimpleName() + ":extrapolateHelixToXPlane for x= " + x );
+ double y = this.getYOnHelixProjection(x);
double s = this.arcLength(x, y);
double z = this.getZ0() + s*this.getTanLambda();
-
+ if(_debug) System.out.println(this.getClass().getSimpleName() + ": y " + y + " s " + s + " z " + z);
return new BasicHep3Vector(x, y, z);
}
@@ -213,9 +254,8 @@
HelicalTrackFit cfit = helfit;
if(cfit==null) {
- SeedTrack st = (SeedTrack) track;
- SeedCandidate seed = st.getSeedCandidate();
- cfit = seed.getHelix();
+ System.out.println(this.getClass().getSimpleName() + ": no helical track fit was provided for interception calculation");
+ System.exit(1);
}
if(_DEBUG) {
System.out.println(this.getClass().getSimpleName() + ": calculateHelixIntercept ---");
@@ -265,7 +305,7 @@
public double getPhi(Hep3Vector position){
// Check if a track has been set
- if(track == null) throw new RuntimeException("Track has not been set!");
+ if(!isTrackSet) throw new RuntimeException("Track has not been set!");
double x = Math.sin(this.getPhi0() - (1/this.getR())*(position.x() - this.getX0()));
double y = Math.cos(this.getPhi0() + (1/this.getR())*(position.y() - this.getY0()));
@@ -287,4 +327,92 @@
return this.getTanLambda()/Math.sqrt(1 + Math.pow(this.getTanLambda(), 2));
}
+
+ /**
+ *
+ */
+ public double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2){
+ return .5*(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1);
+ }
+
+
+ /**
+ *
+ */
+ public boolean sensorContainsTrack(Hep3Vector trackPosition, SiSensor sensor){
+ boolean debug = false;
+ ITransform3D localToGlobal = sensor.getGeometry().getLocalToGlobal();
+
+ Hep3Vector sensorPos = sensor.getGeometry().getPosition();
+ Box sensorSolid = (Box) sensor.getGeometry().getLogicalVolume().getSolid();
+ Polygon3D sensorFace = sensorSolid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
+ if(debug){
+ System.out.println(this.getClass().getSimpleName() + ": Sensor: " + SvtUtils.getInstance().getDescription(sensor));
+ System.out.println(this.getClass().getSimpleName() + ": Track Position: " + trackPosition.toString());
+ }
+
+ List<Point3D> vertices = new ArrayList<Point3D>();
+ for(int index = 0; index < 4; index++){
+ vertices.add(new Point3D());
+ }
+ for(Point3D vertex : sensorFace.getVertices()){
+ if(vertex.y() < 0 && vertex.x() > 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(0, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(0, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if(debug){
+ System.out.println(this.getClass().getSimpleName() + ": Vertex 1 Position: " + vertices.get(0).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 1 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ }
+ else if(vertex.y() > 0 && vertex.x() > 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(1, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(1, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if(debug){
+ System.out.println(this.getClass().getSimpleName() + ": Vertex 2 Position: " + vertices.get(1).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 2 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ }
+ else if(vertex.y() > 0 && vertex.x() < 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(2, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(2, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if(debug){
+ System.out.println(this.getClass().getSimpleName() + ": Vertex 3 Position: " + vertices.get(2).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 3 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ }
+ else if(vertex.y() < 0 && vertex.x() < 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(3, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(3, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ if(debug){
+ System.out.println(this.getClass().getSimpleName() + ": Vertex 4 Position: " + vertices.get(3).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 4 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ }
+ }
+
+ double area1 = this.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z());
+ double area2 = this.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z());
+ double area3 = this.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z());
+ double area4 = this.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());
+
+ if((area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0)) return true;
+
+ return false;
+ }
+
+ public HelicalTrackFit makeHelicalTrackFit(double[] helix_parameters) {
+
+ SymmetricMatrix cov = new SymmetricMatrix(5);
+ double[] chisq = new double[2];
+ int[] ndf = new int[2];
+ Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>();
+ Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
+ HelicalTrackFit fit = new HelicalTrackFit(helix_parameters, cov, chisq, ndf,smap,msmap);
+ return fit;
+ }
+
}