 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps on MAIN | |||
| recon/tracking/FieldMap.java | +175 | added 1.1 | |
| /HPSTrack.java | +338 | -5 | 1.1 -> 1.2 |
| monitoring/HPSCalibrationListener.java | +3 | -1 | 1.3 -> 1.4 |
| +516 | -6 | ||
Read in and use the field map
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; + + } +}
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);
+ }
+ }
+
+
}
diff -u -r1.3 -r1.4 --- HPSCalibrationListener.java 13 Jul 2012 23:05:54 -0000 1.3 +++ HPSCalibrationListener.java 23 Jul 2012 21:45:02 -0000 1.4 @@ -6,12 +6,13 @@
import java.util.TimeZone; import org.lcsim.conditions.ConditionsSet; import org.lcsim.hps.recon.ecal.HPSEcalConditions;
+import org.lcsim.hps.recon.tracking.FieldMap;
import org.lcsim.hps.recon.tracking.HPSSVTCalibrationConstants; /** * * @author Sho Uemura <[log in to unmask]>
- * @version $Id: HPSCalibrationListener.java,v 1.3 2012/07/13 23:05:54 omoreno Exp $
+ * @version $Id: HPSCalibrationListener.java,v 1.4 2012/07/23 21:45:02 mgraham Exp $
*/
public class HPSCalibrationListener implements EtEventListener {
@@ -20,6 +21,7 @@
// Load calibration constants and bad SVT channels found during QA
HPSSVTCalibrationConstants.loadCalibrationConstants(new Date(System.currentTimeMillis()), -1);
HPSEcalConditions.loadPedestals();
+ FieldMap.loadFieldMap();
}
@Override
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1