Author: [log in to unmask]
Date: Mon Jan 4 17:41:58 2016
New Revision: 4079
Log:
Fix bug in correction to perigee track parameters. Refactor old HPS helix class for this purpose. Attempt to restructure to pull out utility methods and try to delete old stuff, work in progress. Add relation between seed and GBL trackin the hold GBL refitter. Save a path length map for the gbl trajectory.
Added:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/HpsHelicalTrackFit.java
- copied, changed from r4077, java/trunk/tracking/src/main/java/org/hps/recon/tracking/HPSTrack.java
Removed:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/HPSTrack.java
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
Copied: java/trunk/tracking/src/main/java/org/hps/recon/tracking/HpsHelicalTrackFit.java (from r4077, java/trunk/tracking/src/main/java/org/hps/recon/tracking/HPSTrack.java)
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/HPSTrack.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/HpsHelicalTrackFit.java Mon Jan 4 17:41:58 2016
@@ -1,97 +1,42 @@
package org.hps.recon.tracking;
-import static org.lcsim.constants.Constants.fieldConversion;
import hep.physics.matrix.SymmetricMatrix;
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
import java.util.HashMap;
-import java.util.List;
import java.util.Map;
-import org.hps.util.Pair;
import org.lcsim.event.MCParticle;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.fit.helicaltrack.MultipleScatter;
-import org.lcsim.geometry.FieldMap;
-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;
/**
- * Class HPSTrack: extension of HelicalTrackFit to include HPS-specific
- * variables other useful
- * things.
+ * Extension of {@link HelicalTrackFit} to include HPS-specific information and utilities.
*
- * @author mgraham created on 6/27/2011
+ * @author mgraham <[log in to unmask]>
+ * @author phansson <[log in to unmask]>
+ *
*/
-// TODO: Check how much of the added behavior and information is really needed here.
-// It might be better to pull out a lot of the analytical functionality into a utility class.
-// FIXME: There are many instance variables marked as unused in my IDE.
-public class HPSTrack extends HelicalTrackFit {
- 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,
- // which is where the tracking returns it's parameters by default
- private Hep3Vector _pDocaY;
- private Hep3Vector _posDocaY;
- // the position & momentum of the track at the intersection of the target (z=0)
- private Hep3Vector _pTarget;
- private Hep3Vector _posTarget;
- // the position & momentum of the track at DOCA to the beam axis (z)
- private Hep3Vector _pDocaZ;
- private Hep3Vector _posDocaZ;
- 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;
- private MCParticle mcParticle;
+public class HpsHelicalTrackFit extends HelicalTrackFit {
- public HPSTrack(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf, Map<HelicalTrackHit, Double> smap, Map<HelicalTrackHit, MultipleScatter> msmap, BeamSpot beam) {
- super(pars, cov, chisq, ndf, smap, msmap);
- _beam = beam;
- calculateParametersAtTarget();
- calculateParametersAtDocaY();
- calculateParametersAtDocaZ();
+ private MCParticle mcParticle = null;
+ private double[] refPoint = {0.,0.};
+
+ public HpsHelicalTrackFit(HelicalTrackFit htf) {
+ super(htf.parameters(), htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
}
- public HPSTrack(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf, Map<HelicalTrackHit, Double> smap, Map<HelicalTrackHit, MultipleScatter> msmap) {
- super(pars, cov, chisq, ndf, smap, msmap);
- calculateParametersAtTarget();
- calculateParametersAtDocaY();
- calculateParametersAtDocaZ();
+ public HpsHelicalTrackFit(double[] parameters, SymmetricMatrix covariance, double[] chiSquared, int[] ndf, Map<HelicalTrackHit, Double> sMap, Map<HelicalTrackHit, MultipleScatter> msMap) {
+ super(parameters, covariance, chiSquared, ndf, sMap, msMap);
}
-
- public HPSTrack(HelicalTrackFit htf, BeamSpot beam) {
- super(htf.parameters(), htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
- _beam = beam;
- calculateParametersAtTarget();
- calculateParametersAtDocaY();
- calculateParametersAtDocaZ();
+
+ public HpsHelicalTrackFit(double[] parameters, SymmetricMatrix covariance, double[] chiSquared, int[] ndf, Map<HelicalTrackHit, Double> sMap, Map<HelicalTrackHit, MultipleScatter> msMap, double[] refPoint) {
+ super(parameters, covariance, chiSquared, ndf, sMap, msMap);
+ this.refPoint = refPoint;
}
-
- public HPSTrack(HelicalTrackFit htf) {
- super(htf.parameters(), htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
- calculateParametersAtTarget();
- calculateParametersAtDocaY();
- calculateParametersAtDocaZ();
- }
-
- public HPSTrack(double[] parameters, SymmetricMatrix covariance, double[] chiSquared, int[] ndf, Map<HelicalTrackHit, Double> sMap, Map<HelicalTrackHit, MultipleScatter> msMap, MCParticle mcParticle) {
- super(parameters, covariance, chiSquared, ndf, sMap, msMap);
-
- // Set the MC particle associated with this fit
- this.mcParticle = mcParticle;
- }
+
/**
* Get map of the the track trajectory within the uniform bfield
@@ -136,260 +81,25 @@
return traj;
}
- private void calculateParametersAtTarget() {
- double pTot = this.p(bField);
- // currently, PathToXPlane only returns a single path (no loopers!)
- List<Double> paths = HelixUtils.PathToXPlane(this, 0.0, 100.0, 1);
- Hep3Vector posTargetTrkSystem = HelixUtils.PointOnHelix(this, paths.get(0));
- Hep3Vector dirTargetTrkSystem = HelixUtils.Direction(this, paths.get(0));
- _posTarget = CoordinateTransformations.transformVectorToDetector(posTargetTrkSystem);
- _pTarget = VecOp.mult(pTot, CoordinateTransformations.transformVectorToDetector(dirTargetTrkSystem));
-
- }
-
- private void calculateParametersAtDocaY() {
- double pTot = this.p(bField);
- Hep3Vector posDocaYTrkSystem = HelixUtils.PointOnHelix(this, 0);
- Hep3Vector dirDocaYTrkSystem = HelixUtils.Direction(this, 0);
- _posDocaY = CoordinateTransformations.transformVectorToDetector(posDocaYTrkSystem);
- _pDocaY = VecOp.mult(pTot, CoordinateTransformations.transformVectorToDetector(dirDocaYTrkSystem));
-
- }
-
- private void calculateParametersAtDocaZ() {
- double pTot = this.p(bField);
- double sAtDocaZ = findPathToDocaZ();
- Hep3Vector posDocaZTrkSystem = HelixUtils.PointOnHelix(this, sAtDocaZ);
- Hep3Vector dirDocaZTrkSystem = HelixUtils.Direction(this, sAtDocaZ);
- _posDocaZ = CoordinateTransformations.transformVectorToDetector(posDocaZTrkSystem);
- _pDocaZ = VecOp.mult(pTot, CoordinateTransformations.transformVectorToDetector(dirDocaZTrkSystem));
- }
-
- public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step) {
- System.out.println("Don't use this anymore! Use the contructor with FieldMap provided");
- Hep3Vector out= new BasicHep3Vector(666,666,666);
- Hep3Vector[] ret={out};
- return ret;
- }
-
- public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step, boolean debugOK) {
- System.out.println("Don't use this anymore! Use the contructor with FieldMap provided");
- Hep3Vector out= new BasicHep3Vector(666,666,666);
- Hep3Vector[] ret={out};
- return ret;
- }
-
- /**
- * Get the position and direction on the track using B-field map for
- * extrapolation
- *
- * @param start = starting z-position of extrapolation
- * @param zFinal = final z-position
- * @param step = step size
- * @param bmap = the B-Field map
- * @return position[0] and direction[1] at Z=zfinal
- */
- public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step, FieldMap bmap) {
- return this.getPositionAtZMap(start, xFinal, step, bmap, false);
- }
-
- public Hep3Vector[] getPositionAtZMap(double start, double xFinal, double step, FieldMap bmap, boolean debugOk) {
-
- double startFringe = start;
-
- _debugForward = false;
- if (xFinal > 900)
- _debugForward = debugOk ? true : false;
- // 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 = bmap.getField(new BasicHep3Vector(rTmp.z(), rTmp.z(), rTmp.x())).z(); // note weird co ordinates for track vs det
- 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;
- }
-
- // Go with finer granularity in the end
- rTmp = _trajectory.getPointAtDistance(0);
- xtmp = rTmp.x();
- pTmp = VecOp.mult(pTot, _trajectory.getUnitTangentAtLength(step));
- pXTmp = pTmp.x();
- step = step / 10.0;
-
- 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 = bmap.getField(new BasicHep3Vector(rTmp.z(), rTmp.z(), rTmp.x())).z(); // note weird co ordinates for track vs det
-
- 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 xfinal (tracking) : x = " + xFinal + "; y = " + pointInTrking.y() + "; z = " + pointInTrking.z());
- Hep3Vector[] out = {CoordinateTransformations.transformVectorToDetector(pointInTrking), CoordinateTransformations.transformVectorToDetector(pTmp)};
-
- return out;
- }
-
- 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
- double s = -30;
- double minDist = 999999;
- double minS = s;
- double dist = 999998;
- // once the distance starts increasing, quit and return
- while (dist < minDist) {
- Hep3Vector pos = HelixUtils.PointOnHelix(this, s);
- dist = pos.y() * pos.y() + pos.z() * pos.z();
- s += step;
- }
- 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);
- }
-
- public Trajectory getTrajectory() {
- return this._trajectory;
- }
/**
* Get the MC Particle associated with the HelicalTrackFit
*
* @return mcParticle :
*/
- public MCParticle getMCParticle() {
+ public MCParticle getMcParticle() {
return this.mcParticle;
}
- /**
- * Set the MC Particle associated with the HelicalTrackFit
- *
- * @param mcParticle :
- */
- public void setMCParticle(MCParticle mcParticle) {
+ public void setMcParticle(MCParticle mcParticle) {
this.mcParticle = mcParticle;
}
+
+ public double[] getRefPoint() {
+ return refPoint;
+ }
+
+
+
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java Mon Jan 4 17:41:58 2016
@@ -6,6 +6,7 @@
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.SpacePoint;
import hep.physics.vec.VecOp;
+
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
@@ -15,10 +16,13 @@
import java.util.List;
import java.util.Map;
import java.util.Set;
+
import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.EventQuality.Quality;
import org.hps.recon.tracking.gbl.HelicalTrackStripGbl;
+
import static org.lcsim.constants.Constants.fieldConversion;
+
import org.lcsim.detector.ITransform3D;
import org.lcsim.detector.solids.Box;
import org.lcsim.detector.solids.Point3D;
@@ -86,6 +90,67 @@
return extrapolateHelixToXPlane(getHTF(track), x);
}
+ /**
+ * Change reference point of helix (following L3 Internal Note 1666.)
+ * @param newRefPoint - The new reference point in XY
+ */
+ public static double[] getParametersAtNewRefPoint(double[] newRefPoint, HpsHelicalTrackFit helicalTrackFit) {
+ return getParametersAtNewRefPoint(newRefPoint, helicalTrackFit.getRefPoint(),helicalTrackFit.parameters());
+ }
+
+ /**
+ * Change reference point of helix (following L3 Internal Note 1666.)
+ * @param newRefPoint - The new reference point in XY
+ */
+ public static double[] getParametersAtNewRefPoint(double[] newRefPoint, double[] __refPoint,
+ double[] parameters) {
+
+ double phi0 = parameters[HelicalTrackFit.phi0Index];
+ double curvature =parameters[HelicalTrackFit.curvatureIndex];
+ double dca = parameters[HelicalTrackFit.dcaIndex];
+ double slope = parameters[HelicalTrackFit.slopeIndex];
+ double z0 = parameters[HelicalTrackFit.z0Index];
+
+ //take care of phi0 range if needed (this matters for dphi below I think)
+ // L3 defines it in the range [-pi,pi]
+ if(phi0 > Math.PI)
+ phi0 -= Math.PI*2;
+
+ double dx = newRefPoint[0] - __refPoint[0];
+ double dy = newRefPoint[1] - __refPoint[1];
+ double sinphi = Math.sin(phi0);
+ double cosphi = Math.cos(phi0);
+ double R = 1.0/curvature;
+
+ // calculate new phi
+ double phinew = Math.atan2( sinphi - dx/(R-dca) , cosphi + dy/(R-dca) );
+
+ // difference in phi
+ // watch out for ambiguity
+ double dphi = phinew - phi0;
+ if (Math.abs( dphi ) > Math.PI)
+ throw new RuntimeException("dphi is large " + dphi + " from phi0 " + phi0
+ + " and phinew " + phinew + " take care of the ambiguity!!??");
+
+ // calculate new dca
+ double dcanew = dca + dx*sinphi - dy*cosphi + (dx*cosphi + dy*sinphi)*Math.tan( dphi/2. );
+
+ // path length from old to new point
+ double s = -1.0*dphi/curvature;
+
+ // new z0
+ double z0new = z0 + s*slope;
+
+ // new array
+ double[] params = new double[5];
+ params[HelicalTrackFit.phi0Index] = phinew;
+ params[HelicalTrackFit.curvatureIndex] = curvature;
+ params[HelicalTrackFit.dcaIndex] = dcanew;
+ params[HelicalTrackFit.slopeIndex] = slope;
+ params[HelicalTrackFit.z0Index] = z0new;
+ return params;
+ }
+
/**
* Extrapolate helix to a position along the x-axis. Re-use HelixUtils.
@@ -801,11 +866,11 @@
/**
- * Transform MCParticle into a Helix object. Note that it produces the helix
- * parameters at nominal x=0 and assumes that there is no field at x<0
- *
- * @param mcp MC particle to be transformed
- * @param org origin to be used for the track
+ * Transform MCParticle into a {@link HelicalTrackFit} object. Note that it produces the {@link HelicalTrackFit}
+ * parameters at nominal reference point at origin and assumes that there is no field at x<0
+ *
+ * @param mcp - MC particle to be transformed
+ * @param origin - origin to be used for the track
* @return {@link HelicalTrackFit} object based on the MC particle
*/
public static HelicalTrackFit getHTF(MCParticle mcp, Hep3Vector origin, double Bz) {
@@ -871,27 +936,7 @@
return htf;
}
- public static StraightLineTrack findSLTAtZ(Track trk1, double zVal, boolean useFringe) {
- SeedTrack s1 = (SeedTrack) trk1;
- HelicalTrackFit htf1 = s1.getSeedCandidate().getHelix();
- HPSTrack hpstrk1 = new HPSTrack(htf1);
- Hep3Vector pos1;
- if (useFringe) {
- // broken because you need ot provide the Field Map to get this...
-// pos1 = hpstrk1.getPositionAtZMap(100.0, zVal, 5.0)[0];
- } else
- pos1 = TrackUtils.extrapolateTrack(trk1, zVal);
- // System.out.printf("%s: Position1 at edge of fringe %s\n",this.getClass().getSimpleName(),pos1.toString());
- Helix traj = (Helix) hpstrk1.getTrajectory();
- if (traj == null) {
- SpacePoint r0 = new SpacePoint(HelixUtils.PointOnHelix(htf1, 0));
- traj = new Helix(r0, htf1.R(), htf1.phi0(), Math.atan(htf1.slope()));
- }
- HelixConverter converter = new HelixConverter(0.);
- StraightLineTrack slt1 = converter.Convert(traj);
- // System.out.printf("%s: straight line track: x0=%f,y0=%f,z0=%f dz/dx=%f dydx=%f targetY=%f targetZ=%f \n",this.getClass().getSimpleName(),slt1.x0(),slt1.y0(),slt1.z0(),slt1.dzdx(),slt1.dydx(),slt1.TargetYZ()[0],slt1.TargetYZ()[1]);
- return slt1;
- }
+
public static MCParticle getMatchedTruthParticle(Track track) {
boolean debug = false;
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java Mon Jan 4 17:41:58 2016
@@ -1,5 +1,12 @@
package org.hps.recon.tracking.gbl;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.Map;
+
+import org.hps.recon.tracking.gbl.FittedGblTrajectory.GBLPOINT;
+import org.hps.recon.tracking.gbl.matrix.SymMatrix;
+import org.hps.recon.tracking.gbl.matrix.Vector;
import org.lcsim.event.Track;
/**
@@ -9,6 +16,7 @@
*
*/
public class FittedGblTrajectory {
+
public enum GBLPOINT {
IP(0), LAST(1), VERTEX(2);
private int numVal;
@@ -34,6 +42,7 @@
}
}
+
public static enum GBLPARIDX {
QOVERP(0),YTPRIME(1),XTPRIME(2),XT(3),YT(4);
private int _value;
@@ -44,18 +53,86 @@
return _value;
}
};
+
private GblTrajectory _traj;
private double _chi2;
private double _lost;
private int _ndf;
private Track _seed = null;
- private GBLTrackData _t = null;
+ private Map<Integer, Double> pathLengthMap = null;
+
+ /**
+ * Default constructor.
+ *
+ * @param traj
+ * @param chi2
+ * @param ndf
+ * @param lost
+ */
public FittedGblTrajectory(GblTrajectory traj, double chi2, int ndf, double lost) {
_traj = traj;
_chi2 = chi2;
_ndf = ndf;
_lost = lost;
}
+
+ /**
+ * Find the index (or label) of the GBL point on the trajectory from the {@link GBLPOINT}.
+ * @param point
+ * @return
+ */
+ public int getPointIndex(GBLPOINT point) {
+ int gblPointIndex;
+ if (point.compareTo(GBLPOINT.IP) == 0)
+ gblPointIndex = 1;
+ else if (point.compareTo(GBLPOINT.LAST) == 0)
+ gblPointIndex = _traj.getNumPoints();
+ else
+ throw new RuntimeException("This GBL point " + point.toString() + "( " + point.name() + ") is not valid");
+ return gblPointIndex;
+ }
+
+
+ /**
+ * Find the corrections and covariance matrix for a particular {@link GBLPOINT}
+ * @param point
+ * @param locPar
+ * @param locCov
+ */
+ public void getResults(GBLPOINT point, Vector locPar, SymMatrix locCov) {
+
+ // find the GBL point index
+ int gblPointIndex = getPointIndex(point);
+
+ // Get the result from the trajectory
+ int ok = _traj.getResults(gblPointIndex, locPar, locCov);
+
+ // check that the fit was ok
+ if( ok != 0)
+ throw new RuntimeException("Trying to extract GBL corrections for fit that failed!?");
+ }
+
+ /**
+ * Find the path length to this point.
+ * @param point - {@link GBLPOINT} point
+ * @return path length
+ */
+ public double getPathLength(GBLPOINT point) {
+ int gblPointIndex = getPointIndex(point);
+ return getPathLength(gblPointIndex);
+ }
+
+ /**
+ * Find the path length to this point.
+ * @param iLabel - GBL point index
+ * @return path length
+ */
+ public double getPathLength(int iLabel) {
+ if( !this.pathLengthMap.containsKey(iLabel) )
+ throw new RuntimeException("This iLabel " + iLabel + " doesn't exists in the path length map.");
+ return this.pathLengthMap.get(iLabel);
+ }
+
public void set_seed(Track seed) {
_seed = seed;
}
@@ -74,5 +151,18 @@
public int get_ndf() {
return _ndf;
}
+
+ public void setPathLengthMap(Map<Integer, Double> pathLengthMap) {
+ this.pathLengthMap = pathLengthMap;
+ }
+
+ public Map<Integer, Double> getPathLengthMap() {
+ if (this.pathLengthMap == null)
+ throw new RuntimeException("No path length map has been set on this trajectory!");
+ return this.pathLengthMap;
+ }
+
+
+
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLFileIO.java Mon Jan 4 17:41:58 2016
@@ -11,8 +11,8 @@
import java.util.logging.Level;
import java.util.logging.Logger;
-import org.hps.recon.tracking.gbl.GBLOutput.ClParams;
-import org.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
+import org.hps.recon.tracking.gbl.GblUtils.ClParams;
+import org.hps.recon.tracking.gbl.GblUtils.PerigeeParams;
import org.hps.svt.alignment.RunAlignment;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
@@ -67,34 +67,34 @@
addLine(String.format("Track perPar (R phi0 slope d0 z0) %.12f %.12f %.12f %.12f %.12f",htf.R(),htf.phi0(),htf.slope(),htf.dca(),htf.z0()));
}
- String getPerTrackParamStr(PerigeeParams perPar) {
+ String getPerTrackParamStr(GblUtils.PerigeeParams perPar) {
return String.format("Track perPar (R theta phi d0 z0) %.12f %.12f %.12f %.12f %.12f",1.0/perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
}
- void printPerTrackParam(PerigeeParams perPar) {
+ void printPerTrackParam(GblUtils.PerigeeParams perPar) {
addLine(this.getPerTrackParamStr(perPar));
}
- String getPerTrackParamTruthStr(PerigeeParams perPar) {
+ String getPerTrackParamTruthStr(GblUtils.PerigeeParams perPar) {
return String.format("Truth perPar (kappa theta phi d0 z0) %.12f %.12f %.12f %.12f %.12f",perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
}
- void printPerTrackParamTruth(PerigeeParams perPar) {
+ void printPerTrackParamTruth(GblUtils.PerigeeParams perPar) {
addLine(this.getPerTrackParamTruthStr(perPar));
}
- String getClTrackParamTruthStr(ClParams perPar) {
+ String getClTrackParamTruthStr(GblUtils.ClParams perPar) {
return String.format("Truth clPar (q/p lambda phi xT yT) %.12f %.12f %.12f %.12f %.12f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
}
- void printClTrackParamTruth(ClParams perPar) {
+ void printClTrackParamTruth(GblUtils.ClParams perPar) {
addLine(this.getClTrackParamTruthStr(perPar));
}
- String getClTrackParamStr(ClParams perPar) {
+ String getClTrackParamStr(GblUtils.ClParams perPar) {
return String.format("Track clPar (q/p lambda phi xT yT) %.12f %.12f %.12f %.12f %.12f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
}
- void printClTrackParam(ClParams perPar) {
+ void printClTrackParam(GblUtils.ClParams perPar) {
addLine(String.format("%s",this.getClTrackParamStr(perPar)));
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java Mon Jan 4 17:41:58 2016
@@ -169,7 +169,7 @@
if (isMC) {
// find the truth particle for this track
- mcp = getMatchedTruthParticle(trk);
+ mcp = TrackUtils.getMatchedTruthParticle(trk);
// check if this is an A' event
for(MCParticle part : mcParticles) {
@@ -222,8 +222,8 @@
// Use the truth helix as the initial track for GBL?
//htf = htfTruth;
// Get perigee parameters to curvilinear frame
- PerigeeParams perPar = new PerigeeParams(htf, bFieldVector.z());
- PerigeeParams perParTruth = new PerigeeParams(htfTruth, bFieldVector.z());
+ GblUtils.PerigeeParams perPar = new GblUtils.PerigeeParams(htf, bFieldVector.z());
+ GblUtils.PerigeeParams perParTruth = new GblUtils.PerigeeParams(htfTruth, bFieldVector.z());
//GBLDATA
gtd.setPerigeeTrackParameters(perPar);
@@ -234,8 +234,8 @@
// Get curvilinear parameters
if (textFile != null) {
- ClParams clPar = new ClParams(htf, bFieldVector.z());
- ClParams clParTruth = new ClParams(htfTruth, bFieldVector.z());
+ GblUtils.ClParams clPar = new GblUtils.ClParams(htf, bFieldVector.z());
+ GblUtils.ClParams clParTruth = new GblUtils.ClParams(htfTruth, bFieldVector.z());
textFile.printClTrackParam(clPar);
textFile.printClTrackParamTruth(clParTruth);
@@ -246,7 +246,7 @@
}
// find the projection from the I,J,K to U,V,T curvilinear coordinates
- Hep3Matrix perToClPrj = getPerToClPrj(htf);
+ Hep3Matrix perToClPrj = GblUtils.getPerToClPrj(htf);
//GBLDATA
for (int row = 0; row < perToClPrj.getNRows(); ++row) {
@@ -780,171 +780,10 @@
}
}
- MCParticle getMatchedTruthParticle(Track track) {
- boolean debug = false;
-
- Map<MCParticle, Integer> particlesOnTrack = new HashMap<MCParticle, Integer>();
-
- if (debug) {
- System.out.printf("getmatched mc particle from %d tracker hits on the track \n", track.getTrackerHits().size());
- }
-
- for (TrackerHit hit : track.getTrackerHits()) {
- List<MCParticle> mcps = ((HelicalTrackHit) hit).getMCParticles();
- if (mcps == null) {
- System.out.printf("%s: warning, this hit (layer %d pos=%s) has no mc particles.\n", this.getClass().getSimpleName(), ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString());
- } else {
- if (debug) {
- System.out.printf("%s: this hit (layer %d pos=%s) has %d mc particles.\n", this.getClass().getSimpleName(), ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString(), mcps.size());
- }
- for (MCParticle mcp : mcps) {
- if (!particlesOnTrack.containsKey(mcp)) {
- particlesOnTrack.put(mcp, 0);
- }
- int c = particlesOnTrack.get(mcp);
- particlesOnTrack.put(mcp, c + 1);
- }
- }
- }
- if (debug) {
- System.out.printf("Track p=[ %f, %f, %f] \n", track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[1]);
- System.out.printf("FOund %d particles\n", particlesOnTrack.size());
- for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
- System.out.printf("%d hits assigned to %d p=%s \n", entry.getValue(), entry.getKey().getPDGID(), entry.getKey().getMomentum().toString());
- }
- }
- Map.Entry<MCParticle, Integer> maxEntry = null;
- for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
- if (maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0) {
- maxEntry = entry; //if ( maxEntry != null ) {
- } // if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue;
- } //}
- //maxEntry = entry;
- if (debug) {
- if (maxEntry != null) {
- System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
- maxEntry.getKey().getPDGID(), maxEntry.getKey().getMomentum().toString(),
- track.getCharge(), track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[2]);
- } else {
- System.out.printf("No truth particle found on this track\n");
- }
- }
- return maxEntry == null ? null : maxEntry.getKey();
- }
-
-// private BasicMatrix getJacPerToCl(HelicalTrackFit htf) {
-// System.out.printf("%s: getJacPerToCl\n", this.getClass().getSimpleName());
-// //use propoerly normalized B-field
-// Hep3Vector Bnorm = VecOp.mult(Constants.fieldConversion, _B);
-// //init jacobian to zero
-// BasicMatrix j = new BasicMatrix(5, 5);
-// initZero(j);
-// double lambda = Math.atan(htf.slope());
-// double q = Math.signum(htf.R());
-// double theta = Math.PI / 2.0 - lambda;
-// Hep3Vector T = HelixUtils.Direction(htf, 0.);
-// Hep3Vector p = VecOp.mult(htf.p(Math.abs(_B.z())), T);
-// double pT = htf.pT(Math.abs(_B.z()));
-// Hep3Vector H = VecOp.mult(1. / (Bnorm.magnitude()), Bnorm);
-// Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
-// Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
-// Hep3Vector U = VecOp.mult(-1, J);
-// Hep3Vector V = VecOp.cross(T, U);
-// double alpha = VecOp.cross(H, T).magnitude();
-// Hep3Vector N = VecOp.mult(1. / alpha, VecOp.cross(H, T));
-// Hep3Vector K = Z;
-// double Q = -Bnorm.magnitude() * q / p.magnitude();
-// double kappa = -1.0 * q * Bnorm.z() / pT;
-//
-// if (this._debug != 0) {
-// System.out.printf("%s: Bnorm=%s mag(Bnorm)=%f\n", this.getClass().getSimpleName(), Bnorm.toString(), Bnorm.magnitude());
-// System.out.printf("%s: p=%s |p|=%f pT=%f\n", this.getClass().getSimpleName(), p.toString(), p.magnitude(), pT);
-// System.out.printf("%s: q=%f\n", this.getClass().getSimpleName(), q);
-// System.out.printf("%s: q/p=%f\n", this.getClass().getSimpleName(), q / p.magnitude());
-// System.out.printf("%s: T=%s\n", this.getClass().getSimpleName(), T.toString());
-// System.out.printf("%s: H=%s\n", this.getClass().getSimpleName(), H.toString());
-// System.out.printf("%s: kappa=%f\n", this.getClass().getSimpleName(), kappa);
-// System.out.printf("%s: alpha=%f Q=%f \n", this.getClass().getSimpleName(), alpha, Q);
-// System.out.printf("%s: J=%s \n", this.getClass().getSimpleName(), J.toString());
-// System.out.printf("%s: V=%s \n", this.getClass().getSimpleName(), V.toString());
-// System.out.printf("%s: N=%s \n", this.getClass().getSimpleName(), N.toString());
-// System.out.printf("%s: TdotJ=%f \n", this.getClass().getSimpleName(), VecOp.dot(T, J));
-// System.out.printf("%s: VdotN=%f \n", this.getClass().getSimpleName(), VecOp.dot(V, N));
-// System.out.printf("%s: TdotK=%f \n", this.getClass().getSimpleName(), VecOp.dot(T, K));
-// System.out.printf("%s: UdotN=%f \n", this.getClass().getSimpleName(), VecOp.dot(U, N));
-// }
-//
-// j.setElement(0, 0, -1.0 * Math.sin(theta) / Bnorm.z());
-//
-// j.setElement(0, 1, q / (p.magnitude() * Math.tan(theta)));
-//
-// j.setElement(1, 1, -1);
-//
-// j.setElement(1, 3, -alpha * Q * VecOp.dot(T, J) * VecOp.dot(V, N));
-//
-// j.setElement(1, 4, -alpha * Q * VecOp.dot(T, K) * VecOp.dot(V, N));
-//
-// j.setElement(2, 2, 1);
-//
-// j.setElement(2, 3, -alpha * Q * VecOp.dot(T, J) * VecOp.dot(U, N) / Math.cos(lambda));
-//
-// j.setElement(2, 4, -alpha * Q * VecOp.dot(T, K) * VecOp.dot(U, N) / Math.cos(lambda));
-//
-// j.setElement(3, 3, -1);
-//
-// j.setElement(4, 4, VecOp.dot(V, K));
-//
-// if (_debug > 0) {
-// System.out.printf("%s: lambda= J(1,1)=%f * theta + J(1,3)=%f * eps + J(1,4)=%f * z0 \n",
-// this.getClass().getSimpleName(),
-// j.e(1, 1), j.e(1, 3), j.e(1, 4));
-//
-// }
-//
-// return j;
-//
-// }
- /**
- * Transform MCParticle into a Helix object. Note that it produces the helix
- * parameters at nominal x=0 and assumes that there is no field at x<0
- *
- * @param mcp MC particle to be transformed
- * @return helix object based on the MC particle
- */
-// private HelicalTrackFit getHTF(MCParticle mcp) {
-// Hep3Vector org = this._hpstrans.transformVectorToTracking(mcp.getOrigin());
-// Hep3Vector p = this._hpstrans.transformVectorToTracking(mcp.getMomentum());
-// // Move to x=0 if needed
-// if(org.x() < 0.) {
-// double dydx = p.y()/p.x();
-// double dzdx = p.z()/p.x();
-// double delta_x = -1. * org.x();
-// double y = delta_x * dydx;
-// double z = delta_x * dzdx;
-// double x = org.x() + delta_x;
-// if( Math.abs(x) > 1e-8) throw new RuntimeException("Error: origin is not zero!");
-// Hep3Vector old = org;
-// org = new BasicHep3Vector(x,y,z);
-// System.out.printf("org %s p %s -> org %s\n", old.toString(),p.toString(),org.toString());
-// } else {
-// org = this._hpstrans.transformVectorToTracking(mcp.getOrigin());
-// }
-//
-//
-//
-// HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1*((int)mcp.getCharge()), -1.0*this._B.z());
-// double par[] = new double[5];
-// par[HelicalTrackFit.dcaIndex] = helixParamCalculator.getDCA();
-// par[HelicalTrackFit.slopeIndex] = helixParamCalculator.getSlopeSZPlane();
-// par[HelicalTrackFit.phi0Index] = helixParamCalculator.getPhi0();
-// par[HelicalTrackFit.curvatureIndex] = 1.0/helixParamCalculator.getRadius();
-// par[HelicalTrackFit.z0Index] = helixParamCalculator.getZ0();
-// SymmetricMatrix cov = new SymmetricMatrix(5);
-// for(int i=0;i<cov.getNRows();++i) cov.setElement(i, i, 1.);
-// HelicalTrackFit htf = new HelicalTrackFit(par, cov, new double[2], new int[2], null, null);
-// return htf;
-// }
- private double truthTrackFitChi2(PerigeeParams perPar, PerigeeParams perParTruth, SymmetricMatrix covariance) {
+
+
+
+ private double truthTrackFitChi2(GblUtils.PerigeeParams perPar, GblUtils.PerigeeParams perParTruth, SymmetricMatrix covariance) {
//re-shuffle the param vector to match the covariance order of parameters
BasicMatrix p = new BasicMatrix(1, 5);
p.setElement(0, 0, perPar.getD0());
@@ -1005,190 +844,6 @@
return Math.sqrt(Math.pow(E1 + E2, 2) - VecOp.add(p1vec, p2vec).magnitudeSquared());
}
- private static BasicMatrix getPerParVector(double kappa, double theta, double phi, double d0, double z0) {
- BasicMatrix perPar = new BasicMatrix(1, 5);
- perPar.setElement(0, 0, kappa);
- perPar.setElement(0, 1, theta);
- perPar.setElement(0, 2, phi);
- perPar.setElement(0, 3, d0);
- perPar.setElement(0, 4, z0);
- return perPar;
- }
-
- private static BasicMatrix getPerParVector(HelicalTrackFit htf, double B) {
- if (htf != null) {
- double kappa = -1.0 * Math.signum(B) / htf.R();
- double theta = Math.PI / 2.0 - Math.atan(htf.slope());
- return getPerParVector(kappa, theta, htf.phi0(), htf.dca(), htf.z0());
- }
- return new BasicMatrix(1, 5);
- }
-
- /**
- *
- * Store perigee track parameters.
- *
- * @author Per Hansson Adrian <[log in to unmask]>
- *
- */
- public static class PerigeeParams {
-
- private final BasicMatrix _params;
-
- public PerigeeParams(HelicalTrackFit htf, double B) {
- _params = getPerParVector(htf, B);
- }
-
- public PerigeeParams(double kappa, double theta, double phi, double d0, double z0) {
- this._params = getPerParVector(kappa, theta, phi, d0, z0);
- }
-
- public BasicMatrix getParams() {
- return _params;
- }
-
- public double getKappa() {
- return _params.e(0, 0);
- }
-
- public double getTheta() {
- return _params.e(0, 1);
- }
-
- public double getPhi() {
- return _params.e(0, 2);
- }
-
- public double getD0() {
- return _params.e(0, 3);
- }
-
- public double getZ0() {
- return _params.e(0, 4);
- }
- }
-
- /**
- * Computes the projection matrix from the perigee XY plane variables dca
- * and z0 into the curvilinear xT,yT,zT frame (U,V,T)
- *
- * @param htf input helix to find the track direction
- * @return 3x3 projection matrix
- */
- static Hep3Matrix getPerToClPrj(HelicalTrackFit htf) {
- Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
- Hep3Vector T = HelixUtils.Direction(htf, 0.);
- Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
- Hep3Vector K = Z;
- Hep3Vector U = VecOp.mult(-1, J);
- Hep3Vector V = VecOp.cross(T, U);
- Hep3Vector I = VecOp.cross(J, K);
-
- BasicHep3Matrix trans = new BasicHep3Matrix();
- trans.setElement(0, 0, VecOp.dot(I, U));
- trans.setElement(0, 1, VecOp.dot(J, U));
- trans.setElement(0, 2, VecOp.dot(K, U));
- trans.setElement(1, 0, VecOp.dot(I, V));
- trans.setElement(1, 1, VecOp.dot(J, V));
- trans.setElement(1, 2, VecOp.dot(K, V));
- trans.setElement(2, 0, VecOp.dot(I, T));
- trans.setElement(2, 1, VecOp.dot(J, T));
- trans.setElement(2, 2, VecOp.dot(K, T));
- return trans;
-
- /*
- Hep3Vector B = new BasicHep3Vector(0, 0, 1); // TODO sign convention?
- Hep3Vector H = VecOp.mult(1 / bfield, B);
- Hep3Vector T = HelixUtils.Direction(helix, 0.);
- Hep3Vector HcrossT = VecOp.cross(H, T);
- double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
- double Q = Math.abs(bfield) * q / p;
- Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
- Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
- Hep3Vector K = Z;
- Hep3Vector U = VecOp.mult(-1, J);
- Hep3Vector V = VecOp.cross(T, U);
- Hep3Vector I = VecOp.cross(J, K);
- Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T)); //-cross(T,H)/alpha = -cross(T,Z) = -J
- double UdotI = VecOp.dot(U, I); // 0,0
- double NdotV = VecOp.dot(N, V); // 1,1?
- double NdotU = VecOp.dot(N, U); // 0,1?
- double TdotI = VecOp.dot(T, I); // 2,0
- double VdotI = VecOp.dot(V, I); // 1,0
- double VdotK = VecOp.dot(V, K); // 1,2
- */
- }
-
-
- /**
- *
- * Store curvilinear track parameters.
- *
- * @author Per Hansson Adrian <[log in to unmask]>
- *
- */
- public static class ClParams {
-
- private BasicMatrix _params = new BasicMatrix(1, 5);
-
- public ClParams(HelicalTrackFit htf, double B) {
-
- if (htf == null) {
- return;
- }
-
- Hep3Matrix perToClPrj = getPerToClPrj(htf);
-
- double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
- double z0 = htf.z0();
- Hep3Vector vecPer = new BasicHep3Vector(0., d0, z0);
- //System.out.printf("%s: vecPer=%s\n",this.getClass().getSimpleName(),vecPer.toString());
-
- Hep3Vector vecCl = VecOp.mult(perToClPrj, vecPer);
- //System.out.printf("%s: vecCl=%s\n",this.getClass().getSimpleName(),vecCl.toString());
- double xT = vecCl.x();
- double yT = vecCl.y();
- //double zT = vecCl.z();
-
- double lambda = Math.atan(htf.slope());
- double q = Math.signum(htf.R());
- double qOverP = q / htf.p(Math.abs(B));
- double phi = htf.phi0();
-
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), qOverP);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), lambda);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), phi);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XT.getValue(), xT);
- _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YT.getValue(), yT);
- }
-
- public BasicMatrix getParams() {
- return _params;
- }
-
- double getQoverP() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
- }
-
- double getLambda() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
- }
-
- double getPhi() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
- }
-
- double getXt() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.XT.getValue());
- }
-
- double getYt() {
- return _params.e(0, FittedGblTrajectory.GBLPARIDX.YT.getValue());
- }
-
- }
-
-
/**
*
* {@link HelicalTrackStripGbl} that explicitly uses the given unit vectors when accessed.
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java Mon Jan 4 17:41:58 2016
@@ -3,7 +3,7 @@
import hep.physics.vec.BasicHep3Matrix;
import hep.physics.vec.Hep3Matrix;
-import org.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
+import org.hps.recon.tracking.gbl.GblUtils.PerigeeParams;
import org.lcsim.event.GenericObject;
/**
@@ -77,7 +77,7 @@
/**
* @param perPar is the perigee parameters that is added to object
*/
- public void setPerigeeTrackParameters(PerigeeParams perPar) {
+ public void setPerigeeTrackParameters(GblUtils.PerigeeParams perPar) {
this.bank_double[GBLDOUBLE.PERKAPPA] = perPar.getKappa();
this.bank_double[GBLDOUBLE.PERTHETA] = perPar.getTheta();
this.bank_double[GBLDOUBLE.PERPHI] = perPar.getPhi();
@@ -85,8 +85,8 @@
this.bank_double[GBLDOUBLE.PERZ0] = perPar.getZ0();
}
- public PerigeeParams getPerigeeTrackParameters() {
- return new PerigeeParams(this.bank_double[GBLDOUBLE.PERKAPPA],
+ public GblUtils.PerigeeParams getPerigeeTrackParameters() {
+ return new GblUtils.PerigeeParams(this.bank_double[GBLDOUBLE.PERKAPPA],
this.bank_double[GBLDOUBLE.PERTHETA],
this.bank_double[GBLDOUBLE.PERPHI],
this.bank_double[GBLDOUBLE.PERD0],
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java Mon Jan 4 17:41:58 2016
@@ -1,23 +1,234 @@
package org.hps.recon.tracking.gbl;
+import java.util.logging.Logger;
+
import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Matrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Matrix;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.MaterialSupervisor;
import org.hps.recon.tracking.MultipleScattering;
+import org.hps.recon.tracking.gbl.matrix.Matrix;
+import org.hps.recon.tracking.gbl.matrix.Vector;
+import org.lcsim.constants.Constants;
import org.lcsim.detector.IDetectorElement;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.recon.tracking.seedtracker.ScatterAngle;
/**
- * A class providing various utilities related to GBL
+ * A class with only static utilities related to GBL
*
* @author Per Hansson Adrian <[log in to unmask]>
*
*/
public class GblUtils {
-
+
+ public static Logger LOGGER = Logger.getLogger(GblUtils.class.getName());
+
+
+
+
+ /**
+ * Private constructor to avoid instantiation.
+ */
private GblUtils() {
}
+
+
+ /**
+ *
+ * Store local curvilinear track parameters.
+ *
+ * @author Per Hansson Adrian <[log in to unmask]>
+ *
+ */
+ public static class ClParams {
+
+ private BasicMatrix _params = new BasicMatrix(1, 5);
+
+ public ClParams(HelicalTrackFit htf, double B) {
+
+ if (htf == null) {
+ return;
+ }
+
+ Hep3Matrix perToClPrj = getPerToClPrj(htf);
+
+ double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
+ double z0 = htf.z0();
+ Hep3Vector vecPer = new BasicHep3Vector(0., d0, z0);
+ //System.out.printf("%s: vecPer=%s\n",this.getClass().getSimpleName(),vecPer.toString());
+
+ Hep3Vector vecCl = VecOp.mult(perToClPrj, vecPer);
+ //System.out.printf("%s: vecCl=%s\n",this.getClass().getSimpleName(),vecCl.toString());
+ double xT = vecCl.x();
+ double yT = vecCl.y();
+ //double zT = vecCl.z();
+
+ double lambda = Math.atan(htf.slope());
+ double q = Math.signum(htf.R());
+ double qOverP = q / htf.p(Math.abs(B));
+ double phi = htf.phi0();
+
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), qOverP);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), lambda);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), phi);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.XT.getValue(), xT);
+ _params.setElement(0, FittedGblTrajectory.GBLPARIDX.YT.getValue(), yT);
+ }
+
+ public BasicMatrix getParams() {
+ return _params;
+ }
+
+ double getQoverP() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
+ }
+
+ double getLambda() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
+ }
+
+ double getPhi() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
+ }
+
+ double getXt() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.XT.getValue());
+ }
+
+ double getYt() {
+ return _params.e(0, FittedGblTrajectory.GBLPARIDX.YT.getValue());
+ }
+
+ }
+
+
+
+
+
+ /**
+ *
+ * Store perigee track parameters.
+ *
+ * @author Per Hansson Adrian <[log in to unmask]>
+ *
+ */
+ public static class PerigeeParams {
+
+ private final BasicMatrix _params;
+
+ public PerigeeParams(HelicalTrackFit htf, double B) {
+ _params = GblUtils.getPerParVector(htf, B);
+ }
+
+ public PerigeeParams(double kappa, double theta, double phi, double d0, double z0) {
+ this._params = GblUtils.getPerParVector(kappa, theta, phi, d0, z0);
+ }
+
+ public BasicMatrix getParams() {
+ return _params;
+ }
+
+ public double getKappa() {
+ return _params.e(0, 0);
+ }
+
+ public double getTheta() {
+ return _params.e(0, 1);
+ }
+
+ public double getPhi() {
+ return _params.e(0, 2);
+ }
+
+ public double getD0() {
+ return _params.e(0, 3);
+ }
+
+ public double getZ0() {
+ return _params.e(0, 4);
+ }
+ }
+
+
+
+
+
+
+ /**
+ * Get corrected perigee parameters.
+ * @param locPar - GBL local curvilinear corrections
+ * @param helicalTrackFit - helix
+ * @param bfield - B-field strength
+ * @return corrected parameters
+ */
+ public static double[] getCorrectedPerigeeParameters(Vector locPar, HelicalTrackFit helicalTrackFit, double bfield) {
+
+
+ // Explicitly assign corrections to local variables
+ double qOverPCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
+ double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
+ double yTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
+ double xTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XT.getValue());
+ double yTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YT.getValue());
+
+
+ // Get helix parameters
+ double qOverP = helicalTrackFit.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helicalTrackFit.slope(), 2)));
+ double d0 = -1.0 * helicalTrackFit.dca(); // correct for different sign convention of d0 in perigee frame
+ double z0 = helicalTrackFit.z0();
+ double phi0 = helicalTrackFit.phi0();
+ double lambda = Math.atan(helicalTrackFit.slope());
+
+ // calculate new d0 and z0
+ Hep3Matrix perToClPrj = GblUtils.getPerToClPrj(helicalTrackFit);
+
+ Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
+ Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
+
+ //d0
+ double d0_corr = corrPer.y();
+ double dca_gbl = -1.0 * (d0 + d0_corr);
+
+ //z0
+ double z0_corr = corrPer.z();
+ double z0_gbl = z0 + z0_corr;
+
+ //calculate new slope
+ double lambda_gbl = lambda + yTPrimeCorr;
+ double slope_gbl = Math.tan(lambda_gbl);
+
+ // calculate new curvature
+ double qOverP_gbl = qOverP + qOverPCorr;
+ double C_gbl = Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl / Math.cos(lambda_gbl);
+
+ //calculate new phi0
+ double phi0_gbl = phi0 + xTPrimeCorr - corrPer.x() * C_gbl;
+
+ LOGGER.info("qOverP=" + qOverP + " qOverPCorr=" + qOverPCorr + " qOverP_gbl=" + qOverP_gbl + " ==> pGbl=" + 1.0 / qOverP_gbl + " C_gbl=" + C_gbl);
+
+ LOGGER.info(String.format("corrected helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", dca_gbl, z0_gbl, C_gbl, slope_gbl, phi0_gbl, Math.abs(1 / qOverP_gbl)));
+
+ double parameters_gbl[] = new double[5];
+ parameters_gbl[HelicalTrackFit.dcaIndex] = dca_gbl;
+ parameters_gbl[HelicalTrackFit.phi0Index] = phi0_gbl;
+ parameters_gbl[HelicalTrackFit.curvatureIndex] = C_gbl;
+ parameters_gbl[HelicalTrackFit.z0Index] = z0_gbl;
+ parameters_gbl[HelicalTrackFit.slopeIndex] = slope_gbl;
+
+ return parameters_gbl;
+
+ }
+
+
+
+
public static BasicMatrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) {
/*
@@ -105,4 +316,172 @@
throw new UnsupportedOperationException("Should not happen. This problem is only solved with the MaterialSupervisor.");
}
}
+
+
+
+
+
+ /**
+ * Calculate the Jacobian from Curvilinear to Perigee frame.
+ * @param helicalTrackFit - original helix
+ * @param helicalTrackFitAtIPCorrected - corrected helix at this point
+ * @param bfield - magnitude of B-field
+ * @return
+ */
+ public static Matrix getCLToPerigeeJacobian(HelicalTrackFit helicalTrackFit, HpsHelicalTrackFit helicalTrackFitAtIPCorrected, double bfield) {
+
+ /*
+ * This part is taken from:
+ // Strandlie, Wittek, NIMA 566, 2006
+ Matrix covariance_gbl = new Matrix(5, 5);
+ //helpers
+ double Bz = -Constants.fieldConversion * Math.abs(bfield); // TODO sign convention and should it be it scaled from Telsa?
+ double p = Math.abs(1 / qOverP_gbl);
+ double q = Math.signum(qOverP_gbl);
+ double tanLambda = Math.tan(lambda_gbl);
+ double cosLambda = Math.cos(lambda_gbl);
+ // Hep3Vector B = new BasicHep3Vector(0, 0, Bz); // TODO sign convention?
+ Hep3Vector H = new BasicHep3Vector(0, 0, 1);
+ Hep3Vector T = HelixUtils.Direction(helix, 0.);
+ Hep3Vector HcrossT = VecOp.cross(H, T);
+ double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
+ double Q = Bz * q / p;
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+ Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
+ Hep3Vector K = Z;
+ Hep3Vector U = VecOp.mult(-1, J);
+ Hep3Vector V = VecOp.cross(T, U);
+ Hep3Vector I = VecOp.cross(J, K);
+ Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T));
+ double UdotI = VecOp.dot(U, I);
+ double NdotV = VecOp.dot(N, V);
+ double NdotU = VecOp.dot(N, U);
+ double TdotI = VecOp.dot(T, I);
+ double VdotI = VecOp.dot(V, I);
+ double VdotK = VecOp.dot(V, K);
+ covariance_gbl.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK / TdotI);
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1);
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), -1 * Bz / cosLambda);
+ // covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1 * q * Bz * tanLambda / (p * cosLambda));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), q * Bz * alpha * Q * tanLambda * UdotI * NdotV / (p * cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), q * Bz * alpha * Q * tanLambda * VdotI * NdotV / (p * cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1 / TdotI);
+ covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1);
+ covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), alpha * Q * UdotI * NdotV / TdotI);
+ covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), alpha * Q * VdotI * NdotV / TdotI);
+
+ covariance_gbl.print(15, 13);
+ */
+
+ // Sho's magic below
+
+ // Use projection matrix
+ //TODO should this not be the corrected helix?
+ Hep3Matrix perToClPrj = getPerToClPrj(helicalTrackFit);
+ Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
+ double C_gbl = helicalTrackFitAtIPCorrected.curvature();
+ double lambda_gbl = Math.atan(helicalTrackFitAtIPCorrected.slope());
+ double qOverP_gbl = helicalTrackFitAtIPCorrected.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helicalTrackFitAtIPCorrected.slope(), 2)));
+
+ Matrix jacobian = new Matrix(5, 5);
+ jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -clToPerPrj.e(1, 0));
+ jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -clToPerPrj.e(1, 1));
+ jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
+ jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(0, 1) * C_gbl);
+ jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), Constants.fieldConversion * Math.abs(bfield) / Math.cos(lambda_gbl));
+ jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl * Math.tan(lambda_gbl) / Math.cos(lambda_gbl));
+ jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), clToPerPrj.e(2, 0));
+ jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(2, 1));
+ jacobian.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Math.pow(Math.cos(lambda_gbl), -2.0));
+
+ return jacobian;
+ }
+
+
+
+
+
+ /**
+ * Computes the projection matrix from the perigee XY plane variables dca
+ * and z0 into the curvilinear xT,yT,zT frame (U,V,T) with reference point (0,0,0)
+ * for the perigee frame.
+ *
+ * @param htf input helix to find the track direction
+ * @return 3x3 projection matrix
+ */
+ static Hep3Matrix getPerToClPrj(HelicalTrackFit htf) {
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+ Hep3Vector T = HelixUtils.Direction(htf, 0.);
+ Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
+ Hep3Vector K = Z;
+ Hep3Vector U = VecOp.mult(-1, J);
+ Hep3Vector V = VecOp.cross(T, U);
+ Hep3Vector I = VecOp.cross(J, K);
+
+ BasicHep3Matrix trans = new BasicHep3Matrix();
+ trans.setElement(0, 0, VecOp.dot(I, U));
+ trans.setElement(0, 1, VecOp.dot(J, U));
+ trans.setElement(0, 2, VecOp.dot(K, U));
+ trans.setElement(1, 0, VecOp.dot(I, V));
+ trans.setElement(1, 1, VecOp.dot(J, V));
+ trans.setElement(1, 2, VecOp.dot(K, V));
+ trans.setElement(2, 0, VecOp.dot(I, T));
+ trans.setElement(2, 1, VecOp.dot(J, T));
+ trans.setElement(2, 2, VecOp.dot(K, T));
+ return trans;
+
+ /*
+ Hep3Vector B = new BasicHep3Vector(0, 0, 1); // TODO sign convention?
+ Hep3Vector H = VecOp.mult(1 / bfield, B);
+ Hep3Vector T = HelixUtils.Direction(helix, 0.);
+ Hep3Vector HcrossT = VecOp.cross(H, T);
+ double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
+ double Q = Math.abs(bfield) * q / p;
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+ Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
+ Hep3Vector K = Z;
+ Hep3Vector U = VecOp.mult(-1, J);
+ Hep3Vector V = VecOp.cross(T, U);
+ Hep3Vector I = VecOp.cross(J, K);
+ Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T)); //-cross(T,H)/alpha = -cross(T,Z) = -J
+ double UdotI = VecOp.dot(U, I); // 0,0
+ double NdotV = VecOp.dot(N, V); // 1,1?
+ double NdotU = VecOp.dot(N, U); // 0,1?
+ double TdotI = VecOp.dot(T, I); // 2,0
+ double VdotI = VecOp.dot(V, I); // 1,0
+ double VdotK = VecOp.dot(V, K); // 1,2
+ */
+ }
+
+
+
+
+
+ private static BasicMatrix getPerParVector(double kappa, double theta, double phi, double d0, double z0) {
+ BasicMatrix perPar = new BasicMatrix(1, 5);
+ perPar.setElement(0, 0, kappa);
+ perPar.setElement(0, 1, theta);
+ perPar.setElement(0, 2, phi);
+ perPar.setElement(0, 3, d0);
+ perPar.setElement(0, 4, z0);
+ return perPar;
+ }
+
+
+
+
+
+ private static BasicMatrix getPerParVector(HelicalTrackFit htf, double B) {
+ if (htf != null) {
+ double kappa = -1.0 * Math.signum(B) / htf.R();
+ double theta = Math.PI / 2.0 - Math.atan(htf.slope());
+ return getPerParVector(kappa, theta, htf.phi0(), htf.dca(), htf.z0());
+ }
+ return new BasicMatrix(1, 5);
+ }
+
+
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java Mon Jan 4 17:41:58 2016
@@ -14,10 +14,12 @@
import java.util.logging.Formatter;
import java.util.logging.Level;
import java.util.logging.Logger;
+
import org.apache.commons.math3.util.Pair;
-
import org.hps.recon.tracking.TrackUtils;
+
import static org.hps.recon.tracking.gbl.MakeGblTracks.makeCorrectedTrack;
+
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
@@ -52,6 +54,8 @@
private final String track2GblTrackRelationName = "TrackToGBLTrack";
private final String gblTrack2StripRelationName = "GBLTrackToStripData";
private final String outputTrackCollectionName = "GBLTracks";
+ private final String trackRelationCollectionName = "MatchedToGBLTrackRelations";
+
private MilleBinary mille;
private String milleBinaryFileName = MilleBinary.DEFAULT_OUTPUT_FILE_NAME;
@@ -184,6 +188,8 @@
LOGGER.info(trackFits.size() + " fitted GBL tracks before adding to event");
List<Track> newTracks = new ArrayList<Track>();
+
+ List<LCRelation> trackRelations = new ArrayList<LCRelation>();
List<GBLKinkData> kinkDataCollection = new ArrayList<GBLKinkData>();
@@ -201,6 +207,10 @@
// Add the track to the list of tracks
newTracks.add(trk.getFirst());
+
+ // Create relation from seed to GBL track
+ trackRelations.add(new BaseLCRelation(fittedTraj.get_seed(), trk.getFirst()));
+
kinkDataCollection.add(trk.getSecond());
kinkDataRelations.add(new BaseLCRelation(trk.getSecond(), trk.getFirst()));
}
@@ -210,6 +220,7 @@
// Put the tracks back into the event and exit
int flag = 1 << LCIOConstants.TRBIT_HITS;
event.put(outputTrackCollectionName, newTracks, Track.class, flag);
+ event.put(trackRelationCollectionName, trackRelations, LCRelation.class, 0);
event.put(GBLKinkData.DATA_COLLECTION, kinkDataCollection, GBLKinkData.class, 0);
event.put(GBLKinkData.DATA_RELATION_COLLECTION, kinkDataRelations, LCRelation.class, 0);
@@ -222,12 +233,15 @@
public static FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac, boolean debug) {
// path length along trajectory
double s = 0.;
+ int iLabel;
+
// jacobian to transport errors between points along the path
Matrix jacPointToPoint = new Matrix(5, 5);
jacPointToPoint.UnitMatrix();
// Vector of the strip clusters used for the GBL fit
List<GblPoint> listOfPoints = new ArrayList<GblPoint>();
+ Map<Integer, Double> pathLengthMap = new HashMap<Integer, Double>();
// Store the projection from local to measurement frame for each strip cluster
Map< Integer, Matrix> proL2m_list = new HashMap<Integer, Matrix>();
@@ -236,6 +250,10 @@
//start trajectory at refence point (s=0) - this point has no measurement
GblPoint ref_point = new GblPoint(jacPointToPoint);
listOfPoints.add(ref_point);
+
+ // save path length to each point
+ iLabel = listOfPoints.size();
+ pathLengthMap.put(iLabel, s);
// Loop over strips
int n_strips = hits.size();
@@ -405,7 +423,10 @@
// Add this GBL point to list that will be used in fit
listOfPoints.add(point);
- int iLabel = listOfPoints.size();
+ iLabel = listOfPoints.size();
+
+ // save path length to each point
+ pathLengthMap.put(iLabel, s);
// Update MS covariance matrix
msCov.set(1, 1, msCov.get(1, 1) + scatErr.get(0) * scatErr.get(0));
@@ -493,8 +514,10 @@
LOGGER.fine("locPar " + aCorrection.toString());
-//
- return new FittedGblTrajectory(traj, dVals[0], iVals[0], dVals[1]);
+ FittedGblTrajectory fittedTraj = new FittedGblTrajectory(traj, dVals[0], iVals[0], dVals[1]);
+ fittedTraj.setPathLengthMap(pathLengthMap);
+
+ return fittedTraj;
}
@Override
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java Mon Jan 4 17:41:58 2016
@@ -5,19 +5,19 @@
import hep.physics.vec.Hep3Matrix;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
+
import java.util.ArrayList;
import java.util.Collection;
-import java.util.Collections;
-import java.util.Comparator;
import java.util.List;
import java.util.logging.Level;
import java.util.logging.Logger;
+
import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.CoordinateTransformations;
+import org.hps.recon.tracking.HpsHelicalTrackFit;
import org.hps.recon.tracking.MultipleScattering;
import org.hps.recon.tracking.TrackType;
import org.hps.recon.tracking.TrackUtils;
-import static org.hps.recon.tracking.gbl.GBLOutput.getPerToClPrj;
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
@@ -73,7 +73,7 @@
}
// Set state at vertex
- Pair<double[], SymmetricMatrix> correctedHelixParams = getGblCorrectedHelixParameters(helix, fittedTraj.get_traj(), bfield, FittedGblTrajectory.GBLPOINT.IP);
+ Pair<double[], SymmetricMatrix> correctedHelixParams = getGblCorrectedHelixParameters(helix, fittedTraj, bfield, FittedGblTrajectory.GBLPOINT.IP);
trk.setTrackParameters(correctedHelixParams.getFirst(), bfield);// hack to set the track charge
trk.getTrackStates().clear();
@@ -81,7 +81,7 @@
trk.getTrackStates().add(stateVertex);
// Set state at last point
- Pair<double[], SymmetricMatrix> correctedHelixParamsLast = getGblCorrectedHelixParameters(helix, fittedTraj.get_traj(), bfield, FittedGblTrajectory.GBLPOINT.LAST);
+ Pair<double[], SymmetricMatrix> correctedHelixParamsLast = getGblCorrectedHelixParameters(helix, fittedTraj, bfield, FittedGblTrajectory.GBLPOINT.LAST);
TrackState stateLast = new BaseTrackState(correctedHelixParamsLast.getFirst(), ref, correctedHelixParamsLast.getSecond().asPackedArray(true), TrackState.AtLastHit, bfield);
trk.getTrackStates().add(stateLast);
@@ -109,135 +109,62 @@
* Compute the updated helix parameters and covariance matrix at a given
* point along the trajectory.
*
- * @param helix - original seed track
- * @param traj - fitted GBL trajectory
- * @param point - the point along the track where the result is computed.
+ * @param htf - original seed track
+ * @param fittedGblTraj - fitted GBL trajectory
+ * @param point - the GBL point along the track where the result is computed.
* @return corrected parameters.
*/
- public static Pair<double[], SymmetricMatrix> getGblCorrectedHelixParameters(HelicalTrackFit helix, GblTrajectory traj, double bfield, FittedGblTrajectory.GBLPOINT point) {
-
- // get seed helix parameters
- double qOverP = helix.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helix.slope(), 2)));
- double d0 = -1.0 * helix.dca(); // correct for different sign convention of d0 in perigee frame
- double z0 = helix.z0();
- double phi0 = helix.phi0();
- double lambda = Math.atan(helix.slope());
-
- LOGGER.info("GblPoint: " + point.toString() + "( " + point.name() + ")");
- LOGGER.info(String.format("original helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", helix.dca(), helix.z0(), helix.curvature(), helix.slope(), helix.phi0(), helix.p(Math.abs(bfield))));
- LOGGER.info("original helix covariance:\n" + helix.covariance());
-
- // get corrections from GBL fit
+ public static Pair<double[], SymmetricMatrix> getGblCorrectedHelixParameters(HelicalTrackFit htf, FittedGblTrajectory fittedGblTraj, double bfield, FittedGblTrajectory.GBLPOINT point) {
+
+ LOGGER.info("Get corrected parameters at GblPoint: " + point.toString() + "( " + point.name() + ")");
+
+ // Get corrections from GBL fit
Vector locPar = new Vector(5);
SymMatrix locCov = new SymMatrix(5);
- int pointIndex;
- if (point.compareTo(FittedGblTrajectory.GBLPOINT.IP) == 0) {
- pointIndex = 1;
- } else if (point.compareTo(FittedGblTrajectory.GBLPOINT.LAST) == 0) {
- pointIndex = traj.getNumPoints();
- } else {
- throw new RuntimeException("This GBLPOINT " + point.toString() + "( " + point.name() + ") is not valid");
- }
-
- traj.getResults(pointIndex, locPar, locCov); // vertex point
-// locCov.print(10, 8);
- double qOverPCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
- double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
- double yTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue());
- double xTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XT.getValue());
- double yTCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.YT.getValue());
-// System.out.println((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr);
- LOGGER.info((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr);
-
- // calculate new d0 and z0
-// Hep3Matrix perToClPrj = traj.get_track_data().getPrjPerToCl();
- Hep3Matrix perToClPrj = getPerToClPrj(helix);
-
- Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
- Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
-
- //d0
- double d0_corr = corrPer.y();
- double dca_gbl = -1.0 * (d0 + d0_corr);
-
- //z0
- double z0_corr = corrPer.z();
- double z0_gbl = z0 + z0_corr;
-
- //calculate new slope
- double lambda_gbl = lambda + yTPrimeCorr;
- double slope_gbl = Math.tan(lambda_gbl);
-
- // calculate new curvature
- double qOverP_gbl = qOverP + qOverPCorr;
-// double pt_gbl = (1.0 / qOverP_gbl) * Math.cos(lambda_gbl);
-// double C_gbl = Constants.fieldConversion * Math.abs(bfield) / pt_gbl;
- double C_gbl = Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl / Math.cos(lambda_gbl);
-
- //calculate new phi0
- double phi0_gbl = phi0 + xTPrimeCorr - corrPer.x() * C_gbl;
-
- LOGGER.info("qOverP=" + qOverP + " qOverPCorr=" + qOverPCorr + " qOverP_gbl=" + qOverP_gbl + " ==> pGbl=" + 1.0 / qOverP_gbl + " C_gbl=" + C_gbl);
-
- LOGGER.info(String.format("corrected helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", dca_gbl, z0_gbl, C_gbl, slope_gbl, phi0_gbl, Math.abs(1 / qOverP_gbl)));
-
- /*
- // Strandlie, Wittek, NIMA 566, 2006
- Matrix covariance_gbl = new Matrix(5, 5);
- //helpers
- double Bz = -Constants.fieldConversion * Math.abs(bfield); // TODO sign convention and should it be it scaled from Telsa?
- double p = Math.abs(1 / qOverP_gbl);
- double q = Math.signum(qOverP_gbl);
- double tanLambda = Math.tan(lambda_gbl);
- double cosLambda = Math.cos(lambda_gbl);
- // Hep3Vector B = new BasicHep3Vector(0, 0, Bz); // TODO sign convention?
- Hep3Vector H = new BasicHep3Vector(0, 0, 1);
- Hep3Vector T = HelixUtils.Direction(helix, 0.);
- Hep3Vector HcrossT = VecOp.cross(H, T);
- double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
- double Q = Bz * q / p;
- Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
- Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
- Hep3Vector K = Z;
- Hep3Vector U = VecOp.mult(-1, J);
- Hep3Vector V = VecOp.cross(T, U);
- Hep3Vector I = VecOp.cross(J, K);
- Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T));
- double UdotI = VecOp.dot(U, I);
- double NdotV = VecOp.dot(N, V);
- double NdotU = VecOp.dot(N, U);
- double TdotI = VecOp.dot(T, I);
- double VdotI = VecOp.dot(V, I);
- double VdotK = VecOp.dot(V, K);
- covariance_gbl.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK / TdotI);
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1);
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), -1 * Bz / cosLambda);
- // covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1 * q * Bz * tanLambda / (p * cosLambda));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), q * Bz * alpha * Q * tanLambda * UdotI * NdotV / (p * cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), q * Bz * alpha * Q * tanLambda * VdotI * NdotV / (p * cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1 / TdotI);
- covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1);
- covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), alpha * Q * UdotI * NdotV / TdotI);
- covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), alpha * Q * VdotI * NdotV / TdotI);
-
- covariance_gbl.print(15, 13);
- */
- // Sho's magic
- Matrix jacobian = new Matrix(5, 5);
- jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -clToPerPrj.e(1, 0));
- jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -clToPerPrj.e(1, 1));
- jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
- jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(0, 1) * C_gbl);
- jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), Constants.fieldConversion * Math.abs(bfield) / Math.cos(lambda_gbl));
- jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl * Math.tan(lambda_gbl) / Math.cos(lambda_gbl));
- jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), clToPerPrj.e(2, 0));
- jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(2, 1));
- jacobian.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Math.pow(Math.cos(lambda_gbl), -2.0));
-
-// jacobian.print(15, 13);
+
+ // Extract the corrections to the track parameters and the covariance matrix from the GBL trajectory
+ fittedGblTraj.getResults(point, locPar, locCov);
+
+ // Use the super class to keep track of reference point of the helix
+ HpsHelicalTrackFit helicalTrackFit = new HpsHelicalTrackFit(htf);
+ double[] refIP = helicalTrackFit.getRefPoint();
+
+ // Calculate new reference point for this point
+ // This is the intersection of the helix with the plane
+ // The trajectory has this information already in the form of a map between GBL point and path length
+ double pathLength = fittedGblTraj.getPathLength(point);
+ Hep3Vector refPointVec = HelixUtils.PointOnHelix(helicalTrackFit, pathLength);
+ double[] refPoint = new double[]{refPointVec.x(), refPointVec.y()};
+
+ LOGGER.info("pathLength " + pathLength + " -> refPointVec " + refPointVec.toString());
+
+ // Propagate the helix to new reference point
+ double[] helixParametersAtPoint = TrackUtils.getParametersAtNewRefPoint(refPoint, helicalTrackFit);
+
+ // Create a new helix with the new parameters and the new reference point
+ HpsHelicalTrackFit helicalTrackFitAtPoint = new HpsHelicalTrackFit(helixParametersAtPoint, helicalTrackFit.covariance(),
+ helicalTrackFit.chisq(), helicalTrackFit.ndf(), helicalTrackFit.PathMap(),
+ helicalTrackFit.ScatterMap(), refPoint);
+
+ // find the corrected perigee track parameters at this point
+ double[] helixParametersAtPointCorrected = GblUtils.getCorrectedPerigeeParameters(locPar, helicalTrackFitAtPoint, bfield);
+
+ // create a new helix
+ HpsHelicalTrackFit helicalTrackFitAtPointCorrected = new HpsHelicalTrackFit(helixParametersAtPointCorrected, helicalTrackFit.covariance(),
+ helicalTrackFit.chisq(), helicalTrackFit.ndf(), helicalTrackFit.PathMap(),
+ helicalTrackFit.ScatterMap(), refPoint);
+
+ // change reference point back to the original one
+ double[] helixParametersAtIPCorrected = TrackUtils.getParametersAtNewRefPoint(refIP, helicalTrackFitAtPointCorrected);
+
+ // create a new helix for the new parameters at the IP reference point
+ HpsHelicalTrackFit helicalTrackFitAtIPCorrected = new HpsHelicalTrackFit(helixParametersAtIPCorrected, helicalTrackFit.covariance(),
+ helicalTrackFit.chisq(), helicalTrackFit.ndf(), helicalTrackFit.PathMap(),
+ helicalTrackFit.ScatterMap(), refIP);
+
+
+ // Calculate the updated covariance
+ Matrix jacobian = GblUtils.getCLToPerigeeJacobian(helicalTrackFit, helicalTrackFitAtIPCorrected, bfield);
Matrix helixCovariance = jacobian.times(locCov.times(jacobian.transpose()));
SymmetricMatrix cov = new SymmetricMatrix(5);
for (int i = 0; i < 5; i++) {
@@ -248,13 +175,75 @@
}
}
LOGGER.info("corrected helix covariance:\n" + cov);
-
+
+ double parameters_gbl[] = helicalTrackFitAtIPCorrected.parameters();
+
+ /*
+
+ // Get seed helix parameters
+ double qOverP = helicalTrackFit.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helicalTrackFit.slope(), 2)));
+ double d0 = -1.0 * helicalTrackFit.dca(); // correct for different sign convention of d0 in perigee frame
+ double z0 = helicalTrackFit.z0();
+ double phi0 = helicalTrackFit.phi0();
+ double lambda = Math.atan(helicalTrackFit.slope());
+
+ // calculate new d0 and z0
+// Hep3Matrix perToClPrj = traj.get_track_data().getPrjPerToCl();
+ Hep3Matrix perToClPrj = getPerToClPrj(helicalTrackFit);
+
+ Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
+ Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
+
+ //d0
+ double d0_corr = corrPer.y();
+ double dca_gbl = -1.0 * (d0 + d0_corr);
+
+ //z0
+ double z0_corr = corrPer.z();
+ double z0_gbl = z0 + z0_corr;
+
+ //calculate new slope
+ double lambda_gbl = lambda + yTPrimeCorr;
+ double slope_gbl = Math.tan(lambda_gbl);
+
+ // calculate new curvature
+ double qOverP_gbl = qOverP + qOverPCorr;
+// double pt_gbl = (1.0 / qOverP_gbl) * Math.cos(lambda_gbl);
+// double C_gbl = Constants.fieldConversion * Math.abs(bfield) / pt_gbl;
+ double C_gbl = Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl / Math.cos(lambda_gbl);
+
+ //calculate new phi0
+ double phi0_gbl = phi0 + xTPrimeCorr - corrPer.x() * C_gbl;
+
+ LOGGER.info("qOverP=" + qOverP + " qOverPCorr=" + qOverPCorr + " qOverP_gbl=" + qOverP_gbl + " ==> pGbl=" + 1.0 / qOverP_gbl + " C_gbl=" + C_gbl);
+
+ LOGGER.info(String.format("corrected helix: d0=%f, z0=%f, omega=%f, tanlambda=%f, phi0=%f, p=%f", dca_gbl, z0_gbl, C_gbl, slope_gbl, phi0_gbl, Math.abs(1 / qOverP_gbl)));
+ */
+
+
+
+// jacobian.print(15, 13);
+// Matrix helixCovariance = jacobian.times(locCov.times(jacobian.transpose()));
+// SymmetricMatrix cov = new SymmetricMatrix(5);
+// for (int i = 0; i < 5; i++) {
+// for (int j = 0; j < 5; j++) {
+// if (i >= j) {
+// cov.setElement(i, j, helixCovariance.get(i, j));
+// }
+// }
+// }
+// LOGGER.info("corrected helix covariance:\n" + cov);
+
+ /*
double parameters_gbl[] = new double[5];
parameters_gbl[HelicalTrackFit.dcaIndex] = dca_gbl;
parameters_gbl[HelicalTrackFit.phi0Index] = phi0_gbl;
parameters_gbl[HelicalTrackFit.curvatureIndex] = C_gbl;
parameters_gbl[HelicalTrackFit.z0Index] = z0_gbl;
parameters_gbl[HelicalTrackFit.slopeIndex] = slope_gbl;
+ */
+
+
return new Pair<double[], SymmetricMatrix>(parameters_gbl, cov);
}
|
