hps-java/src/main/java/org/lcsim/hps/recon/tracking
diff -N FieldMap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FieldMap.java 23 Jul 2012 21:45:02 -0000 1.1
@@ -0,0 +1,175 @@
+package org.lcsim.hps.recon.tracking;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.Reader;
+import java.util.*;
+import org.lcsim.conditions.ConditionsManager;
+import org.lcsim.conditions.ConditionsSet;
+import org.lcsim.hps.monitoring.HPSCalibrationListener;
+import org.lcsim.hps.util.Pair;
+
+/**
+
+ @author Mathew Graham <[log in to unmask]> $Id: FieldMap.java,v 1.14
+ 2012/07/23 23:02:57 mgraham Exp $
+ */
+public class FieldMap {
+
+ // TODO: Change all pairs such that FPGA is the fist value
+ // TODO: Change all map keys to type SiSensor?
+ public static final double DIPOLE_EDGE = 914; // mm
+ public static Map<Pair<Integer/*
+ zbin
+ */, Integer/*
+ xbin
+ */>, Double/*
+ b_y
+ */> fieldMap=new HashMap<Pair<Integer,Integer>,Double>();;
+ public static Map<Pair<Integer/*
+ zbin
+ */, Integer/*
+ zbin
+ */>, Pair<Double/*
+ zpos
+ */, Double/*
+ xpos
+ */>> fieldBins=new HashMap<Pair<Integer,Integer>,Pair<Double,Double>>();;
+ private static boolean fieldMapLoaded = false;
+
+ /**
+ Default Ctor
+ */
+ private FieldMap() {
+ }
+
+ public static void loadFieldMap() {
+ loadFieldMap(null);
+ }
+
+ public static void loadFieldMap(Date date, int runNumber) {
+ loadFieldMap(date);
+ }
+
+ public static void loadFieldMap(Date date) {
+ //write something here to read in constants from calibration file
+ ConditionsManager conditions = ConditionsManager.defaultInstance();
+
+ String filePath = null;
+
+ if (date != null) {
+ ConditionsSet calibSet = conditions.getConditions("FieldMap/bfieldmap");
+ filePath = HPSCalibrationListener.getCalibForDate(calibSet, date);
+ }
+
+ if (filePath == null) {
+ filePath = "FieldMap/bfieldmap.dat";
+ }
+
+ try {
+ Reader baselineReader = conditions.getRawConditions(filePath).getReader();
+ loadBField(baselineReader);
+ } catch (IOException e) {
+ throw new RuntimeException("couldn't get baseline file", e);
+ }
+
+
+
+
+ }
+
+ private static void loadBField(Reader baselineReader) throws IOException {
+
+ BufferedReader bufferedBaseline = new BufferedReader(baselineReader);
+ String line;
+ Double xval=0.0,zval=0.0,bY;
+ Integer iz=-1,ix=0;
+ while (true) {
+ try {
+ line = bufferedBaseline.readLine();
+ } catch (IOException e) {
+ throw new RuntimeException("couldn't parse baseline file", e);
+ }
+ if (line == null) {
+ break;
+ }
+
+ if (line.indexOf("#") != -1) {
+ line = line.substring(0, line.indexOf("#"));
+ }
+
+ StringTokenizer tok = new StringTokenizer(line);
+ System.out.println(line.toString());
+ List<Double> vals = new ArrayList<Double>();
+
+ while ((vals = getNumbersInLine(tok)).size() > 0) {
+ System.out.println(getNumbersInLine(tok).size());
+ for (int i = 0; i < vals.size(); i++) {
+ double val = vals.get(i);
+ if (i == 0 && val >= 0) {
+ xval = 0.0;
+ zval = val * 10.0;//convert from cm to mm
+ ix = 0;
+ iz++;
+ } else {
+ bY = val / 10000.0 * 0.491 / 0.5;//gauss-->tesla and normalize by our nominal bfield (file is for 0.5T)
+ bY = Math.abs(bY);
+ Pair zx = new Pair(zval, xval);
+ Pair izix = new Pair(iz, ix);
+// System.out.println("Putting B = "+bY+" & (Z,X) = "+zval+","+xval+")");
+ fieldMap.put(izix, bY);
+// System.out.println("Making _fieldBins pair = ("+iz+","+ix+")");
+ fieldBins.put(izix, zx);
+ xval += 10.0;
+ ix++;
+ }
+ }
+ }
+ }
+ fieldMapLoaded=true;
+
+ }
+
+ public static boolean fieldMapLoaded(){
+ return fieldMapLoaded;
+ }
+
+ private static List<Double> getNumbersInLine(StringTokenizer tok) throws IOException {
+ List<Double> nums = new ArrayList<Double>();
+ while (tok.hasMoreTokens()) {
+ String tokVal = tok.nextToken();
+// System.out.println(tokVal);
+ nums.add(Double.valueOf(tokVal).doubleValue());
+ }
+// System.out.println("Returning list");
+ return nums;
+ }
+
+ public static double getFieldFromMap(double zval, double xval) {
+ Pair<Integer, Integer> bin = getFieldBin(zval, xval);
+ return fieldMap.get(bin);
+ }
+
+ private static Pair<Integer, Integer> getFieldBin(double zval, double xval) {
+ Integer iz, ix;
+ int nZ = 150;
+ int nX = 25;
+ double zNew = Math.abs(zval - DIPOLE_EDGE / 2);
+ double xNew = Math.abs(xval);
+ for (iz = 0; iz < nZ; iz++) {
+ Pair tmp = fieldBins.get(new Pair(iz, 0));
+ Double zpair = (Double) tmp.getFirstElement();
+ if (zpair > zNew)
+ break;
+ }
+ for (ix = 0; ix < nX; ix++) {
+ Pair tmp = fieldBins.get(new Pair(1, ix));
+ Double xpair = (Double) tmp.getSecondElement();
+ if (xpair > xNew)
+ break;
+ }
+ Pair izix = new Pair(iz, ix);
+ return izix;
+
+ }
+}
hps-java/src/main/java/org/lcsim/hps/recon/tracking
diff -u -r1.1 -r1.2
--- HPSTrack.java 7 Jul 2011 20:57:38 -0000 1.1
+++ HPSTrack.java 23 Jul 2012 21:45:02 -0000 1.2
@@ -4,9 +4,17 @@
*/
package org.lcsim.hps.recon.tracking;
+import java.util.HashMap;
import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
+import java.io.IOException;
+import java.io.StreamTokenizer;
+import java.util.ArrayList;
+import java.util.HashMap;
import java.util.List;
import java.util.Map;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
@@ -15,15 +23,27 @@
import org.lcsim.fit.helicaltrack.MultipleScatter;
import org.lcsim.hps.event.BeamSpot;
import org.lcsim.hps.event.HPSTransformations;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.Line;
+import org.lcsim.util.swim.Trajectory;
+import org.lcsim.spacegeom.SpaceVector;
+import static org.lcsim.constants.Constants.fieldConversion;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.hps.util.Pair;
/**
- * Class HPSTrack: extension of HelicalTrackFit to include HPS-specific variables
- * other useful things.
- * @author mgraham
- * created on 6/27/2011
+ Class HPSTrack: extension of HelicalTrackFit to include HPS-specific variables
+ other useful things.
+
+ @author mgraham created on 6/27/2011
*/
public class HPSTrack extends HelicalTrackFit {
+ public static final double ECAL_FACE = 1450; // mm
+ public static final double DIPOLE_EDGE = 914; // mm
private BeamSpot _beam;
//all of the variables defined below are in the jlab (detector) frame
//this position & momentum are measured at the DOCA to the Y-axis,
@@ -37,7 +57,12 @@
private Hep3Vector _pDocaZ;
private Hep3Vector _posDocaZ;
private HPSTransformations _detToTrk;
- private double bField = 0.5; // make this set-able
+ private double bField = 0.491; // make this set-able
+ private boolean _debug = false;
+ private boolean _debugForward = false;
+ private Trajectory _trajectory;
+ Map<Pair<Integer, Integer>, Double> _fieldMap;
+ Map<Pair<Integer, Integer>, Pair<Double, Double>> _fieldBins;
public HPSTrack(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf,
Map<HelicalTrackHit, Double> smap, Map<HelicalTrackHit, MultipleScatter> msmap,
@@ -58,6 +83,51 @@
calculateParametersAtDocaY();
calculateParametersAtDocaZ();
}
+
+ public HPSTrack(HelicalTrackFit htf) {
+ super(htf.parameters(), htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
+ _detToTrk = new HPSTransformations();
+ calculateParametersAtTarget();
+ calculateParametersAtDocaY();
+ calculateParametersAtDocaZ();
+ }
+
+ public Map<Integer, Double[]> trackTrajectory(double zStart, double zStop, int nSteps) {
+ Map<Integer, Double[]> traj = new HashMap<Integer, Double[]>();
+ double step = (zStop - zStart) / nSteps;
+ Double zVal = zStart;
+ Integer nstep = 0;
+ while (zVal < zStop) {
+ Double[] xyz = {0.0, 0.0, 0.0};
+ double s = HelixUtils.PathToXPlane(this, zVal, 1000.0, 1).get(0);
+ xyz[0] = HelixUtils.PointOnHelix(this, s).y();
+ xyz[1] = HelixUtils.PointOnHelix(this, s).z();
+ xyz[2] = zVal;
+ traj.put(nstep, xyz);
+ zVal += step;
+ nstep++;
+ }
+ return traj;
+ }
+
+ public Map<Integer, Double[]> trackDirection(double zStart, double zStop, int nSteps) {
+ Map<Integer, Double[]> traj = new HashMap<Integer, Double[]>();
+ double step = (zStop - zStart) / nSteps;
+ Double zVal = zStart;
+
+ Integer nstep = 0;
+ while (zVal < zStop) {
+ Double[] xyz = {0.0, 0.0, 0.0};
+ double s = HelixUtils.PathToXPlane(this, zVal, 1000.0, 1).get(0);
+ xyz[0] = HelixUtils.Direction(this, s).y();
+ xyz[1] = HelixUtils.Direction(this, s).z();
+ xyz[2] = zVal;
+ traj.put(nstep, xyz);
+ zVal += step;
+ nstep++;
+ }
+ return traj;
+ }
private void calculateParametersAtTarget() {
double pTot = this.p(bField);
@@ -88,6 +158,253 @@
_pDocaZ = VecOp.mult(pTot, _detToTrk.transformVectorToDetector(dirDocaZTrkSystem));
}
+// public Hep3Vector getPositionAtZ(double xFinal, double fringeHalfWidth, double step) {
+// double startFringeUpstream = -2 * fringeHalfWidth;
+// double stopFringeUpstream = 0;
+// if (_debug)
+// System.out.println(this.toString());
+//
+//
+// }
+ public Hep3Vector getPositionAtZ(double xFinal, double start, double stop, double step) {
+
+ double startFringe = start;
+ double stopFringe = stop;
+ // _debugForward = false;
+ // if (xFinal > 900)
+ // _debugForward = true;
+ // if looking upstream, we'll be propagating backwards
+ if (xFinal < 0) {
+ step = -step;
+ startFringe = stop;
+ stopFringe = start;
+ }
+ double fringeHalfWidth = Math.abs(stopFringe - startFringe) / 2;
+ double fringeMid = (stopFringe + startFringe) / 2;
+ if (_debugForward)
+ System.out.println(this.toString());
+
+ double sStartFringe = HelixUtils.PathToXPlane(this, startFringe, 1000.0, 1).get(0);
+ if (_debugForward)
+ System.out.println("path to end of fringe = " + sStartFringe + "; xFinal = " + xFinal);
+ double xtmp = startFringe;
+ double ytmp = HelixUtils.PointOnHelix(this, sStartFringe).y();
+ double ztmp = HelixUtils.PointOnHelix(this, sStartFringe).z();
+ double Rorig = this.R();
+ double xCtmp = this.xc();
+ double yCtmp = this.yc();
+ double q = Math.signum(this.curvature());
+ double phitmp = getPhi(xtmp, ytmp, xCtmp, yCtmp, q);
+ if (_debugForward) {
+ System.out.println("\nOriginal xtmp = " + xtmp
+ + "; ytmp = " + ytmp
+ + "; ztmp = " + ztmp
+ + "; phitmp = " + phitmp);
+
+ System.out.println("nOriginal Rorig = " + Rorig
+ + "; xCtmp = " + xCtmp
+ + "; yCtmp = " + yCtmp);
+ }
+ if (_debugForward)
+ System.out.println("Original Direction at Fringe: " + HelixUtils.Direction(this, startFringe).toString());
+ double Rtmp = Rorig;
+ // now start stepping through the fringe field
+ Hep3Vector r0Tmp = HelixUtils.PointOnHelix(this, sStartFringe);
+ SpacePoint r0 = new SpacePoint(r0Tmp);
+ double pTot = this.p(bField);
+ Hep3Vector dirOrig = HelixUtils.Direction(this, sStartFringe);
+ Hep3Vector p0 = VecOp.mult(pTot, dirOrig);
+ Hep3Vector dirTmp = dirOrig;
+ SpacePoint rTmp = r0;
+ Hep3Vector pTmp = p0;
+ double pXOrig = p0.x();
+ double pXTmp = pXOrig;
+ double totalS = sStartFringe;
+ if (_debugForward) {
+ double tmpdX = xFinal - startFringe;
+ Hep3Vector fooExt = extrapolateStraight(dirOrig, tmpdX);
+ Hep3Vector fooFinal = VecOp.add(fooExt, r0Tmp);
+ System.out.println("Extrpolating straight back from startFringe = (" + fooFinal.x() + "," + fooFinal.y() + "," + fooFinal.z() + ")");
+
+ }
+ //follow trajectory while: in fringe field, before end point, and we don't have a looper
+ while (Math.signum(step) * xtmp < Math.signum(step) * stopFringe && Math.signum(step) * xtmp < Math.signum(step) * xFinal && Math.signum(pXOrig * pXTmp) > 0) {
+ if (_debugForward) {
+ System.out.println("New step in Fringe Field");
+ System.out.println("rTmp = " + rTmp.toString());
+ System.out.println("pTmp = " + pTmp.toString());
+ System.out.println("OriginalHelix pos = " + HelixUtils.PointOnHelix(this, totalS));
+ System.out.println("OriginalHelix Momentum = " + VecOp.mult(pTot, HelixUtils.Direction(this, totalS)));
+ }
+
+ double fringeFactor = getFringe(Math.signum(step) * (fringeMid - rTmp.x()), fringeHalfWidth);
+// double myBField=bField * fringeFactor;
+ double myBField=FieldMap.getFieldFromMap(rTmp.x(),rTmp.y());
+ if (_debugForward) {
+ System.out.println("rTmp.x() = " + rTmp.x() + " field = " +myBField);
+ }
+ setTrack(pTmp, rTmp, q, myBField);
+ rTmp = _trajectory.getPointAtDistance(step);
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ xtmp = rTmp.x();
+ if (_debugForward) {
+ System.out.println("############## done... #############");
+
+ System.out.println("\n");
+ }
+ totalS += step;
+ }
+ //ok, done with field...extrapolate straight back...
+ Hep3Vector pointInTrking;
+ if (Math.signum(step) * xtmp < Math.signum(step) * xFinal) {
+ //get direction of the track before it hits fringe field
+ double deltaX = xFinal - xtmp;
+ Hep3Vector dir = _trajectory.getUnitTangentAtLength(0);
+// double deltaY = deltaX * dir.y();
+// double deltaZ = deltaX * dir.z();
+ Hep3Vector delta = extrapolateStraight(dir, deltaX);
+// double deltaZ = Math.sqrt(deltaX*deltaX+deltaY*deltaY)* dir.z();
+ pointInTrking = new BasicHep3Vector(xFinal, delta.y() + ytmp, delta.z() + ztmp);
+
+ if (_debugForward) {
+ System.out.println("Pointing straight forward from xtmp = " + xtmp
+ + "; ytmp = " + ytmp
+ + "; ztmp = " + ztmp
+ + "; deltaX= " + deltaX);
+ System.out.println("Directions: " + dir.toString());
+ System.out.println("Position at ECal: x = " + xFinal + "; y = " + pointInTrking.y() + "; z = " + pointInTrking.z());
+
+ }
+ } else { // still in the fringe field...just return the current position
+// pointInTrking = new BasicHep3Vector(xFinal, ytmp, ztmp);
+ pointInTrking = new BasicHep3Vector(rTmp.x(), rTmp.y(), rTmp.z());
+ }
+ return _detToTrk.transformVectorToDetector(pointInTrking);
+ }
+
+
+ public Hep3Vector getPositionAtZMap(double start, double xFinal, double step) {
+
+ double startFringe = start;
+
+ _debugForward = false;
+ if (xFinal > 900)
+ _debugForward = true;
+ // if looking upstream, we'll be propagating backwards
+ if (xFinal < 0) {
+ step = -step;
+ }
+ if (_debugForward)
+ System.out.println(this.toString());
+
+ double sStartFringe = HelixUtils.PathToXPlane(this, startFringe, 1000.0, 1).get(0);
+ if (_debugForward)
+ System.out.println("path to end of fringe = " + sStartFringe + "; xFinal = " + xFinal);
+ double xtmp = startFringe;
+ double ytmp = HelixUtils.PointOnHelix(this, sStartFringe).y();
+ double ztmp = HelixUtils.PointOnHelix(this, sStartFringe).z();
+ double Rorig = this.R();
+ double xCtmp = this.xc();
+ double yCtmp = this.yc();
+ double q = Math.signum(this.curvature());
+ double phitmp = getPhi(xtmp, ytmp, xCtmp, yCtmp, q);
+ if (_debugForward) {
+ System.out.println("\nOriginal xtmp = " + xtmp
+ + "; ytmp = " + ytmp
+ + "; ztmp = " + ztmp
+ + "; phitmp = " + phitmp);
+
+ System.out.println("nOriginal Rorig = " + Rorig
+ + "; xCtmp = " + xCtmp
+ + "; yCtmp = " + yCtmp);
+ }
+ if (_debugForward)
+ System.out.println("Original Direction at Fringe: " + HelixUtils.Direction(this, startFringe).toString());
+ double Rtmp = Rorig;
+ // now start stepping through the fringe field
+ Hep3Vector r0Tmp = HelixUtils.PointOnHelix(this, sStartFringe);
+ SpacePoint r0 = new SpacePoint(r0Tmp);
+ double pTot = this.p(bField);
+ Hep3Vector dirOrig = HelixUtils.Direction(this, sStartFringe);
+ Hep3Vector p0 = VecOp.mult(pTot, dirOrig);
+ Hep3Vector dirTmp = dirOrig;
+ SpacePoint rTmp = r0;
+ Hep3Vector pTmp = p0;
+ double pXOrig = p0.x();
+ double pXTmp = pXOrig;
+ double totalS = sStartFringe;
+ if (_debugForward) {
+ double tmpdX = xFinal - startFringe;
+ Hep3Vector fooExt = extrapolateStraight(dirOrig, tmpdX);
+ Hep3Vector fooFinal = VecOp.add(fooExt, r0Tmp);
+ System.out.println("Extrpolating straight back from startFringe = (" + fooFinal.x() + "," + fooFinal.y() + "," + fooFinal.z() + ")");
+
+ }
+ //follow trajectory while: in fringe field, before end point, and we don't have a looper
+ while (Math.signum(step) * xtmp < Math.signum(step) * xFinal && Math.signum(pXOrig * pXTmp) > 0) {
+ if (_debugForward) {
+ System.out.println("New step in Fringe Field");
+ System.out.println("rTmp = " + rTmp.toString());
+ System.out.println("pTmp = " + pTmp.toString());
+ System.out.println("OriginalHelix pos = " + HelixUtils.PointOnHelix(this, totalS));
+ System.out.println("OriginalHelix Momentum = " + VecOp.mult(pTot, HelixUtils.Direction(this, totalS)));
+ }
+
+ double myBField=FieldMap.getFieldFromMap(rTmp.x(),rTmp.y());
+ if (_debugForward) {
+ System.out.println("rTmp.x() = " + rTmp.x() + " field = " +myBField);
+ }
+ setTrack(pTmp, rTmp, q, myBField);
+ rTmp = _trajectory.getPointAtDistance(step);
+ pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
+ pXTmp = pTmp.x();
+ xtmp = rTmp.x();
+ if (_debugForward) {
+ System.out.println("############## done... #############");
+
+ System.out.println("\n");
+ }
+ totalS += step;
+ }
+ //ok, done with field.
+ Hep3Vector pointInTrking = new BasicHep3Vector(rTmp.x(), rTmp.y(), rTmp.z());
+ if (_debugForward) {
+ System.out.println("Position at ECal: x = " + xFinal + "; y = " + pointInTrking.y() + "; z = " + pointInTrking.z());
+ }
+ return _detToTrk.transformVectorToDetector(pointInTrking);
+ }
+
+ private double getPhi(double x, double y, double xc, double yc, double sign) {
+ // System.out.println("Math.atan2(y - yc, x - xc)="+Math.atan2(y - yc, x - xc));
+ return Math.atan2(y - yc, x - xc) - sign * Math.PI / 2;
+ }
+
+ private Hep3Vector extrapolateStraight(Hep3Vector dir, double deltaX) {
+ double deltaY = deltaX * dir.y();
+ double deltaZ = deltaX * dir.z();
+ return new BasicHep3Vector(deltaX, deltaY, deltaZ);
+ }
+
+ //field that changes linearly from Bmax->0
+ private double getFringe(double x, double halfWidth) {
+// System.out.println("x = " + x + "; halfWidth = " + halfWidth);
+ if (x / halfWidth > 1)
+ return 1;
+ if (x / halfWidth < -1)
+ return 0;
+
+ return (1.0 / 2.0) * (1 + x / halfWidth);
+ }
+
+ private Hep3Vector getDirection(double phi, double sign) {
+ double ux = Math.cos(phi) * this.sth();
+ double uy = Math.sin(phi) * this.sth();
+ double uz = this.cth();
+ // Return the direction unit vector
+ return new BasicHep3Vector(ux, uy, uz);
+ }
+
private double findPathToDocaZ() {
double step = 0.1;//100 um step size
double maxS = 100.0; //go to 10cm
@@ -103,4 +420,20 @@
}
return minS;
}
+
+ private void setTrack(Hep3Vector p0, SpacePoint r0, double q, double B) {
+ SpaceVector p = new CartesianVector(p0.v());
+ double phi = Math.atan2(p.y(), p.x());
+ double lambda = Math.atan2(p.z(), p.rxy());
+ double field = B * fieldConversion;
+
+ if (q != 0 && field != 0) {
+ double radius = p.rxy() / (q * field);
+ _trajectory = new Helix(r0, radius, phi, lambda);
+ } else {
+ _trajectory = new Line(r0, phi, lambda);
+ }
+ }
+
+
}